基于InSAR技術(shù)的地表沉陷變形邊界提取精度研究

韋金龍，趙利安

(1.山西長平煤業(yè)有限責(zé)任公司， 山西 晉城 048006； 2.遼寧工程技術(shù)大學(xué)， 遼寧 阜新 123000)

煤炭開采造成大面積的地表沉陷，其邊界范圍的確定是設(shè)定保護(hù)煤柱和保護(hù)地面設(shè)施的關(guān)鍵，對防止地質(zhì)災(zāi)害的發(fā)生起著重大的意義[1-2]. 傳統(tǒng)測量手段(GPS、水準(zhǔn)測量等)會受環(huán)境與技術(shù)的制約，而且存在監(jiān)測范圍小、工作量大、效率低、時空分辨率低等缺陷[3-4]. 時序InSAR技術(shù)(interferometric synthetic aperture radar，InSAR)具有易操作、成本小、全天候、不受天氣情況影響、監(jiān)測范圍廣、高空間分辨率、測量精度高等優(yōu)點(diǎn)，可以做到以平面為基準(zhǔn)的區(qū)域監(jiān)測[5-8]. 基于傳統(tǒng)D-InSAR技術(shù)的短基線集(small baseline subset，SBAS) 技術(shù)能夠有效地提高觀測數(shù)據(jù)的時間采樣率與空間密度，監(jiān)測一段時間范圍內(nèi)的沉陷情況，獲得研究區(qū)域內(nèi)的沉降信息和演化規(guī)律[9-12].

對于SBAS-InSAR技術(shù)而言，當(dāng)?shù)叵虏傻V產(chǎn)生非常大量級的形變，且地表上往往覆蓋茂盛的植被時，會造成大面積失相干現(xiàn)象，使監(jiān)測結(jié)果受大氣效應(yīng)的嚴(yán)重影響[13-17].為此，選擇利用短基線集(SBAS-InSAR)技術(shù)對采煤沉陷盆地的邊界進(jìn)行提取，根據(jù)實測數(shù)據(jù)和地表動態(tài)沉降預(yù)測，驗證沉陷變形邊界提取的有效性，為形變量級較大的礦區(qū)地表沉降研究提供新的解決思路。

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

SBAS-InSAR(小基線集，Small baseline subset algorithm，SBAS)方法是由Berardino等最先提出，利用同一區(qū)域的多景SAR影像，選取這些影像上滿足一定條件的像素進(jìn)行時間序列分析，與PS-InSAR不同之處在于，SBAS-InSAR技術(shù)需要選擇多組干涉對，且干涉對需滿足時間基線和空間基線達(dá)到一定閾值，并對干涉對進(jìn)行相位空間的多視化處理，盡管降低了空間分辨率，但是最大程度保證了干涉相位的相干性。

假設(shè)覆蓋同一實驗區(qū)按照時間順序t0，t1,…，tn排列的N+1幅SAR影像，自動選擇公共主影像，并將其它SAR影像與該公共主影像進(jìn)行配準(zhǔn)。依據(jù)干涉條件進(jìn)行組合，N+1幅SAR影像可以組合成M幅干涉圖，其中M的大小將滿足以下條件：

(1)

對于從影像tA和主影像tB(tB>tA)時刻獲取的SAR影像，生成的第j幅差分干涉圖，其像素的干涉相位與方位向坐標(biāo)x、距離向坐標(biāo)r存在以下關(guān)系：

δφj(x,r)=φ(tB,x,r)-φ(tA,x,r)≈

(2)

為了將地表形變時間序列表示為具有物理意義的形式，對式(2)中相位表示為獲取的兩個觀測時間內(nèi)的平均相位速度與時間之積，則平均相位速度可以表示為：

(3)

第j幅干涉圖的相位值表示為每一時間段內(nèi)平均速率在主、輔影像的積分，其表達(dá)式為：

(4)

