基于SBAS-InSAR 的采煤區地表形變監測與分析

宋旭斌，路 鑫

（潞安化工集團有限公司 古城煤礦，山西 長治 046000）

0 引 言

煤炭作為我國的主體能源，長期擔負著國家能源安全和經濟高速發展的重任。但大規模、高強度的地下開采往往會破壞地層的原始應力平衡，而且隨著采空區范圍擴大，直接頂在自身重力和上覆巖層荷載作用下破斷，最終引起地表的移動變形[1]。因開采導致的地表沉陷易造成建筑物墻體開裂、交通道路損壞、農田積水等地質災害問題，嚴重制約了礦區城市的可持續發展。因此，對采空區地表進行實時精密形變監測，獲取動態沉降信息勢在必行。

常規開采監測手段主要有精密水準測量和全球導航系統（Global Navigation Satellite System,GNSS），但在實際工程中往往面臨觀測點易丟失、觀測精度易受干擾、成本投入巨大、時空分辨率低等問題，無法準確獲取礦區整體的形變特征。合成孔徑雷達干涉測量(Interferometry Synthetic Aperture Radar,InSAR)是一種新型的對地觀測技術，它僅對同一區域的多景SAR 影像進行復共軛相乘即可提取大范圍、毫米級的地面形變信息[2]，具有全天時、高精度、自動化等優勢。而小基線集（Small Baseline Subsets InSAR,SBAS-InSAR）技術在差分InSAR（Differential InSAR,D-InSAR）基礎上提出，有效克服了后者易受大氣效應、時空失相干影響的缺陷[3]，它通過對SAR 影像上的高相干點進行相位模型建立和相位變化分析，即使在永久散射體PS（Persistent Scatterer）密度較低的礦區也能獲取更為精確的長時序形變結果。

許多學者都曾將SBAS 技術運用于礦區的形變監測，并取得了一系列豐碩的成果。李達等[4]利用SBAS 技術對13 景TerraSAR-X 進行處理，獲取了榆林某礦的地表沉降信息，發現地表下沉值與時間呈線性關系；欒元重等[5]借助19 景Sentinel-1A 影像和SBAS 技術對沉陷漏斗區域內的工作面進行時序監測，并利用水準數據驗證了SBAS 技術的可靠性；柴華彬等[6]提出利用反距離加權插值方法（Inverse Distance Weighted,IDW）將實測水準數據與SBAS 解譯結果相融合，以解決部分區域因失相干導致InSAR 形變值缺失的問題；張香凝等[7]以蒲河礦區為例，結合SBAS 技術和數值模擬共同探究了煤礦區地面沉降的時空變形規律和機制。

山西長治地區煤炭資源豐富，但因地下開采引起的地質災害頻發，亟需開展相關的地表沉陷監測工作。本文在前人研究的基礎上，借助SBAS-InSAR 技術對29 景Sentinel-1A 影像進行解譯處理，以期獲取古城煤礦S1306 工作面在2021年10 月12 日至2022 年10 月7 日期間的時序沉降結果和開采影響范圍。研究成果對礦區開采規劃和地表塌陷的預警防治均具有重要的現實意義。

1 研究區與數據源

1.1 研究區概況

本文研究區位于長治市古城煤礦東南方向一采區，整體地勢呈西高東低，地貌類型以河谷平原為主。研究區分布著大面積農田、村莊房屋以及電力線等重要的基礎設施，且覆蓋有較厚的第四系黃土沉積；地下山西組3 號煤層目前有多個工作面正在回采，其中S1306 工作面走向長約1 982 m，傾向寬約363 m，平均采厚6.4 m，煤層傾角約為0～6°，底板標高-352—-397 m。開采方法以傾斜長壁式綜采法為主，全部垮落法管理頂板。目前，該采區周邊已出現大量明顯的地裂縫及路面損壞情況，亟需對該區域進行地表沉陷監測工作。研究區的地理位置如圖1 所示。

圖1 研究區地理位置Fig.1 Geographical location of the study area

1.2 數據源

