劉金霖，譚志祥，范洪冬

(中國礦業(yè)大學(xué) 環(huán)境與測繪學(xué)院，江蘇 徐州 221116)

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

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

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

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

1 試驗(yàn)原理及流程

1.1 概率積分法

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

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

1.2 DS-InSAR分時(shí)段累加方法

1.2.1 FaSHPS 同質(zhì)點(diǎn)選取原理

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

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

其中，P{＊}為概率；μ(L) 為像元L的數(shù)學(xué)期望；z1-α/2為置信度為α?xí)r對應(yīng)的分位點(diǎn)； Var[A(L)]為像元L幅度的真實(shí)方差；N為樣本數(shù)。若像元后向散射特性穩(wěn)定，變異系數(shù)為常量(0.52)，則式(2)變換為

根據(jù)式(3)，若平均振幅A(L) 在該置信區(qū)間內(nèi)，即可判定兩者為同質(zhì)像元。

1.2.2 特征值分解相位優(yōu)化

由于分布式目標(biāo)像元內(nèi)地物的后散射特性是相同的，且都不占據(jù)像元信號的主導(dǎo)地位，其相位穩(wěn)定性也較差，因此要對相位進(jìn)行優(yōu)化，以達(dá)到去除噪聲相位的目的。試驗(yàn)選用特征值分解法[25-27]進(jìn)行相位優(yōu)化，其通過求解相干矩陣的特征值和特征向量，進(jìn)而確定不同散射信號對應(yīng)的特征值，并找到最大特征值對應(yīng)的特征向量，其所代表的散射體被視為主導(dǎo)散射體，從而實(shí)現(xiàn)對相位的優(yōu)化處理。

定義分布式目標(biāo)的相干性矩陣T的表達(dá)式為

式中，y=[y1,y2,y3,… ,yN]為分布式目標(biāo)的同質(zhì)像元在N景影像上的復(fù)數(shù)觀測量經(jīng)過時(shí)間維振幅歸一化處理后的復(fù)數(shù)觀測值；NSHPs為同質(zhì)點(diǎn)數(shù)量；Ω為全部分布式目標(biāo)同質(zhì)點(diǎn)的集合。

由特征值分解可得：

式中，λi為相干矩陣的特征值，當(dāng)λi越大時(shí)，其對應(yīng)的散射相位越占優(yōu)勢；μi為對應(yīng)的特征向量。

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

式中，λ1為最大特征值；μ1為其對應(yīng)的特征向量，即優(yōu)化后的分布式目標(biāo)相位。

在同質(zhì)像元選取和相位優(yōu)化完成后，通過計(jì)算擬合度，設(shè)置時(shí)間相干性閾值，從而選取相干性較好的分布式目標(biāo)點(diǎn)，進(jìn)行時(shí)序InSAR處理，便可獲得該時(shí)間段內(nèi)研究區(qū)的地表沉降信息。

1.2.3 三次插值法

1959年，DAVIDON首先提出三次插值法，其原理是通過構(gòu)建多項(xiàng)式來實(shí)現(xiàn)插值。

定義一個(gè)在給定區(qū)間[m,n]上的三次多項(xiàng)式為

求解得到：

假設(shè)，＊α為區(qū)間[m,n]上的插值點(diǎn)，則

求解得到插值點(diǎn)＊α為

根據(jù)三次插值原理，對每個(gè)時(shí)間段的地表沉降結(jié)果進(jìn)行插值后累加，即可獲得整個(gè)時(shí)間段內(nèi)研究區(qū)的地表沉降結(jié)果。

1.3 試驗(yàn)流程

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

圖1 試驗(yàn)流程Fig.1 Experimental flow chart

2 研究區(qū)概況

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

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

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

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

表1 工作面參數(shù)Table 1 Parameters of the working face

3 概率積分法預(yù)計(jì)結(jié)果分析

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

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

圖3 研究區(qū)預(yù)計(jì)下沉等值線Fig.3 Anticipated subsidence contour lines in the study area

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

4 DS-InSAR結(jié)果分析

4.1 試驗(yàn)數(shù)據(jù)

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

4.2 DS-InSAR結(jié)果

4.2.1 可靠性分析

為驗(yàn)證DS-InSAR監(jiān)測結(jié)果的可靠性，選取2021年3月至10月村莊內(nèi)23個(gè)監(jiān)測點(diǎn)的水準(zhǔn)實(shí)測數(shù)據(jù)，與DS-InSAR監(jiān)測結(jié)果進(jìn)行對比分析。其水準(zhǔn)點(diǎn)與村莊的位置關(guān)系如圖4所示。DS-InSAR監(jiān)測結(jié)果與水準(zhǔn)實(shí)測數(shù)據(jù)對比結(jié)果如圖5所示。將水準(zhǔn)實(shí)測數(shù)據(jù)視為真值，DS-InSAR監(jiān)測結(jié)果作為測量值，計(jì)算得到DS-InSAR監(jiān)測結(jié)果的最大誤差、平均誤差以及均方根誤差，結(jié)果見表2。

表2 DS-InSAR監(jiān)測結(jié)果精度評價(jià)Table 2 Accuracy evaluation of DS-InSAR monitoring results

圖4 水準(zhǔn)點(diǎn)與村莊位置關(guān)系Fig.4 Relationship between the leveling points and the village

圖5 InSAR與水準(zhǔn)方法監(jiān)測結(jié)果對比Fig.5 Comparison of monitoring results between InSAR and leveling

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

4.2.2 研究區(qū)時(shí)序監(jiān)測結(jié)果分析

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

圖6 各時(shí)間段研究區(qū)地表累計(jì)沉降Fig.6 Cumulative surface subsidence in the study area for each time period

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

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

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

4.2.3 監(jiān)測結(jié)果修正與對比分析

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

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

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

5 結(jié) 論

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

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