基于魯棒M 估計協方差矩陣相位優化的DSInSAR 礦區形變監測研究

曾祥凱，孫鳳娜，陳東興，段文再，滕文龍

（1.山東省煤田地質局第二勘探隊，山東 濟寧 272000；2.陽泉高新技術產業開發區管理委員會，山西 陽泉 045000；3.中國礦業大學環境與測繪學院，江蘇 徐州 221116）

0 引言

煤炭作為我國主體能源，2022 年煤炭消費量占能源消費總量的56.2%，較2021 年上升了0.3 個百分點，在支撐我國經濟社會平穩較快發展進程中持續發揮著能源安全保障的“壓艙石”和“穩定器”作用，即便在當前“雙碳”目標下，其仍將在較長時間內特別是我國能源轉型發展中發揮著兜底保障作用。然而，隨著煤炭資源的不斷開采，礦區滑坡、巖體坍塌、地表塌陷等一系列地質災害問題頻發，嚴重威脅周邊人民的生命安全和財產安全[1]，如2022 年2 月22 日，內蒙古新井煤業露天煤礦重大坍塌事故造成53 人死亡，直接經濟損失達20 430.25 萬元。因此，針對煤礦開展長時序地表形變高精度監測，分析其地表塌陷范圍、時序變化趨勢等對礦區地質災害的預警及防控具有重要的參考價值[2]。

隨著國內外星載合成孔徑雷達（Synthetic Aperture Radar，SAR）系統的發展，干涉SAR（Interferometric SAR，InSAR）技術憑借其全天時、全天候、廣域、高精度的優勢在礦區形變監測領域展現出卓越的應用前景[3]，有效彌補了傳統水準及導航定位測量點狀監測、耗時耗力的短板。尤其是后續發展的多時相InSAR 技術，其通過開展長時序高相干點目標監測，減少了差分InSAR 技術受時空失相干、大氣延遲等干擾因素的影響，提高了形變監測的精度及可靠性[4-5]。LEE 等[6]采用永久性散射體InSAR（Persistent Scatterers InSAR，PSInSAR）技術獲取了韓國最大露天石灰石礦的升軌、降軌邊坡形變速率場，并融合分析了礦區邊坡沉降變化規律，為露天礦區邊坡滑坡預警提供重要參考；PAWLUSZEK 等[7]有效融合了PSInSAR 技術與DInSAR 技術的優勢，進一步減少了大氣延遲誤差的影響，并結合DInSAR 技術獲取了大梯度形變處的形變速率場。然而，由于PSInSAR 技術采用單主影像干涉對構網，集中于選取長時間序列內始終保持高相干性的目標點開展監測分析，致使監測點密度嚴重受限[8]，針對于此，基于多主影像、短時空基線構網的小基線集（Small Baseline Subset InSAR，SBAS）技術得到進一步發展及應用[9]。劉澤洲等[10]以甘肅省、青海省交界處窯街煤礦為例論證分析了SBAS 礦區形變監測的可靠性及有效性，并指出SBAS 技術在礦區地表沉降監測與分析方面具有良好的應用前景；張連蓬等[11]、鄭美楠等[12]分別利用SBAS 技術對徐州市東部礦區及淮南礦區的關閉礦井開展了長時序地表形變監測與分析，獲取了其對應關閉礦井后的地表次生沉陷及閉礦后地表多維變形時空演化規律。為了進一步提高SBAS 技術的形變監測精度，柴華彬等[13]將SBAS 技術水準數據融合，提出融合反距離差值的水準數據SBAS-InSAR 技術，進一步提高了礦區沉降監測的監測精度及可靠性；丁寧等[14]將SBAS 技術與GPS 數據模型相融合，在提高監測精度的基礎上進一步有效反演出了龍首礦露天礦坑的三維形變場信息。盡管SBAS 技術較PSInSAR 技術在監測點密度上得到部分提升，但針對農田、植被覆蓋率較高的礦區環境，其有效監測點仍不能滿足礦區精細化形變監測的需求，可提供豐富監測點密度的分布式散射體InSAR（Distributed Scatterers InSAR，DSInSAR）技術得到廣泛的應用與發展[15]。李毅等[16]將DSInSAR技術應用于沛北礦區監測，有效驗證了DSInSAR 技術在礦區監測方面的有效性及優越性；LI 等[17]通過將空間域濾波與時間維濾波融合處理，有效改善了相位優化性能，進一步提高了DSInSAR 技術在礦區大梯度處的監測效果。盡管國內外學者針對DSInSAR技術在礦區形變監測方面的應用開展了大量研究[18]，但針對農田等非人工地物覆蓋為主的霄云煤礦（圖1）尚未開展精細化的地表形變監測研究。

