陳洋，陶秋香，劉國林，王路遙，王鳳云，王珂

山東科技大學(xué) 測繪與空間信息學(xué)院，山東青島 266590

0 引言

煤炭作為重要的能源基礎(chǔ)，是工業(yè)現(xiàn)代化和社會發(fā)展的有力支撐.但隨著能源需求的日益提升，煤炭開采在為世界經(jīng)濟創(chuàng)造巨大價值的同時，其引發(fā)的災(zāi)害問題也不容忽視：礦區(qū)大規(guī)模持續(xù)開采造成的沉降、塌陷，不僅損害地面基礎(chǔ)設(shè)施，還容易誘發(fā)各種地質(zhì)災(zāi)害，對人民財產(chǎn)安全造成巨大威脅(Yang et al., 2019; Saedpanah and Amanollahi, 2019; Sahoo and Khaoash, 2020).因此，為預(yù)防災(zāi)害發(fā)生，保證安全生產(chǎn)，必須對礦區(qū)進行長期有效的沉降監(jiān)測與防治.

合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar, InSAR)是一種地面沉降監(jiān)測新技術(shù)，可監(jiān)測到厘米級甚至毫米級的地面變化，是探測地震、火山、礦區(qū)等形變的有效工具(Klein et al., 2017; 韓宇飛等, 2010; 申文豪等, 2019).20世紀90年代，該技術(shù)開始用于礦區(qū)地表沉降監(jiān)測，所獲結(jié)果在空間、時間上均連續(xù)，并且可以很好地反映礦區(qū)沉降的結(jié)構(gòu)特征(Raymond and Rudant, 1997).相較于精密水準、GPS等傳統(tǒng)測量手段，InSAR表現(xiàn)出全天候、低成本等的獨特優(yōu)勢，為礦區(qū)地面沉降監(jiān)測提供了新的思路與技術(shù)支持(Yang et al., 2018; Liu et al., 2013).但礦區(qū)開采引起的地面沉降往往呈現(xiàn)沉降速度快、量值大，空間不連續(xù)等特征，InSAR技術(shù)受限于SAR影像的固定特性，在下沉較為嚴重的沉降盆地中心會出現(xiàn)失相干現(xiàn)象，無法得到準確的形變信息，很大程度上限制了該技術(shù)在實際工程中的推廣應(yīng)用(Luo et al., 2019; Carlà et al., 2018; Zhang et al., 2020).

為解決InSAR監(jiān)測礦區(qū)大梯度沉降時出現(xiàn)的失相干問題，國內(nèi)外學(xué)者首先從技術(shù)本身出發(fā)，通過減少對流層延遲、選取長波段數(shù)據(jù)等方法提升InSAR的大梯度沉降監(jiān)測能力，后又結(jié)合實際地質(zhì)模型展開進一步的研究與探索，其中概率積分法是應(yīng)用最為廣泛的一種礦區(qū)沉降預(yù)計模型(Xu et al., 2006; Wang et al., 2020; Ge et al., 2007; 陶秋香等, 2012; Fan et al., 2014).該方法是由我國學(xué)者劉寶琛、廖國華等基于隨機介質(zhì)理論提出的一種沉降預(yù)計方法，利用高精度監(jiān)測點，通過選定適當?shù)臄?shù)學(xué)模型，完成參數(shù)優(yōu)化反演，獲取沉降中心的下沉信息，有助于解決InSAR在礦區(qū)大梯度沉降監(jiān)測中的失相干問題(劉寶琛和廖國華, 1965; 劉寶琛和戴華陽, 2016).但該方法由于受到模型自身特性的限制，在沉降邊緣處收斂過快，尤其模型通常依賴成本高、耗時多的水準數(shù)據(jù)，導(dǎo)致其在礦區(qū)沉降監(jiān)測中仍存在一定的局限性.

綜合InSAR與概率積分法在礦區(qū)地表沉降監(jiān)測中的優(yōu)缺點，本文聯(lián)合InSAR與概率積分兩種監(jiān)測方法，獲取準確的礦區(qū)沉降盆地.該方法首先選取覆蓋礦區(qū)工作面的SAR影像集，利用InSAR技術(shù)，按照時間序列兩兩差分，并對其進行時序累積，獲取礦區(qū)多時段累積沉降監(jiān)測結(jié)果；然后基于相干性條件，選取InSAR沉降盆地邊緣的高相干點，與沉降中心附近少量插值后的水準數(shù)據(jù)作為特征點，通過曲面擬合法對概率積分模型參數(shù)進行優(yōu)化反演，得到概率積分預(yù)計的礦區(qū)地表沉降盆地；進而以相干性為基準，通過形變梯度理論確定兩種監(jiān)測結(jié)果的融合閾值，獲取完整的沉降盆地.論文以山東某礦區(qū)為研究區(qū)，利用該方法獲取礦區(qū)地表沉降信息，結(jié)果表明，該方法結(jié)合了InSAR與概率積分法的優(yōu)勢，既能反映礦區(qū)開采對周邊環(huán)境的影響，又能提取沉降中心準確的地表下沉信息，實現(xiàn)礦區(qū)沉降的高效、準確提取.

