時間序列InSAR技術(shù)在礦區(qū)地表形變監(jiān)測中的應(yīng)用
——以內(nèi)蒙古霍林河露天礦區(qū)為例

王洪明，李如仁，覃怡婷，劉竹青，顧 駿

（1.沈陽建筑大學(xué)交通工程學(xué)院,遼寧 沈陽 110168；2.扎魯特旗扎哈淖爾煤業(yè)有限公司,內(nèi)蒙古通遼 029100；3.沈陽市勘察測繪研究院,遼寧 沈陽 110004）

0 引言

露天礦的開采導(dǎo)致礦區(qū)生態(tài)遭到不同程度的破壞，礦區(qū)及其周邊的地質(zhì)災(zāi)害主要表現(xiàn)形式是地表形變，同時礦區(qū)的形變對城區(qū)地表也有不同程度的影響[1?2]。傳統(tǒng)的監(jiān)測手段（如水準(zhǔn)測量、GPS 等）對于礦區(qū)形變提取有一定的局限性，水準(zhǔn)測量方法過于依賴人工，在連續(xù)性和實時性方面略有不足；GPS 方法技術(shù)成本過高，且對于面域的沉降范圍有較大的限制。合成孔徑雷達可對地表進行全天時、全天候?qū)Φ赜^測[3]。

雷達干涉測量技術(shù)可以對露天礦區(qū)及其周邊地表形變進行監(jiān)測，具有大范圍和高精度等優(yōu)點[3]。InSAR的形變監(jiān)測發(fā)展較迅速，最初的D-InSAR技術(shù)只能獲取地表短時間的形變，后來采用時序InSAR技術(shù)提取長時間跨度的形變[4?6]。傳統(tǒng)的時序InSAR技術(shù)對于礦區(qū)及其周邊地表的形變提取存在明顯的不足，由于露天礦地表的相干性差，故只能提取分布式點目標(biāo)，而城區(qū)地表具有足夠相干性可以提取永久散射體目標(biāo)。劉一霖等[7]利用時序研究開采進程中地表大量級沉陷的完整形變時間序列，得到開采工作面地表形變的時空演變規(guī)律；仝云霄等[8]利用經(jīng)典的D-InSAR技術(shù)對礦區(qū)沉降進行了精細化分析，與最鄰近時間的精密水準(zhǔn)監(jiān)測結(jié)果總體趨勢吻合，為礦區(qū)整體規(guī)劃、災(zāi)害預(yù)警和開采沉陷治理等提供參考；張學(xué)東等[9]以唐山市為例，使用相干目標(biāo)短基線InSAR技術(shù)監(jiān)測礦業(yè)城市地面沉降，查明了其地面沉降量及其空間分布特征，最大年沉降速率達到?46.8 mm/a，主城區(qū)沉降速率普遍低于?11 mm/a；任文靜等[10]、蔣金雄等[11]利用小基線集（SBAS）方法和相干矩陣特征值分解（T-EVD）的DSInSAR技術(shù)分別對礦區(qū)進行地表沉降提取，為礦區(qū)地表形變監(jiān)測提供了良好的應(yīng)用前景；謝文斌等[12]利用時序InSAR技術(shù)對撫順市地表形變信息進行提取，結(jié)合信息熵方法分析地表不均勻形變現(xiàn)狀?？傮w而言，目前在礦區(qū)及其周邊城區(qū)地表形變提取及分析研究較少。隨著露天礦的逐年開采和城市的迅速發(fā)展，礦區(qū)沉陷及礦業(yè)城市地表沉降是亟待探究的復(fù)雜問題。

傳統(tǒng)的時序分析方法誤差較大，也無法整體評估出礦區(qū)形變對其周邊城區(qū)地表形變產(chǎn)生的影響。鑒于此，基于露天礦區(qū)及其周邊城區(qū)地表的相干性差異巨大的特點，文章提出改進的時序InSAR技術(shù)，此方法在點目標(biāo)提取的基礎(chǔ)上，使用狄龍尼三角網(wǎng)將分布式目標(biāo)與永久散射體目標(biāo)進行統(tǒng)一構(gòu)網(wǎng)，進行城市與礦區(qū)的形變一體化反演。以霍林河礦區(qū)及其周邊城區(qū)為試驗案例，獲得了露天礦區(qū)及其周邊城區(qū)地表的形變速率。

1 研究區(qū)與數(shù)據(jù)

