甘肅金川礦區(qū)地表三維形變InSAR監(jiān)測

楊旺， 何毅， 張立峰, 王文輝， 陳有東， 陳毅

(1.蘭州交通大學(xué)測繪與地理信息學(xué)院,蘭州 730070； 2.地理國情監(jiān)測技術(shù)應(yīng)用國家地方聯(lián)合工程研究中心,蘭州 730070； 3.甘肅省地理國情監(jiān)測工程實驗室,蘭州 730070)

0 引言

礦產(chǎn)資源在經(jīng)濟發(fā)展中起到了極大的推動作用，其中，礦業(yè)占了國內(nèi)生產(chǎn)總值的5.5%，礦產(chǎn)資源的開發(fā)利用已成為助推中國經(jīng)濟發(fā)展的重要動力[1]。但是，礦產(chǎn)資源的大規(guī)模開發(fā)不可避免地導(dǎo)致地表形變、塌陷和其他潛在的災(zāi)害隱患。中國地質(zhì)環(huán)境監(jiān)測院調(diào)查結(jié)果顯示，2015年末，中國礦山因地下開采形成的采空區(qū)體積達(dá)到12.8億m3，引發(fā)了礦區(qū)坍塌、塌陷等多種地質(zhì)災(zāi)害，造成大量的財產(chǎn)損失和人員傷亡[2]。因此，有必要監(jiān)測礦區(qū)地表形變，并進(jìn)行多維分析。

目前，已有眾多學(xué)者利用GPS和水準(zhǔn)測量等方式對礦區(qū)的地表形變進(jìn)行了監(jiān)測[3-5]，但這些方法存在監(jiān)測范圍小、效率低、成本高和監(jiān)測點密度低的缺點[6]。合成孔徑雷達(dá)干涉測量(interferometric synthetic aperture Radar,InSAR)是一種利用相位差來獲取地表形變信息的高精度測量技術(shù)[7]，該方法監(jiān)測點的密度高、范圍大且效率高[8]。諸多學(xué)者的研究表明，InSAR技術(shù)是監(jiān)測礦區(qū)地表形變的有效手段[9-12]，例如Li等[13]使用TS-InSAR(time series InSAR)技術(shù)結(jié)合建筑物破壞水平經(jīng)驗?zāi)Ｐ蛯茨系V區(qū)的居民住宅進(jìn)行了穩(wěn)定性評價； Ma等[14]利用SBAS-InSAR(small baseline subset InSAR)技術(shù)對神東煤田的開采形變狀況進(jìn)行了調(diào)查和分析?？偨Y(jié)上述研究并結(jié)合文獻(xiàn)綜述獲知，雷達(dá)數(shù)據(jù)的合理選擇對礦區(qū)地表形變的監(jiān)測至關(guān)重要。目前主流的雷達(dá)衛(wèi)星有C波段的Sentinel-1和GF-3，L波段的ALOS，X波段的TerraSAR-X，COSMO-SkyMed和TanDEM-X等。由于Sentinel-1A數(shù)據(jù)開放，且具備影像幅寬大、重返周期短、成像質(zhì)量高等優(yōu)勢，在礦區(qū)地表形變監(jiān)測中得到了廣泛應(yīng)用[15-16]。然而，常規(guī)InSAR技術(shù)只能獲取視線向(line-of-sight,LOS)的形變，聯(lián)合多個軌道模式的影像不僅可以求解地表的三維形變，而且可以緩解由雷達(dá)側(cè)視引起的透視收縮、疊掩和陰影[17]。