1 基本原理

1.1 InSAR形變監(jiān)測及最大形變梯度

實際觀測中，InSAR獲取的相位與雷達參數(shù)、天線位置、天線入射角及地面目標高程緊密相關(guān)，干涉過程中不可避免會受到各種誤差因素的影響，導(dǎo)致干涉相位信息中不僅包含高程和形變信息，還存在其他干擾項(Massonnet and Feigl, 1998)，即

φ=φflat+φtopo+φdef+φatm+φnoise,

(1)

由式(1)可以看出，InSAR獲取的相位信息主要包括：平地相位φflat，地形相位φtopo，地表形變相位φdef，電離層延遲和對流層延遲引起的大氣相位φatm以及噪聲相位φnoise.顯然，要想獲取地表形變相位，必須將干擾相位從干涉圖中分離.

當去除各類干擾后，其形變相位在干涉圖中表現(xiàn)為一系列干涉條紋，且干涉條紋所包含的相位信息并不是絕對相位，而是各個像元之間的相對相位.為此，Massonnet和Feigl(1998)提出形變梯度的概念，用以代表任意兩點間的相對位移，并給出了 InSAR 可監(jiān)測的最大形變梯度為

(2)

式中，dmax表示InSAR理論上能夠監(jiān)測到的形變梯度最大值，μ表示單個SAR影像像元的尺寸，λ表示雷達波長.由式(2)可見，dmax理論上與雷達本身參數(shù)有關(guān)，雷達波長越長、單位像元越小，能夠監(jiān)測到的形變梯度越大.以多視處理后的C波段Sentinel(分辨率為20 m)和L波段 ALOS(分辨率為10 m)數(shù)據(jù)為例，經(jīng)過1∶4多視處理的Sentinel影像，理論上能監(jiān)測到1.4×10-3的形變梯度，而經(jīng)過1∶3多視處理的ALOS影像，理論上能監(jiān)測到11.5×10-3的形變梯度.

式(2)是假定理想條件下，InSAR能監(jiān)測到的最大形變梯度.然而在實際情況中，InSAR不可避免地會受到軌道誤差、時空失相干、大氣延遲效應(yīng)等的影響，最大形變梯度監(jiān)測值往往比理論值要小.針對這一現(xiàn)象，Baran等(2005)分析研究了相干性對InSAR技術(shù)形變梯度的影響，給出了以影像相干性為基準的形變梯度閾值表達式，即

Dmax=dmax+0.002(γ-1),

(3)

式中，Dmax為實際InSAR最大可監(jiān)測形變梯度，γ為影像相干性.由式(3)可以看出，InSAR可監(jiān)測最大形變梯度與影像相干性呈線性正相關(guān)，相干性越高，實際可監(jiān)測最大形變梯度越大.對于SAR影像，一般存在某一相干閾值，使得Dmax=0.例如，對Sentinel影像，相干系數(shù)γ=0.3時，Dmax=0，即對于Sentinel數(shù)據(jù)來說，當影像相干性小于0.3時，就無法監(jiān)測到地表形變.通常情況下，礦區(qū)常常具有植被覆蓋率高、短時間內(nèi)形變量值大的特點，為了更為有效地提取形變信息，在礦區(qū)利用InSAR技術(shù)進行地表沉降監(jiān)測往往對相干性有著更高的要求.

1.2 概率積分模型

開采沉降研究中的巖體往往成因復(fù)雜，由于各類地質(zhì)作用，使得巖體生成過程中被各種斷層、節(jié)理等不均勻介質(zhì)切割，產(chǎn)生一系列次生裂隙及非連續(xù)面，而礦山開采會再次擾動巖體的原生結(jié)構(gòu)，呈現(xiàn)出更加明顯的非連續(xù)性，因此非連續(xù)介質(zhì)模型更適合開采沉降的應(yīng)用研究.20世紀50年代，Litwiniszyn(1956)基于非連續(xù)介質(zhì)模型，首次提出隨機介質(zhì)理論，并將其用于巖層移動規(guī)律的研究，認為隨機介質(zhì)的顆粒介質(zhì)模型可用來描述由于礦山開采引發(fā)的地下巖層及地表移動規(guī)律.20世紀60年代，我國學(xué)者劉寶琛和廖國華(1965)在此基礎(chǔ)上提出概率積分法，該方法將礦山巖層移動看作一個服從正態(tài)分布的隨機現(xiàn)象，結(jié)合概率密度函數(shù)模擬地表沉降盆地的形態(tài).對于走向長為D3、傾向長為D1的工作面，其開采引發(fā)的地面任意點(x,y)的沉降值可表示為

