融合DS-InSAR與PIM的煤礦開采影響范圍確定方法

劉金霖，譚志祥，范洪冬

(中國礦業大學 環境與測繪學院，江蘇 徐州 221116)

煤礦開采會導致地面沉降，對沉陷區內的房屋、鐵路、公路等工建(構)筑物造成一定程度的損害。我國部分煤礦資源位于村莊下方，此類煤礦開采時一定程度上會導致村莊房屋損毀；無開采作業時，很多房屋因質量等原因也會造成房屋出現裂縫等問題，因此煤礦企業僅希望承擔因煤炭開采引起的損毀責任，然而損害原因不清經常造成礦村糾紛，影響礦區和諧發展，因此開采損害技術鑒定是一項非常必要的工作[1]。

開采損害技術鑒定包括確定煤礦的開采影響邊界和影響程度兩個方面的內容。其中，開采損害影響程度可根據相關規范要求進行判定，一般將建筑物的破壞情況劃分為Ⅰ、Ⅱ、Ⅲ、Ⅳ等級。然而，確定開采影響的邊界比較復雜[2]。煤礦開采損害鑒定離不開對開采沉陷影響因素、監測手段、沉陷規律、沉陷預計等的研究。蔣旭梓[3]等對濟寧井工礦區不同開采條件對沉陷階段的影響進行了研究；孫國慶[4]等對頂板下沉因素進行了分析；康新亮[5]等對多種地質下開采沉陷規律分析進行了研究；張向陽[6]等對沉陷區內路基沉降變形規律進行了研究；蔣春梅[7]等對地面下沉規律進行了分析；孫慶先[8]等提出了物探技術和開采沉陷學理論相結合的開采損害技術鑒定方法；文獻[9-10]等對InSAR測量技術進行了研究；鄧喀中[11]等利用D-InSAR技術對開采沉陷進行了監測研究；郭山川[12]等利用DInSAR技術對黃土高原礦區地表形變進行了監測；李柱[13]等利用DS-InSAR技術對烏達煤田火區長時序地表形變進行了監測；張文龍[14]等結合PS-InSAR和DS-InSAR技術對韓城礦區采空區地面塌陷進行了識別；楊雋[15]等將SBAS-InSAR技術對形變特征的監測用于礦區損害鑒定分析；李學良[16]對煤礦老采空區覆巖移動變形監測方法進行了分析；吳群英[17]等對荒漠化礦區的生態環境進行了監測；陳毓[18]對地質災害風險評估與適應性評價進行了研究；張官進[19]等對礦區的開采沉陷進行了預測；譚志祥[20]等提出在準確掌握采礦資料時，鑒定方法有地表移動變形預計、移動角計算、相似材料模擬以及根據房屋和地表裂縫特征判別等方法；李曉斌[21]等對高強度開采地表損傷程度分類進行了研究。

上述研究為開采影響邊界研究提供了有益借鑒，但針對多個工作面較長時間開采導致的房屋損害鑒定問題，傳統方法為理論預計結果，缺少實測數據支撐，導致礦村雙方對邊界位置存在較大爭議，而常規InSAR技術由于時空失相干，無法獲取長時間的地表沉降信息。

基于此，筆者提出了一種基于DS-InSAR的分時段累加方法，實現對華東地區某村莊2016年7月至2021年10月期間的地表沉降監測，同時結合概率積分法，綜合分析劃定了該煤礦3個工作面的開采影響邊界，得到了礦村雙方的認可，解決了礦村矛盾，為開采損害鑒定提供了新的方法。

1 試驗原理及流程

1.1 概率積分法

1965年，我國學者劉寶琛、廖國華根據波蘭學者J.LITWINISYN的隨機介質理論提出了概率積分法，并在我國廣泛應用[22]。根據概率積分法，煤礦開采后引起地面某點的下沉值為