研究區(qū)霍林河煤田地處大興安嶺南段脊部，是一個山間盆地，四周為中低山巒所環(huán)抱（圖1），海拔標(biāo)高1 100 ～1 350 m。盆地內(nèi)部地勢較平坦，東北和西南兩端為低丘陵，中間區(qū)段是開擴的平原，海拔標(biāo)高930 ～980 m。研究區(qū)內(nèi)部包括中部與北部兩個露天礦和霍林郭勒市市區(qū)。為了較好覆蓋研究區(qū)，采用29 景Sentinel 1 降軌數(shù)據(jù)及90 m 分辨率Version4 SRTM 數(shù)據(jù)，Sentinel 1 數(shù)據(jù)具體參數(shù)如表1。

表1 研究區(qū)SAR 數(shù)據(jù)參數(shù)表Table 1 SAR data parameter table in the research area

圖1 研究區(qū)Google earth 影像圖Fig.1 Google Earth image of the research area

2 時序InSAR 處理

由于露天礦區(qū)及其城區(qū)周邊地表相干性差異巨大，為保持礦區(qū)周邊具備一定的相干性，采用改進的時序InSAR技術(shù)，該方法通過嚴格篩選差分干涉數(shù)據(jù)，對差分干涉數(shù)據(jù)提取分布式散射體與永久散射體，作為形變分析的目標(biāo)點，通過進行回歸分析與模型精化，得到研究區(qū)地表形變速率。改進的時序InSAR技術(shù)既繼承了傳統(tǒng)時序InSAR 的技術(shù)優(yōu)勢，又很好的解決了目標(biāo)研究區(qū)域相干性差異性較大的問題，具體技術(shù)流程如圖2所示。

圖2 改進時序InSAR技術(shù)流程圖Fig.2 Flow chart of improved sequence InSAR technical

2.1 預(yù)處理與差分干涉圖生成

通過GAMMA 預(yù)先設(shè)置時間基線為30 d，空間基線為臨界基線的30% 進行干涉生成，為了更好的篩選干涉數(shù)據(jù)[11?13]，選用SRTM 數(shù)據(jù)除去差分干涉對原始的地形相位，通過對差分干涉數(shù)據(jù)的易發(fā)生形變區(qū)域進行目視分析比較剔除干涉效果差的像對，最后選擇34 個差分干涉結(jié)果進行后續(xù)試驗處理，且為了保證低相干區(qū)域結(jié)果的準(zhǔn)確性對9 個相對組成的子集進行3D 解纏，低相干地區(qū)參照3D 解纏區(qū)域進行時序處理。篩選出部分的差分干涉結(jié)果如圖3所示，6 張差分干涉圖在礦區(qū)位置都有明顯的形變漏斗產(chǎn)生。

圖3 差分干涉結(jié)果圖Fig.3 Differential interference result diagram

2.2 干涉點目標(biāo)選取

由于礦區(qū)及其周邊相干性較低，導(dǎo)致干涉結(jié)果較差，需要對礦區(qū)及其周邊分布式目標(biāo)點進行提取，文章提出經(jīng)過參數(shù)優(yōu)化的分布式散射點選取方法對研究區(qū)進行選點，主要包括以下三步：

（1）將同一區(qū)域分干涉結(jié)果像素數(shù)據(jù)進行時間維度平均，采用24×24 窗口進行計算，將其中心區(qū)視為均值即為參考值，鄰域像素值作為待估計值，為了減小誤差，考慮空間臨近關(guān)系，只把臨近點像素作為待估計值[14]。

（2）運用瑞麗分布函數(shù)獲得置信域區(qū)間，函數(shù)如式（1）：

式中：γref——參考像素均值；

γarb——任意臨近像素均值；

z1?α/2——標(biāo)準(zhǔn)正態(tài)分位點，服從標(biāo)準(zhǔn)正態(tài)分布，為了增加同質(zhì)點樣本數(shù)量，α預(yù)先設(shè)置值為50%；

P——概率函數(shù)模型；

L——多視視數(shù)，得到置信域區(qū)間并獲得同質(zhì)點[15?16]。

（3）為進一步提高樣本選取精度，在小樣本數(shù)據(jù)情況下，傳統(tǒng)的高斯假設(shè)不能夠逼近真實的分布，可以根據(jù)SAR 強度服從指數(shù)分布的性質(zhì)，其累加服從伽馬分布，可以得到更加精確的置信域區(qū)間。

式中，σ代表SAR 強度指數(shù)分布的期望，根據(jù)伽馬分布性質(zhì)可知其分位點為 2/α，在有N 個樣本的情況下，則可得置信域區(qū)間為：

將所有樣本均值帶入計算，提高同質(zhì)樣本的選取精度，進而得到分布式目標(biāo)點。分布式目標(biāo)（DS 點）在干涉體系中占據(jù)大比例，包括裸漏的巖石、低矮植被等，對礦區(qū)選取DS 點可以彌補永久散射體數(shù)量不足的缺點，避免失相干，進而提高干涉處理精度。