(4)

其中：

(5)

式(4)、(5)均在工作面坐標系下，其中W0表示最大下沉值，W0(x)、W0(y)分別表示待求點在走向、傾向主斷面上投影點的下沉值，u表示概率積分的開采單元，m為煤層厚度，H0為平均開采深度，H1、H2分別為下山、上山開采深度，q為下沉系數(shù)，tanβ為主要影響角正切，tanβ1、tanβ2分別表示下山、上山主要影響角正切，S3、S4分別表示左右邊界的拐點偏移距，S1、S2分別表示下山，上山拐點偏移距，θ為開采影響傳播角.

1.3 InSAR與概率積分法聯(lián)合的礦區(qū)地表沉降監(jiān)測方法

InSAR技術(shù)是獲取礦區(qū)沉降的高效手段，能夠準確獲取沉降的邊緣信息，但由于SAR影像自身特性的限制及各種噪聲誤差的干擾，該技術(shù)在礦區(qū)沉降量大的中心區(qū)無法探測到準確的沉降信息，監(jiān)測結(jié)果與實際情況相差較大.概率積分法在沉降盆地中心具有較高的監(jiān)測精度，但開采煤層容易受外力影響產(chǎn)生裂縫，加大煤柱上覆巖層的變形，導(dǎo)致概率積分法預(yù)計的沉降盆地邊緣收斂過快，在沉降邊緣的監(jiān)測精度降低(王正帥和鄧喀中, 2012).綜合InSAR與概率積分法兩種方法的沉降監(jiān)測優(yōu)勢，可以獲取準確完整的礦區(qū)地表沉降信息.其具體步驟如下：

(1)選取覆蓋礦區(qū)的N+1景SAR影像，按照時間序列對其進行兩兩差分處理，得到N個干涉對，并通過時序累積，獲取礦區(qū)時序形變信息dInSAR.

(2)收集礦區(qū)實際資料及工作面水準監(jiān)測數(shù)據(jù)，考慮到水準點高程值在各個監(jiān)測時段內(nèi)近似呈線性變化，對其進行時間域上的分段線性插值，實現(xiàn)InSAR與水準在時間上的統(tǒng)一，預(yù)測出SAR成像日期下，各水準點的實際沉降信息，即

(6)

式中，dlevling表示水準點在SAR成像日期下的沉降信息，Tlevling表示水準點觀測日期，TInSAR、T′InSAR分別表示距水準點觀測日期最近的前后兩個SAR成像日期，dInSAR、d′InSAR為這兩個SAR成像日期下該水準點的InSAR沉降監(jiān)測結(jié)果.

(3)選取InSAR沉降盆地邊緣的高相干點，與沉降中心附近少量插值后的水準數(shù)據(jù)作為特征點，并將其投影至工作面坐標系.假設(shè)坐標為(x,y)的特征點的沉降為D(x,y)，則按照曲面擬合法可將其表示為一個與特征點坐標(x,y)及概率積分預(yù)計參數(shù)B有關(guān)的函數(shù)，即

D(x,y)=f[(x,y);B],

(7)

用DZk(k=1,2,…,n)表示各特征點對應(yīng)的實際沉降，并使曲面擬合值D(x,y)與實際監(jiān)測結(jié)果的偏差滿足平方和最小，即

(8)

式中，V為各特征點實際監(jiān)測值相對于最小二乘擬合值的偏差.從而即可得到礦區(qū)概率積分參數(shù)的反演結(jié)果為

B=[q;tanβ;tanβ1;tanβ2;S3;S4;S1;S2;θ].

(9)

(4)根據(jù)反演的概率積分參數(shù)，按照式(4)、(5)，建立礦區(qū)概率積分沉降模型，并將沉降結(jié)果轉(zhuǎn)換至地理坐標系下.進一步利用克里金插值法，對區(qū)域化沉降進行最優(yōu)估計，獲取連續(xù)的概率積分沉降結(jié)果dPIM.

