基于SBAS-InSAR的榆神礦區(qū)地表動態(tài)沉陷特征研究

李小濤，張孝榮，黃景偲，湯伏全

（1.陜西一八五煤田地質(zhì)有限公司，陜西 榆林 719000； 2.西安科技大學測繪科學與技術學院，陜西 西安 710000）

D-InSAR作為一種新型的地面觀測技術，已逐漸成為雷達遙感技術的熱點研究方向[1]，通過InSAR技術可以獲得cm甚至mm的高精度表面變形信息[2-3]，因此該技術在地震、冰川漂移、火山噴發(fā)、城區(qū)地面塌陷、山體滑坡、泥石流等災害監(jiān)測中發(fā)揮著重要作 用[4]，但在礦區(qū)大范圍開采沉陷監(jiān)測中應用相對較少。本文選取2018-01-06～2019-01-13的31期Sentinel-1A雷達影像，基于小基線集（SBAS-InSAR）技術對榆神礦區(qū)地表沉降進行動態(tài)監(jiān)測，基于變形序列結果分析地下開采工作面與地表變形的量化關系，實現(xiàn)了榆神礦區(qū)大范圍的遙感沉降監(jiān)測。

1 關鍵技術

1.1 SBAS-InSAR原理

小基線集（SBAS-InSAR）技術由Berardino于2002年最先提出[5]，該方法的基本思想是為了有效抑制時空失相干及大氣延遲帶來的誤差，將獲取的 SAR 影像按照時間基線及空間基線閾值分成若干小基線集組合，對小基線組合內(nèi)的干涉圖進行差分干涉處理，并利用奇異值分解法（SVD）計算出研究區(qū)的時序形變量。

設有按時間序列t0,t1,…,tn獲取N+1幅單視復數(shù)影像，首先將他每年以任意影像為主影像進行配準，然后設定垂直基線閾值，將垂直基線小于該閾值的SAR影像歸為一組，共分為L組。對每組內(nèi)的影像進行差分干涉處理，最終L組影像共得到M幅差分干涉圖，假設N為奇數(shù)，則差分干涉圖的個數(shù)M可以表示為：

以t0時刻為初始時刻，對于任意時刻ti(i=1,…,N)相對于初始時刻的差分相位φ(ti)為未知參數(shù)，觀測量為數(shù)據(jù)處理獲取的差分干涉相位δφ(tk)(k=1,…,m)如果所有的差分干涉圖都正確解纏，而且差分干涉相位被校正到某個穩(wěn)定或形變信息已知的高相干像元(x,r)上，該像元則為參考像元。參考像元解纏后的差分干涉相位為[6]：

式中，φ(tb,x,r)和φ(ta,x,r)分別為d(tb,x,r)和d(ta,x,r)對應的解纏相位值；φn為隨機噪聲；λ為波長；d(tb,x,r)和d(ta,x,r)分別為tb和ta時間內(nèi)得累計LOS向形變量；v(x)為像元x處的形變速率；Ti為主輔圖像的時間間隔。

對任意干涉組合均能表示上式的線性方程組，矩陣表示如下式：

式中，B為系數(shù)矩陣，每行對應一幅干涉圖，每列對應一個時間上的SAR圖像，影像分別為1和-1，其余列為0。如果m≥n，且B的秩是N，為滿秩矩陣，則最小二乘解得形變速率v的估計值：

為了保證干涉圖的高相干性，通常選取時間基線和空間基線都比較小的干涉組合，使得原有的SAR影像被分成L個獨立的干涉組合，由于系數(shù)矩陣B的秩 L＜N，導致上式秩虧，方程組的解不唯一，因此需要對B進行奇異值分解。

式中，U是m×m階酉矩陣；∑是半正定m×n階對角矩陣，其對角線上的元素∑i為B的奇異值；而VT，即V的共軛轉(zhuǎn)置，是n×n階酉矩陣。因此形變速率的估值v為：

將形變速率積分即可獲得每個SAR影像的累計地表形變相位。

1.2 下沉量提取

由于InSAR技術只能監(jiān)測雷達衛(wèi)星視線向的一維位移量[7-8]，如要進一步分析井下工作面開采與地面下沉范圍之間的關系，需要獲取地面真實的三維位移量。以往許多學者針對確定礦區(qū)下沉移動盆地邊界的問題時，一般采用忽略水平位移的疊加影響，直接將LOS形變量投影到下沉方向作為下沉量的方法。

本文基于開采沉陷的對稱原理（如圖1所示），相對于各工作面對應地表下沉盆地中心對稱的兩點P1和P3之間下沉量相等，水平位移大小相等方向相反的特征，對P1和P3點進行相加就能抵消掉水平位移，從而分別提取出該點的真實下沉量，如公式（8）[9]。

圖1 工作面地表對稱點的水平位移和下沉示意圖