式中，L1，L2分別為工作面走向長度和傾向寬度；1l，2l分別為沿走向、傾向計算的開采寬度；m為煤層厚度；S0，Sdown，Sup分別為走向、傾向下山和上山拐點偏移距；Wmax為充分采動時地表最大下沉值；α為煤層傾角；q為下沉系數；θ為開采影響傳播角；β為主要影響傳播角；r為主要影響半徑；H為采深。

1.2 DS-InSAR分時段累加方法

1.2.1 FaSHPS 同質點選取原理

同質像元識別是指識別在地表形變監測區域內SAR影像數據集中具有相同變化特征的像元，包括參數檢驗和非參數假設檢驗的兩種識別方法。筆者選用參數檢驗中的快速同質點選取方法(FaSHPS)[23-24]，它是將假設檢驗問題轉化為置信區間估計問題，比較兩個像元是否服從相同的分布函數，從而判斷它們是否為同質像元。

根據中心極限定理，在SAR影像數目足夠的情況下，像元平均振幅A(L) 可以近似服從高斯分布，區間估計可表示為

其中，P{＊}為概率；μ(L) 為像元L的數學期望；z1-α/2為置信度為α時對應的分位點； Var[A(L)]為像元L幅度的真實方差；N為樣本數。若像元后向散射特性穩定，變異系數為常量(0.52)，則式(2)變換為

根據式(3)，若平均振幅A(L) 在該置信區間內，即可判定兩者為同質像元。

1.2.2 特征值分解相位優化

由于分布式目標像元內地物的后散射特性是相同的，且都不占據像元信號的主導地位，其相位穩定性也較差，因此要對相位進行優化，以達到去除噪聲相位的目的。試驗選用特征值分解法[25-27]進行相位優化，其通過求解相干矩陣的特征值和特征向量，進而確定不同散射信號對應的特征值，并找到最大特征值對應的特征向量，其所代表的散射體被視為主導散射體，從而實現對相位的優化處理。

定義分布式目標的相干性矩陣T的表達式為

式中，y=[y1,y2,y3,… ,yN]為分布式目標的同質像元在N景影像上的復數觀測量經過時間維振幅歸一化處理后的復數觀測值；NSHPs為同質點數量；Ω為全部分布式目標同質點的集合。

由特征值分解可得：

式中，λi為相干矩陣的特征值，當λi越大時，其對應的散射相位越占優勢；μi為對應的特征向量。

因散射體間相互獨立，則相干性矩陣可分解為主散射體信號Tsignal和噪聲信號Tnoise，即

式中，λ1為最大特征值；μ1為其對應的特征向量，即優化后的分布式目標相位。

在同質像元選取和相位優化完成后，通過計算擬合度，設置時間相干性閾值，從而選取相干性較好的分布式目標點，進行時序InSAR處理，便可獲得該時間段內研究區的地表沉降信息。

1.2.3 三次插值法

1959年，DAVIDON首先提出三次插值法，其原理是通過構建多項式來實現插值。

定義一個在給定區間[m,n]上的三次多項式為

求解得到：

假設，＊α為區間[m,n]上的插值點，則

求解得到插值點＊α為

根據三次插值原理，對每個時間段的地表沉降結果進行插值后累加，即可獲得整個時間段內研究區的地表沉降結果。

1.3 試驗流程

試驗的具體處理流程如圖1所示。①根據已知地質資料，選取合適的預計參數，利用概率積分法預計鑒定區的沉降情況；②利用快速同質點選取算法和特征值分解相位優化方法獲得DS點；③經過相位解纏即可提取地表形變情況，對每個時間段的沉降結果進行插值疊加后可獲取鑒定區長時序的地表形變信息，對PIM預計結果進行修正、驗證；④確定多個工作面的開采影響邊界，完成開采損害鑒定。

圖1 試驗流程Fig.1 Experimental flow chart

2 研究區概況

研究區域位于華東地區中部，其地表屬淮河沖積平原，地形平坦，海拔標高+23.52～+27.58 m，地勢西高東低，地表水系密布。需鑒定的村莊北側分布有3個工作面，分別為121302，121303，121304工作面，工作面位置如圖2所示。