金川礦區(qū)是我國最大的鎳礦開采基地之一，其中龍首礦作為金川鎳礦的主力礦山，由于連續(xù)開采，2016年出現(xiàn)閉環(huán)形裂縫，就此停止采礦。2018年首次采用無底柱分段崩落法繼續(xù)開采，該方法具有機械化程度高、開采強度大以及采礦成本相對較低等特點，但引起的地表形變量級較大[18-19]。地表形變是龍首礦嚴(yán)重的次生災(zāi)害，對礦山安全生產(chǎn)活動造成了巨大威脅。因此，有必要探究復(fù)采以來龍首礦區(qū)地表形變的多維時空特征。基于此，丁寧等[20]利用SBAS-InSAR技術(shù)對龍首礦露天礦坑2009—2011年和2015—2016年的地表形變進(jìn)行了監(jiān)測，分析了三維形變特征，但該研究存在監(jiān)測時段不連續(xù)以及最新研究資料欠缺的疏漏，并缺乏對金川礦區(qū)整體形變的監(jiān)測分析； 李強[18]聯(lián)合SBAS-InSAR、無人機傾斜攝影測量和水準(zhǔn)測量對龍首礦工業(yè)化試驗區(qū)的地表形變特征進(jìn)行了研究。仍存在監(jiān)測范圍小、數(shù)據(jù)一致性差、形變機理揭示不夠精細(xì)的缺陷[20]。因此，本文基于升降軌模式下3個軌道(升軌128、降軌33和135)的Sentinel-1A數(shù)據(jù)，利用SBAS-InSAR技術(shù)并結(jié)合先驗條件的最小二乘迭代法，獲取了金川礦區(qū)管轄范圍內(nèi)所有生產(chǎn)區(qū)2018—2020年間的地表三維形變時空特征，并就其形變機理展開了全面剖析。

1 研究區(qū)概況及數(shù)據(jù)源

1.1 研究區(qū)概況

金川礦區(qū)位于龍首山斷裂帶，北接阿拉善臺塊，南鄰北祁連山縫合帶[20](圖1)。該礦是金昌市最早發(fā)現(xiàn)并開采利用的老礦山[21]，擁有我國最大的鎳儲蓄量，也是我國鉑金屬的主要產(chǎn)地。該礦體呈NW-SE走向分布，長約6 000 m，寬約300 m，中間部分厚度約1 000 m[22]。礦區(qū)被一系列SW-NE走向的沖斷層切割，將侵入礦體分成若干段，甘肅省地質(zhì)調(diào)查院將其劃分為Ⅰ，Ⅱ，Ⅲ，Ⅳ礦區(qū)[23]，其中Ⅰ和Ⅲ礦區(qū)由金川公司龍首礦開采和管理，Ⅱ礦區(qū)由西二礦開采和管理，Ⅳ礦區(qū)由東三礦開采和管理[19](圖1(b))。金川礦區(qū)土地覆蓋以沙化土地為主，包含砂、黏土和礫石等，植被覆蓋率較低，后向散射強度較高，開采區(qū)地形起伏程度相對較小，適用 InSAR進(jìn)行地表形變監(jiān)測。

(a) 金昌市地形特征及影像覆蓋范圍(b) 金川礦區(qū)采礦區(qū)分布

1.2 數(shù)據(jù)源

Sentinel-1A衛(wèi)星是歐洲航天局和歐盟委員會推出的哥白尼計劃中的首個空間對地觀測衛(wèi)星，搭載C波段雷達(dá)傳感器[24]，重訪周期為12 d，包含4種工作模式： 條帶模式(strip map,SM)，干涉寬幅模式(interferometric wide,IW)，超寬幅模式(extra-wide swath,EW)以及波束采集模式(wave mode acquisitions,WM)[25]。其中IW模式采用TOPSAR(Terrain observation with progressive scans SAR)技術(shù)成像，幅寬為240 km,空間分辨率為5 m×20 m，影像質(zhì)量較高，噪聲更易處理。此外，相同極化(VV)比交叉極化對地物的介電常數(shù)和表面粗糙度更敏感[26]，因此，本研究選擇IW模式、VV極化的雷達(dá)影像。