本文在獲取地表下沉量時，均根據(jù)公式（8）求出采礦區(qū)任意地表點的下沉量。該下沉量剔除了水平位移的影響，能夠獲取真實的三維形變量。

2 研究區(qū)監(jiān)測實驗

2.1 研究區(qū)概況

榆神礦區(qū)跨陜西省榆林市榆陽區(qū)和神木市，是國家14個大型煤炭基地中陜北基地的主力礦區(qū)之一，總面積約5 265 km2，已探明儲量733.5億噸。礦區(qū)位于鄂爾多斯高原東南部，毛烏素沙漠的東南緣與陜北黃土高原北部的接壤地帶，總地勢西北高、東南低。榆神礦區(qū)西北部以陜蒙邊界為界，西南部以榆神礦區(qū)北邊界為界，東北部以神府礦區(qū)西邊界為界。礦區(qū)地理坐標為109 08′24″E～110 27′59″E，38 19′33″N～39 11′23″N。

2.2 實驗數(shù)據(jù)及數(shù)據(jù)處理

本文的實驗數(shù)據(jù)選取2018-01-06～2019-01-01的31期Sentinel-1A雷達影像，其工作模式為條帶模式，幅寬為250 km，距離向分辨率為5 m，方位向分辨率為20 m，時間分辨率約為12 d[10-11]。較高的時間、空間分辨率保證了影像對的相干性和變形監(jiān)測精度。

DEM數(shù)據(jù)選用美國地質(zhì)調(diào)查局網(wǎng)站公布的SRTM DEM數(shù)據(jù)，并且選擇覆蓋范圍為（38 N～39 N，109 E～110 E）的DEM。

本實驗基于小基線原理，選取時間基線不超過48 d，垂直基線不超過200 m，30幅影像得到了78個干涉良好的干涉像對。對選取的78個干涉像對分別進行差分干涉，去除地形相位，濾波處理；再利用干涉處理得到的干涉圖，依據(jù)相干系數(shù)法，對高相干點進行相位解纏；然后利用奇異值分解（SVD）算法， 計算出高相干點的線性形變速率及高程殘差，通過時空濾波及形變速率的積分累積運算，得到時序累積形變結果。

3 金雞灘礦地表動態(tài)沉陷與地下開采關系分析

本論文研究基于GAMMA軟件進行SBAS技術處理，繪制金雞灘106工作面從2018-01-06～2019-01-01期間的地表時序形變圖2，該圖中的黑色虛線框為井下工作面開采面。

從圖2可以看出2018-03-31之前地表并未出現(xiàn)地表形變，到2018-05-06，地面開始出現(xiàn)了形變，且往后形變范圍持續(xù)向東北方向延長，呈長條狀。根據(jù)動態(tài)變形情況可推斷106工作面的開采影響時間是2018-04，經(jīng)實地調(diào)查礦區(qū)的開采資料，該工作面確實于2018年四月中旬開始開采，與監(jiān)測到的地表形變情況相符。

圖2 金雞灘106工作面地表的時序形變圖

為了揭示監(jiān)測到的地表動態(tài)沉降特征與煤礦地下開采工作面推進情況的關系，根據(jù)實際采礦資料，該變形區(qū)不同時段對應的井下工作面推進距離增量與變形區(qū)域前沿位置距工作面開切眼的水平距離之間的相對關系如表1所示。

表1 變形區(qū)前沿擴展與工作面實際推進距離的關系

從表1可見，地表變形區(qū)域的前沿位置擴展的距離與相應時間段地下工作面的推進距離之比值始終在1∶1左右，其動態(tài)特征符合開采沉陷的基本規(guī)律，也與該觀測站的實測數(shù)據(jù)一致。說明隨著地下工作面的不斷推進，地表變形區(qū)向前擴展的速度與地下工作面推進的速度基本相同，變形區(qū)域隨著地下工作面的推進而同步發(fā)展。

4 榆神礦區(qū)地表動態(tài)沉陷監(jiān)測結果分析

基于以上對金雞灘106工作面地表形變與地下工作面開采關系的研究基礎，可以將SBAS-InSAR技術和本文的下沉量提取方法應用到大范圍變形區(qū)域的時序形變分析。首先可利用干涉圖堆疊技術搜索出整個榆神礦區(qū)的形變區(qū)域分布，如圖3 為整個榆神礦區(qū)在2018年3月到2019年1月間的形變累積速率圖。然后根據(jù)就搜索出的形變區(qū)分塊分別進行時序InSAR數(shù)據(jù)處理，并提取出地表下沉量。

圖3 榆神礦區(qū)2018年InSAR監(jiān)測變形速率分布

從累積形變結果圖中總共監(jiān)測到20個變形速率大于0.1m/year的變形區(qū)域，分別位于杭來灣、榆樹灣、曹家灘、小保當、柳巷、千樹塔等煤礦開采區(qū)，其中 9號形變區(qū)為前文的金雞灘礦區(qū)的106工作面。