通過時序技術(shù)分別提取永久散射體點目標(biāo)，將分布式點目標(biāo)與永久散射體目標(biāo)進行疊加（圖4），采用Delaunay 三角網(wǎng)法對兩種相干目標(biāo)點統(tǒng)一構(gòu)網(wǎng)。差分干涉圖中,相鄰目標(biāo)點相位差表示為：

圖4 干涉點目標(biāo)分布圖Fig.4 Image of interference point target distribution

式中：Δυ——相鄰點之間的形變速率；

Δε——高程誤差；

Δω——殘余相位（包括大氣延遲相位，非線性形變相位和噪聲相位）；

B⊥——干涉對的空間基線；

T——干涉對的時間基線；

λ——波長；

R——傳感器到目標(biāo)的距離；

θ——雷達入射角。

2.3 干涉點目標(biāo)形變處理

對目標(biāo)點進行第一次回歸分析時可以得到其高程誤差、線性形變速率以及殘余相位，由于第一次回歸分析主要是為了估算初始的形變速率范圍與數(shù)據(jù)處理效率，并沒有顧及大氣的與非線性形變相位的影響，因此在回歸分析時，需要通過設(shè)置相位標(biāo)準(zhǔn)差小于1.5 rad剔除質(zhì)量較低的點，通過兩次迭代可以使高程誤差的趨近于0，標(biāo)準(zhǔn)差顯著降低[17]。在慮性形變的條件下，形變速率增量以及高程誤差增量通過整體相位相干系數(shù)最大化模型來進行模型優(yōu)化，公式為：

其中，γ越大，表示模型估計與差分相位越相似誤差越小，以 γ作為權(quán)重，使用帶權(quán)的最小二乘方法，從已知參考點解算出稀疏格網(wǎng)上目標(biāo)點的時序線性沉降速率 υ和高程誤差ε。將解算的每個目標(biāo)點的線性形變和高程誤差從初始差分干涉圖中減去,即得到殘余相位 Δω，其包含了大氣擾動相位、非線性形變相位以及噪聲相位。

再對點目標(biāo)進行第二次回歸分析，目的是使得非線性形變相位于大氣延遲相位分離，因為大氣延遲相位在1 km2范圍內(nèi)具有較強的空間相關(guān)性，在空間域上表現(xiàn)為平滑的低頻信號，在時間域上呈現(xiàn)高頻信號可被當(dāng)作隨機噪聲處理[16?18]。因此根據(jù)其各自的特征進行合適的空間濾波，可將非線性形變相位和大氣延遲相位分離，得到精確的非線性形變量。最后，根據(jù)計算的線性形變量和非線性形變量，得到形變相位進行相位并將其轉(zhuǎn)到為形變量，視線向形變分解為垂直向形變。

3 結(jié)果與分析

對時序InSAR 結(jié)果進行區(qū)域裁剪，去除失相干區(qū)域，得到礦區(qū)及其周邊城區(qū)的形變速率圖與部分形變累積量如圖5、圖6，可得到在2019年全年中部露天礦區(qū)靠近城市側(cè)的最大沉降速率達到569 mm/a，北部露天礦區(qū)遠離城市側(cè)最大沉降速率達到630 mm/a，中部露天礦附近城區(qū)地表最大沉降速率達到23 mm/a，北部露天礦附近城區(qū)地表最大沉降速率達到26 mm/a。由于露天礦區(qū)周邊地質(zhì)條件不穩(wěn)定，露天礦旁道路年形變速率最大達到71 mm/a，城區(qū)的西側(cè)沉降速率明顯高于東側(cè)。

圖5 地表年形變速率圖Fig.5 Annual surface deformation rate map

圖6 形變累積量圖Fig.6 Deformation cumulant diagram

3.1 地表形變及影響機理分析

礦區(qū)與城區(qū)之間地質(zhì)條件不夠穩(wěn)定，最大沉降速率達到71 mm/a，平均沉降速率達到28 mm/a。試驗區(qū)在1—3月平均沉降量5.3 mm、4—6月平均沉降量14.7 mm、7—9月平均沉降量16.2 mm、10—12月平均沉降量12.5 mm，其中4—6月和7—9月兩個時間階段沉降累積量較大。三個月為一組時間段，其平均沉降量、最大沉降量見表2。

表2 研究區(qū)累積沉降量表Table 2 Cumulative settlement scale

