礦山開采后地表沉降邊界的D-InSAR 監測與分析*

趙 恒，李亞星，楊可明*，竇 智

（1.淮北雙龍礦業有限責任公司地測科，安徽 淮北235000；2.中國礦業大學（北京）地球科學與測繪工程學院，北京100083）

地下煤層開采會對其上覆的巖層結構產生破壞，進而引起地面塌陷與沉降變形，煤礦區的采動嚴重影響威脅了采區地表的土地、水體、建（構）筑物等的安全。為了保障煤礦區生產和生活安全的同時盡可能多回收煤炭資源，需要科學地確定開采沉降范圍，以便明確資源開采位置和界定相關責任，降低煤礦企業的開采風險。

礦山開采引發地表沉降影響范圍的常規分析主要是依據開采工作面的實際地質采礦條件、地表巖移觀測數據以及巖移規律等資料，運用實測數據分析、巖層移動角幾何圖解、概率積分法（Probability Integral Method，PIM）預計和數值模擬計算等方法[1-4]，但這些常規方法存在較多且明顯的不足，如：巖移觀測站實測數據分析與幾何圖解主要體現在開采工作面的走向與傾向方向上，觀測點數量與圖解范圍相當有限；PIM 與數值模擬方法雖然能大范圍確定地表沉降邊界，但不能反映不同的開采工作面存在不同的上覆巖層巖移與物理參數而引發實際沉降邊界的不同變化，導致所確定的邊界與實際位置存在一定甚至較大的偏差。而微波遙感的合成孔徑雷達干涉測量（Interferometric Synthetic Aperture Radar，InSAR）技術并不需要礦山開采的地質采礦資料，而是基于側視雷達電磁波的地面反射回波相位干涉處理技術來獲取地表的沉降信息，常用技術有1989 年Gabriel 等首次提出的差分干涉測量（Differential InSAR，D-InSAR）[5]、2000 年Ferretti正式提出的永久散射體干涉測量（Permanent Scatterer InSAR，PS-InSAR）[6-7]和2002 年Berardino 提出的小基線集干涉測量（Small Baseline Subset InSAR，SBAS-InSAR）[8]等技術。InSAR 技術由于可以克服因地下開采條件復雜、地表移動變形，以及易受地下開采、覆巖巖性及結構、地下水等多種因素影響而難以準確界定地下開采所形成影響邊界的局限，所以現今已成為礦山開采形變監測[9-10]、地質災害分析與環境影響評價[11-13]等開采沉陷研究的前沿技術。

D-InSAR 是目前較為先進的一種地表形變監測方法，與常規方法相比在監測地表微小形變方面有著顯著的優勢，它不僅能以厘米量級甚至毫米量級監測地表形變，且能對大面積測區進行快速監測，從而得到采區主要形變區域與范圍[14-15]。本文將結合安徽省淮北礦區某煤礦井下311 工作面煤層開采下保護地面重要建（構）筑物運營安全等實際需求，利用D-InSAR 干涉測量技術監測煤層工作面采動影響的地表沉降邊界，以下沉10mm 的位置為確定邊界的依據，并基于現場實測數據、SBAS-InSAR 技術分析結果和概率積分法變形預測等對所監測到的地表沉降邊界進行驗證，進而能為煤礦安全開采和最大程度回收地下煤炭資源提供科學的理論依據與技術保障。

1 理論與方法

1.1 D-InSAR 技術原理

D-InSAR 技術以其大范圍、低成本、速度快、高精度等優勢，已被廣泛應用于各種地表形變監測中。如圖1 所示，D-InSAR 的基本思想從兩幅SAR 影像的干涉相位中分離出地表形變所引起的相位分量，從而對影像區域的地表沿雷達視線向地表形變進行制圖。干涉相位是兩幅配準過后的SAR 影像對應像元的相位差，其組成可以表示如下[16]：

圖1 平地-高程相位函數模型

