SBAS-InSAR在礦區地表形變監測中的應用與分析

王道順，李偉*，孫建，鄭國棟

(1.山東省魯南地質工程勘察院(山東省地質礦產勘查開發局第二地質大隊)，山東 兗州 272100；2.高分辨率對地觀測系統濟寧數據與應用中心，山東 兗州 272100；3.自然資源部采煤沉陷區綜合治理工程技術創新中心，山東 兗州 272100)

0 引言

煤炭資源作為我國的重要能源之一，極大地推動了國家的經濟和社會發展。隨著開采強度的日益增加，由此引起的生態環境問題更加凸顯[1-3]。煤炭開采造成的地表沉陷屬于典型地質災害之一，往往導致礦區及一定范圍內的地表產生不同程度的沉降，對地面房屋設施、道路、農田水利基礎設施等建(構)筑物造成不同程度的損毀，為經濟發展和社會穩定帶來重大隱患，因此對煤炭開采引起的地表沉降進行監測具有重要意義。

以往的礦區地表沉降監測多采用布設巖移觀測站的方式，利用水準測量、GNSS測量等傳統方法進行監測，此類方法監測精度可觀，但作業成本較高，獲取的監測數據量有限，不能完全反映礦區地表沉降特征，難以實現大范圍的地表沉降監測。合成孔徑雷達差分干涉測量(D-InSAR)通過微波遙感主動測量的方式，可以克服傳統監測方法的不足，且具有全天時、全天候作業的優點，可獲取大范圍、高精度的空間延續性地表形變信息。已有學者將其應用于礦區的地表形變監測研究[4-9]，但該方法使用數據量較少，在長時間序列的監測過程中容易受到時空失相干、大氣延遲等因素的影響，降低其監測精度。在此基礎上發展而來的時序InSAR技術彌補了上述不足，極大地提高了形變監測精度，已廣泛應用于礦區形變監測、城市沉降監測等方面，其中以永久散射體技術(Permanent Scatterer InSAR，PS-InSAR)和短基線集方法(Small Baseline Subset InSAR，SBAS-InSAR)應用最為普遍[10-15]。

PS-InSAR技術已成功應用于礦區的地表沉降監測，為了獲得在時間序列上散射特性穩定的PS點，PS-InSAR往往需要20景以上的數據量參與計算。該方法基于大量PS點的迭代回歸計算，導致運算效率偏低[16]。礦區多位于郊區、農田等低相干地區，若所用數據量偏少、時間間隔較長，利用PS-InSAR技術難以獲取足夠的PS點，導致無法完整揭示礦區沉降規律。本文在分析SBAS-InSAR技術方法的基礎上，以濟寧北部礦區為研究區域(圖1)，采用C波段Sentinel-1A雷達數據和SBAS-InSAR時序方法開展礦區地表沉降監測和時序分析，獲取監測范圍內的地面沉降時空分布特征，可為采煤沉陷區的調查和治理提供技術支撐。

圖1 研究區域影像圖

1 SBAS-InSAR原理

SBAS-InSAR由Berardino等人[17]于2002年提出，不同于PS-InSAR的單主影像，該方法是一種基于多主影像的InSAR時間序列方法。它通過短基線原則，將大量SAR數據組合為具有多個主影像的干涉子集，每個子集內的干涉對基線長度均低于臨界基線值，時間基線也盡可能短，集合間的SAR影像基線距大，通過這種方式克服了時間和空間上的失相關，因此SBAS-InSAR可以通過較少的數據量來獲取較可靠的監測結果，且其監測對象為分布式散射體，該方法更適合于礦區的地表形變監測。

基于短基線集原則，SBAS-InSAR具有多個主影像，但仍需選取其中一景影像作為公共主影像進行配準；組成各干涉子集后，利用外部參考數字高程模型(Digital Elevation Model，DEM)數據模擬并去除每個干涉相對的地形相位，然后生成時間序列差分干涉圖集；相位解纏后得到每個相干目標的相位信息，包括形變相位、大氣延遲相位、軌道誤差相位等信息，各誤差相位可在時間序列上采用濾波方法或多項式模型予以去除；由于SBAS-InSAR具有多個主影像，各干涉子集在聯合求解時容易出現方程秩虧現象，因此引入奇異值分解方法，利用最小二乘原理得到地表時間序列形變信息。

2 研究區域與數據源

2.1 研究區域

濟寧市位于山東省魯西南地區，礦產資源豐富，是魯西煤田基地的重要組成部分，是全國重要的煤炭基地之一[18-20]。本次實驗研究區域位于濟寧市北部地區，覆蓋范圍大致為E116°28′01″～116°44′34″、N35°29′40″～35°42′35″，研究區域內礦區集中分布，主要有義橋煤礦、唐陽煤礦、義能煤礦、魯西煤礦、新驛煤礦等礦區(圖1)。

2.2 數據源