本研究選用3個軌道(升軌128，降軌33和135)的Sentinel-1A數(shù)據(jù)反演金川礦區(qū)時序三維形變，下載獲取了覆蓋研究區(qū)128軌道(2018年1月2日—2020年10月6日)81景、33軌道(2018年 1月7日—2020年10月23日)81景和135軌道(2018年3月27日—2020年10月11日)71景影像，影像參數(shù)如表1所示[19]。此外，本文利用數(shù)字高程模型(digital elevation model，DEM)去除地形相位并對形變結(jié)果進(jìn)行地理編碼，DEM選用美國國家航空航天局、美國國家測繪局和德國、意大利航天局聯(lián)合繪制的SRTM(shuttle Radar topography mission)地形數(shù)據(jù)，空間分辨率為90 m[26]。為了驗證本文形變監(jiān)測結(jié)果的有效性和準(zhǔn)確性，選取部分水準(zhǔn)測量數(shù)據(jù)[18]與研究區(qū)的時序InSAR結(jié)果做了檢驗。

2 研究方法

2.1 SBAS-InSAR技術(shù)

SBAS-InSAR是由Berardino等[25]和Lanari等[27]提出的時間序列InSAR分析技術(shù)，該技術(shù)通過將獲取的數(shù)據(jù)進(jìn)行分組配對，得到一系列短的空間基線差分干涉圖集，利用奇異值分解法(singular value decomposition, SVD)將獨立的數(shù)據(jù)集進(jìn)行連接并生成相干目標(biāo)的平均形變速率和時間序列累計沉降值[28]。SBAS-InSAR技術(shù)的優(yōu)勢是充分利用了多時相的SAR數(shù)據(jù)，提高了形變測量的時間采樣率和空間覆蓋率[28]。本文數(shù)據(jù)處理中，SBAS-InSAR流程如圖2所示。具體流程如下:

1)干涉對組合和連接圖生成。將獲取的數(shù)據(jù)按時間和空間閾值進(jìn)行分組配對，綜合考慮時相數(shù)、影像質(zhì)量和計算機的處理能力設(shè)置閾值，本研究中3個軌道時間基線閾值均設(shè)為240 d，空間基線閾值設(shè)為45%。超級主影像由軟件自動選取，原則為時間基線和最小，影像日期分別為： 2018年7月13日(升軌128)、2018年12月28(降軌33)、2019年3月22日(降軌135)，生成時空基線圖如圖3所示。

圖2 SBAS-InSAR處理流程

(a) 128軌道時間基線 (b) 128軌道空間基線

(c) 33軌道時間基線 (d) 33軌道空間基線

(e) 135軌道時間基線 (f) 135軌道空間基線

2)干涉工作流處理，包括干涉圖生成、去平、濾波和相位解纏。首先，利用復(fù)數(shù)像對共軛相乘得到干涉圖，為了提高原始影像的信噪比和相干性，多視比設(shè)置為4∶1； 其次，利用參考DEM去除地形相位； 接著選用Goldstein自適應(yīng)濾波算法進(jìn)行干涉相位濾波； 最后，利用最小費用流解纏方法獲得差分干涉序列的解纏相位圖，解纏相干性閾值設(shè)置為0.2。

3)軌道精煉和重去平。各軌道數(shù)據(jù)分別選取30個地面控制點(ground control points，GCP)，且滿足GCP相干性良好、遠(yuǎn)離形變區(qū)，用殘余相位優(yōu)化方法來估算和去除殘余的固定相位和解纏后還存在的相位坡道。

4)SBAS-InSAR反演。首先，通過選擇和細(xì)化GCP來估算相位坡道[29]； 其次，利用線性模型估計位移速率和殘余地形； 最后，采用SVD方法獲取時間序列地表形變。在估計平均形變速率后，通過利用時間域高通濾波和空間域低通濾波去除大氣相位延遲。

5)反演結(jié)果進(jìn)行地理編碼。DEM作為坐標(biāo)參考，生成坐標(biāo)系(GCS-WGS-84)下的LOS形變和垂直向形變。

2.2 三維形變構(gòu)建方法

2.2.1 InSAR觀測方程

SAR衛(wèi)星多采用近極地軌道，雷達(dá)波束照射到同一地面區(qū)域時飛行狀態(tài)有2種: 一種是由南向北，稱為升軌; 一種是由北向南，稱為降軌(圖4 (a))。由于升降軌使用同一衛(wèi)星平臺，成像時各項物理參數(shù)幾乎未發(fā)生變化，所以影像質(zhì)量一致性較好。升降軌成像最大的區(qū)別就是雷達(dá)波束入射角和衛(wèi)星飛行的方位角不同，利用這一特征，可以構(gòu)建形變監(jiān)測結(jié)果的三維分解模型[30]。InSAR觀測結(jié)果是將真實形變投影到LOS方向上，LOS向的形變又可以分解到垂直、東西、南北向上，其幾何關(guān)系如圖4(b)表示，其關(guān)系式為：