將各差分干涉圖進(jìn)行組合，寫成矩陣形式為：

Bv=δφ

(5)

式中，B表示M×N的系數(shù)矩陣。則vT可以表示為：

(6)

當(dāng)M≥N時，即系數(shù)矩陣滿秩，利用最小二乘法即可求解出φ的估計值：

(7)

當(dāng)M≤N時，系數(shù)矩陣B容易出現(xiàn)秩虧現(xiàn)象，根據(jù)奇異值分解(SVD)方法計算得到B的廣義逆矩陣，從而獲得形變速率v的最小范數(shù)解，最終以此方法對每一時間段內(nèi)的速度進(jìn)行積分，然后進(jìn)行累加就可以獲取整個時間范圍內(nèi)的形變量。

SBAS-InSAR技術(shù)方法的主要步驟如下：

1) 配準(zhǔn)參考影像的選取與影像配準(zhǔn)。SBAS-InSAR數(shù)據(jù)處理過程中，需要選擇滿足一定時間和空間基線閾值要求的影像進(jìn)行組合，因此可以選取SAR影像集中的任意一景作為配準(zhǔn)的參考影像。然而，受時間失相干、幾何失相干以及多普勒質(zhì)心失相干等因素的影響，選擇不同的配準(zhǔn)參考影像，得到的配準(zhǔn)精度存在差異，從而會影響差分干涉相位的質(zhì)量。為了減少配準(zhǔn)誤差的影響，需要依據(jù)最優(yōu)準(zhǔn)則選取SAR影像序列中的一景作為配準(zhǔn)參考影像，以保證配準(zhǔn)精度最優(yōu)。

總之，大致在“文化大革命”及其以前十年間，是抗戰(zhàn)勝利紀(jì)念活動高潮的回落時期，活動的規(guī)模和規(guī)格也大受影響。從當(dāng)時國內(nèi)國際局勢來看，可以確信的原因包括蘇聯(lián)政局的動蕩、中蘇關(guān)系的逐步惡化、中國國內(nèi)工作重心的轉(zhuǎn)變乃至“文化大革命”的動亂等。但紀(jì)念活動仍得以延續(xù)，主要方式有：召開中央一級的紀(jì)念大會，逢五逢十由《人民日報》發(fā)表社論，以東北地區(qū)為主的地方性群眾集會、瞻仰抗戰(zhàn)紀(jì)念設(shè)施等。

2) 干涉對選取。將配準(zhǔn)后的N+1幅SAR影像自由組合，得到N×(N+2)/2個干涉對，計算所有干涉對的時間基線和空間基線，根據(jù)給定的時間基線閾值和空間基線閾值，從所有可能組成的干涉對中剔除不符合閾值的部分干涉對，其余干涉對繼續(xù)用于后續(xù)的數(shù)據(jù)處理。

3) 差分干涉處理。使用DEM進(jìn)行差分處理，可以消除平地相位和地形相位的影響，進(jìn)而得到每個干涉對的差分干涉相位，并根據(jù)數(shù)據(jù)處理的需要，對差分干涉相位進(jìn)行多視以提高干涉相位的相干性。此時差分干涉相位包含形變、DEM誤差、大氣誤差、軌道誤差和隨機(jī)噪聲等成分。

4) SVD形變參數(shù)求解。由于SBAS干涉對的組網(wǎng)中，解纏的差分干涉相位是兩景影像間的相對各干涉對主影像的變化量，需要采用SVD分解的方法，將這些相對變化量轉(zhuǎn)換到單一參考時間的相對變化量，從而得到原始形變序列。此時的形變序列還包含大氣誤差、DEM殘差和軌道誤差等。

5) 形變序列恢復(fù)。由于原始形變序列受大氣誤差和軌道誤差的影響，只有去除此類誤差才能恢復(fù)出準(zhǔn)確的形變序列。同PS-InSAR的處理策略，根據(jù)大氣延遲誤差和軌道誤差在時空域上的不同信號特征，通過時間高通濾波和空間低通濾波將其分離開來，獲得準(zhǔn)確的形變序列。