式中，φint表示干涉相位；φtop、φflat、φorb、φatm和φnoise分別表示高程相位、平地相位、軌道誤差相位、大氣相位和相位噪聲；φdef表示形變相位。軌道誤差可以通過精密定軌和模型擬合去除，噪聲一般通過濾波來抑制。礦區沉降空間范圍較小，大氣相位常常可以忽略不計。高程相位和平地相位可通過函數模型來去除。

圖1 中S1、S2分別表示主輔影像成像時的天線位置。B 為空間基線，B||和B⊥分別表示平行和垂直于雷達入射方向的基線分量，稱為平行基線和垂直基線。α 為空間基線與水平方向的夾角，θ 為雷達入射角，Δθ 為高程引起的雷達入射角變化。R1、R2表示天線與目標點P 之間的斜距。R 表示P 點等斜距投影到參考橢球面的R1距離（R=R1）。h 為目標點P 的高程。根據圖中幾何關系和泰勒公式，有：

式中，λ 為波長。因此可以得出平地相位和高程相位的表達式：

形變相位φdef與地表視線向形變dlos之間的關系如下：

1.2 地表沉降邊界分析方法

D-InSAR 是在InSAR 的基礎上發展起來的監測地表微小形變的技術，在監測工礦區地表形變方面具有連續空間覆蓋，高度自動化和高精度的優勢。一般常用二軌法D-InSAR 處理流程，即采用外部的數字高程模型（Digital Elevation Model，DEM）數據模擬并去除地形相位。二軌法D-InSAR 的基本處理流程如圖2 所示。

圖2 二軌法D-InSAR 基本處理流程

1.3 概率積分法地表移動變形預計思想

地表移動與變形計算的概率積分法（PIM）是指把巖層移動看作服從統計規律的隨機過程，從而將開采引起的地表下沉剖面表示成概率密度函數積分公式的預計方法。已有大量理論推導及實例驗證表明，當地下開采工作面的頂板管理方法采用自然垮落時，概率積分法得到的地表沉降盆地穩定后的預計結果與實測數據具有較高的擬合優度和精度，是用于礦區開采沉陷預計的一種行之有效且廣泛使用的方法。

2 結果與分析

2.1 研究區工作面概況與遙感數據源

2.1.1 311 工作面概述

圖3 為研究區煤礦的311 工作面及其地表走向觀測線的沉降監測點位示意圖。311 工作面位于該礦井Ⅱ1 采區上山西翼，工作面可采走向長550m、面長240m。該工作面整體呈一單斜構造，煤層傾角8°～18°，平均傾角為10°，煤層厚度為2.6m～4.2m，平均采厚為3.0m。該工作面標高-160m～-280m，地面標高為+31m～+32.5m，2018 年2 月9 日開始回采，于2019 年7 月27 日開采結束。

圖3 井下開采的311 工作面及地表走向觀測線點位圖

2.1.2 遙感數據源

所使用的SAR 數據為Sentinel-1A（哨兵-1A）IW 成像模式SLC 數據，此SAR 數據為C 波段數據，空間分辨率為5m×20m（距離向×方位向），其他詳細參數見表1 所示。根據311 工作面開采的時間節點，獲取了2018 年2月8 日至2019 年9 月7 日間的46 景哨兵1A 影像數據。

表1 所用SAR 數據主要參數

2.2 基于D-InSAR 技術的開采沉降邊界分析