DLOS=DUcosθ-DNsinθcos(α-3π/2)-DEsinθsin(α-3π/2)，

(1)

式中：DLOS為真實形變投影到LOS向形變值；DU，DN和DE分別表示垂直向、南北向和東西向形變(約定垂直向上、北向、東向為正)；θ為SAR傳感器入射角;α為衛(wèi)星飛行方位角。

(a) 衛(wèi)星飛行軌道示意圖 (b) 地表形變?nèi)S幾何關(guān)系

2.2.2 結(jié)合先驗條件的最小二乘迭代法

結(jié)合先驗條件的最小二乘迭代法[30]是利用已知條件或者假設(shè)可信已知條件作為多余觀測，迭代使用最小二乘法求解矛盾方程獲得變量近似解的間接求解方法。Sentinel-1A升降軌道均為近極地軌道，三維形變分解中南北向?qū)π巫冃畔O不敏感[31]，則可以將南北向位移值向極小值逼近作為初始先驗條件，結(jié)合式(1)進(jìn)行最小二乘解算，依據(jù)LOS形變中三維分量的貢獻(xiàn)程度和理論敏感度的高低，依次作為先驗約束條件迭代使用最小二乘求解獲得三維分解值，具體步驟如下：

1)獲得多個軌道時序InSAR結(jié)果之后，根據(jù)式(1)，聯(lián)立三維求解方程組，下標(biāo)表示3個不同軌道，矩陣形式表示為：

(2)

首先設(shè)DN為0，增加多余觀測，此時未知數(shù)為DE和DU，而3個軌道可構(gòu)成3個觀測方程，利用最小二乘法可以求得DE和DU，由于DU分解敏感度較高，認(rèn)為DU數(shù)值是可靠的。此步驟將DN設(shè)為0，因為南北向分解敏感度較低，設(shè)為0可作為條件約束，但不表示南北向位移為0，在前兩步約束求解之后可依據(jù)DE和DU求解DN的近似值。

2) 獲得DU之后，將其代入式(2)構(gòu)成多余觀測，2個未知數(shù)為DN和DE，再次運用最小二乘法求解DN和DE，為了提高DU求解精度，只將DE作為可靠值代入下一步。

3) 設(shè)式(1)—(2)中求得的DU和DE視為真值，此時未知數(shù)只有一個即DN。3個條帶的方程均可求出DN，此時構(gòu)成2個多余觀測，再次運用最小二乘求解得到DN。

3 金川礦區(qū)地表三維形變特征

3.1 地表三維形變速率特征

利用SBAS-InSAR和三維形變解算后，得到三維形變速率(圖5)。圖5(a)為垂直向形變速率，圖中龍首礦、西二礦、東三礦開采區(qū)出現(xiàn)了3個明顯的形變區(qū)A，B和C，從西北至東南方向依次分布，其中A區(qū)形成一個長軸約為2.6 km，短軸約為0.9 km，面積約為1.1 km2的橢圓狀形變區(qū)，垂直向最大形變速率為-132.4 mm/a； B區(qū)形成一個形變面積約為1.8 km2的圓形形變區(qū)，垂直向最大形變速率為-162.8 mm/a； C區(qū)形成一個面積約為0.71 km2的半圓形形變區(qū)，垂直向最大形變速率為-52 mm/a。圖5(b)為東西向形變速率，圖中3個形變區(qū)的西部均向東發(fā)生少量形變，最大形變速率為20.5 mm/a； 東部向西收縮，最大形變速率約為-23.8 mm/a。圖5(c)為南北向形變速率，3個沉降區(qū)西部均向南發(fā)生形變，最大形變速率約為-32.4 mm/a； 東部均向北發(fā)生形變，最大形變速率約為59.2 mm/a。綜上，B形變區(qū)整體形變速率最大，A區(qū)次之，C區(qū)由于東三礦建成投產(chǎn)時間較短，開采力度較小，整體形變速率最小。