(5)統(tǒng)計N個干涉對的相干情況，計算其平均相干系數(shù)γavg，根據(jù)式(3)得到該相干條件下InSAR可監(jiān)測到的最大形變梯度Dmax.基于形變梯度理論，考慮到Sentinel-1A影像的分辨率是20×20 m，可以計算得到N個干涉對能監(jiān)測到的最大形變量dmax|InSAR(Massonnet and Feigl, 1998)，即

dmax|InSAR=μ·Dmax·N,

(10)

將dmax|InSAR作為InSAR與概率積分法沉降結(jié)果融合的閾值，從而得到完整的礦區(qū)地表沉降盆地，即

(11)

其數(shù)據(jù)處理流程如圖1所示.

圖1 InSAR與概率積分法聯(lián)合的礦區(qū)沉降監(jiān)測數(shù)據(jù)處理流程

2 研究區(qū)及數(shù)據(jù)源

2.1 研究區(qū)概況

試驗選取山東某礦區(qū)為研究區(qū)域，其地理位置及范圍如圖2a所示.礦區(qū)地屬黃河沖積平原，地形平坦，地面標高+40.53～+45.89 m，自然地形坡度為0.01%.礦區(qū)內(nèi)交通便利，與鄰近各縣、鄉(xiāng)均有公路相通，且有兗新鐵路、220國道在其南部通過.區(qū)內(nèi)河渠眾多，主要有宋金河、鄄鄆新河、鄆巨河等，水源豐富，地表以植被為主.區(qū)內(nèi)村莊眾多，人口稠密，因此對于礦區(qū)來說，一方面要保證開采安全，另一方面，又要考慮到周邊建筑、道路等地物的保護問題.同時兼顧礦區(qū)安全開采和地面維護，完成高效率、低損耗的生產(chǎn)，迫切需要有效的礦區(qū)地面沉降監(jiān)測.

圖2b中白色矩形表示該礦區(qū)1301工作面的開采范圍，該工作面上覆地表主要以植被為主，西北部為已開采的1300工作面老采空區(qū)，且有宋金河、鄆城新河自其上方穿過.工作面平均煤厚6.8 m，自2016年10月開始開采，傾斜寬223.4 m，計劃開采長度2265.8 m，至2018年3月，工作面推進1385.6 m，平均開采深度為938 m，地面標高為+41.8～+44 m.圖2(b)中綠色標記點表示工作面上方的水準點，沿傾向線自西向東布設(shè)有H1～H69，共69個水準點，沿走向線自南向北布設(shè)有Z1～Z92，共92個水準點，測量工作自2016年10月16日開始，至2018年3月24日結(jié)束，共19期.

圖2 研究區(qū)地理范圍

2.2 數(shù)據(jù)源

哨兵1號(Sentinel-1)衛(wèi)星于2014年4月成功發(fā)射，載有C波段合成孔徑雷達.該衛(wèi)星在軌運行高度為693 km，數(shù)據(jù)更新周期為12天，是一個可實現(xiàn)時間、空間連續(xù)觀測的雷達成像系統(tǒng)，實現(xiàn)了單、雙極化等若干種不同的極化方式.可提供連續(xù)影像(白天、夜晚和各種天氣)，以其大范圍、多模式、多應(yīng)用的特點為更多用戶提供了數(shù)據(jù)服務(wù).為獲取礦區(qū)的地面沉降情況，本試驗選取覆蓋礦區(qū)的21景C波段、VV極化的Sentinel-1A升軌數(shù)據(jù)(成像中心入射角約為38.92°)，時間跨度為2016年10月16日—2018年3月4日.具體參數(shù)如表1所示.

表1 SAR影像參數(shù)

InSAR技術(shù)需要通過外部數(shù)字高程模型(DEM)模擬地形相位，并將其從SAR干涉相位中移除，以提取地表沉降相位和沉降值.由于研究區(qū)域地形平坦，實驗選取3*3弧秒的SRTM3 DEM，可以滿足礦區(qū)地面沉降監(jiān)測的精度需求.SRTM3 DEM數(shù)據(jù)的采集間隔為3弧秒，地面分辨率為90 m，平均精度為16 m.

3 聯(lián)合法提取的礦區(qū)地表沉降結(jié)果及分析

3.1 InSAR礦區(qū)沉降監(jiān)測

將表1中的20個干涉對按照時間序列進行兩兩差分處理，考慮到干涉相位不可避免會受到多種噪聲的影響，利用自適應(yīng)濾波方法進行噪聲去除，得到濾波后增強的差分干涉相位及相應(yīng)的相干系數(shù)圖，如圖3、圖4所示.

