基于時(shí)序InSAR技術(shù)的甘肅金川銅鎳礦區(qū)地表沉降研究

李強(qiáng)，王衛(wèi)紅，楊寧，熊高翔 ,雷歡歡，李凱

(1.四川省冶金地質(zhì)勘查院，成都 610051；2.國家遙感中心綿陽科技城分部，四川 綿陽 621010；3.西南科技大學(xué)環(huán)境與資源學(xué)院，四川 綿陽 621010)

0 引言

甘肅省金川銅鎳礦區(qū)位于阿拉善臺(tái)塊的南部邊緣隆起區(qū)，南鄰北祁連山加里東地槽邊緣拗陷帶。礦區(qū)生產(chǎn)采用無底柱分段崩落式采礦方法，在采礦爆破時(shí)會(huì)加速地表的沉陷，形成地裂縫。同時(shí)地下開采會(huì)引起巖層及地表產(chǎn)生移動(dòng)和變形，誘發(fā)各種地質(zhì)災(zāi)害[1]。因此，掌握地表變形規(guī)律，能夠有效保證人民群眾的生命財(cái)產(chǎn)安全、國家經(jīng)濟(jì)的可持續(xù)發(fā)展以及社會(huì)的穩(wěn)定[2]，對(duì)礦區(qū)災(zāi)害預(yù)警、開采沉陷控制和治理有重要指導(dǎo)意義。

礦區(qū)地表沉陷是一個(gè)長(zhǎng)期過程，具有地形地質(zhì)條件復(fù)雜、沉降速度快等特點(diǎn)，所以利用GPS測(cè)量、大地水準(zhǔn)測(cè)量和全站儀測(cè)量等傳統(tǒng)的地表形變監(jiān)測(cè)方法監(jiān)測(cè)礦區(qū)沉陷區(qū)面臨諸多難題。第一，工業(yè)法采礦多采用放炮方式，導(dǎo)致水準(zhǔn)點(diǎn)穩(wěn)定性存在問題，加之采礦是一個(gè)長(zhǎng)期過程，傳統(tǒng)測(cè)量周期長(zhǎng)，工作量大，耗時(shí)耗力[3]；第二，傳統(tǒng)監(jiān)測(cè)方式只是在局部得到沉降信息，不利于體現(xiàn)整個(gè)礦區(qū)的沉降規(guī)律[4]；第三，傳統(tǒng)測(cè)量具有地域局限性，需要預(yù)先知道沉降位置才能布置觀測(cè)點(diǎn)。合成孔徑雷達(dá)干涉(Interferometric Synthetic Aperture Radar)技術(shù)，即InSAR測(cè)量技術(shù)是近30年來迅速發(fā)展起來的一種空間對(duì)地觀測(cè)技術(shù)，基于合成孔徑雷達(dá)復(fù)數(shù)影像的相位信息獲取地表變化信息，可用于監(jiān)測(cè)厘米級(jí)或更微小的地球表面形變[5]。相比可見光和紅外遙感等傳統(tǒng)測(cè)量技術(shù)，合成孔徑雷達(dá)干涉技術(shù)具有覆蓋范圍廣、精度高、不受天氣狀況影響、空間采樣密度高等優(yōu)勢(shì)。目前，InSAR技術(shù)已廣泛應(yīng)用于區(qū)域地表形變監(jiān)測(cè)[6]、隱患點(diǎn)探測(cè)[7]、冰川凍融監(jiān)測(cè)[8]、洪水淹沒區(qū)變化監(jiān)測(cè)[9]等多個(gè)領(lǐng)域。