本文從歐空局（European Space Agency,ESA）哥白尼開放數據訪問中心（https://scihub.copernicus.eu/dhus/#/home）收集了29 景Sentinel-1A 升軌單視復數影像，時間跨度為2021 年10 月12 日至2022 年10 月7 日，工作模式均為干涉寬幅（Interferometric Wide-swath，IW）模式，該模式基于遞進式地形掃描(Terrain Observation with Progressive Scans SAR，TOPSAR)技術，能夠將3 個子條帶合并成1 景高質量的SAR 雷達數據[11]，其地面分辨率可達5 m×20 m。實驗影像的其他參數見表1。另外，在數據處理過程中，引入了POD精密軌道文件用來去除衛星軌道的系統性誤差，并利用30 m 分辨率的STRM1 DEM 數據消除了干涉圖中地形相位的影響。

表1 Sentinel-1A 基本參數Table 1 Basic parameters of Sentinel-1A

2 技術原理與數據處理流程

2.1 SBAS-InSAR 原理

小基線集技術最早由國外學者Berardino 等人在2002 年提出，其基本原理如下[12-15]：假設選取（t0,t1,t2,…,tN）時段內覆蓋同一區域的N+1 幅SAR影像，基于設定的時空基線閾值，組合成M 幅多視差分干涉對，其中M 的取值范圍滿足：

若在tA和tB（tB>tA）時刻對2 幅SAR 影像進行干涉處理，并生成了第j 幅已解纏的差分干涉圖，在去除地形相位、大氣延遲相位和噪聲誤差之后，任意像元位置（x,r）（x 和r 分別為方位向和距離向坐標）處的差分干涉相位可表示為：

式中：λ 為雷達波長；φ(tB,x,r)和φ(tA,x,r)分別表示2 幅影像上的相位值；d(tB,x,r)和d(tA,x,r)分別為tA和tB時刻相對于初始時刻t0在雷達視線向(Line of Sight,LOS)上的形變，一般假設t0時刻的形變相位值為0；vj為tA至tB時段內的平均形變速率。

將用于干涉處理的M 對影像數據按照時序進行排列，則所有干涉圖相位可表示為如下矩陣形式：

式中：B 為1 個M×N 的系數矩陣；V 為速度矢量。

由于矩陣B 的秩小于N，根據最小二乘法無法得到唯一的結果，因此需要采用奇異值分解法（Singular Value Decomposition,SVD）聯合多個基線集，來解決法方程的秩虧問題，并得到速度V 的廣義逆解。之后將各個時間段內的形變速率在時間域上進行積分，即可獲取整個監測時段內地表在LOS 向上的累積時序形變結果。

2.2 SBAS-InSAR 數據處理流程

SBAS-InSAR 數據處理關鍵步驟和參數設置如下：①選取2021 年11 月29 日的影像作為主影像，其余從影像與之配準，并設置空間基線百分比閾值為45%，時間基線閾值為60 d，多視比為4∶1（距離向∶方位向），最后共生成104 個干涉對；②利用最小費用流(Minimum Cost Flow,MCF)方法進行相位解纏，Goldstein 濾波技術去除干涉圖噪聲，并逐個檢查調整不理想的像對；③參考干涉處理步驟生成的相干圖和解纏相位圖，選取20～30 個穩定且分布均勻的地面控制點（Ground Control Point,GCP）完成軌道精煉，去除殘余軌道誤差相位；④選擇線性模型（Linear Model）進行二次解纏，并估算平均位移速率和殘余地形相位；⑤根據大氣延遲相位、形變相位和噪聲相位表現出的不同的時空特性，利用時間域的高通濾波（365 d）和空間域的低通濾波（1 200 m）分離并剔除大氣延遲誤差，得到覆蓋整個觀測時間的地表形變時間序列；⑥將最終解譯結果從斜距坐標系轉換至WGS-84 坐標系，即完成地理編碼。本文的數據處理流程如圖2所示。

圖2 SBAS 數據處理流程Fig.2 Data dealing process of SBAS

3 解譯結果與分析

3.1 空間演化分析

利用SBAS-InSAR 技術對覆蓋S1306 工作面的29 景影像進行解譯處理，獲取了監測時段內研究區LOS 向的時序累積形變結果。之后忽略水平位移的影響，利用公式dv=dLOS/cosθ 將解譯結果從視線向投影至垂直向，其中θ 為衛星入射角。另外，解譯結果中的正值表示地面抬升，負值表示地面下沉。

S1306 工作面自2022 年1 月開始回采，目前僅開采到距離始采線約850 m 的位置。但圖3 顯示，在監測時段內，該工作面周邊地表已經發生了明顯沉陷。

圖3 S1306 工作面時序形變Fig.3 Time series deformation of No.S1306 Face