圖2 工作面與村莊位置Fig.2 Position of the working face and the village

由圖2中工作面與村莊位置可知，工作面自北向南推進，采用后退式傾向長壁一次采全高采煤法采煤，全部垮落法管理頂板。

3個工作面的具體參數見表1，鑒定區表土層厚度580 m，工作面上覆巖層以砂巖、泥巖為主，巖性中硬偏軟。經過筆者現場勘察和調研，村莊部分房屋受到不同程度的損壞，村民認為是工作面開采導致，但是企業認為是由于村民在村莊北側大量取土引起房屋部分地面及墻體產生裂痕。因此，需要對該村莊進行開采損害鑒定，以判定村礦之間的責任。

表1 工作面參數Table 1 Parameters of the working face

3 概率積分法預計結果分析

文獻[28-30]中給出了該礦首采工作面及附近煤礦10多個地質采礦條件相似工作面的巖層移動參數，結合鑒定區工作面具體情況，通過綜合分析，確定采用的概率積分法預計參數為：下沉系數q=0.99，水平移動系數b=0.35，主要影響傳播角正切tanβ=1.8，開采影響傳播角θ=86°，拐點偏移距s=0.03H。

基于上述參數，首先利用中國礦業大學開采損害及防護研究所研制的礦區開采沉陷預測分析系統(MSAS)，對開采后引起的村莊地表沉陷情況進行了計算，其下沉等值線如圖3所示。

圖3 研究區預計下沉等值線Fig.3 Anticipated subsidence contour lines in the study area

由圖3中的研究區預計下沉等值線結果可知，3個工作面的開采使村莊部分地表產生了沉降，影響最大的位于村莊東北部，開采影響邊界位于村莊南部，附近煤礦的開采對村莊影響范圍較廣。礦村雙方認為概率積分法是理論預計模型，不能反映實際地表沉降情況，對開采影響邊界具體位置存在較大爭議。因此筆者利用DS-InSAR分時段累加方法，獲取研究區地表沉降結果，對概率積分法預計結果進行修正和驗證。

4 DS-InSAR結果分析

4.1 試驗數據

研究使用歐空局發射的Sentinel-1A衛星獲取的試驗數據，從中選取153景覆蓋研究區域的衛星升軌影像，時間跨度為2016年7月至2021年10月。試驗采用了90 m分辨率的航天飛機雷達測圖計劃數字高程模型作為外部DEM數據，以實現地形相位去除的目的。為了降低時空失相干的影響，試驗將影像分為2016年7月至2018年6月，2018年6月至2020年6月和2020年6月至2021年10月3個時間段，采用多主影像策略，對數據進行DS-InSAR技術處理。

4.2 DS-InSAR結果

4.2.1 可靠性分析

為驗證DS-InSAR監測結果的可靠性，選取2021年3月至10月村莊內23個監測點的水準實測數據，與DS-InSAR監測結果進行對比分析。其水準點與村莊的位置關系如圖4所示。DS-InSAR監測結果與水準實測數據對比結果如圖5所示。將水準實測數據視為真值，DS-InSAR監測結果作為測量值，計算得到DS-InSAR監測結果的最大誤差、平均誤差以及均方根誤差，結果見表2。

表2 DS-InSAR監測結果精度評價Table 2 Accuracy evaluation of DS-InSAR monitoring results

圖4 水準點與村莊位置關系Fig.4 Relationship between the leveling points and the village

圖5 InSAR與水準方法監測結果對比Fig.5 Comparison of monitoring results between InSAR and leveling

結合圖5和表2分析可知，DS-InSAR獲取的地表沉降趨勢基本與水準數據一致，誤差較小，證明DSInSAR技術在該研究區的監測結果具有較高的可靠性。

4.2.2 研究區時序監測結果分析