對(duì)于InSAR技術(shù)在礦山監(jiān)測(cè)方面的應(yīng)用，國內(nèi)外學(xué)者展開了大量研究，并取得高精度形變的觀測(cè)結(jié)果。1993年，Massonnet等[10]利用ERS-1衛(wèi)星的SAR圖像，對(duì)加利福尼亞州蘭德斯地震產(chǎn)生的運(yùn)動(dòng)進(jìn)行監(jiān)測(cè)研究，監(jiān)測(cè)精度高于以前的空間成像技術(shù)，精度優(yōu)于3 cm；2005年，Tomás等[11]在用D-InSAR技術(shù)監(jiān)測(cè)1993—1995年期間塞古拉河周邊地面沉降時(shí)發(fā)現(xiàn)，某些區(qū)域的沉降達(dá)8 cm，對(duì)預(yù)測(cè)和治理該區(qū)域地面沉降有一定指導(dǎo)作用；2010年，Herrera等[12]利用相干像素技術(shù)(CPT)繪制和監(jiān)測(cè)La Union礦區(qū)的地面運(yùn)動(dòng)，證明此方法可用于廢棄礦區(qū)大部分地面運(yùn)動(dòng)；2014年，Saygin等[13]分析了D-InSAR的缺陷，利用PS-InSAR技術(shù)對(duì)土耳其西北部成交地區(qū)采煤沉陷區(qū)域進(jìn)行監(jiān)測(cè)，充分說明了PS-InSAR技術(shù)在礦區(qū)沉降監(jiān)測(cè)具有指導(dǎo)性意義；2017年，朱建軍等[14]著重對(duì)InSAR監(jiān)測(cè)方法進(jìn)行了分類，分析了PSInSAR、SBAS-InSAR及D-InSAR，提出改進(jìn)的干涉測(cè)量技術(shù)可以用于礦區(qū)的沉降監(jiān)測(cè)；2018年，李達(dá)等[15]利用SBAS-InSAR技術(shù)對(duì)陜西省榆林市某礦區(qū)西北部進(jìn)行地表沉降監(jiān)測(cè)，證明了SBAS-InSAR技術(shù)在礦山沉降監(jiān)測(cè)的可靠性，為礦區(qū)地表沉降監(jiān)測(cè)提供新手段；2019年，黃長(zhǎng)軍等[16]采用持續(xù)散射干涉測(cè)量法(PSI)，避免了傳統(tǒng)干涉測(cè)量的限制，有效揭示了山東葛亭煤礦開采沉陷的演變過程，證明PSI法是測(cè)量采礦區(qū)沉陷的有效方法，可為礦山災(zāi)害提供預(yù)警。

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

金川銅鎳礦床由4個(gè)礦區(qū)組成，龍首礦是主力礦山之一，20世紀(jì)60年代建成投產(chǎn)，在生產(chǎn)鎳銅的同時(shí)回收鉑族金屬，是我國鉑族金屬的主要生產(chǎn)基地。金川銅鎳硫化物礦床的金屬礦體產(chǎn)于超基性巖體中，部分礦體產(chǎn)于含礦巖體底板與圍巖接觸帶附近。在礦體下方存在斷裂構(gòu)造，由于開采強(qiáng)度較大，故產(chǎn)生的地質(zhì)災(zāi)害分布廣。龍首礦西二采區(qū)位于Ⅰ礦區(qū)以西，礦內(nèi)斷裂構(gòu)造發(fā)育，其中F8斷裂規(guī)模最大，角礫巖層厚20～40 m，在4行勘探線直接與礦體接觸，地質(zhì)環(huán)境復(fù)雜。受斷裂影響，在2016年出現(xiàn)明顯的地面開裂現(xiàn)象，地面裂隙主要分布在5行—7行這一區(qū)段，呈環(huán)形閉合環(huán)，直接威脅礦山的生產(chǎn)安全，地下的采礦作業(yè)因此停采。在停采后，為監(jiān)測(cè)其地面沉陷情況，礦區(qū)管理人員對(duì)沉陷區(qū)進(jìn)行維護(hù)，在裂縫附近安裝控制點(diǎn)進(jìn)行重點(diǎn)監(jiān)測(cè)，直到2018年年底沉陷區(qū)趨于穩(wěn)定后，選用無底柱分段崩落法進(jìn)行試驗(yàn)性復(fù)采。在此之前，金川礦區(qū)的地表沉降監(jiān)測(cè)研究通常采用傳統(tǒng)變形測(cè)量方法，故采用InSAR技術(shù)對(duì)金川礦區(qū)進(jìn)行地表監(jiān)測(cè)具有必要性。

本次研究工作利用2018年1月—2018年11月的Sentinel-1A降軌影像和同期的升軌影像對(duì)西二采區(qū)進(jìn)行沉降監(jiān)測(cè)。所有影像均為寬幅單視復(fù)數(shù)影像數(shù)據(jù)(slc)，影像覆蓋范圍如圖1所示，斜距向和方位向分辨率為3.7 m和22.7 m，單景影像覆蓋范圍約為251 km×251 km。本文采用了C波段升軌影像共26景，并采用地面分辨率為30 m的SRTM DEM數(shù)據(jù)來消除干涉相位中的地形相位。

圖1 金昌市區(qū)域Sentinel-1A雷達(dá)衛(wèi)星覆蓋示意及區(qū)位圖

2 數(shù)據(jù)處理

2.1 SBAS-InSAR和時(shí)序D-InSAR基本原理

