馬 飛，李日印，劉小鵬，侯景怡，崔喆森，張垚杰

（1.長治學(xué)院計(jì)算機(jī)系，山西長治 046011；2.自然資源部第二地形測量隊(duì)，陜西西安 710054）

大規(guī)模開發(fā)利用礦產(chǎn)資源會直接影響礦區(qū)安全[1]。例如：2023 年2 月22 日，內(nèi)蒙古新井煤業(yè)有限公司露天煤礦發(fā)生大面積邊坡失穩(wěn)事故，造成人員財(cái)產(chǎn)損失。因此，對礦區(qū)地表變形進(jìn)行監(jiān)測具有重要的現(xiàn)實(shí)意義[2-4]。

差分干涉合成孔徑雷達(dá)(DInSAR)、小基線集In-SAR(SBAS-InSAR)和永久散射體InSAR(PS-InSAR)技術(shù)都是監(jiān)測采礦沉陷、凍土沉陷、建筑物沉陷和滑坡的常用方法[5-7]，具有空間覆蓋面大、成本低、精度高、全天候觀測能力等優(yōu)點(diǎn)，而常規(guī)測量技術(shù)（如水準(zhǔn)測量[8]、三維激光掃描[9]和全球?qū)Ш叫l(wèi)星系統(tǒng)[10]等）無法做到這一點(diǎn)。為了滿足開采沉陷長期監(jiān)測的需要，削弱大氣相位的影響，時序InSAR 提出了一種簡單有效的方法，即通過對一系列展開的干涉圖按時間跨度加權(quán)和平均來估計(jì)平均視線開采沉陷速度[11-12]。小基線集InSAR (SBAS-InSAR) 可以很好地估計(jì)變形結(jié)果[13-15]，在許多開采沉陷區(qū)得到了應(yīng)用?，F(xiàn)對內(nèi)蒙古新井煤礦應(yīng)用SBAS-InSAR技術(shù)，識別該礦區(qū)地表沉陷范圍，對沉陷區(qū)下沉數(shù)值進(jìn)行解算，分析災(zāi)前微小形變規(guī)律，以預(yù)防同類事故的再次發(fā)生。

1 SBAS-InSAR技術(shù)原理

SBAS-InSAR 通過設(shè)置特定的時間基線閾值來合并小的基線干擾對，根據(jù)最小參數(shù)準(zhǔn)則執(zhí)行最優(yōu)的形變相位序列，以限制時間和空間消相干。假設(shè)已經(jīng)采集了在時間序列(t0，t1，…，tn)相同區(qū)域中獲得的（n+1）個SAR 圖像。設(shè)置特定的時間基線閾值后，可以形成m個干涉對，且m滿足：

SBAS-InSAR 在進(jìn)行參數(shù)反演之前，需要將所有圖像與超級主影像配準(zhǔn)，并通過多視處理、地形差分和相位解纏來生成m幅差分干涉圖圖像?？紤]在時刻ta和t（b假設(shè)ta

式中：φ(tb,x)和φ(ta,x)分別為相位形變值；λ為雷達(dá)波長；φn為隨機(jī)噪聲；d(tb,x)和d(ta,x)分別為衛(wèi)星視線方向的形變量；v(x)為像元的形變速率；Ti為干涉對的時間基線長度。

這樣可以組成m個方程組：

式中：A為m×n的系數(shù)矩陣，每行對應(yīng)一幅干涉圖，每列對應(yīng)一個時間上的SAR 圖像，主影像所在列+1，輔圖像所在列-1，其余列為0。

如果m≥n，且A的秩是n，則最小二乘解得φ的估計(jì)值：

由式（4）可知，SBAS 技術(shù)的關(guān)鍵在于求時間序列干涉圖間形變速率的最小范數(shù)最小二乘解，可以利用時序干涉圖中解算出的形變相位值，作為時間域的約束條件解，簡化形變速率解算的復(fù)雜度，達(dá)到反演形變速率的目的。

假設(shè)V為形變速率，P為形變參數(shù)，以下約束條件被滿足：

將式（5）代入BV=δφ得：

在干涉圖中任意像元（x,r）的解算模型轉(zhuǎn)換為：

式中：LOS方向形變的平均加速度變化率?、相位平均速率、平均加速度用下式即可解得：

2 研究區(qū)概況及實(shí)驗(yàn)結(jié)果

