王志崗, 周文韜

(1.四川科技職工大學(xué) 應(yīng)急管理系，四川 成都 610000；2.西南科技大學(xué) 環(huán)境與資源學(xué)院，四川 綿陽 621010；3.國家遙感中心綿陽科技城分部，四川 綿陽 621010)

礦產(chǎn)資源的開采一方面可以促進(jìn)我國經(jīng)濟(jì)建設(shè)和社會(huì)發(fā)展，同時(shí)也會(huì)破壞礦區(qū)地下巖層原有的應(yīng)力結(jié)構(gòu)，從而造成一系列地質(zhì)災(zāi)害[1-4]。2016年3月，甘肅省金昌市金川西二礦1 610 m水平面5～7行范圍的膠結(jié)充填體發(fā)生大面積垮塌并貫通至地表，從而導(dǎo)致停產(chǎn)[5]。從2019年起，該區(qū)域開始使用無底柱分段崩落法進(jìn)行采礦，但該方法對地表影響較大[6]。截至2021年7月，該區(qū)域地表已出現(xiàn)3個(gè)大小不一的塌陷坑，對附近礦區(qū)人員的生命財(cái)產(chǎn)、建筑設(shè)施等具有一定的安全隱患。因此，開展地表形變監(jiān)測工作，對后續(xù)工作開展具有重要的意義。

不同于傳統(tǒng)大地測量方法，合成孔徑雷達(dá)干涉測量(Interferometric Synthetic Aperture Radar, InSAR)技術(shù)自1989年首次提出以來，就因其具有空間分辨率高、覆蓋范圍廣等優(yōu)點(diǎn)在礦區(qū)地表形變監(jiān)測中受到廣泛關(guān)注[7-10]。經(jīng)過國內(nèi)外學(xué)者的研究，一些新的InSAR技術(shù)逐漸出現(xiàn)，例如合成孔徑雷達(dá)差分干涉測量(Differential InSAR, D-InSAR)技術(shù)、永久散射體合成孔徑雷達(dá)差分干涉測量(Persistent Scatterer InSAR, PS-InSAR)技術(shù)、短基線集合成孔徑雷達(dá)差分干涉測量(Small Baseline Subset InSAR, SBAS-InSAR)技術(shù)和人工角反射器差分干涉測量(Corner Reflector InSAR, CR-InSAR)技術(shù)等[11-15]。然而，上述技術(shù)僅能監(jiān)測研究區(qū)地表沿SAR衛(wèi)星視線向(Line of Sight, LOS)形變，而礦區(qū)復(fù)雜的地形信息難以靠單一LOS向形變獲得[16-18]。

因此，本文將利用SBAS-InSAR技術(shù)和InSAR幾何分解原理，基于升降軌Sentinel-1A數(shù)據(jù)解算甘肅省金昌市金川西二礦地表二維形變場，對該區(qū)域地表形變情況進(jìn)行分析，為礦區(qū)開采沉陷監(jiān)測和預(yù)警工作提供思路和方法。

1 研究區(qū)概況及數(shù)據(jù)基礎(chǔ)

1.1 研究區(qū)概況

金昌市(101°4′35″～102°43′40″E，37°47′10″～ 39°00′30″N)位于我國西北部，是甘肅省下轄地級市，也是我國重要的硫化鎳銅礦資源地之一，被譽(yù)為“祖國的鎳都”[19]。本文研究區(qū)選為金昌市金川西二礦，位于金昌市區(qū)西南部，如圖1所示。該區(qū)域降雨量少，日照充足，平均海拔達(dá)1 700 m左右，地下水較少，生態(tài)環(huán)境總體偏弱。礦區(qū)主要以超基性巖型礦體以及圍巖中接觸交代的礦體構(gòu)成。由于地下長期采礦和F8斷層不斷活動(dòng)，研究區(qū)地表形成較大塌陷區(qū)，對地下開采產(chǎn)生極大的影響。

1.2 數(shù)據(jù)基礎(chǔ)

本文采用由歐洲航天局(European Space Agency, ESA)針對哥白尼全球觀測計(jì)劃研制的Sentinel-1衛(wèi)星星座的Sentinel-1A數(shù)據(jù)作為實(shí)驗(yàn)數(shù)據(jù)。自2014年4月發(fā)射以來，該衛(wèi)星可以為用戶提供免費(fèi)的雷達(dá)影像，并能夠?qū)θ蚍秶牡乇硇巫?、海洋、冰川等進(jìn)行監(jiān)測[20-22]。

