基于時序InSAR技術的新井煤礦地表沉降監測

馬 飛，李日印，劉小鵬，侯景怡，崔喆森，張垚杰

（1.長治學院計算機系，山西長治 046011；2.自然資源部第二地形測量隊，陜西西安 710054）

大規模開發利用礦產資源會直接影響礦區安全[1]。例如：2023 年2 月22 日，內蒙古新井煤業有限公司露天煤礦發生大面積邊坡失穩事故，造成人員財產損失。因此，對礦區地表變形進行監測具有重要的現實意義[2-4]。

差分干涉合成孔徑雷達(DInSAR)、小基線集In-SAR(SBAS-InSAR)和永久散射體InSAR(PS-InSAR)技術都是監測采礦沉陷、凍土沉陷、建筑物沉陷和滑坡的常用方法[5-7]，具有空間覆蓋面大、成本低、精度高、全天候觀測能力等優點，而常規測量技術（如水準測量[8]、三維激光掃描[9]和全球導航衛星系統[10]等）無法做到這一點。為了滿足開采沉陷長期監測的需要，削弱大氣相位的影響，時序InSAR 提出了一種簡單有效的方法，即通過對一系列展開的干涉圖按時間跨度加權和平均來估計平均視線開采沉陷速度[11-12]。小基線集InSAR (SBAS-InSAR) 可以很好地估計變形結果[13-15]，在許多開采沉陷區得到了應用。現對內蒙古新井煤礦應用SBAS-InSAR技術，識別該礦區地表沉陷范圍，對沉陷區下沉數值進行解算，分析災前微小形變規律，以預防同類事故的再次發生。

1 SBAS-InSAR技術原理

SBAS-InSAR 通過設置特定的時間基線閾值來合并小的基線干擾對，根據最小參數準則執行最優的形變相位序列，以限制時間和空間消相干。假設已經采集了在時間序列(t0，t1，…，tn)相同區域中獲得的（n+1）個SAR 圖像。設置特定的時間基線閾值后，可以形成m個干涉對，且m滿足：

SBAS-InSAR 在進行參數反演之前，需要將所有圖像與超級主影像配準，并通過多視處理、地形差分和相位解纏來生成m幅差分干涉圖圖像。考慮在時刻ta和t（b假設ta

式中：φ(tb,x)和φ(ta,x)分別為相位形變值；λ為雷達波長；φn為隨機噪聲；d(tb,x)和d(ta,x)分別為衛星視線方向的形變量；v(x)為像元的形變速率；Ti為干涉對的時間基線長度。

這樣可以組成m個方程組：

式中：A為m×n的系數矩陣，每行對應一幅干涉圖，每列對應一個時間上的SAR 圖像，主影像所在列+1，輔圖像所在列-1，其余列為0。

如果m≥n，且A的秩是n，則最小二乘解得φ的估計值：

由式（4）可知，SBAS 技術的關鍵在于求時間序列干涉圖間形變速率的最小范數最小二乘解，可以利用時序干涉圖中解算出的形變相位值，作為時間域的約束條件解，簡化形變速率解算的復雜度，達到反演形變速率的目的。

假設V為形變速率，P為形變參數，以下約束條件被滿足：

將式（5）代入BV=δφ得：

在干涉圖中任意像元（x,r）的解算模型轉換為：

式中：LOS方向形變的平均加速度變化率?、相位平均速率、平均加速度用下式即可解得：

2 研究區概況及實驗結果

2.1 研究區概況

圖1為內蒙古新井煤業露天礦區地理位置，該礦區位于賀蘭山南段，礦區面積1.344 8 km2，地勢北高南低，海拔高程1 159.03～1 490.70 m。礦區為低山丘陵區，屬高原侵蝕性丘陵地貌，地形切割一般，基巖裸露，植被稀疏，為荒漠地區。氣候屬大陸性干旱荒漠氣候，日溫差大，夏季炎熱，冬季嚴寒，氣候干燥，雨量稀少，年最高氣溫38℃，最低-28℃，平均8℃～9℃，年降水量為146～ 198 mm，年蒸發量最高2 000 mm以上。地表植被屬于荒漠化草原向草原化荒漠過渡地帶。地表多沙質化、礫石化和有龜裂結皮。

圖1 新井煤業露天煤礦地理位置

通過分析新井煤業露天礦區地表植被和地形地貌特征，該礦區地表屬于荒漠地區，植被稀少，適用C波段數據，但地面建筑物、工礦構筑物、大型巖石等永久散射體物較少，適合采用SBAS-InSAR對礦區地表沉降情況進行監測。

2.2 SAR數據選取及預處理

