D-InSAR 技術與SBAS-InSAR 技術在礦區地面沉降監測中的應用研究

王 鵬，程海強，王 翔

（1.山東省三河口礦業有限責任公司地質測量科，山東 濟寧 277600；2.棗莊礦業（集團）有限責任公司地質測量部，山東 棗莊 277000；3.山東建筑大學測繪地理信息學院，山東 濟南 250101）

0 引言

煤炭資源作為我國儲量豐富的礦產資源，在推動我國經濟快速發展的同時，也導致礦區范圍內的地面發生不同程度的沉降，對其上覆房屋、橋梁以及水利設施等建（構）筑物造成不同程度的損毀，因此，對礦區沉降情況進行監測具有重要意義[1-2]。傳統地面沉降監測是利用水準測量、GNSS 測量等方法進行，不僅工作量大、成本高，而且無法獲取大范圍、全沉降盆地的沉降信息[3]。InSAR 技術是利用主動微波成像技術獲取覆蓋研究區重復軌道的SAR 影像，結合差分干涉測量技術獲取地面沉降信息，具有全天時、全天候、大范圍、高精度等優點，克服了傳統監測方法的不足[4]。但D-InSAR 技術監測精度受時空去相干和大氣延遲等因素影響，進行長時間序列地面沉降監測時并不能提供可靠的沉降信息。為克服相關問題，InSAR 技術通過對長時間跨度內多景SAR 影像聯合處理，在提高地面沉降監測精度的同時，提供長時間序列的地面沉降監測結果，使得該技術在礦區沉降監測中得到廣泛應用[5-6]。目前，已有不少學者基于不同的InSAR技術對礦區沉降情況進行了監測與分析[7]。冷紅偉等[8]采用一種融合多視和全分辨率的D-InSAR 方法對濟寧礦區下沉盆地進行了沉降監測，結果表明2017 年2 月17 日—3 月11 日期間，研究區內最大沉降出現在濟寧市新驛煤礦，最大沉降量達到-74 mm，最大沉降速率達到-3.36 mm/d。張樹衡等[9]集成D-InSAR 技術和SBAS-InSAR 技術對淮南市丁集煤礦進行了地面沉降監測，有效解決了SBAS-InSAR 技術在礦區沉降監測中沉降漏斗邊緣精度低和部分區域像元空洞的問題，獲得了連續高精度的礦區地面沉降結果。

為驗證D-InSAR 技術和SBAS-InSAR 技術在礦區地面沉降中的監測能力，本文以濟寧市某礦區為研究區域，收集了2022 年11 月2 日—2023 年8 月29 日期間共26 景Sentinel-1A 影像，分別采用D-InSAR 技術和SBAS-InSAR 技術進行數據處理，獲取了該時段內研究區地面沉降監測結果，并結合研究區開采工作面內的實測水準數據對其精度進行驗證，對比分析了兩種技術在礦區地面沉降監測中的應用能力。

1 原理與方法

1.1 D-InSAR 技術基本原理

D-InSAR 技術通過兩次重復軌道觀測獲取同一地區不同時刻的兩幅SAR 影像，結合外部DEM 數據去除干涉圖中的地形相位而得到差分干涉圖，進而獲取研究區范圍內地面沉降監測結果情況[10]，主要原理見式（1）。

式中：φflat為平地相位；φtop為地形引起的相位；φdef為視線方向的形變引起的相位；φatm為大氣延遲相位；φnoise為噪聲相位；k為整周模糊度。其中，φflat可根據基線進行估算[11]，φtop可根據已有的DEM 估算，若忽略大氣相位和噪聲相位的影響，通過將估算出的φflat和 φtop去除，則可估算出地面沉降后的干涉相位 φdef，見式（2）。

式中：R1為主影像斜距；R2為輔影像斜距；Δr為沉降導致的地面點位移；等式右側第一部分-（R1-R2）為形變前的干涉相位；第二部分Δr為由地面沉降引起的形變相位。

1.2 SBAS-InSAR 技術基本原理

SBAS-InSAR 技術通過設定時間-空間基線閾值，將滿足閾值要求的時間序列SAR 影像數據分為多個具有較短時空基線的相互獨立的干涉子集，以保證干涉對的相干性，進而得到更準確的地面沉降監測結果[12]。

假設SAR 衛星在研究區內獲得了不同時段的N+1 幅SAR 圖像，利用設定的時間和空間閾值可以獲得k幅差分干涉圖，則第k幅影像在方位-距離像元坐標系（x,r）中的差分干涉相位 δφk（x，y）見式（3）。