對于每個時間段，計算其他SAR影像相對于第一景影像的沉降值，然后將其從衛星的視線方向(Line of Sight，LOS)轉換為垂直方向，即可得到研究區域在該時間段內的累積地表沉降值。筆者將累積沉降圖與地下工作面的位置進行疊加，以便更加直觀展示出沉降區域的形成與工作面開采之間的關系。2016年7月至2018年6月，2018年6月至2020年6月和2020年6月至2021年10月3個時間段的研究區地表累積沉降信息如圖6所示。

圖6 各時間段研究區地表累計沉降Fig.6 Cumulative surface subsidence in the study area for each time period

由圖6(a)可知，2016年7月至2018年6月研究區地表出現了較為明顯的沉降區域，其位置基本與121303工作面和121304工作面相對應。由于工作面區域地表沉降較大、中央積水，DS-InSAR無法監測到工作面中心的沉降值，可監測到的最大沉降值為782 mm，位于121303工作面和121304工作面邊緣，村莊內監測到的最大沉降值為309 mm，位于村莊北部。由此可見，121303工作面和121304工作面的開采對村莊中部和北部均產生一定影響。由圖6(b)可知，2018年6月至2020年6月的DS-InSAR監測結果中也出現了較為明顯的沉降區域，其位置基本與121302工作面相對應，監測到的最大沉降值為804 mm。該時間段內村莊內依舊受到工作面開采的影響，地表產生了一定沉降。根據圖6(c)可知，2020年6月至2021年10月，121303，121304和121302三個工作面均開采結束，但是工作面附近仍有殘余沉降，雖然其沉降值明顯減小，但仍對村莊東北部產生了一定影響。

為了更好地展示2016年7月至2021年10月期間，121303，121304和121302工作面的開采對村莊產生的沉降影響，筆者將3個時間段的DS-InSAR監測結果，通過插值處理后，對同一地理位置的沉降值進行疊加，獲得了累計沉降情況，如圖7所示。根據DS-InSAR疊加沉降結果，3個工作面的開采對地表產生最大沉降值約為1 496 mm，位于工作面邊界，對村莊影響較大，村莊地表產生了較為明顯的下沉，計算得到的最大沉降值約為576 mm，位于村莊東北部。

圖7 2016年7月至2021年10月研究區地表累計沉降Fig.7 Cumulative surface subsidence of the study area from July 2016 to October 2021

4.2.3 監測結果修正與對比分析

由于概率積分法預計結果在采空區邊界上方收斂較快[31]，導致其預計得到的下沉邊界較實際要小，筆者根據DS-InSAR技術獲取的實際地表沉降結果，繪制出研究區地表下沉等值線，據此對概率積分法預計得到的村莊下沉等值線進行向外偏移修正，最終劃定出3個工作面對村莊的開采影響邊界，位置如圖8所示?？梢?，開采影響邊界以內的房屋等建筑均受到3個工作面開采的影響。

圖8 研究區DS-InSAR下沉等值線Fig.8 Subsidence contour lines in the study area using DS-InSAR

通常開采沉陷影響邊界由巖層移動的邊界角確定，但是邊界角受采深、采寬、表土層厚度、覆巖巖性、采動程度等諸多因素的影響，通常很難準確選取，用類比法獲得的邊界角往往存在較大誤差，從而導致其確定的開采沉陷邊界與實際相差較大；由上述研究結果可知，InSAR技術在確定開采區域邊界上方小變形時結果較準確，在條件許可情況下可采用InSAR技術結合PIM方法確定開采沉陷影響邊界。

5 結 論

(1) 利用DS-InSAR對時間跨度長達5 a的153景Sentinel-1A數據進行分段處理，并將分段沉降監測結果進行疊加，獲取了研究區域長時序大變形的地表沉降信息，有效減小了時空失相干對InSAR監測結果的影響，提升了InSAR技術在長時間地表形變監測方面的能力。

(2) 通過DS-InSAR對概率積分法的預計結果進行驗證及修正，劃定了開采影響邊界，得到了礦村雙方的認可，為開采損害技術鑒定提供了新的參考，豐富了開采沉陷學領域的學科內容。