圖1 研究區概況Fig.1 General situation of the study area

本文在現有研究成果的基礎上，提出基于魯棒M 估計協方差矩陣特征值分解相位優化的DSInSAR技術，并對霄云煤礦開展地表精細化形變監測研究，對其地表沉陷規律開展進一步分析，既拓寬了DSInSAR 技術在復雜礦區環境下的形變監測應用范疇，也進一步揭示了霄云煤礦的長時序地表沉陷規律。

1 研究區概況及實驗數據源

1.1 研究區概況

霄云煤礦位于濟寧市金鄉縣霄云鎮境內，煤礦井田東西長8.1 km，南北寬2.9 km，面積23.42 km2，東北與魚臺縣接壤，主井井口距金鄉縣城20 km，距魚臺縣城25 km，如圖1（a）所示。井田內地勢平坦，為黃淮沖積平原，但煤田內部構造復雜，褶曲不發育，斷層比較發育，局部地段揭露巖漿巖侵入體，以斷裂構造為主，且均為正斷層，并發育次一級波幅較小的波狀褶皺，斷裂構造的方向性明顯，為北傾的單斜構造。

霄云煤礦于2012 年8 月建成投產，地質儲量10 120.2 萬t，可采儲量4 791.5 萬t，其含煤地層主要為二疊系下統山西組及石炭系上統太原組，主要開采煤層為3 煤層，可采范圍內平均厚度3.49 m，傾角一般在10°～21°之間，開采上下限為-430～-1 500 m。全礦井3 煤層共劃分了4 個采區，其中，一采區位于霄云煤礦西南部，二采區位于霄云煤礦東部，三采區和四采區位于霄云煤礦西北部。

結合可持續發展大數據國際研究中心公布的全球30 m 精細土地覆蓋數據集[19]，對研究區所在研究時段（2022 年）的土地覆蓋類型分布情況進行簡要分析，如圖1（b）所示。研究區內主要包含耕地、草地、水體、裸地、建設用地五種地物類型，以耕地和建設用地為主，占比分別達到約82.0%和16.9%。霄云煤礦內耕地、建設用地的地物類型占比分別約82.1%、14.6%。研究區內較高的耕地地物占比嚴重限制了長時間序列下保持高相干性的PS 點識別，進而制約了形變信息的有效反演，因此，本文采用可提供豐富點密度的DSInSAR 技術開展實驗分析。

1.2 實驗數據源

本文實驗采用歐洲空間局向全球開源提供的Sentinel-1A SAR 衛星數據，選取覆蓋霄云煤礦研究區2021 年11 月至2022 年12 月的34 景數據，時間間隔為12 d，實驗數據為Sentinel-1A 衛星IW 成像模式、升軌、VV 極化方式。此外，實驗采用TanDEM-X 計劃中TerraSAR-X/TanDEM-X 編隊衛星生產的90 m分辨率WorldDEMTM產品為外部數字高程模型，以便更精確地去除地形相位。

鑒于霄云礦區耕地廣泛分布的地理環境，實驗采用最鄰近干涉對的構網方式，以避免長時間基線、空間基線引起的時空失相干因素影響，基于此構網的干涉對時空基線信息如圖2 所示。

圖2 SAR 影像時空基線組合Fig.2 Spatiotemporal baseline combination of SAR images

2 實驗方法及原理