(c) 南北向形變速率(d) 三維分解幾何條件

3.2 地表三維形變時序累計形變量特征

為分析礦區(qū)形變過程，解算2018年1月—2020年10月6個時段的三維累積形變值(圖6)，便可知曉該過程。從圖6中可以看出，3個形變區(qū)的形變范圍逐漸擴大，相對于2018年1月，垂直向累積形變值在逐漸增大，研究時段內(nèi)的最大累積形變值為-408.9 m。此外，東西向和南北向的累積形變值也在逐漸增大，其中，東西向形變范圍為-85.4～152.8 mm，南北向形變范圍為-60.6～54.0 mm。金川礦區(qū)垂直向的形變趨勢主要表現(xiàn)為地面沉降，西二礦開采區(qū)的累積形變值最大，龍首礦次之，東三礦開采區(qū)的累積形變值最小。水平位移趨勢均表現(xiàn)為形變區(qū)西北部和東南部向中部收縮。

(a) 2018.1—2018.9(垂直向)(b) 2018.1—2019.3(垂直向)(c) 2018.1—2019.8(垂直向)

(d) 2018.1—2020.1(垂直向)(e) 2018.1—2020.6(垂直向)(f) 2018.1—2020.10(垂直向)

(g) 2018.1—2018.9(東西向)(h) 2018.1—2019.3(東西向)(i) 2018.1—2019.8(東西向)

(j) 2018.1—2020.1(東西向)(k) 2018.1—2020.6(東西向)(l) 2018.1—2020.10(東西向)

(m) 2018.1—2018.9(南北向)(n) 2018.1—2019.3(南北向)(o) 2018.1—2019.8(南北向)

(p) 2018.1—2020.1(南北向)(q) 2018.1—2020.6(南北向)(r) 2018.1—2020.10(南北向)

為進(jìn)一步分析形變區(qū)的時序特征，求取圖5(a)中剖面線a，b，c和d在2018年1月—2020年10月期間6個時段的剖面圖(圖7)。

(a) 剖面線a (b) 剖面線b

(c) 剖面線c (d) 剖面線d

剖面線a的坐標(biāo)方位角為124°，全長8 233 m，呈NW-SE走向同時穿過A，B和C形變較為顯著的區(qū)域，在1 796 m，4 321 m和6 282 m處出現(xiàn)3個沉降漏斗； 剖面線b，c和d坐標(biāo)方位角均為214°，全長分別為2 703 m，3 920 m以及3 274 m，在b，c和d上1 451 m，1 952 m和1 528 m處分別出現(xiàn)沉降漏斗。由剖面線趨勢可以發(fā)現(xiàn)，3個形變區(qū)中沿NW-SE向的剖面形變曲線基本符合Peck公式[32]，NE-SW向的剖面圖顯示： 形變區(qū)A，B和C的東北邊界出現(xiàn)驟降趨勢，可能原因是受到自西向南開采活動的影響且受到斷層的控制。圖7中剖面線橫坐標(biāo)取特定值時，各時相累積形變曲線的間距近似相等，可推斷3個形變區(qū)在近勻速沉降。為驗證該論斷，在形變區(qū)均勻抽樣選取測試點，位置分布如圖5(a)所示。繪制測試點累積形變曲線見圖8，結(jié)果顯示，50個測試點的形變曲線表現(xiàn)為均近直線形變趨勢。因此，可認(rèn)為金川礦區(qū)的3個形變區(qū)，整體在勻速沉降，形變速率伴隨空間位置的不同而發(fā)生變化。

圖8 測試點形變曲線

4 討論

4.1 多軌道相關(guān)性分析