圖3 濾波后的差分干涉圖示例

圖4 相干系數(shù)圖示例

由圖3可見，礦區(qū)在監(jiān)測初期，干涉條紋較少，之后隨著時間推移，矩形區(qū)域內(nèi)干涉條紋變得密集，因而判定該區(qū)域在研究時段內(nèi)可能發(fā)生地面沉降.結(jié)合圖4可知，第18、19個干涉對在監(jiān)測時段內(nèi)，相干性較好，相應(yīng)的差分干涉圖中顯現(xiàn)出清晰完整的干涉條紋；相比而言，第9、12個干涉對噪聲較大，出現(xiàn)嚴重的失相干，這是由于夏季植被覆蓋率較高或沉降梯度過大，超出影像的實際監(jiān)測能力導(dǎo)致的，對比此時的差分干涉圖，在沉降區(qū)域，干涉條紋雜亂，甚至未形成完整的條紋形狀.

將濾波增強后的差分干涉相位進一步進行解纏處理，并將解纏后的真實相位轉(zhuǎn)換成雷達視線向上的地表形變值，并通過投影、地理編碼得到地理坐標系下垂直向上的真實形變量.為了實現(xiàn)研究區(qū)域時序沉降盆地的監(jiān)測，假設(shè)起始監(jiān)測日期2016年10月16日無形變發(fā)生，將InSAR監(jiān)測結(jié)果進行時序累加，得到2016年10月16日—2018年3月4日期間各成像時刻的累積沉降結(jié)果，如圖5所示.

圖5 InSAR累積沉降結(jié)果

由圖5分析可以得到，1301工作面自2016年10月開始開采，至2016年12月，InSAR監(jiān)測到工作面起采線附近開始出現(xiàn)漏斗狀的沉降盆地，說明此時工作面的開采活動已對上覆地表產(chǎn)生影響.但由于開采時間較短，沉降量較小，截至2016年12月3日，InSAR監(jiān)測到的最大累積沉降為22 mm.隨后，由于開采工作面的持續(xù)推進，礦區(qū)沉降范圍及沉降量相對于監(jiān)測初期而言，明顯增大，至2017年6月1日，InSAR監(jiān)測到的累積沉降達到225 mm.隨著工作面持續(xù)向北推進，監(jiān)測結(jié)果顯示沉降中心開始沿工作面走向逐步向北移動，至2017年11月28日，監(jiān)測到的沉降中心距離起采線約160 m，最大累積沉降為420 mm.由此可見，InSAR監(jiān)測到的沉降位置、范圍及空間變化趨勢與礦區(qū)實際開采工作面相符，但在沉降中心處的大梯度沉降監(jiān)測能力明顯不足.隨著工作面的持續(xù)開采，礦區(qū)沉降繼續(xù)加大，至2018年3月4日，水準監(jiān)測到工作面最大累積沉降達到2813 mm，而InSAR監(jiān)測到的最大沉降僅為548 mm，與礦區(qū)實際存在較大差異.

3.2 特征點選取

為了構(gòu)建概率積分模型，進一步提取InSAR沉降結(jié)果中的高相干特征點.由InSAR監(jiān)測最大形變梯度式(3)可知，對于Sentinel-1A SAR影像來說，當相干性系數(shù)小于0.3時，InSAR監(jiān)測結(jié)果是不可信的.在實際情況中，由于礦區(qū)短時間內(nèi)沉降過大，其監(jiān)測精度對相干性有更高的要求.因此，根據(jù)圖4礦區(qū)的相干情況，這里將0.5作為相干性閾值，認為相干性大于0.5的特征點為高相干點，提取InSAR沉降盆地邊緣的13個高相干點，并根據(jù)概率積分原理，選取臨近礦區(qū)最大沉降點和拐點的5個水準點，共18個特征點，其疊加在工作面的分布情況如圖6所示.表2給出了各點的相干性系數(shù)及沉降信息.

表2 特征點相干系數(shù)及沉降信息

圖6 特征點分布圖

3.3 概率積分沉降預(yù)計

結(jié)合選取的18個特征點和概率積分模型，利用曲面擬合法，反演得到該礦區(qū)的概率積分參數(shù)為

B=[q;tanβ;tanβ1;tanβ2;S3;S4;S1;S2;θ]

=[1.3;1.2;1.9;1.6;47.5;47.5;34.6;34.6;84.6°]

