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

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

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

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

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

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

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

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

這樣可以組成m個(gè)方程組：

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

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

由式（4）可知，SBAS 技術(shù)的關(guān)鍵在于求時(shí)間序列干涉圖間形變速率的最小范數(shù)最小二乘解，可以利用時(shí)序干涉圖中解算出的形變相位值，作為時(shí)間域的約束條件解，簡(jiǎn)化形變速率解算的復(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，地勢(shì)北高南低，海拔高程1 159.03～1 490.70 m。礦區(qū)為低山丘陵區(qū)，屬高原侵蝕性丘陵地貌，地形切割一般，基巖裸露，植被稀疏，為荒漠地區(qū)。氣候?qū)俅箨懶愿珊祷哪畾夂?，日溫差大，夏季炎熱，冬季?yán)寒，氣候干燥，雨量稀少，年最高氣溫38℃，最低-28℃，平均8℃～9℃，年降水量為146～ 198 mm，年蒸發(fā)量最高2 000 mm以上。地表植被屬于荒漠化草原向草原化荒漠過(guò)渡地帶。地表多沙質(zhì)化、礫石化和有龜裂結(jié)皮。

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

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

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

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

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

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

2.3.1 短基線集干涉圖組合

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

圖2 SAR干涉對(duì)時(shí)空基線圖

2.3.2 干涉圖生成

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

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

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

圖4 為新井礦業(yè)露天礦區(qū)形變速率圖，紅色代表形變方向與衛(wèi)星視線方向一致，可解算為下沉值，藍(lán)色代表形變方向與衛(wèi)星視線方向相反，可解算為抬升值。圖中可見(jiàn)，該礦區(qū)大部分區(qū)域相對(duì)穩(wěn)定，只有在此次滑坡事故中的核心區(qū)發(fā)生劇烈形變，整體形變速率非常大，沉降中心位置最大形變速率達(dá)到-160～ -209 mm/a，說(shuō)明此處滑坡2022 年開(kāi)始已經(jīng)出現(xiàn)不穩(wěn)定狀態(tài)。另外，在該區(qū)域的東南側(cè)也出現(xiàn)一個(gè)小范圍的形變區(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ù)值提取出來(lái)并進(jìn)行處理。該區(qū)域共解算得到2 475 個(gè)形變點(diǎn)，圖4 顯示該區(qū)域自邊坡的坡頂向坡底發(fā)生形變量由大到小，2022年1月至11 月坡頂累計(jì)形變量超過(guò)170 mm，坡底累計(jì)形變量超過(guò)150 mm。

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

圖5 形變速率剖面圖

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

3 結(jié)論

(1)利用SBAS-InSAR 方法可以有效識(shí)別出露天采礦區(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 個(gè)形變點(diǎn)發(fā)現(xiàn)該區(qū)域最大形變量超過(guò)160 mm，沉降速率超過(guò)200 mm/a；在核心災(zāi)變區(qū)東南側(cè)亦有一個(gè)緩慢形變區(qū)，累積形變量稍小，但需要有關(guān)部門關(guān)注，避免發(fā)生新的地質(zhì)災(zāi)害。

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

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

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