圖9給出了3個軌道垂直向形變速率的相關(guān)性密度散點圖。從圖9中可以看出，升軌128與降軌33和135軌道散點圖的相關(guān)系數(shù)分別為0.68和0.81，由于衛(wèi)星飛行的方向和視角不同，雖然導(dǎo)致了散點圖左右兩側(cè)的偏離度增加，但卻形成了三維形變分解的幾何條件。由于衛(wèi)星飛行方向視角接近，降軌33和135軌道的密度散點圖表現(xiàn)出較高的相關(guān)性，相關(guān)系數(shù)為0.96，相關(guān)直線的斜率接近0.5。綜上，3個不同軌道既構(gòu)成了立體觀測條件，其監(jiān)測結(jié)果又體現(xiàn)出較好的相關(guān)性，說明本文SBAS-InSAR處理結(jié)果可信度較好，三維形變分解方法較為合理。為了進(jìn)一步證明時序InSAR的監(jiān)測性能，選取實測數(shù)據(jù)對本文的監(jiān)測結(jié)果進(jìn)行驗證。

(a) 33軌道與128軌道(b) 135軌道與128軌道(c) 135軌道與33軌道

4.2 三維分解結(jié)果對比驗證

根據(jù)2019年6月和7月的水準(zhǔn)實測數(shù)據(jù)[18]，選擇點位清晰、完整性較好的水準(zhǔn)點1和2位置如圖5(a)所示，用于驗證本文三維形變的結(jié)果。由于水準(zhǔn)實測數(shù)據(jù)記錄了所在點位的當(dāng)月垂直向累積形變值，因此僅對垂直向分解結(jié)果進(jìn)行驗證。提取水準(zhǔn)點位置20 m半徑圓形緩沖區(qū)內(nèi)同時期的SBAS-InSAR時序垂向形變矢量監(jiān)測點，求其平均，將結(jié)果與水準(zhǔn)監(jiān)測值進(jìn)行比較。從圖10可以看出，2個水準(zhǔn)點的分解結(jié)果與水準(zhǔn)測量結(jié)果的趨勢基本相同，最大誤差的絕對值分別為1.3 mm和0.7 mm。說明，SBAS-InSAR的監(jiān)測結(jié)果與水準(zhǔn)測量結(jié)果具有較高的一致性，證明本文實驗參數(shù)的設(shè)置較為合理，表明本文中的三維解算結(jié)果總體上是可靠的。

(a) 水準(zhǔn)點1對比驗證(b) 水準(zhǔn)點2對比驗證

圖10 垂向分解結(jié)果與水準(zhǔn)數(shù)據(jù)的驗證

上述通過水準(zhǔn)數(shù)據(jù)驗證了基于先驗條件的最小二乘迭代法進(jìn)行地表形變?nèi)S分解的可靠性，為凸顯本文方法的優(yōu)勢，將本文方法與常用三維分解算法進(jìn)行對比分析。

聯(lián)合GPS數(shù)據(jù)和InSAR監(jiān)測數(shù)據(jù)的三維形變求解方法可以高效高精度地求解地表三維形變，但是這種方法對數(shù)據(jù)有較高的要求，主要為： 時間域統(tǒng)一、空間域統(tǒng)一、數(shù)據(jù)融合權(quán)重的統(tǒng)一[33]。本文方法采用基于同源數(shù)據(jù)多視角觀測的形變?nèi)S分解策略，擺脫了對數(shù)據(jù)統(tǒng)一性的要求，并且在精度控制上不斷進(jìn)步[34-36]。直接求解法和基于先驗條件的最小二乘迭代法均為基于同源數(shù)據(jù)多視角觀測的形變?nèi)S分解策略，直接求解法是對式(1)直接聯(lián)立方程組進(jìn)行求解。將利用直接求解法計算獲得的三維形變速率與本文方法計算的結(jié)果進(jìn)行對比分析，圖11為直接求解法和基于先驗條件的最小二乘法求解的2019年三維形變速率。從圖11中可以看出，直接求解法解算的三維形變速率東西向和南北向均沒本文方法效果好，尤其南北向的三維形變速率精度較差，這可能是直接求解法中由于InSAR觀測對南北向形變極不敏感，解算結(jié)果必然存在誤差放大現(xiàn)象，誤差最大可達(dá)數(shù)十厘米[37]，獲取的形變并不能代表真實的運動情況，因此，本文方法嘗試在有效降低形變解算誤差的同時，重建了更為真實可靠的三維形變場。此外，本文將圖10中水準(zhǔn)數(shù)據(jù)形變量換算至形變速率，并與直接求解法和本文方法進(jìn)行垂直向?qū)Ρ闰炞C(表2)，結(jié)果顯示： 直接求解法解算結(jié)果整體偏小，最大誤差達(dá)到33.8 mm/a； 本文方法中最大誤差為6.2 mm/a,結(jié)果較為理想，證明了本文方法的可靠性和可行性。依據(jù)精準(zhǔn)的三維形變監(jiān)測結(jié)果，可以預(yù)防并防治相關(guān)災(zāi)害的發(fā)生。然而，形變機理的科學(xué)分析，是減緩形變演化行之有效的方法，鑒于此，就金川礦區(qū)2018年初—2020年末形變發(fā)生的機理進(jìn)行討論。