實(shí)驗(yàn)分別選取Sentinel-1A升軌25景和降軌27景影像，統(tǒng)一時(shí)間跨度為2020.08.19～2021.06.27，實(shí)驗(yàn)數(shù)據(jù)參數(shù)如表1所示。實(shí)驗(yàn)處理過程統(tǒng)一時(shí)間基線閾值120 d，臨界基線閾值45%。其中升軌數(shù)據(jù)以2020.12.17期影像為主影像，生成187個(gè)干涉對；降軌數(shù)據(jù)以2020.08.19期影像為主影像，生成215個(gè)干涉對。同時(shí)，在數(shù)據(jù)處理過程中引入外部數(shù)字高程模型(Digital Elevation Model, DEM)數(shù)據(jù)和精密軌道數(shù)據(jù)來消除地形起伏和軌道誤差帶來的影響。

2 技術(shù)與方法

2.1 SBAS-InSAR技術(shù)原理

在過去幾十年里，InSAR技術(shù)因具有監(jiān)測范圍廣、精度高等優(yōu)點(diǎn)而被廣泛應(yīng)用。但隨著技術(shù)的進(jìn)步，人們發(fā)現(xiàn)傳統(tǒng)InSAR同樣伴隨大氣效應(yīng)、時(shí)空失相干等誤差的缺點(diǎn)[23-24]。SBAS-InSAR技術(shù)由Berardino和Lanari于二十世紀(jì)初首次提出[25]，該方法可以在一定程度上克服傳統(tǒng)InSAR遇到的難題，有效改善時(shí)空失相干和大氣延遲等問題。

SBAS-InSAR技術(shù)的主要原理如下：假設(shè)有t0,t1，…,tn景按時(shí)間序列排序的N+1幅SLC (Single Look Complex)影像，以任意一景影像為超級主影像對所有影像進(jìn)行配準(zhǔn)和重采樣，通過設(shè)定時(shí)間基線閾值和臨界基線閾值可得到M個(gè)干涉對如下：

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

(1)

然后利用精密軌道數(shù)據(jù)和外部DEM數(shù)據(jù)去除平地效應(yīng)和地形效應(yīng)，得到M個(gè)差分干涉圖。若有時(shí)間tA>tB，以t0為參考時(shí)間，則第j幅干涉圖在x處的干涉相位值為：

(2)

式中，λ為雷達(dá)中心波長；φ、μ為干涉相位。則LOS向形變D可表示為：

D(tB,x,r)-D(tA,x,r)=vi(tB-tA)

(3)

式中，vi為從tA至tB的LOS向平均形變速率。由上式可轉(zhuǎn)化為：

Bv=δφ

(4)

式中，B為系數(shù)矩陣。

在SBAS-InSAR技術(shù)解算過程中，由于多主影像模式容易導(dǎo)致矩陣B秩虧，因此采用奇異值分解法(Singular Value Decomposition, SVD)求解，使得計(jì)算有唯一解[26]。該方法的技術(shù)流程如圖2所示。

2.2 基于升降軌InSAR獲取礦區(qū)地表二維形變場

根據(jù)2.1節(jié)的技術(shù)方法，我們僅能獲取單一軌道沿衛(wèi)星LOS向的一維形變，LOS向形變的分解關(guān)系如圖3所示。在復(fù)雜的礦區(qū)地表形變過程中，LOS向的一維形變往往不能準(zhǔn)確反映真實(shí)形變特征。當(dāng)結(jié)合多軌道SAR衛(wèi)星數(shù)據(jù)聯(lián)合解算時(shí)，則可能獲取研究區(qū)的多維形變場。

由圖3可知，地面某點(diǎn)在東西向、南北向和垂直向的位移是LOS向形變分量。當(dāng)同時(shí)考慮地面點(diǎn)水平和垂直位移時(shí)，則有下式為：

DLOS=-sin(θ)cos(α)DE+sin(θ)sin(α)DN+DVcos(θ)

(5)

式中，θ為衛(wèi)星入射角；α為衛(wèi)星方位角。

圖3 InSAR幾何分解原理