在形變監測的初期階段，即2021 年10 月12日至2021 年12 月23 日（圖3a～圖3c），研究區地表總體較為穩定，零星區域發生小幅沉降，量值約-13 mm，產生該現象的原因可歸結為工作面所處的農田區域相干性較差，導致SBAS 在反演形變的過程中發生解纏誤差，影響了結果的準確性；從2022 年1 月16 日起（圖3d），工作面的始采線位置開始出現明顯下沉，最大沉降量約-32.6 mm，符合工作面開始回采的時間；隨著工作面持續推進，地表沉降范圍逐漸朝西南方向（即工作面開采方向）擴張；到了2022 年5 月16 日（圖3h），始采線附近地表已經形成了一個較為明顯的沉降盆地，盆地中心可探測到的最大下沉值約-178.6 mm；隨著采動影響逐漸傳遞至工作面中部位置（圖3i～圖3k），始采線區域的下沉速度逐漸趨于平緩，但仍處在下沉狀態；直至2022 年10 月7 日（圖3m），此時S1306 工作面的開采影響邊界稍超過實際開采進度位置，但總體與實際情況相符，說明SBAS 技術在礦區形變監測中具有一定的可靠性和準確性。

3.2 時間演化分析

為了進一步探究研究區地表在采動影響下的時序形變過程，在S1306 工作面的走向和傾向上共選取8 個特征點（編號點A～點H），并以2021 年10 月12 日影像為參考基準，生成各個點的時序累積沉降曲線。具體點位分布如圖4 所示。

圖4 特征點位及其時序累積沉降值Fig.4 The location of feature points and their time-series cumulative subsidence values

由圖4 可知，在2021 年10 月12 日至2022 年1 月16 日期間，走向線上4 個特征點沒有發生明顯形變，由于點A 和點B 分別位于始采線和沉陷盆地中心，易受開采擾動，這一時段內兩點均出現小幅度下沉，最大沉降值約-10 mm；而點C 和點D 分布在工作面中部和停采線的位置，點位表現出反復微弱抬升的跡象，可將這種情況其視為SAR圖像噪聲帶來的解算誤差。

2022 年1 月16 日之后，點A 和點B 開始持續加速下沉，且分別于2022 年3 月17 日和2022 年5 月16 日，其形變速率才逐漸放緩，但形變量仍在持續增加，這是由于采動作用對始采線區域地表的影響即將達到“臨界點”，該區域的形變經歷了“開始下沉—劇烈下沉—緩慢停止下沉”的過程，一定程度上符合常規開采沉陷規律。點C 和點D在2022 年8 月8 日之前，下沉曲線吻合度較高，但之后點C 的沉降量陡增至-87.6 mm，說明工作面中部區域此時開始受采動的強烈影響，與實際工程進度一致。同時，發生迅速沉降的原因還可能是由于點C 處在土質松軟的農田，在夏季強降水的滲透下，土層的抗剪強度會大大降低，從而引起地表的濕化崩解和垮落[12]。

傾向線上的4 個特征點均勻分布在工作面的兩側，在2022 年1 月16 日之前，由于點E 和點H遠離開采中心，點F 和點G 位于沉陷盆地邊緣，因此后者呈現小幅持續下沉的態勢，前者沉降與起伏交替，受采動影響較小；之后點F 和點G 開始加速下沉，但點F 的沉降幅度與速率要大于點G。同樣，點E 和點H 于2022 年4 月10 日也開始經歷一個緩慢下沉的過程，但總體上點E 的沉降量級要高于點H。

因此結合圖3 的時序形變結果，可以發現沿開采方向右側的地表形變值要大于左側區域，即始采線附近的沉陷盆地實際上有朝西北方向發育的趨勢，后期需要對東李高村附近地表進行重點監測。

4 結 論

（1）在采動作用下，S1306 工作面周邊地表自2022 年1 月16 日開始出現下沉，之后在始采線附近迅速產生了一個沉陷盆地，且沉陷盆地有朝西北方向發育的趨勢。

（2）目前采動影響僅波及至工作面的中部區域，與實際開采到達位置基本吻合，說明SBAS-InSAR 技術在礦區形變監測中具有良好的適用性。

（3）除了利用時序InSAR 技術獲取不同時段的地表沉陷情況，仍可以考慮利用UAV 激光雷達掃描技術或建立地面移動觀測站來對局部重點區域進行詳細監測，真正發揮“空-天-地”協同監測技術在開采沉陷領域的重要作用。