SBAS-InSAR時序形變分析方法數(shù)據(jù)處理流程見圖1. 該算法也可以通過迭代計算獲取精確的DEM，同時還可以從原始形變序列中分離出大氣誤差和軌道誤差，最終獲得準(zhǔn)確的形變時間序列。

圖1 SBAS-InSAR時序形變分析方法數(shù)據(jù)處理流程圖

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

2.1 研究區(qū)概況

以山西三交河煤礦為例，煤礦位于山區(qū)地形，煤田采區(qū)邊界溝谷縱橫，山勢磅礴，具有復(fù)雜的地形。煤礦的地勢為中間高，左右兩側(cè)低，地勢最高的點(diǎn)位于煤礦西南角，有1 469 m的海拔標(biāo)高；最低點(diǎn)位于煤礦東緣的山谷底，有975 m的海拔標(biāo)高，最大相對高差為498 m. 三交河礦區(qū)10-203工作面的傾斜長度為241 m，走向長度為1 415 m. 10-203工作面位于中部的10號煤層，煤厚1.6～2.5 m，平均厚2.2 m，煤層底板標(biāo)高1 306～1 402 m，地面標(biāo)高980～1 033 m. 煤層向東傾斜，傾角0°～10°，平均傾角為5°. 研究區(qū)工作面地貌遙感影像見圖2.

圖2 工作面地貌遙感影像圖

2.2 數(shù)據(jù)源

試驗使用的影像數(shù)據(jù)為2018年5月—2018年12月的18景和2019年1月—2020年4月的36景降軌雙極化(HH/HV)Sentinel-1數(shù)據(jù)，其中列舉2018年5月—12月的影像信息為例。影像的下載來自歐空局網(wǎng)站，成像模式為IW條帶模式，分辨率為5 m×20 m，并采用Sentinel-1衛(wèi)星精密軌道數(shù)據(jù)(POD精密定軌星歷數(shù)據(jù))作為輔助。2018年5月—12月影像數(shù)據(jù)的具體時間見表1.

DEM 數(shù)據(jù)的下載來自美國宇航局(NASA)網(wǎng)站，分辨率為30 m×30 m，用來模擬地形相位，剔除干涉圖的地形相位成分。

3 結(jié)果與分析

3.1 試驗結(jié)果

共進(jìn)行了兩組試驗，首次監(jiān)測使用SBAS-InSAR技術(shù)，通過對18景Sentinel-1數(shù)據(jù)的處理，結(jié)合ArcGIS軟件，得到2018年5月—2018年12月地面沉降平均速率圖(圖3)和地面沉降累積形變量圖(圖4).

表1 Sentinel-1數(shù)據(jù)表

圖3 地面沉降平均速率圖

圖4 地面沉降累積形變量圖

從圖3可以看出，該煤礦在2018年5月—12月監(jiān)測的時間段內(nèi)呈明顯沉陷盆地狀，沉陷邊界形變速率達(dá)-0.14 m/yr. 從圖4可以看出，在InSAR可監(jiān)測范圍內(nèi)的研究區(qū)最大形變量累積超過7 cm. 在礦區(qū)中央位置，由于煤炭的開采造成礦區(qū)中央位置產(chǎn)生大量級形變，從而造成失相干現(xiàn)象，導(dǎo)致InSAR監(jiān)測產(chǎn)生空白區(qū)域。通過和礦區(qū)開采工作面位置的對比，可將此范圍定為三交河礦區(qū)的沉陷變形邊界。

第二次也采用SBAS-InSAR技術(shù)對試驗區(qū)域進(jìn)行監(jiān)測，通過對2019年01月01日—2020年04月25日共36期影像進(jìn)行了每月沉降量的監(jiān)測，監(jiān)測結(jié)果見圖5，其中選取12幅去平和濾波后的干涉圖進(jìn)行展示，見圖6.