令S=[Sx,Sy,Sz]，則Sx=-sin(θ)cos(α),Sy=sin(θ)sin(α)和Sz=cos(θ)分別為三個(gè)方向的投影矢量。但通常情況下，SAR衛(wèi)星采用極軌法飛行(即近南北向飛行)，導(dǎo)致InSAR技術(shù)對南北向形變監(jiān)測并不敏感，因此可忽略南北向形變，利用升軌和降軌監(jiān)測結(jié)果解算礦區(qū)地表二維形變[27]，即：

(6)

式中，θ1和θ2分別為升軌和降軌衛(wèi)星入射角；α1和α2分別為升軌和降軌衛(wèi)星方位角。

3 結(jié)果與分析

基于2.1節(jié)的SBAS-InSAR技術(shù)原理，本文首先獲取了Sentinel-1A數(shù)據(jù)在升軌和降軌LOS向形變速率圖，如圖4所示(其中圖4(a)為升軌形變速率，圖4(b)為降軌形變速率，以抬升為正、下沉為負(fù)方向)?？梢钥闯觯h(yuǎn)離礦體的區(qū)域大部分趨于穩(wěn)定狀態(tài)，無明顯形變特征。由圖4(a)可知，升軌監(jiān)測結(jié)果的形變區(qū)域主要集中在礦體西南部區(qū)域。而在圖4(b)的降軌監(jiān)測結(jié)果中，形變區(qū)域主要集中在南部以及西南部區(qū)域，礦體西部有局部抬升趨勢。造成這種差異則是受軌道影響，同時(shí)在礦體西部區(qū)域有一個(gè)明顯的山坡，從而導(dǎo)致了這種差異。整體來看，研究區(qū)在升軌結(jié)果中最大年平均形變速率達(dá)到-126.37 mm/y，降軌最大年平均形變速率為-84.03 mm/y，整體呈下沉趨勢，且沉降中心匯集于礦體西南部。

圖4 升降軌SBAS-InSAR形變速率圖

然而，不同軌道的監(jiān)測結(jié)果具有一定差異，不能反映地表的真實(shí)形變情況。因此本文采用2.2節(jié)的方法，解算了研究區(qū)地表二維形變場，如圖5所示(其中圖5(a)為東西向形變速率，圖5(b)為垂直向形變速率，以向東移動(dòng)為正、向西移動(dòng)為負(fù)方向)。

由圖5(a)可知，研究區(qū)地表西部區(qū)域向東移動(dòng)，最大年平均形變速率達(dá)到103.73 mm/y；東部區(qū)域向西移動(dòng)，且最大年平均形變速率為-32.40 mm/y，最后整體匯集于落礦區(qū)域并趨于穩(wěn)定。由圖5(b)可知，研究區(qū)整體呈下沉趨勢，最大年平均下沉速率達(dá)到-135.88 mm/y，最小為-2.47 mm/y，下沉區(qū)域集中于礦體西南部。

圖5 礦區(qū)地表東西向和垂直向形變速率圖

根據(jù)形變結(jié)果來看，礦區(qū)形變區(qū)域整體位于落礦區(qū)域上方及南部區(qū)域，不是落礦區(qū)域的正上方。經(jīng)實(shí)地考察與分析，總結(jié)了造成這種現(xiàn)象的原因，主要有以下兩點(diǎn)：

(1)該區(qū)域西南部存在F8斷層[6]，且地形呈西南高、東北低的特點(diǎn)，因此導(dǎo)致研究區(qū)西部和西南部的形變速率更大；

(2)地下開采影響反應(yīng)至地表需要一定的過程，具有滯后性，因此形變區(qū)域與落礦區(qū)域具有一定偏差?？傮w上來說，研究區(qū)的地表二維形變結(jié)果符合開采沉陷規(guī)律。

4 結(jié) 論

本文利用Sentinel-1A的升降軌SAR影像，采用SBAS-InSAR技術(shù)首先得到了甘肅省金昌市金川西二礦地表LOS向年平均形變速率。通過對升降軌結(jié)果差異的分析，結(jié)合InSAR幾何分解原理獲取了研究區(qū)地表二維形變場。結(jié)果顯示，該區(qū)域在2020.08.19～2021.06.27期間，地表東西向年平均形變速率達(dá)到103.73 mm/y，垂直向年平均形變速率達(dá)到-135.88 mm/y，且在地表形成塌陷坑。此外，通過考察與分析，形變結(jié)果基本符合開采沉陷規(guī)律，驗(yàn)證了方法的可靠性。