SBAS-InSAR(小基線集)技術(shù)是對(duì)于長(zhǎng)時(shí)間序列上的一組SAR復(fù)數(shù)圖像根據(jù)一定的基線約束條件進(jìn)行分組，通過控制空間基線的長(zhǎng)度來提高干涉圖的相干性，對(duì)差分干涉圖進(jìn)行多視處理降低噪聲，提取高相干性單元，使用奇異值分解法求得影像序列間地表形變速率的最小范數(shù)最小二乘解[17]，SBAS詳細(xì)技術(shù)原理見文獻(xiàn)[18]。

假設(shè)按時(shí)間順序(t1，…，tn)在同一地點(diǎn)獲取了N+1幅SAR影像，同時(shí)每一短基線子集至少有2幅影像，可生成M個(gè)干涉圖，M滿足以下不等式

(N+1)/2≤M≤[N(N+1)]/2

(1)

假設(shè)j(j=0，1，…，N)干涉圖是由tA、tB時(shí)間獲取的影像生成(tB>tA)，去除平地效應(yīng)與地形相位影響后, 干涉圖j中任意一個(gè)像素的干涉相位可表示為

Δφj=φ(tB)-φ(tA)=φdef，j+φatm，j+φflat，j+φnoise，j

(2)

式中，φ(tB)、φ(tA)分別表示其時(shí)刻相位值；φdef，j表示視線向(LOS)時(shí)刻的形變相位；φfiat，j表示殘余地形相位誤差；φatm，j表示大氣相位誤差；φnoise，j表示噪聲誤差。因此，相對(duì)基準(zhǔn)時(shí)間t0，每幅SAR影像對(duì)應(yīng)的差分相位時(shí)間序列φ與每個(gè)差分干涉圖相位時(shí)間序列Δφ存在以下關(guān)系

Δφ=Aφ

(3)

式中，A為M*N系數(shù)矩陣，每個(gè)差分干涉圖相位時(shí)間序列Δφ=[Δφ(t0)，Δφ(t1)，…，Δφ(tM)]T，每幅SAR影像對(duì)應(yīng)的差分相位時(shí)間序列φ=[φ(t0)，φ(t1)，…，φ(tN-1)]T。

再通過矩陣的奇異值分解法(SVD)求出最小范數(shù)意義上的最小二乘解，所以觀測(cè)時(shí)刻形變量：

(4)

式中，A+=(ATA)-1AT。

合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)(D-InSAR)是干涉測(cè)量技術(shù)和合成孔徑雷達(dá)技術(shù)的結(jié)合，其技術(shù)的核心思想是利用重復(fù)軌道觀測(cè)得到干涉相位，通過差分處理去除2次觀測(cè)相位中的共有量[19](地形相位、平地相位、大氣相位以及噪聲相位)，傳感器視線方向上的沉降值Δr：

Δr=(λ/4π)φdef

(5)

式中，λ為雷達(dá)傳感器工作波長(zhǎng)；φdef為有2次觀測(cè)期間沿著雷達(dá)實(shí)現(xiàn)方向地表形變引起的形變位移。

將公式(5)求出的Δr轉(zhuǎn)化成垂直方向上的沉降值Δh：

Δh=Δr/cosθ

(6)

式中，θ為傳感器入射角。

2.2 數(shù)據(jù)處理

基于SAR Scape 軟件實(shí)現(xiàn)SBAS和時(shí)序D-InSAR的處理。時(shí)序D-InSAR技術(shù)通過影像配準(zhǔn)、去平地效應(yīng)、地理編碼、DEM相位差分等一系列步驟得到形變圖。圖2所示的SBAS技術(shù)流程主要包括小基線組合選擇生成連接圖，生成差分干涉圖，Delaunay MCF解纏，選取GCP點(diǎn)，線性估算位移速率和DEM校正系數(shù)，大氣相位屏的估計(jì)和減法，總形變量估計(jì)等。

圖2 SBAS技術(shù)處理流程圖

3 結(jié)果與分析

研究中使用SBAS技術(shù)和二軌法時(shí)序D-InSAR技術(shù)獲取研究區(qū)2018年1月—11月的地表沉降速率。為了更直觀地展示研究區(qū)的地表形變情況，在envi中進(jìn)行地理編碼，將形變結(jié)果疊加到谷歌影像得到礦區(qū)沉降圖(圖3)，圖3a為SBAS技術(shù)沉降結(jié)果圖，圖3b為時(shí)序D-InSAR技術(shù)沉降結(jié)果圖，將2018年1月作為形變區(qū)域參考，圖中1區(qū)為西二采區(qū)5行—7行，2區(qū)為露天采場(chǎng)老礦坑。為了便于觀察沉降變化趨勢(shì)，在1區(qū)和2區(qū)選擇了沉降漏斗中心做時(shí)序分析，結(jié)果見圖4，分別經(jīng)過D-InSAR與SBAS技術(shù)計(jì)算，得到其形變速率圖(圖5，圖6)，由于1區(qū)相較于2區(qū)沉降較小，故分離出1區(qū)單獨(dú)出圖(即圖5b和圖6b)。以上結(jié)果表明，1區(qū)和2區(qū)中均有沉降漏斗出現(xiàn)，1區(qū)不明顯，2區(qū)則很明顯。

