基于InSAR技術(shù)對(duì)汕頭潮陽(yáng)區(qū)地質(zhì)災(zāi)害體識(shí)別分析

方 剛，葉 波，紀(jì)振付，張明亮，李建穎，李曉淳

(廣東省有色金屬地質(zhì)局九三一隊(duì)，廣東 汕頭 515041)

合成孔徑雷達(dá)干涉測(cè)量(InSAR)技術(shù)是利用多次重復(fù)觀測(cè)的合成孔徑雷達(dá)獲取地表微小變形，利用具有一定視角差的兩部天線(或一部天線重復(fù)經(jīng)過(guò))來(lái)獲取同一地區(qū)兩幅具有相干性的SAR單視復(fù)數(shù)影像，并根據(jù)其干涉相位數(shù)據(jù)提取地表的高程信息。InSAR技術(shù)可以快速地提供被監(jiān)測(cè)區(qū)域內(nèi)大范圍、準(zhǔn)確的地形與地表形變信息[1]，便于揭示地學(xué)現(xiàn)象的時(shí)空變化規(guī)律，在地質(zhì)災(zāi)害監(jiān)測(cè)與識(shí)別中發(fā)揮了重要作用。

1 數(shù)據(jù)源獲取

本次InSAR地質(zhì)災(zāi)害識(shí)別范圍為廣東省汕頭市潮陽(yáng)區(qū)全區(qū)[2]，總面積666.63km2，重點(diǎn)調(diào)查區(qū)面積103.47km2。選用SAR影像數(shù)據(jù)源為Sentinel-1A雷達(dá)衛(wèi)星數(shù)據(jù)[3]，該衛(wèi)星由歐空局于2014年4月發(fā)射的一顆雷達(dá)衛(wèi)星，與2016年4月發(fā)射的Sentinel-1B衛(wèi)星組成雙星伴飛模式，重訪周期從12天縮短到6天。IW模式采用了最新的TOPs成像技術(shù)，能夠解決寬幅成像時(shí)出現(xiàn)的scalloping效應(yīng)并增強(qiáng)成像輻射性。Sentinel-1衛(wèi)星配備了新一代C波段合成孔徑雷達(dá)系統(tǒng)(波長(zhǎng)5.6cm)，相較于L波段，波長(zhǎng)更短，形變探測(cè)能力更強(qiáng)。

在歐空局官方網(wǎng)站上下載覆蓋研究區(qū)C波段Sentinel-1A衛(wèi)星IW模式下level-1A級(jí)影像，分辨率均為5m×20m，覆蓋研究區(qū)Sentinel-1升軌影像，時(shí)間跨度分別為2017年1月1日至2022年8月15日共168景。采用DEM數(shù)據(jù)主要是30m分辨率的Aster GDEM，其標(biāo)稱高程精度10m～25m，處理Sentinel-1A數(shù)據(jù)前，下載與之對(duì)應(yīng)的精密軌道數(shù)據(jù)并對(duì)軌道信息進(jìn)行修正，有效去除因軌道誤差造成的系統(tǒng)性誤差。

2 數(shù)據(jù)可視性分析

SAR地形可視性取決于SAR采集數(shù)據(jù)的幾何參數(shù)、視線(LOS：light of sight)方向的角度和地形坡度和坡向，根據(jù)這些參數(shù)對(duì)研究區(qū)疊掩和陰影區(qū)分和提取。建立SAR地形可視性分析流程(圖1)，對(duì)Sentinel-1影像數(shù)據(jù)進(jìn)行地形可視性分析(圖2)，結(jié)果表明所選用的影像陰影和疊掩分別較少，僅在山區(qū)少量分布，陰影區(qū)域占0.2%，疊掩區(qū)域占0.01%。

圖1 SAR 地形可視性分析流程圖

3 InSAR數(shù)據(jù)處理流程

3.1 D-InSAR處理流程

該技術(shù)是InSAR形變監(jiān)測(cè)的基本處理方式，利用兩景高分辨率SAR影像可有效進(jìn)行坡體短期厘米級(jí)變形特征識(shí)別，包括滑移面積、形變量信息。處理流程見(jiàn)圖3。判別一定范圍內(nèi)潛在的滑移坡體，要求地質(zhì)災(zāi)害體表面保持良好的相干性[4-6]。兩景高分辨率SAR影像時(shí)間間隔一般1～2個(gè)月，適合短期內(nèi)厘米級(jí)滑坡變形監(jiān)測(cè)。其缺點(diǎn)是對(duì)于長(zhǎng)期毫米級(jí)緩慢形變的滑坡體，差分干涉SAR厘米級(jí)形變監(jiān)測(cè)精度難以識(shí)別，單一監(jiān)測(cè)角度也不足以克服SAR影像疊掩和陰影影響，難以準(zhǔn)確判定長(zhǎng)期毫米級(jí)緩慢形變滑坡的整體范圍、邊界、強(qiáng)度。D-InSAR干涉過(guò)程中生成干涉圖(圖4、圖5)是識(shí)別形變區(qū)域的重要工具，對(duì)所有干涉圖瀏覽，篩選出多處具有明顯干涉條紋的點(diǎn)位。