將這些參數(shù)代入式(4)，計算得到礦區(qū)1301工作面的概率積分沉降預(yù)計盆地，如圖7所示.為獲取連續(xù)的概率積分地面沉降，考慮到Sentinel-1A影像分辨率為20×20 m，在礦區(qū)1301工作面沉降盆地處，以20 m為間隔共選取18396個特征點，計算其沉降值，并將結(jié)果轉(zhuǎn)換至地理坐標系下，進一步利用克里金插值法，對區(qū)域化沉降進行最優(yōu)估計.圖8給出了插值后空間連續(xù)的概率積分礦區(qū)沉降結(jié)果.

由圖7、圖8可以看出，研究時段內(nèi)，概率積分法探測到1301工作面上方呈現(xiàn)明顯的漏斗狀沉降盆地，最大沉降為2528 mm.但該方法測得工作面周邊幾乎無沉降發(fā)生，沉降盆地在邊緣處收斂過快，其原因是在實際情況下，礦區(qū)沉降往往存在復(fù)雜地形、地下水流失及其他采空區(qū)開采等外部因素的影響，而概率積分模型僅考慮開采活動引發(fā)的地面沉降，導(dǎo)致其監(jiān)測到的地面沉降集中在工作面上方.

圖7 概率積分礦區(qū)沉降預(yù)計盆地立體圖

圖8 連續(xù)的概率積分沉降結(jié)果

3.4 聯(lián)合法提取完整的礦區(qū)地表沉降盆地

對圖4所示矩形區(qū)域內(nèi)的相干系數(shù)結(jié)果進行均值濾波處理，并統(tǒng)計20個干涉對中沉降區(qū)域的相干情況，統(tǒng)計結(jié)果如表3所示.分析表3可知，2016年10月16日—2018年3月4日期間，礦區(qū)沉降區(qū)域的相干性為0.39～0.64.其中夏季4—8月，為農(nóng)作物生長時節(jié)，植被覆蓋率高，獲取的SAR影像相干性普遍較差，秋冬季節(jié)SAR影像質(zhì)量相對較好.

根據(jù)表3各干涉對的相干結(jié)果，計算得到研究時段內(nèi)20個干涉對的平均相干系數(shù)γavg為0.502.結(jié)合式(3)、式(10)，考慮到Sentinel-1A的影像分辨率為20×20 m，進一步計算得到此時InSAR能監(jiān)測到的最大形變量dmax|InSAR為162 mm.以此為閾值，對InSAR沉降邊緣和概率積分沉降中心結(jié)果進行融合，得到礦區(qū)完整的沉降盆地，其結(jié)果如圖9所示.

表3 礦區(qū)20個干涉對的相干系數(shù)統(tǒng)計

由圖9可以看出，InSAR與概率積分法通過融合，可以得到礦區(qū)完整的沉降結(jié)果，但在融合的邊界卻出現(xiàn)結(jié)果的不連續(xù)現(xiàn)象.為解決這一問題，本文在融合邊界的兩側(cè)尋找一定區(qū)域內(nèi)的InSAR與概率積分法預(yù)計的沉降數(shù)據(jù)，對其進行加權(quán)平均.

圖9 聯(lián)合法沉降監(jiān)測結(jié)果

由InSAR監(jiān)測結(jié)果可知，20個干涉對的平均相干系數(shù)γavg為0.502，但處于夏季的干涉對干涉效果差，相干性低，其相干系數(shù)低于平均相干系數(shù)，監(jiān)測精度較低.因此，在加權(quán)平均中，選擇相干性大于平均相干系數(shù)γavg的第1、2、3、4、5、17、18、19個干涉對，計算其平均相干系數(shù)為0.58.在該相干條件下，結(jié)合式(3)、式(10)，考慮到Sentinel-1A的影像分辨率為20×20 m，計算得到此時InSAR能監(jiān)測到的最大形變量為90 mm.以此作為加權(quán)平均的下臨界閾值，當沉降量小于90 mm時，以InSAR監(jiān)測值作為最終沉降結(jié)果.

對于概率積分法而言，拐點通常位于沉降盆地主斷面上沉降曲線曲率為零的分界點，根據(jù)概率積分預(yù)計的沉降結(jié)果，以工作面走向起采線與終采線處概率積分沉降結(jié)果的平均值917 mm作為加權(quán)平均的上臨界閾值.當沉降量大于917 mm時，以概率積分預(yù)計的沉降值作為最終沉降結(jié)果.

當沉降量介于90 mm到917 mm時，通過對InSAR與概率積分預(yù)計的沉降值進行加權(quán)平均，確定最終沉降結(jié)果.首先選擇沉降量介于90～917 mm的水準點，共48個，計算其InSAR、概率積分法沉降結(jié)果與水準的絕對差值，如圖10所示.然后計算48個水準點InSAR監(jiān)測結(jié)果、概率積分法預(yù)計結(jié)果的累積均方差，分別為186 mm、208 mm.根據(jù)