(a) 本文方法(垂直向)(b) 本文方法(東西向)(c) 本文方法(南北向)

(d) 直接求解法(垂直向)(e) 直接求解法(東西向)(f) 直接求解法(南北向)

表2 垂直向形變速率對比驗證

4.3 金川礦區(qū)地表形變機理

在采礦區(qū)發(fā)生地表形變的關(guān)鍵因素一般涉及到機械化開采導(dǎo)致的地表應(yīng)力變化、斷層的空間分布和巖性的結(jié)構(gòu)強度[19]，因此，將從地表形變與人機采動、研究區(qū)斷層分布和巖性構(gòu)造3方面分析礦區(qū)地表形變機理。

1)礦區(qū)地表形變與人機采動的關(guān)系。2018年復(fù)采后，采用無底柱分段崩落法采礦[19]，礦物從地下被開采出以后，礦房內(nèi)部和周圍的應(yīng)力狀態(tài)失穩(wěn)，為了重新達(dá)到應(yīng)力平衡，在采礦區(qū)上方覆蓋的巖層以及地表地物重力的作用下，巖體的變形和破壞主要表現(xiàn)為坍塌、斷裂、褶皺等形式，繼而引發(fā)地表形變[37]，形成以礦脈為中心的沉降漏斗。隨著開采工作的有序進(jìn)行，金川礦區(qū)的形變趨于穩(wěn)定。因此，2018—2020年金川礦區(qū)地表形變的原因之一可能是人機采動引起。

2)礦區(qū)地表形變與斷層的關(guān)系。如圖12所示，金川礦區(qū)北部被F1斷裂帶分割，下盤巖系主要包括侏羅紀(jì)—新近紀(jì)陸相碎屑沉積巖，上盤主要包括元古代變質(zhì)巖和古生代沉積巖[23]，巖層曾被多代超基性巖脈侵入，形成了龍首山礦體[38-39]。隨著采礦進(jìn)行，受F1斷層控制，上盤未發(fā)生明顯形變； 與礦體相間分布的F16，F(xiàn)17，F(xiàn)23與F1斷層相交，傾角變化范圍為60°～80°，將巖體分割為4個礦區(qū)，金川礦區(qū)形變區(qū)域也相應(yīng)地表現(xiàn)出間斷性差異，龍首礦區(qū)北鄰F8斷層和F6斷層，形變區(qū)在此斷層邊界上出現(xiàn)間斷，斷層北部形變量整體減小，在剖面線b中1 000 m附近形變量驟增處即為F6斷層切割； 西南部受到F3斷層控制出現(xiàn)近似弧形形變邊界，整個開采區(qū)在斷層的作用下形變僅局限在近似橢圓范圍內(nèi)。西二礦開采區(qū)除北部分布有F16斷層外周圍再無主要斷層存在，故開采區(qū)以礦脈為中心不斷向周圍擴張形成近似圓形形變區(qū)，但北部仍然受到斷層限制擴張范圍有限，且預(yù)測在開采工作不斷進(jìn)行過程中除北部受到斷層控制外，形變區(qū)還會向外繼續(xù)延伸。東三礦開采區(qū)受到F17和F23斷層的控制形變區(qū)范圍基本局限在2個斷層中間形成呈半圓形分布。通過以上分析可推斷，金川礦區(qū)斷層的分布控制該區(qū)域地表形變。

