基于SBAS-InSAR的榆神礦區地表動態沉陷特征研究

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

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

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

1 關鍵技術

1.1 SBAS-InSAR原理

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

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

以t0時刻為初始時刻，對于任意時刻ti(i=1,…,N)相對于初始時刻的差分相位φ(ti)為未知參數，觀測量為數據處理獲取的差分干涉相位δφ(tk)(k=1,…,m)如果所有的差分干涉圖都正確解纏，而且差分干涉相位被校正到某個穩定或形變信息已知的高相干像元(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時間內得累計LOS向形變量；v(x)為像元x處的形變速率；Ti為主輔圖像的時間間隔。

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

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

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

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

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

1.2 下沉量提取

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

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

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

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

2 研究區監測實驗

2.1 研究區概況

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

2.2 實驗數據及數據處理

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

DEM數據選用美國地質調查局網站公布的SRTM DEM數據，并且選擇覆蓋范圍為（38 N～39 N，109 E～110 E）的DEM。

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

3 金雞灘礦地表動態沉陷與地下開采關系分析

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

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

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

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

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

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

4 榆神礦區地表動態沉陷監測結果分析

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

圖3 榆神礦區2018年InSAR監測變形速率分布

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

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

表2 各變形區累計變形面積統計/km2

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

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

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

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

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

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

5 結 語

1）基于時序D-InSAR的煤礦區地表沉降監測相比于傳統測量而言，具有全天時、全天候、覆蓋廣、監測結果更直觀等諸多優勢。可以從時間和空間上反映出井下工作面開采后地表的變形范圍和動態發展趨勢，能夠實現礦區大范圍、長時序的地表變形區域探測，該技術在西部煤礦大規模開采地表沉陷遙感監測中具有普適性。

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