圖4 部分全域干涉圖

圖5 部分局部干涉圖

3.2 SBAS-InSAR處理流程

SBAS-InSAR技術(shù)是基于D-InSAR技術(shù)基礎(chǔ)上的發(fā)展，其核心是時(shí)空基線閾值設(shè)置，對(duì)于試驗(yàn)區(qū)而言，由于高相干點(diǎn)較少，因此除了需要確定時(shí)空基線閾值外，適宜的空間范圍也尤為重要。通過(guò)對(duì)調(diào)查區(qū)SBAS-InSAR處理[7-10]，反演地表形變速率。具體流程(圖6)。

圖6 SBAS-InSAR數(shù)據(jù)處理工作流程圖

3.2.1 數(shù)據(jù)預(yù)處理

(1)SAR主影像選擇和像對(duì)組合。a)計(jì)算所有影像像對(duì)間時(shí)間和空間基線，生成時(shí)間和空間基線分布圖；b)采用時(shí)間和空間基線均滿足給定閾值的像對(duì)組合生成差分干涉圖集。(2)所有SAR影像對(duì)一景影像進(jìn)行配準(zhǔn)、裁剪，并組合生成時(shí)間序列干涉圖集。a)選擇非夏季、時(shí)空基線盡量居中的影像作為配準(zhǔn)參考影像，所有影像對(duì)其進(jìn)行配準(zhǔn)；b)將所有數(shù)據(jù)裁剪成一致的區(qū)域；c)對(duì)所有配準(zhǔn)好的干涉像對(duì)，按時(shí)間和空間基線限制條件，選擇像對(duì)組合。逐像元計(jì)算干涉相位，生成時(shí)間序涉圖集；(3)將DEM與配準(zhǔn)參考影像進(jìn)行配準(zhǔn)，將DEM范圍裁剪成與配準(zhǔn)參考影像一致區(qū)域；(4)將所有主輔影像前置濾波，計(jì)算干涉相位，生成干涉圖；(5)相干系數(shù)計(jì)算；(6)相干點(diǎn)目標(biāo)選取，對(duì)時(shí)間序列干涉圖集的像元進(jìn)行相干點(diǎn)目標(biāo)篩選。a)相干點(diǎn)目標(biāo)選?。籦)相干點(diǎn)目標(biāo)干涉相位序列生成。將滿足上述條件要求的輔影像與主影像進(jìn)行相位干涉處理，提取相干點(diǎn)目標(biāo)干涉相位序列圖。

3.2.2 差分干涉計(jì)算

(1)平地和地形相位去除，對(duì)相干目標(biāo)點(diǎn)組成的干涉圖，進(jìn)行平地和地形相位去除；(2)差分干涉圖濾波；(3)相位解纏。

3.2.3 時(shí)間/空間域變形估算

對(duì)干涉圖的差分干涉相位應(yīng)進(jìn)行時(shí)間和空間域線性變形相位估計(jì)，如有要求，還應(yīng)進(jìn)行非線性變形相位估計(jì)，去除大氣、噪聲等殘余相位，得到點(diǎn)目標(biāo)時(shí)間序列變形相位。

SBAS-InSAR計(jì)算步驟：a)相鄰點(diǎn)間參數(shù)估計(jì)；b)線性變形相位和殘余高程計(jì)算；c)殘余相位低通濾波。從差分干涉相位中減去步驟：a)中兩項(xiàng)相位分量后得到殘余相位，并進(jìn)行空間域低通濾波得到濾波后的殘余相位；奇異值分解處理是根據(jù)短基線像對(duì)組合關(guān)系，對(duì)步驟b)得到的濾波后殘余相位進(jìn)行奇異值分解(SVD)處理，求解每個(gè)影像對(duì)應(yīng)時(shí)刻的大氣相位和非線性變形相位；大氣相位和非線性變形相位計(jì)算。對(duì)奇異值分解得到的大氣相位和非線性變形相位進(jìn)行空間域高通濾波，得到大氣相位，并對(duì)濾波后的相位序列進(jìn)行時(shí)域低通濾波，得到非線性變形相位；時(shí)間序列變形相位計(jì)算是將b)步驟中線性變形相位和大氣相位和非線性變形相位計(jì)算中非線性變形相位相加，結(jié)合時(shí)間基線參數(shù)，得到每個(gè)相干點(diǎn)目標(biāo)的時(shí)間序列變形相位。