D-InSAR 處理的具體技術路線如圖4 所示，所涉及的主要技術內容包括：（1）影像預處理。包括數據格式轉換和目標區裁剪。主要是將原始影像數據處理為標準的SLC 數據，然后根據研究區范圍的大小，對影像進行目標區裁剪。（2）基線估算。估算結果包括時間基線、空間基線、多普勒偏移、一個相位變化周期代表的高程變化等信息，主要目的是評估像對干涉質量。（3）生成干涉圖。首先根據影像的精密星歷計算影像中心偏移值，并設置合適的窗口大小以進行自動配準，精度可以達到1 個像元的水平；然后將主輔影像與已有的外部基準DEM 進行配準，再經過采樣處理之后使配準精度達到0.001 個像元的水平；在此基礎上根據配準好的主輔影像的相位信息進行干涉，生成干涉條紋圖。（4）濾波和相干性計算。濾波處理主要目的是去除大氣延遲和平地效應所引起的誤差相位，同時得到相干系數圖。（5）相位解纏。使影像相位信息與地形信息相對應。（6）控制點選擇。控制點主要用于軌道精煉和重去平，常在去平后的干涉圖上進行選擇，選擇應遵循以下原則：避免在有地形相位未去除和變化的區域選擇控制點；選擇相干性高的區域；控制點應分布于整個范圍內；避免解纏錯誤的區域，如相位孤島等。（7）軌道精煉和重去平。目的是估算和去除殘余的恒定相位和解纏后還存在的相位坡道。（8）相位轉形變和地理編碼。解纏后的相位經過轉換就可得到相對于解纏起始點的視線向形變信息。地理編碼的目的是將所有的結果編碼到所需的坐標系下。

圖4 D-InSAR 技術路線圖

基于獲取的46 景哨兵1A 影像組成45 個干涉對，利用疊加二軌法D-InSAR 技術疊加了2018 年2 月8 日至2019 年9 月7 日的累積沉降，將每對D-InSAR 得到的沉降監測結果按時序疊加即可獲得研究區內的時序沉降，以此可確定311 工作面開采下沉影響邊界，如圖5 所示。311 工作面中心區域出現下沉梯度過大導致的地表形變失相關，因此監測到的最大下沉為-302mm，其原因是因為淹水或下沉梯度過大區域，InSAR 監測失相干，無法監測實際最大下沉量。

2.3 驗證分析與下沉邊界的確定

以下沉值-10mm 作為開采沉陷盆地邊界，聯合D-InSAR、SBAS-InSAR 技術以及概率積分預計法，以多手段、多角度來綜合權衡確定311 工作面南側實際下沉邊界。

2.3.1 基于D-InSAR 技術確定的沉降邊界位置

如圖5 所示，利用D-InSAR 技術監測到的下沉邊界大約在II311 工作面走向觀測線上的D013 點南側、距II311 工作面采空區邊緣的切眼位置平均距離為165m。

圖5 D-InSAR 技術確定的311 工作面實際下沉邊界

2.3.2 基于SBAS-InSAR 技術確定的沉降邊界位置

基于2018 年2 月8 日至2020 年9 月7 日的46 景哨兵1A 影像，再利用SBAS-InSAR 技術對311 工作面開采影響的下沉邊界進行監測，監測到的下沉邊界大約在II311 工作面走向觀測線上的D013 點、距II311 工作面采空區邊緣的切眼位置平均距離為160m。

2.3.3 基于概率積分法確定的沉降邊界位置

如圖6 所示，可看出基于概率積分預計模型所獲得的下沉邊界在II311 工作面走向觀測線上的D013 點、距II311 工作面采空區邊緣的切眼位置平均距離為160m。

2.3.4 基于InSAR 及概率積分法預測的下沉邊界綜合確定

通過綜合上述D-InSAR、SBAS-InSAR 技術與概率積分預計三種方法分析311 工作面實際下沉邊界，得到如圖6 所示多方法確定的311 工作面實際下沉邊界位置，與實際觀測線上的點位相符合。

圖6 多方法綜合的311 工作面南側實際下沉邊界對比圖

3 結束語

工礦區地表形變監測過程復雜，其破壞及危害性不可避免。傳統的變形監測技術費時費力，難以得到大范圍內的整體形變結果，隨著微波遙感的快速發展與應用，D-InSAR 和SBAS-InSAR 等技術已廣泛應用于地表沉降動態監測。本研究主要是利用D-InSAR 技術對淮北礦區某煤礦311 工作面開采的采動影響邊界進行確定研究，并通過概率積分法預計與SBAS-InSAR 干涉測量技術對D-InSAR 監測結果進行驗證，提取到311 工作面開采后的沉降邊界，同時根據311 工作面走向觀測線的觀測點位，在走向觀測線上指定了沉降邊界的具體位置，即在D013 測點處，該測點至311 工作面采空區切眼巷道的平面距離為162m。