表1為覆蓋該礦區的Sentinel-1A數據列表，包括衛星拍攝時間、影像模式、極化方式和軌道方向等信息。根據歐空局(ESA)網站查詢哨兵-1A 數據情況發現，該衛星2022年能夠覆蓋研究區的影像只有15景，觀測時間是2022年1月到2022年11月。根據干涉處理實驗需要，下載了該區域干涉寬幅模式數據，單幅影像幅寬為250 km，距離向分辨率為4 m，方位向分辨率為13 m，多視比例設置為3∶1。同時下載到所有影像的精密軌道數據，用于干涉數據處理過程中的圖像配準、去平地相位等流程，提高數據處理精度。另外，為提高數據處理效率，對原始SAR影像進行了裁切。

表1 覆蓋礦區的Sentinel-1A 數據列表

2.3 時序InSAR 技術數據處理流程

2.3.1 短基線集干涉圖組合

圖2 為干涉對時空基線圖。將選取的15 景影像按照如下策略進行自由組合，形成干涉對：時間基線ΔT≤100 d，空間基線ΔB≤200 m，以20220401 為超級主影像。最終得到43個干涉對。

圖2 SAR干涉對時空基線圖

2.3.2 干涉圖生成

圖3 是在新井煤業露天礦區生成的以20220401為主影像的差分干涉圖，該圖像為SAR 斜距坐標系統。圖中可以看到，礦區中心位置形變條紋非常清晰，隨著時間的推移，礦區中心位置形變量在不斷變化；時間基線越長，干涉條紋周期變化越多；在圖(i)、(j)等干涉對的中心位置出現干涉條紋疊加和失相干。

圖3 以20220401為主影像生成的的干涉圖

2.3.3 SBAS時間序列形變和形變速率結果

圖4 為新井礦業露天礦區形變速率圖，紅色代表形變方向與衛星視線方向一致，可解算為下沉值，藍色代表形變方向與衛星視線方向相反，可解算為抬升值。圖中可見，該礦區大部分區域相對穩定，只有在此次滑坡事故中的核心區發生劇烈形變，整體形變速率非常大，沉降中心位置最大形變速率達到-160～ -209 mm/a，說明此處滑坡2022 年開始已經出現不穩定狀態。另外，在該區域的東南側也出現一個小范圍的形變區，該區域的形變速率在-60～ -110 mm/a，需要引起相關部門的高度關注，避免發生新的地質災害。

圖4 災變核心區形變速率圖

為進一步分析新井煤業露天礦區核心區災前形變情況，將該形變區域數值提取出來并進行處理。該區域共解算得到2 475 個形變點，圖4 顯示該區域自邊坡的坡頂向坡底發生形變量由大到小，2022年1月至11 月坡頂累計形變量超過170 mm，坡底累計形變量超過150 mm。

圖5 為邊坡在2022 年的形變速率剖面圖，橫坐標為距離邊坡坡頂GCP1 號點的距離，縱坐標為形變速率值。為進一步分析該礦區邊坡失穩前的形變規律，在邊坡選取了四個點命名為GCP1～4進行分析，分別分析其剖面線和形變時序圖。圖中可知，邊坡頂部下沉速率遠遠高于坡底位置，自坡頂向坡底形變速率逐漸下降，形變速率從-180～-140 mm/a。由于該區域邊坡角度較大，坡頂區松散層不斷下沉，經過一段時間的累積，需要采取任何進行邊坡加固。

圖5 形變速率剖面圖

圖6 為新井煤業露天礦區災變核心區內邊坡形變點的時序下沉變化情況。由圖可知，在邊坡頂部的GCP1 號點累積形變量最大，在底部的GCP4 號點累積形變量稍小，這與圖5 中的剖面圖一致；該邊坡從2022 年2 月開始出現下沉，2022 年3 月至6 月累積形變趨勢開始增大，這可能與該地區前期轉暖凍融解凍有關，也可能與礦區露天開采活動有關。

3 結論

(1)利用SBAS-InSAR 方法可以有效識別出露天采礦區地表變形，這是一種快速、高效、低廉的技術，可以為礦區地表邊坡失穩、地面塌陷、地面沉降、地裂縫、滑坡、崩塌等地質災害提供預警。

(2)內蒙古新井煤業露天采礦區核心災變區在2022 年1-11 月持續發生下沉形變。分析SBAS-In-SAR 方法獲取的2 475 個形變點發現該區域最大形變量超過160 mm，沉降速率超過200 mm/a；在核心災變區東南側亦有一個緩慢形變區，累積形變量稍小，但需要有關部門關注，避免發生新的地質災害。

(3)內蒙古新井煤業露天采礦區失穩邊坡在災前一年間持續發生下沉，且邊坡坡頂形變速率高于坡底速率，坡頂形變速率超過180 mm/a，頂部持續微小的下沉形變是導致邊坡應力失衡的重要原因。

(4)由于缺少地表監測點，故SBAS-InSAR 的形變監測值無法得到驗證，但該技術獲得的時間序列累積變形能更直觀地反映沉陷盆地的時間序列變化趨勢，在探測變形區域范圍和形態方面具有優勢，有助于定性確定形變類型。

致謝：本實驗的Sentinel-1A 數據和SRTM 外部DEM數據均來自歐空局網站，在此表示感謝。