D-InSAR 技術(shù)與SBAS-InSAR 技術(shù)在礦區(qū)地面沉降監(jiān)測中的應(yīng)用研究

王 鵬，程海強(qiáng)，王 翔

（1.山東省三河口礦業(yè)有限責(zé)任公司地質(zhì)測量科，山東 濟(jì)寧 277600；2.棗莊礦業(yè)（集團(tuán)）有限責(zé)任公司地質(zhì)測量部，山東 棗莊 277000；3.山東建筑大學(xué)測繪地理信息學(xué)院，山東 濟(jì)南 250101）

0 引言

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

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

1 原理與方法

1.1 D-InSAR 技術(shù)基本原理

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

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

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

1.2 SBAS-InSAR 技術(shù)基本原理

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

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

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

則相位為式（5）。

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

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

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

圖1 數(shù)據(jù)處理流程圖Fig.1 Flow chart of data processing

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

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

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

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

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

3 數(shù)據(jù)處理

3.1 D-InSAR 技術(shù)

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

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

3.2 SBAS-InSAR 技術(shù)

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

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

4 結(jié)果分析

4.1 沉降監(jiān)測結(jié)果對比分析

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

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

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

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

4.2 沉降監(jiān)測結(jié)果精度驗(yàn)證

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

圖8 水準(zhǔn)點(diǎn)累積沉降量折線圖Fig.8 Line graph of cumulative subsidence at the level

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

5 結(jié)論

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

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

2）兩種技術(shù)監(jiān)測到的沉降區(qū)域分布基本一致，D-InSAR 技術(shù)容易受到誤差影響，得到的沉降圖形在空間上較為離散，沉降量波動(dòng)較大；而SBAS-InSAR技術(shù)得到的沉降結(jié)果在空間上平滑連續(xù)，沉降量相對穩(wěn)定。

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