不同于擁有主導散射體可保持長時間序列下高相干性的PS 點，DS 點分辨單元內并無主導散射體，致使其易受時空失相干因素的影響，因此，提升相位信噪比的相位優化處理是DSInSAR 技術的關鍵步驟。為保證相位優化在低相干區域展示出更好的優化效果，保障后續相位信息的精準解譯，采用基于魯棒M估計協方差矩陣的特征值分解相位優化方法。隨后，在獲取優化相位信息的基礎上，進一步開展相位解譯處理，首先，對優化后相位開展三維相位解纏處理；然后，依據相位成分的時空特征，通過時間維及空間維的高低通濾波，實現對解纏相位的相位干擾成分分離，進而獲取最終的形變信息。

2.1 基于魯棒M 估計協方差矩陣

假設SAR 圖像中任一像素x滿足經典的零均值復數圓高斯分布（Complex Circular Gaussian，CCG），則在該像素處獲取的時間序列數據滿足多元CCG 分布，其對應的概率密度函數可表示為式（1）。

式中：N為時間序列SAR 影像的個數；C 為協方差矩陣；det（）表示矩陣行列式運算。由式（1）可知其概率密度函數主要由其協方差矩陣表征，在給定具有相似統計特征像素集合 Ω的基礎上，協方差矩陣的最大似然估計量可表示為式（2）。

式中：M為集合 Ω中像素的個數；（）H為共軛轉置運算。在DSInSAR 技術的相位優化處理過程中，通常將上述作為樣本協方差矩陣開展后續處理。然而，該最大似然估計量對鄰域像素集中像素是否滿足CCG 同分布較為敏感，樣本中一個單一的異常值便會使估計量產生較大偏差，因此，該最大似然估計量對異常值不具備一定的魯棒性。

鑒于此，一些研究學者從同質像元的選取算法出發，以提高同質像素識別的準確性，減少異質像素的干擾。本文在先進研究成果的基礎上，采用較為先進的置信區間假設檢驗（Hypothesis Test Of Confidence Interval，HTCI）算法識別同質像素集[20]，其將局部方差最小的候選窗口內的強度均值作為參考像素的估值，進而確定準確的置信區間。在對變量x進行對數變換后y=lnx，其L視數下對應的置信區間可表示為式（3）。

式中：μy為對數瑞利分布對應概率密度函數的期望值；α為顯著性水平；z1-α/2為標準正態分布的 1-α/2百分位。通過上述置信區間，將鄰域窗口中所有樣本均值yˉ滿足該區間的連通像素識別為具有相似分布特征的同質像元，進而確定同質像素集 Ω。然而，按此概率法確定同質像元的方式難免會引入小概率的異質像素，進而影響式（2）最大似然估計量的可靠性。

為避免非高斯散射機制的異質像素或樣本選擇異常等引起的非魯棒估計，本文采用魯棒性較好的M 估計協方差矩陣[21]。針對時序SAR 像素的多元變量x，其協方差矩陣的M 估計量可表示為式（4）。

其可進一步通過迭代重加權的方式進行求解估計，其初始值可采用最大似然估計量，第k+1 次的迭代估計量可表示為式（5）。

式中，w（）為實值加權函數。由于SAR 影像中非高斯像元的分布可建模為復數t分布，其加權函數可進一步推導為式（6）。

為了保證整體估計具有較好的魯棒性，通選設置v＜5。為抑制嚴重的長尾分布的影響，可假設v=0，進而權值函數可等價表示為式（7）。

式中，‖‖為二范數運算。此種估計方法可以避免冗余的迭代運算，且對異常值具有較好的魯棒性。

2.2 基于魯棒M 估計協方差矩陣特征值分解相位優化

FORNARO 等[22]在CAESAR（Component extrAction and sElection SAR）技術框架中首次提出了基于特征值分解的相位優化策略，其主要通過對樣本協方差矩陣開展特征值分解，估計出多元散射機制對應的相位分量，通過提取最大特征值對應的特征向量來有效分離出主導散射機制，實現相位的優化估計處理。通過上述特征值分解相位優化方法的核心思想可知，樣本協方差矩陣是其方法的基礎及核心所在，因此，其協方差矩陣的非魯棒性、非平穩性估計將會對后續的優化相位估計量產生最直接的影響。鑒于此，本文采用上述提到的基于魯棒M 估計協方差矩陣。

