InSAR與概率積分法聯合的礦區地表沉降精細化監測方法

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

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

0 引言

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

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

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

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

1 基本原理

1.1 InSAR形變監測及最大形變梯度

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

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

(1)

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

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

(2)

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

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

Dmax=dmax+0.002(γ-1),

(3)

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

1.2 概率積分模型

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

(4)

其中：

(5)

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

1.3 InSAR與概率積分法聯合的礦區地表沉降監測方法

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

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

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

(6)

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

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

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

(7)

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

(8)

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

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

(9)

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

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

dmax|InSAR=μ·Dmax·N,

(10)

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

(11)

其數據處理流程如圖1所示.

圖1 InSAR與概率積分法聯合的礦區沉降監測數據處理流程

2 研究區及數據源

2.1 研究區概況

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

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

圖2 研究區地理范圍

2.2 數據源

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

表1 SAR影像參數

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

3 聯合法提取的礦區地表沉降結果及分析

3.1 InSAR礦區沉降監測

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

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

圖4 相干系數圖示例

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

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

圖5 InSAR累積沉降結果

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

3.2 特征點選取

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

表2 特征點相干系數及沉降信息

圖6 特征點分布圖

3.3 概率積分沉降預計

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

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°]

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

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

圖7 概率積分礦區沉降預計盆地立體圖

圖8 連續的概率積分沉降結果

3.4 聯合法提取完整的礦區地表沉降盆地

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

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

表3 礦區20個干涉對的相干系數統計

由圖9可以看出，InSAR與概率積分法通過融合，可以得到礦區完整的沉降結果，但在融合的邊界卻出現結果的不連續現象.為解決這一問題，本文在融合邊界的兩側尋找一定區域內的InSAR與概率積分法預計的沉降數據，對其進行加權平均.

圖9 聯合法沉降監測結果

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

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

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

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

(12)

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

(13)

得到加權平均后完整的礦區沉降盆地，如圖11所示.

圖11 加權融合后的聯合法沉降監測結果

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

4 聯合法監測結果精度評估

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

表4 誤差統計

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

圖12 監測結果對比圖

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

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

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

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

5 結論

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

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

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

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

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