為了統(tǒng)計出現(xiàn)的變形區(qū)范圍的動態(tài)變化，利用SBAS-InSAR技術繼續(xù)處理31期SAR影像，將榆神礦區(qū)地表變形區(qū)域分成9個區(qū)塊，分別統(tǒng)計9個變形區(qū)在2018年度不同時段累積的變形區(qū)面積統(tǒng)計如表2所示。

表2 各變形區(qū)累計變形面積統(tǒng)計/km2

從表2可以看出，位于金雞灘礦區(qū)的2號、9號形變區(qū)，杭來灣的3號形變區(qū)、榆樹灣的5號6號形變區(qū)的形變范圍均發(fā)生了明顯的變化，結合圖3中各個形變區(qū)的形態(tài)及表2中形變區(qū)面積變化，可進一步推斷地下工作面的開采動態(tài)情況。

其中，1號和2號形變區(qū)位于金雞灘煤礦區(qū)域，根據(jù)其變形區(qū)的西南端窄，東北端寬的形態(tài)可推斷該地區(qū)開采前進方向是東偏北方向，并且從2018年1月至2019年1月的形變范圍基本無變化，推斷該地下工作面于2018年初停止開采。2號變形區(qū)呈長條狀，2018年4月到10月間形變范圍逐漸增大，根據(jù)其西南端寬而東北端窄，且10月9日到11月2日間的變形區(qū)域基本保持不變的表現(xiàn)特征，可以推斷該地下工作面沿西偏南方向進行，并且于10月份的地下開采工作出現(xiàn)了暫停。結合采礦資料與上述對地下工作面開采情況的推斷完全符合。9號在該監(jiān)測時間范圍內(nèi)其形變面積逐漸增大。

3號形變區(qū)位于杭來灣煤礦，該形變區(qū)呈相鄰兩個長條狀的形態(tài)，分成左上角和右下角兩塊，從表2中的數(shù)據(jù)顯示3號變形區(qū)面積也在逐漸變大，并且從形變圖3上可以看出該形變區(qū)中心出現(xiàn)了明顯的失相干現(xiàn)象，說明在該時間段內(nèi)地下正在處于開采工作狀態(tài)，并且該地下工作面的開采造成了其對應地表出現(xiàn)了大梯度變形。第4變形區(qū)域呈圓狀，且一年內(nèi)形變范圍沒有變化，可推斷該區(qū)域地下處于停采狀態(tài)且存在殘余變形。

5、6號變形區(qū)域位于榆樹灣煤礦，其中5號變形區(qū)域面積基本沒有變化，推斷在2018年初地下開采工作逐漸停止，地表處于殘余變形。6號變形區(qū)在2018年5月前形變面積逐漸增大，但在5月以后變形區(qū)范圍不再變化，推斷該區(qū)地下工作面開采在5月后停止，地表屬于殘余變形。經(jīng)實地調(diào)查該礦6工作面在5月停采，與監(jiān)測到的形變情況相符。

7、8號變形區(qū)域位于柳巷礦區(qū)，其中7號變形區(qū)面積和位移量基本沒有變化，可推斷該時間段內(nèi)地下開采工作已停止，地表處于殘余變形階段。8號區(qū)域于2018年5月期間地表出現(xiàn)變形，呈長條形狀，7月以后變形區(qū)范圍沒有隨時間擴展，推斷該區(qū)域于5月期間開始地下開采工作，然而在七月后地下停止開采，屬于小面積的地下開采。

在沉陷盆地中央下沉梯度較大的情況下，由于SAR影像的波長有限，導致時序累積D-InSAR監(jiān)測結果不能真實地反映沉陷區(qū)中央的絕對變形量。但是利用其微小形變結果能夠較為準確地反演地下的開采時間以及工作面的開采位置及其范圍情況。

5 結 語

1）基于時序D-InSAR的煤礦區(qū)地表沉降監(jiān)測相比于傳統(tǒng)測量而言，具有全天時、全天候、覆蓋廣、監(jiān)測結果更直觀等諸多優(yōu)勢?？梢詮臅r間和空間上反映出井下工作面開采后地表的變形范圍和動態(tài)發(fā)展趨勢，能夠?qū)崿F(xiàn)礦區(qū)大范圍、長時序的地表變形區(qū)域探測，該技術在西部煤礦大規(guī)模開采地表沉陷遙感監(jiān)測中具有普適性。

2）D-InSAR技術對于煤礦地表沉陷區(qū)邊緣的小梯度變形監(jiān)測非常有效，通過地表變形信息可揭示動態(tài)變形邊界與井下工作面推進邊界之間的量化關系。實驗表明，2018年1月至2019年1月榆神礦區(qū)地表動態(tài)變形區(qū)域的范圍和分布特征與地下工作面開采推進范圍和方向呈顯著的線性相關。