基于魯棒M 估計協方差矩陣，其特征值分解可表示為式（9）。

式中：λ為矩陣分解的特征值；u 為特征值對應的特征向量。由式（9）分解后的每個特征值可理解為分辨單元內不同的散射機制，其對應的特征向量即為不同散射機制對應的相位分量。因此，令λ1＞λ2＞λ3＞···＞λN，最大特征值 λ1即為分辨單元內的主導散射機制，將其視為主導信號，其對應的特征向量u1即為相位的主要信息源，其余特征值即為分辨單元內的子散射機制，其對應的特征向量可視為噪聲信號。鑒于此，上述最大特征值對應的特征向量u1即為相位的優化估計量。

近些年，一些相位優化通用函數模型[23]被推導并得到大量的研究，上述基于魯棒M 估計協方差矩陣特征值分解的相位優化同樣也可以推導為通用函數模型的表達式。前述最大特征值對應特征向量的優化估計量可進一步表示為常規的最大化問題解算，即式（10）。

由上述推導可知，基于魯棒M 估計協方差矩陣特征值分解的相位優化策略同樣可表征為相位優化通用函數模型，僅權重wmn定義為即可。

經由上述處理獲取相位優化估計量后，進一步根據優化結果的時域相干性指標篩選DS 候選點，并將其對應的DS 候選點相位替換為優化后相位，同時將其與PS 點對應的原始相位融合，實現PS 候選點與DS 候選點的監測點融合處理，基于此，開展常規的相位穩定性分析、三維相位解纏等相位信息解譯處理。

3 監測結果分析

3.1 DSInSAR 監測結果

為了驗證基于魯棒M 估計協方差矩陣特征值分解相位優化的DSInSAR 形變監測效果，實驗采用常規SBAS 技術及基于常規最大似然估計量協方差矩陣特征值分解相位優化的DSInSAR 監測結果開展對比分析。針對兩種DSInSAR 技術中不同樣本協方差矩陣開展的特征值分解的相位優化結果展開分析，實驗采用時域相干性值作為定量化評價指標，定義為式（12）。

圖3 時域相干性值Fig.3 Temporal coherence value

由圖3 可知，基于魯棒M 估計協方差矩陣特征值分解對應相位優化結果的時域相干性值略優于基于最大似然估計協方差矩陣特征值分解對應的相位優化結果，其對應更多的高相干值像素，表明基于魯棒M 估計協方差矩陣特征值分解的相位優化方法較常規最大似然估計協方差矩陣特征值分解相位優化方法展現出更好的優化效果。本文實驗設置篩選DS 候選點的時域相干性閾值為0.4，基于魯棒M 估計協方差矩陣對應優化結果可篩選出4 142 065 個像素，像素占比約0.80%，而常規優化結果僅能篩選出3 624 099 個像素，像素占比約0.70%，對比提升了517 966 個DS 候選點，其更多的DS 候選點將更有助于后續獲取更高密度的監測點，更有利于形變場信息的豐富及精準解譯。

基于上述優化相位，分別對其開展相位信息解譯處理，基于魯棒M 估計協方差矩陣特征值分解相位優化的DSInSAR（DSInSAR-M）、基于常規最大似然估計協方差矩陣特征值分解相位優化的DSInSAR（DSInSAR-MLE）及SBAS 技術獲取的雷達視線向形變速率結果如圖4 所示。

由圖4 可知，即便在耕地地物分布為主的礦區復雜監測背景下，DSInSAR 技術依然能夠獲取較高密度的監測點信息，而SBAS 技術獲取的監測點主要覆蓋建設用地所對應的地物類別，在主要覆蓋耕地地物的研究區霄云煤礦內僅識別出極少的觀測點信息，極大地限制了該研究區形變場信息的有效反演。具體而言，SBAS 技術僅篩選出約206 819 個可靠監測目標點，DSInSAR-MLE 方法及DSInSAR-M 方法則分別篩選出2 159 739 個可靠監測目標點和2 563 576個可靠監測目標點，本文提出的DSInSAR-M 方法在監測點密度上較SBAS 技術及DSInSAR-MLE 方法分別提升了約11.4 倍、0.2 倍，進一步表明了本文技術較常規時序InSAR 技術在非建設用地區域監測的有效性。值得說明的是，高密度的監測點將更有助于解譯出豐富、可靠的形變場信息，尤其是圖4 局部放大所示的礦區地表沉降的中心區域。