4—6月形變量分析，中部露天礦最大沉降量達到57 mm，北部露天礦最大形變量達到54 mm，城區(qū)地表最大沉降量為4 mm。通過對降雨量的分析，霍林郭勒市4—6月降雨量與往年同期相比減少約30%，土地解凍是礦區(qū)沉降量劇增的主要因素。中部礦區(qū)西側(cè)的地表形變量為正值，是由于礦區(qū)西側(cè)地質(zhì)條件不穩(wěn)定和排土場堆棄排土所致。7—9月形變量分析，中部露天礦最大沉降量達到61 mm，北部露天礦最大沉降量達到59 mm，城市地表沉降量最大值達到12 mm。由于霍林郭勒市7—9月的降雨量與往年同期相比增加約50%，雨水可能導(dǎo)致地表沉降量增大，北部露天礦南側(cè)地表形變量為正值，由于北部礦區(qū)開采排土場積土導(dǎo)致。對比上述兩時間段形變量，可以看出中部與北部露天礦的南側(cè)形變量均從由負值變?yōu)檎?，均為排土堆砌所致，可以從礦區(qū)開采得到驗證。采用Kenda 系數(shù)法相關(guān)系數(shù)驗證雨季時兩者之間關(guān)系，經(jīng)計算得系數(shù)為0.600，在夏季與秋季侵蝕性降雨與地表形變速率有正向的影響關(guān)系，可知在雨季降雨侵蝕力對露天礦坑周邊地表環(huán)境變化產(chǎn)生較大影響，圖7 為侵蝕性降雨量與監(jiān)測點形變量關(guān)系圖。

由于礦區(qū)形變復(fù)雜，將像元值數(shù)目低于總數(shù)5%的剔除，并將形變速率圖與強度影像進行疊加得到圖8，可以看出，城區(qū)地表形變速率與其距離露天礦的距離成反比，即城區(qū)離露天礦越近地表的形變速率值越大，從客觀上也符合礦業(yè)城市地表沉降的規(guī)律，城市并未出現(xiàn)大規(guī)模的不規(guī)則沉降，礦區(qū)與城市之間道路沉降值較大均在32 mm/a 以上，靠近露天礦一側(cè)的道路沉降值遠大于另一側(cè)，道路傾斜嚴重。

圖8 基于強度影像的地表形變速率圖Fig.8 Surface deformation rate map based on intensity image

3.2 GPS 數(shù)據(jù)驗證與分析

為了驗證改進得時序InSAR技術(shù)監(jiān)測結(jié)果的精度，采用同期GPS 結(jié)果進行比對。在圖7 上提取6 個點并與同樣的GPS 點監(jiān)測結(jié)果進行比較（表3），最小互差0.76，說明該InSAR 監(jiān)測方法是可行與可靠的。

圖7 侵蝕性降雨量與監(jiān)測點形變速率圖Fig.7 Chart of erosive rainfall and deformation rate at monitoring points

表3 GPS 監(jiān)測結(jié)果與同期InSAR 監(jiān)測結(jié)果比較Table 3 Comparison of GPS monitoring results with InSAR monitoring results in the same period

4 結(jié)論

（1）考慮到露天礦區(qū)地表的低相干性，提出了一種改進的時序InSAR技術(shù)，采用24×24 窗口同質(zhì)樣本選取方法提取分布式目標(biāo)，最大化的識別礦區(qū)的相干點目標(biāo)，對于霍林郭勒市城區(qū)采用經(jīng)典的永久散射體干涉測量方法提取永久散射體，對兩種干涉點目標(biāo)進行統(tǒng)一的構(gòu)網(wǎng)分析和形變反演，與GPS 數(shù)據(jù)進行對比驗證，證明改進的時序InSAR技術(shù)方法具有較好的監(jiān)測效果，對礦業(yè)城市地表形變監(jiān)測提出了一種新的InSAR 監(jiān)測方法。

（2）將霍林郭勒市的月降雨量數(shù)據(jù)和干涉點月沉降量關(guān)聯(lián)分析，Kenda 系數(shù)達到0.736 證明二者具有較高相關(guān)性，可知侵蝕性降雨對于露天礦區(qū)及其周邊地表沉降具有正向的影響。

（3）通過本方法得到露天礦區(qū)地表沉降速率最大可達346 mm/a，城區(qū)地表形變速率最大可達31 mm/a，露天礦和城區(qū)之間道路沉降量巨大，且公路兩側(cè)形變速率值差異巨大，道路出現(xiàn)明顯的傾斜。剔除形變速率極值，可知城區(qū)地表形變速率從客觀上符合礦業(yè)城市地表沉降的規(guī)律，即城區(qū)形變速率值與其到露天礦的距離成反比，可得城區(qū)并未出現(xiàn)大規(guī)模的不規(guī)則沉降，礦區(qū)與城市之間道路沉降值較大均在32 mm/a 以上，靠近露天礦一側(cè)的道路沉降值遠大于另一側(cè)，道路傾斜嚴重。

感謝歐洲空間局( ESA)為本文提供的SAR 數(shù)據(jù)和國家科學(xué)數(shù)據(jù)服務(wù)平臺提供SRTM-Version4 數(shù)據(jù)。