圖12 金川礦區(qū)地質(zhì)構(gòu)造與累積形變疊加簡圖[23]

3)礦區(qū)地表形變與巖性的關(guān)系。金川礦區(qū)巖性復(fù)雜，區(qū)域差異性較大，龍首礦和西二礦開采區(qū)東北部主要為大理巖巖組、片麻巖巖組，主要由白云石和方解石組成，含有少量滑石、石英、方美石和云母等組成，顆粒均勻細(xì)小，質(zhì)地均勻，硬度、結(jié)構(gòu)強度較小，在應(yīng)力和侵蝕作用下易于發(fā)生形變[28]； 西南部主要為片麻巖和復(fù)片麻巖巖組，主要由長石、石英和云母等組成，此外還包括角閃石、輝石等，顆粒大小不均，質(zhì)地不均，結(jié)構(gòu)強度整體較強。除巖層結(jié)構(gòu)外，受大陸性溫帶干旱氣候和地表巖性的影響，金川礦區(qū)地表形成較厚的一層砂、礫石、黏土和黃土混合松散堆積層，該層整體密度較低，質(zhì)松軟、結(jié)構(gòu)強度較差，堆積層長期受到雨水侵蝕和滲透，降低了黃土、黏土、砂的黏聚強度和基質(zhì)吸力，破壞了堆積層的穩(wěn)定狀態(tài)，砂土覆蓋層軟化，濕陷性增加[29]，致使礦區(qū)人工建筑基礎(chǔ)發(fā)生坍塌，進(jìn)而導(dǎo)致覆蓋層表面發(fā)生形變。金川礦區(qū)地表形變會隨著人工開采繼續(xù)進(jìn)行，且西部開采區(qū)的巖體穩(wěn)定性大于東部開采區(qū)[39]，隨著鎂鐵超鎂鐵侵入巖和超基性巖的多期侵入，節(jié)理裂隙發(fā)育破碎，東部開采區(qū)的沉降也將會更加明顯。因此，推斷巖性特征也是金川礦區(qū)地表形變的控制性因素。

5 結(jié)論

本文采用SBAS-InSAR技術(shù)分別利用3個軌道(升軌128、降軌33和135)的Sentinel-1A數(shù)據(jù)獲取金川礦區(qū)地表形變速率和形變量，并利用結(jié)合先驗條件的最小二乘迭代法解算礦區(qū)三維地表形變，對3個形變區(qū)空間變化規(guī)律和地表形變機理進(jìn)行了分析。主要結(jié)論如下：

1)通過3個軌道垂直向累計形變值相關(guān)性分析實現(xiàn)解算結(jié)果的內(nèi)部互檢，并使用水準(zhǔn)實測數(shù)據(jù)對三維形變分解后垂直向累計形變值進(jìn)行驗證，保證了地面形變監(jiān)測的可靠性。

2)金川礦區(qū)出現(xiàn)3個形變區(qū)，其中B區(qū)形變范圍、整體形變速率和形變量均為最大； A區(qū)次之； C區(qū)最小。形變方式以沉降為主，最大年平均垂直向形變速率達(dá)到-162.8 mm/a，最大垂直向累計形變量達(dá)到-408.9 mm。此外，形變區(qū)也同時發(fā)生了水平面上的位移，主要表現(xiàn)為3個形變區(qū)西北部和東南部向形變區(qū)中部收縮。

3)地質(zhì)斷層和巖性構(gòu)造是礦區(qū)形變的控制性因素，F(xiàn)1斷層限制了礦區(qū)地表形變向東北向延伸，F(xiàn)16，F(xiàn)17和F23斷層的存在使得開采區(qū)形變出現(xiàn)間斷性差異。人機采動是礦區(qū)地表形變的主要因素，無底柱分段崩落法開采造成了開采區(qū)在30個月內(nèi)近似均勻形變。