式中：φ（tb,x,r）、φ（ta,x,r）為tb時刻、ta時刻的相位；d（ta,x,r）、d（tb,x,r）為tb時刻、ta時刻視線向累計沉降量。去除地形殘余相位[13]、大氣延遲相位[14]以及各種噪聲相位[15]后，地面沉降的平均速率vt可以表示為式（4）。

則相位為式（5）。

式中：EJ為主影像獲取時間；SJ為從影像獲取時間；vk為k時刻像元形變速率。

由此，定義A（j,k）=tk-tk-1，且A是一個m×n的秩虧矩陣，得到矩陣方程，見式（6）。

通過奇異值分解法和最小二乘法[16]可求出平均沉降速率相位值，進而可以獲取地面沉降量。本文數據處理流程如圖1 所示。

圖1 數據處理流程圖Fig.1 Flow chart of data processing

2 研究區概況及數據源

濟寧市位于魯西南地區，處在黃淮海平原與魯中南山地交接地帶，地理坐標介于東經115°54′～117°06′，北緯34°25′～35°55′，地勢東高西低，地形以平原洼地為主。濟寧市礦產資源極其豐富，已發現和探明的礦產儲量達70 多種，以煤炭為主，含煤面積占全市面積的45%，經探測，濟寧市煤炭儲量達260 億t，占山東省煤炭總儲量的50%，主要分布在兗州區、曲阜市和微山縣等地區，是全國重點建設的14 個億噸級大型煤炭基地之一[17]。本文研究區地理位置如圖2 所示，該采礦工作面所處區域周邊斷層、褶曲一般發育，巖層產狀變化不大，地質構造條件中等。工作面賦存標高-244.8～-308.3 m，平均采厚5.1 m；煤層傾角4°～12°，平均傾角7°左右；走向長1 070～1 120 m，傾向長102～226 m，面積22.69 萬m2。工作面采煤方法為傾向長臂后退式采煤法，頂板管理方法為全部冒落法。

圖2 研究區地理位置圖Fig.2 Geographic location of the study area

本文選取覆蓋研究區的26 景VV 極化Sentinel-1A數據開展實驗分析，SAR 影像時間跨度為2022 年11 月2 日到2023 年8 月29 日，具體參數見表1。為去除地形效應的影響，本文收集了覆蓋研究區范圍的SRTM DEM 數據，該產品為美國航天航空局“航天飛機雷達地形測繪任務”獲取的產品數據，空間分辨率為30 m，具有高精度、高覆蓋率等優點。此外，為了準確獲取SAR 影像軌道信息，本文收集了可以提供高精度的衛星位置、衛星速度信息及衛星姿態數據等信息的POD 軌道數據，該數據是目前常用的SAR 影像定軌數據，可用來處理數據中矯正衛星運動誤差和大氣延遲誤差，精密軌道數據日期與Sentinel-1A 數據日期一一對應。

表1 Sentinel-1A 數據IW 成像模式基本參數Table 1 Basic parameters of IW imaging mode of Sentinel-1A data

3 數據處理

3.1 D-InSAR 技術

在D-InSAR 技術數據處理中，每間隔約36 d 選擇一幅影像，在2022 年11 月2 日—2023 年8 月29日期間內共有9 幅影像。將相鄰的兩幅影像組合為干涉對，依次對所獲得的8 個干涉對進行干涉處理。SAR 影像配準直接影響干涉圖的質量，本文先基于軌道信息對影像進行粗配準，再基于頻譜差異法進行精配準，然后通過不斷迭代，最終達到0.001 像元的配準精度[18]。研究區地表主要被魚塘、農田和村莊覆蓋，當雷達波束照射到地面時，在表面光滑的水體區域會產生類似鏡面反射的現象，而在植被覆蓋區域會受到植被結構的影響發生多次的散射和反射，導致這些地物特征區域的目標后向散射系數較低，回波信號弱，相干性差且干涉相位噪聲嚴重。為了抑制這些噪聲對干涉相位的影響，以10∶2 的比例對SAR 影像進行多視處理[19]，并采用自適應濾波方法[20]對差分干涉圖進行濾波處理。進一步通過設定相干性閾值，掩膜掉部分相干性較差的區域（主要為水體區域）。在這一過程中，雖然礦區沉降中心的大幅度地表形變也導致該區域的相干性降低，但是為了獲取全面的礦區形變結果，沒有對這一主要研究區域進行掩膜。采用最小費用流法[21]進行相位解纏，獲取解纏后的差分干涉相位圖并將雷達視線方向的解纏相位轉換為形變量，通過地理編碼將獲取的形變結果轉換到地理坐標系下，將8 個干涉對獲取的形變結果導入ArcGIS 軟件中，用柵格計算器將其進行疊加，獲得相應時間段內的累積沉降量。D-InSAR技術獲取的截至2023 年8 月29 日的累積沉降結果如圖3 所示。