根據(jù)監(jiān)測成果圖(圖5)與之前半年的沉降累計量圖(圖4)的比對，可以清晰地看到監(jiān)測區(qū)域內(nèi)存在一些沉降盆地，有較大形變發(fā)生。在采區(qū)工作面中央，分析由于形變量級較大導(dǎo)致失相干現(xiàn)象的產(chǎn)生，InSAR技術(shù)未能提取大量級形變點(diǎn)，但是通過InSAR監(jiān)測可以發(fā)現(xiàn)，形變區(qū)域邊界勾畫較為明顯，且可以清晰看出沉降區(qū)域的發(fā)展，該區(qū)域沉降發(fā)展與采礦規(guī)律高度吻合，呈現(xiàn)“S”型生長曲線模式。從圖5可以看出，沉陷邊界每月沉降量達(dá)0.014 m，累加半年的沉降量符合圖4的沉降累計量0.07 m的數(shù)值，也可以驗證InSAR監(jiān)測沉陷邊界的有效性。

圖5 各月沉降量圖

圖6 短基線差分干涉圖

3.2 精度分析

將收集到的2018年5月—2020年4月水準(zhǔn)數(shù)據(jù)與監(jiān)測值進(jìn)行比較，結(jié)果見表2. 其中A1、A5、B1和B4這4個點(diǎn)為研究區(qū)域四周的邊界點(diǎn)，觀測線及邊界點(diǎn)位置見圖7. 9個監(jiān)測點(diǎn)標(biāo)準(zhǔn)差的最大誤差為5.755 mm，最小為0.475 mm，標(biāo)準(zhǔn)差均小于10 mm. 可以確定，SBAS-InSAR技術(shù)通過監(jiān)測沉陷變形邊界中的地表形態(tài)變化能夠很好地達(dá)到一個亞厘米數(shù)量級的精度,進(jìn)而驗證了沉陷變形邊界提取的準(zhǔn)確性。產(chǎn)生這種誤差的主要原因是實際試驗所得到的PS點(diǎn)位置和實際水準(zhǔn)點(diǎn)的位置之間有所偏差,所以采取了通過計算求取各個水準(zhǔn)點(diǎn)周圍足夠多PS點(diǎn)的平均值和監(jiān)測數(shù)據(jù)相互對比的方法來分析。

表2 監(jiān)測結(jié)果對比表

圖7 觀測線及邊界點(diǎn)位置圖

4 結(jié)束語

針對采煤沉陷變形邊界變形監(jiān)測精度的問題，采用SBAS-InSAR技術(shù)精確提取采煤地表沉陷變形邊界。山西三交河煤礦在2018年5月—2018年12月時間段InSAR數(shù)據(jù)處理結(jié)果表明，地表沉陷盆地變形邊界明顯，沉陷變形邊界形變速率達(dá)-0.14 m/yr，在InSAR可監(jiān)測范圍內(nèi)的研究區(qū)最大形變量累積超過7 cm. 在采區(qū)工作面地表中央位置，由于煤炭開采造成大量級形變，產(chǎn)生失相干現(xiàn)象，導(dǎo)致InSAR監(jiān)測的采煤工作面地表范圍只有部分PS點(diǎn)。將2018年5月—2020年4月觀測線實測水準(zhǔn)數(shù)據(jù)與預(yù)測值進(jìn)行比較，9個預(yù)測點(diǎn)標(biāo)準(zhǔn)差的最大誤差為5.755 mm，最小為0.475 mm，標(biāo)準(zhǔn)差均小于10 mm. 實驗表明，SBAS-InSAR技術(shù)提取地表沉陷變形邊界變形值達(dá)到亞厘米數(shù)量級的精度，為精確劃定采煤沉陷區(qū)范圍提供了新的技術(shù)方法。