結合上述監測結果可知，研究區內共包含四個顯著形變場。其中，有兩個顯著形變場分布在霄云煤礦內（圖4 上方局部顯示），兩種DSInSAR 技術估計的形變結果較為相似，僅在形變場的中心區域有較輕微的差異，獲取的局部監測點分別為19 626 個和25 506 個，其監測到的最大沉降速率均高達530 mm/a；而SBAS 技術在該局部區域僅識別出1 158 個監測點，監測到的最大沉降速率僅約272 mm/a，主要是因為沉降中心區域監測點信息較少，致使了形變速率呈現出上述形變低估及形變信息缺失的現象。另外兩個顯著形變場分布在單縣境內（圖4 左下方局部顯示），SBAS 技 術、DSInSAR-MLE 方法及DSInSAR-M 方法獲取的局部監測點分別為3 768 個、36 682 個、49 789 個，對應可監測到的最大沉降速率分別約為540 mm/a、517 mm/a、559 mm/a，其更高的監測點密度將有助于后續三維相位解纏獲取更準確的解纏相位。

3.2 形變結果分析

基于DSInSAR-M 獲取的形變監測結果，實驗提取了平均間隔36 d 的霄云煤礦雷達視線向形變時間序列結果，如圖5 所示。由圖5 可知，霄云煤礦在實驗研究時段共包含三個沉降形變場，即形變場A、形變場B、形變場C（標記于圖5（l）），其中，形變場A、形變場B 在圖4 形變速率結果中已被識別，而形變場C 因形變速率較小而并未能有效識別，其主要在2022 年7 月之后逐漸開始產生緩慢形變。實驗結合各形變場中心處近似最大沉降點對應的時序形變量開展變化趨勢分析，如圖6 所示。形變場A 從研究時段2021 年11 月7 日起，地表便以近似顯著的線性大形變速率趨勢發生沉降變化，直至2022 年6 月23 日，228 d 內形變場內特征點處產生了約368 mm的沉降量，隨后地表開始趨于緩慢沉降，在后續的180 d 內僅發生了約50 mm 的沉降量，整個研究時段累計產生了最大約418 mm 的沉降量。形變場B 在2021 年11 月7 日 至2022 年5 月18 日，地表變化較為顯著，下沉值達到174 mm，隨后地表沉降趨勢開始逐漸平穩，最大沉降量約234 mm。形變場C 在2022年7 月之前，地表沉降變化非常微弱，下沉值僅約50 mm，而后地表開始產生快速的沉降變化，短短168 d 下沉值便高達約240 mm，最大沉降量達到293 mm，且其形變場范圍開始逐漸擴大。

圖5 形變時間序列Fig.5 Time series of deformation

圖6 最大沉降點處形變時間序列Fig.6 Time series of deformation in maximum subsidence points

實驗分別以局部形變場A 及形變場C 為例，結合其工作面分布開展形變特征分析，局部結果如圖7 所示。在實驗研究時段，局部形變場A 區域僅1315 工作面處于開采階段，其余早期已開采工作面距該時段間隔較久，對應產生的地表殘余變形相對較少，即形變場A 主要受1315 工作面開采影響，該形變場以工作面為中心呈現出近似圓形的形變分布特征，且形變范圍較大。實驗分別沿1315 工作面的走向及切線方向劃定兩個剖面線，如圖7（a）中黑線所示，并提取出剖面線對應時序累積形變量開展分析，如圖7（b）和圖7（c）所示，其形變量均展現出顯著的拋物線趨勢，且隨著日期的變化，均呈現出先快速沉降后緩慢沉降的時序變化趨勢，近似以2022 年6月23 日為分段日期，這與上述最大沉降點對應時序變化分析的結論基本保持一致，其走向及切線向剖面線的最大沉降量分別約為372 mm 和333 mm，結合上述最大特征點的時序變化量可知，實驗劃定的兩處剖面線并未涉及該形變場處的最大沉降點。局部形變場C 主要受1310 工作面和1314 工作面這兩個正在開采的工作面影響，由于兩個開采工作面的走向長度不同，其形變場呈現出不規則的分布特征，結合工作面對應走向剖面線的時序累積形變量分析，在2022 年7 月29 日前兩個剖面線均以較為緩慢的趨勢產生地表沉降，在此之后，均以較大的形變速率開始迅速沉降，1310 工作面剖面線的沉降變化程度比1314 工作面剖面線更加迅速，兩個剖面線的最大沉降量分別達到了286 mm 和153 mm，其最大沉降量的差異主要是由于1314 工作面的開采時間晚于1310工作面所致，與前述特征點時序變化趨勢分析的結論保持一致。

