基于InSAR技術對汕頭潮陽區地質災害體識別分析

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

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

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

1 數據源獲取

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

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

2 數據可視性分析

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

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

3 InSAR數據處理流程

3.1 D-InSAR處理流程

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

圖4 部分全域干涉圖

圖5 部分局部干涉圖

3.2 SBAS-InSAR處理流程

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

圖6 SBAS-InSAR數據處理工作流程圖

3.2.1 數據預處理

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

3.2.2 差分干涉計算

(1)平地和地形相位去除，對相干目標點組成的干涉圖，進行平地和地形相位去除；(2)差分干涉圖濾波；(3)相位解纏。

3.2.3 時間/空間域變形估算

對干涉圖的差分干涉相位應進行時間和空間域線性變形相位估計，如有要求，還應進行非線性變形相位估計，去除大氣、噪聲等殘余相位，得到點目標時間序列變形相位。

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

3.3 關鍵參數設置

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

圖7 影像時間連接圖

3.4 干涉工作流

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

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

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

3.5 軌道精煉和重去平

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

3.6 反演與地理編碼

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

4 解算分析及結果

4.1 解算分析

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

4.2 解算結果

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

本次工作識別出4處斜坡變形，其中2處崩塌，光學影像中可見變形破壞跡象；2處滑坡，變形區與地形較為吻合，地表存在一定變形破壞跡象；其余點位推測均為軟土受地表水與地下水影響形成的區域地表沉降。其中1處存在持續變形跡象，2017年1月至2021年8月存在持續變形，但存在一定季節性特征，約5.5年累計形變量為116mm；崩塌區季節變形明顯，在夏季變形速度增大，曲線波動較大，說明與降雨相關性較高。

5 結 論

本次工作主要采用Sentinel-1進行地表形變探測，由于Sentinel-1為C波段20m分辨率，其分辨率較低，對工程建設開挖導致的地表形變探測效果較好；對于小型滑坡探測能力存在一定限制，而且山區植被茂盛，C波段衛星對植被覆蓋區無法有效性探測。建議采用L波段ALOS-2與機載LiDAR相結合對小型地質災害隱患點進行探測，提高準確性。