圖3 D-InSAR 累積沉降圖Fig.3 Cumulative subsidence results from D-InSAR

3.2 SBAS-InSAR 技術

選 取在2022 年11 月2 日—2023 年8 月9 日 期間內的共26 幅SAR 影像，時間間隔12 d。設定時空基線閾值，空間基線閾值最大為310 m，時間基線閾值為36 d，進行干涉對組合，共生成47 個干涉對。在SBAS-InSAR 技術數據處理過程中，與D-InSAR 技術對解纏相位的疊加計算不同，SBAS-InSAR 的關鍵步驟是對解纏相位進行時間序列分析，從干涉相位中剔除高程殘差和大氣相位。首先，需要提取每個高相干點對應的解纏相位、圖像信息和高程等信息；然后，利用SVD 方法進行小基線集分析，解算出高程改正量和形變序列；最后，采用時空域濾波的方法從形變序列相位中分離出大氣延遲相位，即可獲得最終的時間序列形變相位。SBAS-InSAR 技術獲取的累積沉降結果如圖4 所示。

圖4 SBAS-InSAR 累積沉降圖Fig.4 Cumulative subsidence results from SBAS-InSAR

4 結果分析

4.1 沉降監測結果對比分析

為了對兩種技術在礦區沉降監測應用中的可靠性和監測能力進行驗證分析。以工作面區域為研究對象（圖2 中右圖框線所示），以D-InSAR 技術數據處理差分干涉時間為準，獲取SBAS-InSAR 技術對應時間的累積沉降量，結果如圖5 所示（每組子圖中左側為D-InSAR 結果，右側為SBAS-InSAR 結果）。由圖5 可知，沉降主要集中在該采礦工作面北側，東經117°04′20″、北緯34°51′50″周圍，初期2022 年12月20 日沉降量較小且均勻分散，隨著采礦作業的進行，沉降量不斷增大，沉降中心逐步顯現，至2023 年8 月29 日，從沉降中心向周圍擴散且沉降量逐漸減少的沉降漏斗基本形成。在這一過程中，D-InSAR 技術和SBAS-InSAR 技術的監測結果在沉降中心分布及沉降范圍方面吻合很好。但是受限于D-InSAR 技術分期獨立處理數據，并通過簡單累加的方式獲取累積沉降量的方法，每個獨立干涉對所包含的誤差都會體現在最終的累積沉降結果中，致使其累積沉降結果中存在較多誤差，在圖5 中體現為沉降區域的空間連續性較差。而SBAS-InSAR 技術得益于其時間序列處理能夠較好去除大氣延遲等誤差的特點，得到的累積沉降結果受到的誤差影響更小，在空間上的連續性更好，這有益于沉降區域的讀取及其對時空演變特征的分析。

進一步對兩種技術獲取的研究區地面沉降信息進行統計分析（圖6 和圖7）。由圖6 可知，當監測時間較短，小于半年時，SBAS-InSAR 技術監測得到的最大累積沉降量基本大于D-InSAR 技術監測得到的結果；當監測時間更長時，D-InSAR 技術監測得到的最大累積沉降量明顯超過了SBAS-InSAR 方法監測到的結果。截至2023 年8 月29 日，D-InSAR 技術監測得到的研究區內最大累積沉降量為69 mm，SBASInSAR 技術監測得到的最大沉降量為59 mm。由圖7 可知，在監測初期（2022 年12 月20 日），研究區內沒有監測到超過10 mm 的沉降量，在隨后的9 個月中，沉降超過10 mm 的區域面積不斷增大，截至2023 年8 月29 日，D-InSAR 技術和SBAS-InSAR 技術監測到累積沉降量超過10 mm 的面積分別為0.60 km2和0.87 km2，其中，沉降量大于30 mm 的面積分別為0.10 km2和0.17 km2。在整個監測過程中，SBAS-InSAR 技術識別的沉降量大于10 mm 的區域面積隨時間推移而穩定增大，而D-InSAR 技術卻存在一定的波動。此外，SBAS-InSAR 技術監測結果中，受植被和水體等地物影響相干性降低而未識別的區域面積為0.13 km2，而D-InSAR 技術則由開始的0.22 km2逐漸增加至0.32 km2。這是由于SBAS-InSAR 技術數據處理中通過小基線集技術獲取了更多不同時空基線的干涉對，在后續的時間序列解算中這些在時間維度上有一定重疊的干涉對中的相位信息能夠相互補充，彌補了D-InSAR 技術中某一時間間隔內僅依靠單個干涉對獲取地表形變相位的不足，得到了更豐富的地面沉降信息。