圖10 部分水準點的誤差情況

(12)

計算得到PPIM=0.8PInSAR.這里取PPIM=0.44，PInSAR=0.56，則由

(13)

得到加權(quán)平均后完整的礦區(qū)沉降盆地，如圖11所示.

圖11 加權(quán)融合后的聯(lián)合法沉降監(jiān)測結(jié)果

由圖11可知，2016年10月16日—2018年3月4日期間，聯(lián)合法監(jiān)測到礦區(qū)出現(xiàn)地面沉降，且沉降盆地位置與礦區(qū)1301工作面高度吻合.研究時段內(nèi)，該方法監(jiān)測到工作面周邊存在小量級沉降，這是因為工作面附近存在一系列小斷層，開采工作的持續(xù)推進，擾動了其上覆巖層，產(chǎn)生微小形變，對周邊環(huán)境造成不同程度的影響，與礦區(qū)實際開采進度及地質(zhì)情況相符，說明該方法在沉降盆地邊緣的監(jiān)測效果優(yōu)于概率積分法.由聯(lián)合法在沉降盆地中心的監(jiān)測結(jié)果分析可以得到，截至2018年3月4日，該方法監(jiān)測到1301工作面的沉降中心位于距離起采線約600 m處，沉降中心向東部偏移，這主要是因為1301工作面存在約10°的地面傾斜，沉降中心向下山方向偏移，與礦區(qū)實際情況相符，其監(jiān)測到的最大累積沉降達到2528 mm，在沉降中心的監(jiān)測能力明顯優(yōu)于InSAR手段.

4 聯(lián)合法監(jiān)測結(jié)果精度評估

為進一步驗證聯(lián)合法監(jiān)測結(jié)果的精度，提取礦區(qū)1301工作面傾向線上H1～H69，走向線上Z1～Z92，共161個水準點的聯(lián)合法監(jiān)測結(jié)果、InSAR累積沉降、概率積分法預(yù)計沉降，并收集實際水準測量數(shù)據(jù)，對三種監(jiān)測方法的測量結(jié)果進行定量對比分析，圖12給出了對比結(jié)果圖.以InSAR能監(jiān)測到的最大形變量dmax|InSAR=162 mm為區(qū)分沉降盆地邊緣和中心的閾值，分別統(tǒng)計三種監(jiān)測方法在沉降盆地邊緣和中心的誤差，表4給出了三種監(jiān)測方法的誤差對比情況.

表4 誤差統(tǒng)計

由圖12、表4三種監(jiān)測方法的對比結(jié)果，分析可以得到：

圖12 監(jiān)測結(jié)果對比圖

(1)2016年10月16日至2018年3月4日期間，聯(lián)合法、InSAR和概率積分法均探測到1301工作面上方出現(xiàn)地面沉降，監(jiān)測得到的沉降盆地位置、沉降曲線形態(tài)及變化趨勢基本一致，且與水準監(jiān)測結(jié)果相符.

(2)通過與概率積分法結(jié)果對比，聯(lián)合法在沉降邊緣處與水準吻合程度更高.分析圖12(c)、(e)中二者在傾向線上的監(jiān)測結(jié)果可以得到，在沉降邊緣處，由于1301工作面西部存在小斷層，受其影響，導(dǎo)致沉降范圍向西部擴大，聯(lián)合法監(jiān)測到這一沉降趨勢，在H1～H23范圍內(nèi)，其測得的沉降曲線與水準吻合較好.對比此時概率積分法的監(jiān)測結(jié)果，發(fā)現(xiàn)該方法探測到的沉降范圍略小，沉降曲線在邊緣處收斂較快，分析是由于概率積分法在反演過程中并未考慮外部條件對沉降的影響，造成該方法對于受到斷層影響而引起的沉降不敏感.分析圖12(d)、(f)中走向線上二者的沉降曲線同樣發(fā)現(xiàn)，在沉降盆地的邊緣區(qū)域，聯(lián)合法可以獲得與實際更相符的監(jiān)測結(jié)果.在沉降邊緣處，聯(lián)合法監(jiān)測結(jié)果的平均絕對誤差為30 mm，均方根誤差為40 mm，最大絕對誤差為104 mm，相比概率積分法而言，三種誤差分別提升53%、69%和84%.