3.3 關(guān)鍵參數(shù)設(shè)置

本次試驗(yàn)采用SBAS-InSAR(小基線集)技術(shù)數(shù)據(jù)處理，其核心是時(shí)空基線閾值的設(shè)置，對(duì)于試驗(yàn)區(qū)而言，由于高相干點(diǎn)較少，因此除了需要確定時(shí)空基線閾值外，適宜的空間范圍也尤為重要。依據(jù)多次試驗(yàn)結(jié)果，最終選擇時(shí)空基線分別為24d和2m，超級(jí)主影像為2017/1/25影像，共634組干涉對(duì)。連接圖如圖7所示，圖中方形點(diǎn)代表超級(jí)主影像，圓形點(diǎn)為其他影像。小基線干涉對(duì)組合圖能清晰地反映超級(jí)主影像與其它影像的相互位置以及圖像參與運(yùn)算的情況。由于本研究采用的SAR影像時(shí)序較多，因此干涉對(duì)組合情況較好，每個(gè)時(shí)相都與其它時(shí)相建立了連接，且連接較為均勻。

圖7 影像時(shí)間連接圖

3.4 干涉工作流

InSAR受時(shí)間過(guò)長(zhǎng)而失相干、幾何位置由于DEM影像或衛(wèi)星影像錯(cuò)位失相干、多普勒效應(yīng)造成的中心頻率失相干和由于溫度產(chǎn)生熱噪聲等失相干影響。干涉配準(zhǔn)后獲得的干涉圖表面含有隨機(jī)噪聲，以灰白色圓點(diǎn)或者直接失相干而失去影像，嚴(yán)重時(shí)會(huì)干擾條紋條帶分布或失去干涉條紋。需要對(duì)干涉圖濾波處理，主要目的是減少相位解纏中粗差。

本次研究數(shù)據(jù)量較多，選用Goldstein濾波方法處理。干涉工作流是對(duì)所有配對(duì)的干涉像對(duì)干涉處理，包括相干性生成、去平、濾波和相位解纏等過(guò)程，最終所有數(shù)據(jù)對(duì)都會(huì)被配準(zhǔn)到超級(jí)主影像上，以便為軌道精煉與重去平、第一次反演和第二次反演做準(zhǔn)備。干涉工作流處理過(guò)程中需設(shè)置多視的視數(shù)，可在SAR影像文件中查看相應(yīng)參數(shù)。估算多視的視數(shù)，若輸入的數(shù)據(jù)分辨率不同，程序會(huì)對(duì)每個(gè)像對(duì)多視處理，得到與超級(jí)主影像一致的分辨率。多視可以增加干涉圖的信噪比，提供可靠的相干性，也可提高運(yùn)算速度。此處設(shè)置為4：1。制圖分辨率為20m，解纏方法為Minimum Cost Flow，解纏相關(guān)系數(shù)閾值為0.2，濾波方法為Goldstein。高程精度閾值為20，速率精度閾值為32。

選取2021年不同時(shí)間段對(duì)相同時(shí)間基線的干涉對(duì)比分析(圖8)。結(jié)果表明，研究區(qū)內(nèi)冬季相干性較高，干涉條紋鮮亮，可以識(shí)別到明顯形變區(qū)域；而植被較為茂盛區(qū)域，不管在任何時(shí)間，相干性都較低(0.2～0.3)，這部分區(qū)域噪聲較多，也是Sentinel-1在植被茂盛區(qū)地表形變識(shí)別的不足。此外，在解譯區(qū)東側(cè)與南側(cè)存在水田區(qū)域，沒(méi)有任何相干性。

3.5 軌道精煉和重去平

在研究區(qū)選擇相對(duì)穩(wěn)定區(qū)域，遠(yuǎn)離形變區(qū)的控制點(diǎn)。在軌道精煉和重去平過(guò)程中輸入GCP點(diǎn)，估算和消除殘余的固定相位以及解纏后還存在的相位殘差。GCP點(diǎn)既可選擇在干涉結(jié)果較好像對(duì)的解纏結(jié)果圖上或者干涉圖上，刺點(diǎn)生成GCP文件，也可選擇外部GCP點(diǎn)。在做地表形變監(jiān)測(cè)時(shí)，不能確定刺點(diǎn)位置為穩(wěn)定區(qū)域，尤其是研究區(qū)面積大、地形復(fù)雜，選擇GCP點(diǎn)位置存在極大的不確定性。本次在做地表形變測(cè)量時(shí)，對(duì)比多次遙感影像，將穩(wěn)定點(diǎn)位置(如穩(wěn)定建筑物、路橋等構(gòu)筑物，不存在形變條紋)作為GCP點(diǎn)，以此降低選點(diǎn)不當(dāng)帶來(lái)的誤差。共選用182處GCP點(diǎn)，保證GCP均勻分布(圖9)，控制點(diǎn)標(biāo)準(zhǔn)差小于1。