本文使用覆蓋研究區域2019年9月19日—2020年3月29日間的17景Sentinel-1A C波段單視復數(SLC)數據，該數據為干涉寬幅(IW)成像模式，極化方式VV，空間分辨率5m×20m(方位向×距離向)，具體信息見表1。Sentinel-1A衛星在軌運行高度693km，重訪周期12d，若與Sentinel-1B衛星組網觀測，可將重訪周期縮短至6d，可有效提高監測區域的時間相干性，此外由于Sentinel-1A衛星采用先進的軌道控制技術，其干涉對空間基線最大為162m，多數都在100m范圍內[21]，遠遠小于其空間臨界基線值，因此可以忽略其干涉對空間基線的差異。外部參考DEM選擇由日本太空發展署(JAXA)免費分發的ALOS World 3D數據，空間分辨率30m，用于模擬去除地形相位和地理編碼。

表1 研究區Sentinel-1A影像參數

3 礦區形變監測與分析

3.1 短基線集干涉組合

實驗選擇2019年11月30日的影像作為公共主影像，其余影像作為從影像與主影像進行配準并采樣到主影像像素空間。設置空間基線閾值100m，時間基線閾值60d，共得到55組干涉相對，其中最大空間基線87.47m，最小空間基線4.11m，時空基線連接情況見圖2。

圖2 時空基線圖

3.2 干涉圖生成與相位解纏

將得到的55組像對進行差分干涉處理，得到時間序列干涉圖集，同時生成對應的相干系數圖，由于設置的閾值相對嚴格，多數差分干涉圖干涉效果較好，可以看到明顯的形變干涉條紋。圖3是2019年11月30日—2019年12月24日內生成的差分干涉結果(SAR坐標系)，圖3a為相干系數圖，顏色越亮代表相干性越高；圖3b為差分干涉圖，可以看出，監測區域內有多處明顯的形變干涉條紋。相位解纏選擇最小費用流(MCF)解纏算法，相干性閾值設為0.2，低于相干閾值的像元進行掩膜處理，不參與計算，然后逐個干涉圖進行相位解纏，圖3c為對應的相位解纏圖，黑色代表相干性過低而沒有解纏結果的地方。

3.3 軌道精煉和重去平

對得到的55組干涉像對結果進行目視檢查，剔除相干性過低或解纏錯誤的像對。選取一組干涉效果良好的像對作為參考影像，選擇遠離形變區域且相干性較高的若干穩定地面點作為控制點，使用多項式模型進行軌道精煉，去除殘余的恒定相位和相位坡道。

3.4 獲取監測結果

上述步驟完成后反演區域形變速率，同時對干涉圖進行二次解纏優化。在第一次形變速率反演的基礎上，利用大氣相位在時間上的高通濾波和空間上的低通濾波進行去除，獲取時間序列上的位移結果(圖4)。

a—相干系數圖；b—差分干涉圖；c—解纏相位圖圖3 干涉對相干系數圖、差分干涉圖、解纏相位圖

圖4 監測區域年均沉降速率圖

3.5 結果分析

從監測結果可以看出，監測區域整體基本處于穩定狀態，濟寧北部礦區內明顯分布有6處不同沉降程度的區域(A—F)，分別位于義橋煤礦、唐陽煤礦、義能煤礦、新驛煤礦、魯西煤礦、何崗煤礦，其中C區域沉降最大，年平均沉降速率達到-242mm/y，呈現出雙沉降漏斗特征；A、B兩處沉降區域相距僅640m，已呈現接連成片沉降趨勢，沉降范圍最大；D區域最大沉降速率達到-85.7mm/y；E，F兩處區域沉降范圍最小，E區域最大沉降速率為-70.1mm/y，F區域最大沉降速率達到-91.6mm/y。部分沉降區域中心不同程度的出現了“監測空值”，主要是礦區沉降過大，超出了InSAR可探測的最大形變梯度，致使沉降中心出現失相干，但其沉降邊緣范圍探測不受影響。

4 結論

(1)基于2019年9月19日—2020年3月29日的17景Sentinel-1A數據，運用SBAS-InSAR技術獲取了濟寧北部礦區的年均沉降速率和沉降漏斗分布。

(2)濟寧北部礦區內分布有6處明顯沉降區域，分布于義橋煤礦、唐陽煤礦、義能煤礦、新驛煤礦、魯西煤礦和何崗煤礦，形變范圍與礦區實際空間分布位置一致，其中義能煤礦沉降特征最為明顯，年均沉降速率達到-242mm/y。

(3)受形變梯度影響，部分沉降區域中心出現“監測空值”，但沉降范圍仍可以完整探測。相比較于C波段數據，L波段數據波長更長，在低相干區域更容易穿透地表植被覆蓋，得到更高質量的干涉結果和更大中心沉降量。下一步工作將收集該區域L波段數據，開展更高質量的InSAR礦區地表形變監測，為礦區的地面沉降防治提供決策依據。