3.1 西二采區(qū)5—7行地表沉降分析

由圖3至圖6可以看出，2018年1月—11月西二采區(qū)整體較為穩(wěn)定，使用時(shí)序D-InSAR方法獲得此區(qū)域最大沉降量為11 mm，年平均最大沉降速率為-14 mm/a，有明顯沉降漏斗的部分，沉降量范圍為3 ～11 mm。SBAS方法獲取最大沉降量為13 mm，年平均最大沉降速率-14 mm/a，沉降量范圍為1～14 mm。在停采之后，2種方法監(jiān)測(cè)結(jié)果顯示，在2018年此區(qū)域仍然存在緩慢的形變，依據(jù)我國形變監(jiān)測(cè)基準(zhǔn)[20]以及本礦區(qū)地質(zhì)情況，沉降速率<30 mm/a時(shí)，處于安全范圍，說明5—7行區(qū)域呈平穩(wěn)態(tài)勢(shì)。由于形變較小，2種方法監(jiān)測(cè)結(jié)果基本一致。

圖3 累積沉降量結(jié)果圖

圖4 累積沉降量曲線

圖5 采用SBAS技術(shù)獲得的沉降速率圖

圖6 采用D-InSAR技術(shù)獲得的沉降速率圖

3.2 老礦坑地表沉降分析

由圖3可以看出，老礦坑?xùn)|南部有明顯沉降，相較于5行—7行，主要是邊坡附近受雨水侵蝕等導(dǎo)致沉降并伴隨滑坡，基于時(shí)序D-InSAR方法獲得其最大沉降量在73 mm，最大沉降速率為-83 mm/a；SBAS方法獲取其最大沉降量在78 mm，而最大沉降速率為-85 mm/a，兩者監(jiān)測(cè)結(jié)果均顯示老礦坑沉降量大。從圖5和圖6可以看出，時(shí)序D-InSAR方法與SBAS方法得到的沉降速率圖明顯不同，其大沉降位置在SBAS方法所出得出的結(jié)果更為精確。同時(shí)試驗(yàn)發(fā)現(xiàn)，老礦坑西北部有小部分區(qū)域存在輕微抬升的現(xiàn)象，一年抬升量約為20 mm，現(xiàn)場(chǎng)踏勘結(jié)果顯示邊坡有輕微隆起，而西南部有滑坡痕跡，老礦坑北部有輕微沉降。從2區(qū)結(jié)果可以看出，2種方法均能監(jiān)測(cè)到老礦坑的沉降現(xiàn)象，但由于時(shí)空基線以及大氣效應(yīng)等誤差[21]，致使時(shí)序D-InSAR監(jiān)測(cè)結(jié)果差于SBAS監(jiān)測(cè)結(jié)果，而SBAS方法能減輕時(shí)空失相關(guān)、大氣延遲等影響[22]能獲取更精確的沉降值以及沉降位置。

4 結(jié)論

(1)經(jīng)過2種方法實(shí)驗(yàn)結(jié)果對(duì)比分析，依據(jù)當(dāng)?shù)氐刭|(zhì)情況以及國家標(biāo)準(zhǔn)，西二采區(qū)5行—7行地表總體呈平穩(wěn)態(tài)勢(shì)，在過去一年最大沉降量為11～13 mm，最大沉降速率為-14 mm/a；老礦坑區(qū)域出現(xiàn)了滑坡現(xiàn)象，最大沉降量在75～78 mm，最大沉降速率為-85 mm/a。

(2)停止采礦后，原采礦區(qū)域地面沉降明顯減緩，基本處于穩(wěn)定狀態(tài)，印證了時(shí)序InSAR在礦山地表沉降監(jiān)測(cè)的可行性。

(3)雖然2種方法都能得到沉陷區(qū)地表沉降規(guī)律，但由于SBAS技術(shù)能夠減輕時(shí)空失相關(guān)、大氣延遲等影響，且從2區(qū)域可以看出SBAS監(jiān)測(cè)結(jié)果相比時(shí)序D-InSAR監(jiān)測(cè)結(jié)果更為理想，故SBAS-InSAR方法更能在后期大幅度沉降出現(xiàn)時(shí)，精確地反映開采沉陷區(qū)沉降速率與沉降位置。