(3)通過與InSAR結(jié)果對比，聯(lián)合法在沉降中心處的大梯度監(jiān)測能力明顯提升.圖12(a,e)二者的監(jiān)測結(jié)果顯示，水準監(jiān)測到傾向線上32點處出現(xiàn)2813 mm的最大沉降.InSAR雖也探測到此處為沉降中心，但監(jiān)測到的最大沉降點與水準略有差別，出現(xiàn)在29點處.監(jiān)測結(jié)果與水準相比，最大絕對誤差達到2445 mm，這主要是由于礦區(qū)沉降中心沉降梯度過大，超出Sentinel-1A影像的實際監(jiān)測范圍，加上地表植被、噪聲等降低了SAR影像質(zhì)量，導(dǎo)致其在沉降中心處無法探測出真實的沉降信息.對比聯(lián)合法獲得的沉降結(jié)果可以看出，在沉降中心處，該方法監(jiān)測結(jié)果與水準吻合較好，其監(jiān)測得到的最大沉降點也與水準結(jié)果一致，監(jiān)測到的最大累積沉降為2510 mm，與水準相差303 mm.圖12(b,f)為走向線上兩種技術(shù)在沉降盆地中心的監(jiān)測結(jié)果，表現(xiàn)出與傾向線上相同的特征，InSAR與聯(lián)合法監(jiān)測到的最大沉降分別為469 mm和2510 mm，與水準的最大絕對誤差分別為2396 mm和706 mm.在沉降中心處，聯(lián)合法監(jiān)測結(jié)果的平均絕對誤差為198 mm，均方根誤差為258 mm，最大絕對誤差為706 mm，最大沉降點處的絕對誤差為303 mm，相比InSAR而言，四種誤差分別提升81%、80%、71%和88%.

(4)根據(jù)三種監(jiān)測方法的誤差統(tǒng)計分析可知，礦區(qū)傾向線、走向線上所有水準點的聯(lián)合法監(jiān)測結(jié)果的平均絕對誤差為121 mm，均方根誤差為192 mm，最大絕對誤差為706 mm，最大沉降點處的絕對誤差為303 mm.相比InSAR和概率積分法而言，整體監(jiān)測精度有所提升，四種誤差相較于InSAR分別提升79%、80%、71%和88%，平均絕對誤差與均方根誤差相較于概率積分法分別提升27%和17%.

5 結(jié)論

本文針對InSAR和概率積分法目前各自監(jiān)測礦區(qū)地表沉降存在的不足，基于影像相干性及InSAR形變梯度理論，提出了兩種技術(shù)聯(lián)合進行礦區(qū)沉降監(jiān)測的方法.以山東某礦區(qū)1301工作面為研究區(qū)，利用該方法進行沉降監(jiān)測，獲取了研究時段內(nèi)礦區(qū)沉降盆地.通過與同期水準實測數(shù)據(jù)、InSAR累積沉降和概率積分監(jiān)測結(jié)果進行對比分析，得出如下主要結(jié)論：

(1)聯(lián)合法相比概率積分法而言，在沉降盆地邊緣的監(jiān)測精度有所提升，可以在減少水準觀測的基礎(chǔ)上，很好地避免模型本身的缺陷.該方法在沉降邊緣處與水準的平均絕對誤差、均方根誤差和最大絕對誤差相較于概率積分法分別提升53%、69%和84%，可以解決概率積分法在礦區(qū)地表沉降監(jiān)測中收斂過快的問題，監(jiān)測結(jié)果與水準更為吻合.

(2)聯(lián)合法相比InSAR技術(shù)而言，在沉降盆地中心的監(jiān)測能力明顯提升.該方法在沉降中心處的平均絕對誤差、均方根誤差、最大絕對誤差和最大沉降點處的絕對誤差相較于InSAR分別提升81%、80%、71%和88%，能夠有效彌補InSAR在監(jiān)測礦區(qū)大梯度地表沉降時的不足.

(3)聯(lián)合法結(jié)合InSAR與概率積分法各自的優(yōu)點，可以準確探測出礦區(qū)沉降盆地的位置、范圍及沉降變化趨勢，反映礦區(qū)開采對周邊環(huán)境的影響，且監(jiān)測得到的沉降中心符合礦區(qū)開采的地表移動規(guī)律，能夠提取沉降中心準確的下沉信息.相比概率積分法和InSAR技術(shù)而言，該方法的整體監(jiān)測精度有所提升，可以獲得與實際情況更為吻合的地表沉降.

致謝感謝歐洲航天局(European Space Agency, ESA)提供免費的Sentinel-1A SAR數(shù)據(jù)和美國國家航空航天局(National Aeronautics and Space Administration, NASA)提供SRTM3 DEM數(shù)據(jù).