圖6 研究區最大累積沉降量柱狀圖Fig.6 Histogram of maximum cumulative subsidence in the study area

圖7 研究區各階段沉降面積變化柱狀圖Fig.7 Histogram of subsidence area change in the study area for different stages

4.2 沉降監測結果精度驗證

為了驗證兩種技術獲取的研究區內累積沉降量精度，本文對研究區開采情況進行了調查。該工作面于2022 年11 月開始開采，并在2022 年11 月4 日—2023 年8 月27 日期間進行了8 次水準測量，獲取了整個開采過程中主要區域的地面沉降情況。結合精密水準數據，分別對兩種技術獲取的地面沉降結果進行精度分析，繪制了8 個水準點（圖5 中依次為B01、B04、B08、B12、B16、B20、B24、B28），由D-InSAR、SBAS-InSAR 和水準測量三種方法獲取的累積形變量折線圖，如圖8 所示。

圖8 水準點累積沉降量折線圖Fig.8 Line graph of cumulative subsidence at the level

由圖8 可知，B01 號水準點、B04 號水準點、B08號水準點、B12 號水準點和B28 號水準點位于沉降漏斗邊緣，監測期間內最大累積沉降量不超過50 mm。兩種技術監測到的地面沉降趨勢與水準監測結果總體吻合，其中，SBAS-InSAR 技術監測得到的沉降曲線在不同時刻變化相對平緩，但在沉降量級較大區域的變化特征表現不明顯；D-InSAR 技術監測得到的沉降曲線變化幅度較大，但其對大量級沉降的探測能力比SBAS-InSAR 技術更敏感。這一現象與兩種不同技術形變量獲取的方法有關，其中，SBAS-InSAR技術是通過時間序列解算方法去除大氣效應等誤差影響，提高地面沉降解算結果的精度；而D-InSAR技術的疊加解算方法能夠更為明顯地反映真實地面沉降的形變特征，但在數據處理過程中往往無法消除如大氣效應等誤差的影響，進而影響其形變監測精度[22]。B16 號水準點、B20 號水準點和B24 號水準點位于沉降漏斗中心，監測期間內最大累積沉降量均超過300 mm，最大達到1 822 mm。根據InSAR 技術可監測形變梯度理論[23-24]，受雷達數據分辨率和波長等因素的影響，當地面沉降量過大時（超過InSAR技術可監測形變梯度范圍），兩種技術獲取的沉降監測結果不可靠。

5 結論

本文基于2022 年11 月2 日—2023 年8 月29 日期間的26 景Sentinel-1 影像數據，分別通過D-InSAR和SBAS-InSAR 兩種技術獲取了該階段內濟寧市某礦區的地面沉降信息，并結合煤礦開采工作面的水準測量數據對兩種技術得到的監測結果進行分析，得出結論如下所述。

1）在長期的礦區地面沉降監測過程中，D-InSAR技術監測到的累積沉降量大于SBAS-InSAR 技術，但識別到的沉降面積小于SBAS-InSAR 技術。在研究區內，D-InSAR 技術和SBAS-InSAR 技術監測到的最大累積沉降量分別為69 mm 和59 mm，識別出沉降量大于10 mm 的區域面積分別為0.60 km2和0.87 km2，其中，沉降量大于30 mm 的區域面積分別為0.10 km2和0.17 km2。

2）兩種技術監測到的沉降區域分布基本一致，D-InSAR 技術容易受到誤差影響，得到的沉降圖形在空間上較為離散，沉降量波動較大；而SBAS-InSAR技術得到的沉降結果在空間上平滑連續，沉降量相對穩定。

3）由水準測量數據驗證可得，兩種技術得到的監測結果與真實沉降情況基本吻合，D-InSAR 技術對較大沉降量的監測更為敏感。在沉降漏斗中心地面沉降梯度較大的區域，地面沉降量超出了InSAR 技術可監測范圍，兩種技術的監測精度都大幅下降。