3.6 反演與地理編碼

SBAS反演核心就是第一次反演，主要估算研究區(qū)形變速率和殘余地形。第一次反演相關(guān)系數(shù)閾值0.2，解纏相干系數(shù)閾值0.2，解纏方法為最小費(fèi)用流法。第二次反演是獲得時(shí)間序列上形變量。利用之前軌道精煉和重去平步驟中GCPs點(diǎn)文件移除恒定的相位和斜坡相位，根據(jù)第一次反演結(jié)果進(jìn)行大氣濾波處理，最后得到時(shí)間序列位移結(jié)果。第二次反演時(shí)相關(guān)系數(shù)閾值0.3，大于該閾值參與計(jì)算。軌道精煉方法選用多項(xiàng)式精煉方法，且精煉殘余相位多項(xiàng)式次數(shù)設(shè)置為3。

4 解算分析及結(jié)果

4.1 解算分析

利用時(shí)序InSAR結(jié)果可得整個(gè)研究區(qū)域視線向的年均形變速率，Sentinel-1A升軌衛(wèi)星影像年均形變速率圖(圖10)，正值表示靠近雷達(dá)方向，負(fù)值表示遠(yuǎn)離雷達(dá)方向。潮陽(yáng)區(qū)整體形變量較小，大部分區(qū)域形變速率處于-5mm/a～20mm/a，處于平穩(wěn)狀態(tài)，局部區(qū)域存在形變，這和主城區(qū)建成區(qū)范圍廣、人類工程活動(dòng)有顯著關(guān)系[11]。雷達(dá)信號(hào)對(duì)水體不敏感，干涉效果差，表現(xiàn)為無(wú)地表形變，植被和建設(shè)用地均可表現(xiàn)出抬升和沉降現(xiàn)象。山區(qū)部分區(qū)域植被茂盛，造成失相干現(xiàn)象嚴(yán)重。沉降高值區(qū)域主要分布于東北部和南西部城鎮(zhèn)建設(shè)區(qū)域以及地表起伏較大區(qū)域，沿道路分布。從地質(zhì)災(zāi)害識(shí)別的角度出發(fā)，居民地及其周邊地形起伏區(qū)出現(xiàn)明顯地表形變異常，在遙感解譯時(shí)可以重點(diǎn)考慮由于切坡建房引起的斜坡地質(zhì)災(zāi)害。

4.2 解算結(jié)果

據(jù)2017年1月1日至2022年8月15日168景Sentinel-1影像D-InSAR與SBAS-InSAR處理，按照年形變速率大于20mm/a為強(qiáng)變形區(qū)，年形變速率10mm/a～20mm/a為中等強(qiáng)度變形區(qū)，年形變速率小于10mm/a為弱變形區(qū)。共解譯73處疑似地表形變點(diǎn)。大致分為地表開(kāi)挖(工程施工)、地表侵蝕與地表沉降與斜坡變形三類，主要以地表沉降為主；沉降區(qū)多為平原區(qū)河流、建成區(qū)周圍，軟土分布廣泛，受地表水與地下水影響嚴(yán)重，這類變形往往在地表不出現(xiàn)明顯變形，而以區(qū)域變形為主。

本次工作識(shí)別出4處斜坡變形，其中2處崩塌，光學(xué)影像中可見(jiàn)變形破壞跡象；2處滑坡，變形區(qū)與地形較為吻合，地表存在一定變形破壞跡象；其余點(diǎn)位推測(cè)均為軟土受地表水與地下水影響形成的區(qū)域地表沉降。其中1處存在持續(xù)變形跡象，2017年1月至2021年8月存在持續(xù)變形，但存在一定季節(jié)性特征，約5.5年累計(jì)形變量為116mm；崩塌區(qū)季節(jié)變形明顯，在夏季變形速度增大，曲線波動(dòng)較大，說(shuō)明與降雨相關(guān)性較高。

5 結(jié) 論

本次工作主要采用Sentinel-1進(jìn)行地表形變探測(cè)，由于Sentinel-1為C波段20m分辨率，其分辨率較低，對(duì)工程建設(shè)開(kāi)挖導(dǎo)致的地表形變探測(cè)效果較好；對(duì)于小型滑坡探測(cè)能力存在一定限制，而且山區(qū)植被茂盛，C波段衛(wèi)星對(duì)植被覆蓋區(qū)無(wú)法有效性探測(cè)。建議采用L波段ALOS-2與機(jī)載LiDAR相結(jié)合對(duì)小型地質(zhì)災(zāi)害隱患點(diǎn)進(jìn)行探測(cè)，提高準(zhǔn)確性。