圖7 局部形變特征Fig.7 The local deformation characteristics

3.3 監測精度分析

為進一步驗證本文提出的DSInSAR-M 方法的可靠性，采用圖1（a）所示的23 個水準點開展實驗對比分析。通過插值處理實現InSAR 監測結果與水準數據的時空基準統一，設置水準點對應小范圍緩沖區內的InSAR 監測點為同名點，SBAS 技術、DSInSAR-MLE 方法及DSInSAR-M 方法對應累積沉降量與水準數據的對比結果如圖8 所示。

圖8 InSAR 結果與水準數據對比Fig.8 Comparison between InSAR results and leveling data

由圖8 可知，在形變量較小的點處，三種方法對應的InSAR 方法監測結果與水準數據均具有相對較高的一致性，但在大梯度形變處，SBAS 技術基本無有效同名點，即使存在同名點也與水準數據保持著顯著的差異，最大誤差高達218 mm，其主要是由于稀疏的監測點信息增加了后續相位解纏的難度，嚴重制約了形變信息的解譯精度；兩種DSInSAR 方法監測結果雖與水準數據保持著相似的形變趨勢，但在大梯度形變處也存在輕微的形變低估現象，具體而言，DSInSAR-MLE 方法與DSInSAR-M 方法分別識別出了20 個水準同名點、22 個水準同名點，其最大誤差分別為81 mm、52 mm，SBAS 技術、DSInSAR-MLE方法、DSInSAR-M 方法監測結果與水準數據的絕對平均誤差分別為51 mm、25 mm、17 mm，對應的均方根誤差分別為84 mm、32 mm、22 mm，表明本文提出的DSInSAR-M 方法較傳統監測方法具有相對最高的監測精度，進而證實了本文提出的DSInSAR-M 方法的有效性及可靠性。

4 結論

針對以耕地等自然地表地物為主的霄云煤礦，為克服復雜監測環境對InSAR 形變解譯的影響，本文采用基于魯棒M 估計協方差矩陣特征值分解的DSInSAR 技術，對霄云煤礦2022 年34 景Sentinel-1A影像開展形變監測處理和分析。所得結論如下所述。

1）與常規的SBAS 技術、基于最大似然估計矩陣特征值分解的DSInSAR 技術相比，本文采用的基于魯棒M 估計協方差矩陣特征值分解的DSInSAR技術可以獲取更高密度的監測點，較DSInSAR-MLE方法提升了約0.2 倍，較SBAS 技術則提升了約11.4 倍，有效提升了形變場信息反演的豐富度。

2）在實驗研究時段，霄云煤礦內主要包含3 個顯著的形變場，其中西南部的兩個形變場均呈現出先快速沉降后趨于平穩的變化趨勢，東部形變場呈現出先緩慢變化后急劇沉降的變化趨勢，研究時段內累積形變量最大的為最西部形變場，其最大下沉值高達418 mm。

3）通過與水準數據對比分析，InSAR 技術監測結果與水準數據在低沉降量對應點處具有較高的一致性，但在高沉降量對應點處均呈現出形變低估現象，其中本文提出的DSInSAR-M 方法對應均方根誤差最小，約22 mm，較常規DSInSAR-MLE 方法及SBAS技術具有更高的監測精度。