2.1 研究區(qū)概況

圖1為內(nèi)蒙古新井煤業(yè)露天礦區(qū)地理位置，該礦區(qū)位于賀蘭山南段，礦區(qū)面積1.344 8 km2，地勢北高南低，海拔高程1 159.03～1 490.70 m。礦區(qū)為低山丘陵區(qū)，屬高原侵蝕性丘陵地貌，地形切割一般，基巖裸露，植被稀疏，為荒漠地區(qū)。氣候?qū)俅箨懶愿珊祷哪畾夂?，日溫差大，夏季炎熱，冬季?yán)寒，氣候干燥，雨量稀少，年最高氣溫38℃，最低-28℃，平均8℃～9℃，年降水量為146～ 198 mm，年蒸發(fā)量最高2 000 mm以上。地表植被屬于荒漠化草原向草原化荒漠過渡地帶。地表多沙質(zhì)化、礫石化和有龜裂結(jié)皮。

圖1 新井煤業(yè)露天煤礦地理位置

通過分析新井煤業(yè)露天礦區(qū)地表植被和地形地貌特征，該礦區(qū)地表屬于荒漠地區(qū)，植被稀少，適用C波段數(shù)據(jù)，但地面建筑物、工礦構(gòu)筑物、大型巖石等永久散射體物較少，適合采用SBAS-InSAR對礦區(qū)地表沉降情況進(jìn)行監(jiān)測。

2.2 SAR數(shù)據(jù)選取及預(yù)處理

表1為覆蓋該礦區(qū)的Sentinel-1A數(shù)據(jù)列表，包括衛(wèi)星拍攝時間、影像模式、極化方式和軌道方向等信息。根據(jù)歐空局(ESA)網(wǎng)站查詢哨兵-1A 數(shù)據(jù)情況發(fā)現(xiàn)，該衛(wèi)星2022年能夠覆蓋研究區(qū)的影像只有15景，觀測時間是2022年1月到2022年11月。根據(jù)干涉處理實(shí)驗(yàn)需要，下載了該區(qū)域干涉寬幅模式數(shù)據(jù)，單幅影像幅寬為250 km，距離向分辨率為4 m，方位向分辨率為13 m，多視比例設(shè)置為3∶1。同時下載到所有影像的精密軌道數(shù)據(jù)，用于干涉數(shù)據(jù)處理過程中的圖像配準(zhǔn)、去平地相位等流程，提高數(shù)據(jù)處理精度。另外，為提高數(shù)據(jù)處理效率，對原始SAR影像進(jìn)行了裁切。

表1 覆蓋礦區(qū)的Sentinel-1A 數(shù)據(jù)列表

2.3 時序InSAR 技術(shù)數(shù)據(jù)處理流程

2.3.1 短基線集干涉圖組合

圖2 為干涉對時空基線圖。將選取的15 景影像按照如下策略進(jìn)行自由組合，形成干涉對：時間基線ΔT≤100 d，空間基線ΔB≤200 m，以20220401 為超級主影像。最終得到43個干涉對。

圖2 SAR干涉對時空基線圖

2.3.2 干涉圖生成

圖3 是在新井煤業(yè)露天礦區(qū)生成的以20220401為主影像的差分干涉圖，該圖像為SAR 斜距坐標(biāo)系統(tǒng)。圖中可以看到，礦區(qū)中心位置形變條紋非常清晰，隨著時間的推移，礦區(qū)中心位置形變量在不斷變化；時間基線越長，干涉條紋周期變化越多；在圖(i)、(j)等干涉對的中心位置出現(xiàn)干涉條紋疊加和失相干。

圖3 以20220401為主影像生成的的干涉圖

2.3.3 SBAS時間序列形變和形變速率結(jié)果

圖4 為新井礦業(yè)露天礦區(qū)形變速率圖，紅色代表形變方向與衛(wèi)星視線方向一致，可解算為下沉值，藍(lán)色代表形變方向與衛(wèi)星視線方向相反，可解算為抬升值。圖中可見，該礦區(qū)大部分區(qū)域相對穩(wěn)定，只有在此次滑坡事故中的核心區(qū)發(fā)生劇烈形變，整體形變速率非常大，沉降中心位置最大形變速率達(dá)到-160～ -209 mm/a，說明此處滑坡2022 年開始已經(jīng)出現(xiàn)不穩(wěn)定狀態(tài)。另外，在該區(qū)域的東南側(cè)也出現(xiàn)一個小范圍的形變區(qū)，該區(qū)域的形變速率在-60～ -110 mm/a，需要引起相關(guān)部門的高度關(guān)注，避免發(fā)生新的地質(zhì)災(zāi)害。

圖4 災(zāi)變核心區(qū)形變速率圖

為進(jìn)一步分析新井煤業(yè)露天礦區(qū)核心區(qū)災(zāi)前形變情況，將該形變區(qū)域數(shù)值提取出來并進(jìn)行處理。該區(qū)域共解算得到2 475 個形變點(diǎn)，圖4 顯示該區(qū)域自邊坡的坡頂向坡底發(fā)生形變量由大到小，2022年1月至11 月坡頂累計(jì)形變量超過170 mm，坡底累計(jì)形變量超過150 mm。

圖5 為邊坡在2022 年的形變速率剖面圖，橫坐標(biāo)為距離邊坡坡頂GCP1 號點(diǎn)的距離，縱坐標(biāo)為形變速率值。為進(jìn)一步分析該礦區(qū)邊坡失穩(wěn)前的形變規(guī)律，在邊坡選取了四個點(diǎn)命名為GCP1～4進(jìn)行分析，分別分析其剖面線和形變時序圖。圖中可知，邊坡頂部下沉速率遠(yuǎn)遠(yuǎn)高于坡底位置，自坡頂向坡底形變速率逐漸下降，形變速率從-180～-140 mm/a。由于該區(qū)域邊坡角度較大，坡頂區(qū)松散層不斷下沉，經(jīng)過一段時間的累積，需要采取任何進(jìn)行邊坡加固。

圖5 形變速率剖面圖

圖6 為新井煤業(yè)露天礦區(qū)災(zāi)變核心區(qū)內(nèi)邊坡形變點(diǎn)的時序下沉變化情況。由圖可知，在邊坡頂部的GCP1 號點(diǎn)累積形變量最大，在底部的GCP4 號點(diǎn)累積形變量稍小，這與圖5 中的剖面圖一致；該邊坡從2022 年2 月開始出現(xiàn)下沉，2022 年3 月至6 月累積形變趨勢開始增大，這可能與該地區(qū)前期轉(zhuǎn)暖凍融解凍有關(guān)，也可能與礦區(qū)露天開采活動有關(guān)。

3 結(jié)論

(1)利用SBAS-InSAR 方法可以有效識別出露天采礦區(qū)地表變形，這是一種快速、高效、低廉的技術(shù)，可以為礦區(qū)地表邊坡失穩(wěn)、地面塌陷、地面沉降、地裂縫、滑坡、崩塌等地質(zhì)災(zāi)害提供預(yù)警。

(2)內(nèi)蒙古新井煤業(yè)露天采礦區(qū)核心災(zāi)變區(qū)在2022 年1-11 月持續(xù)發(fā)生下沉形變。分析SBAS-In-SAR 方法獲取的2 475 個形變點(diǎn)發(fā)現(xiàn)該區(qū)域最大形變量超過160 mm，沉降速率超過200 mm/a；在核心災(zāi)變區(qū)東南側(cè)亦有一個緩慢形變區(qū)，累積形變量稍小，但需要有關(guān)部門關(guān)注，避免發(fā)生新的地質(zhì)災(zāi)害。

(3)內(nèi)蒙古新井煤業(yè)露天采礦區(qū)失穩(wěn)邊坡在災(zāi)前一年間持續(xù)發(fā)生下沉，且邊坡坡頂形變速率高于坡底速率，坡頂形變速率超過180 mm/a，頂部持續(xù)微小的下沉形變是導(dǎo)致邊坡應(yīng)力失衡的重要原因。

(4)由于缺少地表監(jiān)測點(diǎn)，故SBAS-InSAR 的形變監(jiān)測值無法得到驗(yàn)證，但該技術(shù)獲得的時間序列累積變形能更直觀地反映沉陷盆地的時間序列變化趨勢，在探測變形區(qū)域范圍和形態(tài)方面具有優(yōu)勢，有助于定性確定形變類型。

致謝：本實(shí)驗(yàn)的Sentinel-1A 數(shù)據(jù)和SRTM 外部DEM數(shù)據(jù)均來自歐空局網(wǎng)站，在此表示感謝。