基于升降軌Sentinel-1A SBAS InSAR的西南山區滑坡隱患早期識別研究

劉排英 賀少帥 王鵬生 孫 亮

1 石家莊鐵路職業技術學院測繪工程系，石家莊市四水廠路18號，050041 2 中國航天空氣動力技術研究院，北京市云崗西路17號，100074

滑坡具有隱蔽性強、破壞程度大、發生頻率高和分布廣泛等特點。貴州省位于中國西南地區，地貌類型主要為高原、山地、丘陵和盆地[1]，氣候溫暖濕潤、多云多雨、植被茂密，這會增大滑坡的隱蔽性和多發性，也會增加光學遙感成像和解譯滑坡的難度。《2020年度貴州省地質災害防治方案》提出“空天地”一體化的地災隱患調查，使用InSAR開展廣域地表形變監測，識別、提取出更多地質災害隱患點。

衛星InSAR利用微波頻段的相位信息探測雷達視線向cm級甚至mm級微小地表形變，具有廣覆蓋、全天候、全天時、高重訪和主動遙感的優勢。由于傳統D-InSAR受時空去相干和大氣擾動等限制，以PS InSAR[2]和SBAS InSAR[3]為代表的時序InSAR應運而生，可以在很大程度上解決D-InSAR的應用限制。在山區環境中，分布式散射體(distributed scatterers，DS)數量遠大于永久散射體數量，因此SBAS InSAR比PS InSAR更適用于山區場景的滑坡形變監測，能夠顯著減小SAR影像的數量需求。

2014年以來，中高分辨率Sentinel-1A數據免費向用戶開放，極大推進了InSAR在我國的普及和應用。針對重大地質災害隱患的早期識別，我國學者相繼提出空-天-地“三查”技術體系和形態-形變-形式“三形觀測”的技術思路，強調InSAR在地質災害監測中的普查作用，使得InSAR在地質災害調查中逐漸得到認可[4-5]。本文以貴州省最近一次滑塌體量大且致災嚴重的雞場鎮滑坡為研究背景，獲取2018-07～2019-07覆蓋貴州省水城縣雞場鎮滑坡發生前的31期升軌和30期降軌Sentinel-1A數據，結合30 m分辨率的SRTM DEM數據，使用SBAS InSAR開展滑坡發生前的地表形變監測和隱患早期識別研究。該研究可分析和評估SBAS InSAR在中國西南山區地表形變監測中的應用效果，為貴州省以及中國西南山區滑坡隱患調查與早期識別提供技術參考。

1 研究區域及數據

貴州省位于中國西南地區，地貌類型以高原、山地、丘陵和盆地為主，其中山地和丘陵占92.5%[1]。貴州省年平均降雨量超過1 000 mm，主要集中在夏季。凹凸不平的地貌類型和時域集中的強降雨特征易引發滑坡、泥石流等地質災害。2019-07-23水城縣雞場鎮發生大型滑坡(圖1白色輪廓區)，該滑坡體屬于構造侵蝕、剝蝕型中山地貌，滑坡在高海拔地帶發生，演變為沿兩條溝道分布的長滑移距離型泥石流，整個滑坡區域長約1 300 m，前后緣高差約460 m，體積約為191.2×104m3，屬于典型的高位遠程滑坡，滑坡體和碎屑流的主滑方向約為26°[6]。

圖1 研究區及Sentinel-1A覆蓋范圍Fig.1 Study area and Sentinel-1A coverage

研究區范圍為104°35′～104°49′E、26°12′～26°27′N(圖1紅框)。實驗采用歐空局Sentinel-1A漸進式地形觀測模式干涉寬幅成像模式VV極化SAR數據集，右視側視角約為39°。覆蓋雞場鎮滑坡的Sentinel-1A數據有升軌第128軌、降軌第164軌，升軌成像時段為2018-07-25～2019-07-20，共31期；降軌成像時段為2018-07-27～2019-07-22，共30期，均為災前影像。SAR數據集空間覆蓋范圍如圖1藍框所示。輔助數據有POD精密定軌星歷數據和分辨率為30 m的SRTM DEM數據，其中POD數據用于精化軌道坐標，SRTM DEM數據用于SAR數據精配準、去地形相位和地理編碼等。

2 研究方法

本研究使用FaSHPS算法[7]統計、提取同質像元作為DS候選點，設定干涉對集的平均相干閾值篩選同質像元作為DS點。由于研究區植被較多，去相干情況較為嚴重，因此使用像對干涉相干性最優且連通的原則生成像對干涉連接圖[8]。使用StaMPS軟件進行SBAS InSAR處理，反演地表形變場。整體技術流程如圖2所示。

圖2 總體技術流程Fig.2 Overall technical flow chart

貴州省水城縣地形起伏較大、植被茂密，跨季影像相互配準難度較大。本文使用POD精密定軌星歷數據精化軌道坐標，計算輔影像相對于主影像的初始偏移量，使精度達到像元級。在此基礎上，使用強度互相關算法從主影像、輔影像中分區塊提取若干同名像點，對每個區塊同名像點方位向和距離向的偏移量進行最小二乘擬合，求得配準多項式系數，配準精度達到0.1像元。配準多項式系數確定后，使用三次卷積內插法將輔影像重采樣到主影像空間上。

假設共有成像日期為[t0,t1,…，tN]的N+1景SAR影像，并且所有SAR影像均已精配準。通過設定時間、空間基線閾值或使用干涉相干性最優且連通的原則，得到M個干涉像對，則干涉像對數M與影像數N+1滿足關系式：

(1)

影像tA和影像tB生成的第i幅干涉圖(tA

為消除Sentinel-1A數據在山地、丘陵等地區的陰影、疊掩現象，本文對研究區成像時段相對一致的升軌和降軌數據分別進行獨立處理，獲取視線向形變速率(正值表示地表向靠近衛星方向發生形變，負值表示向遠離衛星方向發生形變)。根據InSAR成像幾何可知，視線向形變量在水平方向只能分解出東西方向分量，無南北方向分量，這也是衛星InSAR成像幾何的固有缺陷。水平形變量VE-W和垂直形變量VV表達式為：

VE-W=VLOS·sinα，VV=VLOS·cosα

(2)

式中，VLOS為雷達視線向形變速率，α為降軌和升軌入射角。

為了更好地分析坡體形變特征，將VLOS投影至斜坡方向：

(3)

式中，γ為雷達衛星視向線方向和斜坡坡度方向間夾角，詳細公式見文獻[9]。

3 數據處理與結果分析

3.1 SBAS InSAR處理

使用POD精密定軌星歷數據精化Sentinel-1A軌道坐標，生成單視復數(single look complex, SLC)影像，并按照研究區范圍裁剪SLC數據。使用干涉相干性最優且像對連通的原則分別對31期(ID:0～30)升軌影像和30期(ID:0～29)降軌影像進行干涉像對連接，分別生成64個干涉對和68個干涉對(圖3)。

圖3 干涉像對Fig.3 Interference image pairs

對研究區SLC數據進行SBAS InSAR處理：1)選擇升軌、降軌第1期(ID:0)影像作為超級主影像，多視視數設為方位向×距離向=1×4，使用SRTM DEM進行主輔影像精配準、差分干涉處理，計算相干圖、幅度圖，分辨率分別為13.9 m×14.5 m(升軌)和14.0 m×14.5 m(降軌)；2)使用Goldstein算法對差分干涉圖進行濾波處理；3)設定相干系數閾值為0.2進行掩膜，使用MCF算法對差分干涉圖進行相位解纏；4)設定平均相干閾值為0.41，升軌共提取26 038個DS點，降軌共提取50 049個DS點；5)在DS點上建立和解算觀測方程，使用時空濾波器去除大氣相位和高程殘差相位，使用SVD算法解算DS點的時間序列形變；6)使用SRTM DEM對反演的地表形變和幅度圖進行地理編碼。

3.2 結果分析與討論

升軌、降軌Sentinel-1A影像SBAS InSAR反演的2018-07～2019-07研究區地表雷達視線向平均形變速率結果如圖4所示。DS點空間分布存在3種情況：1)西向陽坡DS點多數存在于升軌Sentinel-1A影像上，少量存在于降軌影像上；2)東向陽坡DS點多數存在于降軌Sentinel-1A影像上，少量存在于升軌影像上；3)地勢較為平坦的DS點同時存在于升軌和降軌Sentinel-1A影像上。上述現象是由衛星雷達右視側視成像造成，可通過融合升降軌時序InSAR形變場提升地勢起伏較大山區、丘陵的形變監測能力。

雞場鎮滑坡體(圖4紫色區域)在發生滑移前植被覆蓋度高，DS點較稀疏，上沿及周邊DS點形變速率約為0，說明雞場鎮滑坡發生前滑坡體上沿及周邊區域并未出現大型地表形變，表現為穩定狀態。查詢雞場鎮滑坡發生前兩周天氣狀況可知(表1)，雞場鎮滑塌體在滑坡發生當天有強降雨，在滑坡發生前兩周內有11 d出現不同程度降雨。因此初步判斷雞場鎮滑坡是連續降雨、強降雨誘發的突發性滑坡，已超出Sentinel-1A的12 d重訪周期監測能力，這與Wang等[1]的研究結論一致。

圖4 研究區Sentinel-1A視線向形變結果Fig.4 Sentinel-1A line-of-sight deformation results in the study area

表1 水城縣歷史天氣

在研究區內劃定出A、B、C、D、E五個形變區，主要分布于北盤江兩側斜坡體上(圖4黑框)。將各形變區升降軌視線向形變沿豎直方向、東西水平方向分解，沿坡向合成，結果如圖5所示。豎直形變量向上(抬升)為正值，向下(沉降)為負值；水平形變量向西為正值，向東為負值。各形變區沿坡向形變量及方向示意圖中箭頭線段長度與坡向形變量成正比。由圖5等高線分布情況可知：1)山坡坡度相對較陡處坡向形變量相對較大；2)地勢平坦地區DS點密度明顯大于地勢陡峭處DS點密度；3)地勢平坦地區主要表現為豎直形變，地勢陡峭處同時存在豎直形變和水平形變。5個形變區的潛在威脅目標如圖5所示：1)形變區A。形變分布在斜坡體、采礦區或加工廠，形變可能與地下采礦和礦物加工的抽排水有關，主要威脅對象有工廠、居民住房、學校和道路；2)形變區B。形變分布在斜坡體，形變可能由斜坡體失穩造成，主要威脅對象有居民住房和道路；3)形變區C。形變分布在斜坡體、加工廠，部分斜坡體具有明顯的老滑坡痕跡，形變可能由礦物加工的抽排水或斜坡失穩引起，主要威脅對象有居民住房、工廠、道路和河流；4)形變區D。形變分布在斜坡體，形變可能由斜坡體失穩造成，主要威脅對象有居民住房和道路；5)形變區E為露天采礦區，附近有煤炭倉儲中心、煤炭加工廠和發電廠等，形變可能與采礦和礦物加工的抽排水有關，主要威脅對象有居民住房、道路和附近工廠。

4 結 語

使用SBAS InSAR對升軌、降軌Sentinel-1A數據進行處理，提取貴州省水城縣研究區5個明顯形變區。結合Google Earth推斷引起地表形變的主要原因有斜坡體失穩、地下/露天采礦和礦物加工的抽排水等，潛在威脅對象包括居民住房、學校、工廠和道路等。雞場鎮滑坡發生前升軌、降軌Sentinel-1A的SBAS InSAR形變場并未出現明顯形變信息，且在滑坡發生當天及前兩周內有11 d出現不同程度降雨，因此初步判斷雞場鎮滑坡是由外部因素(如降雨)觸發的突發性滑坡。

SBAS InSAR能有效應用于山區滑坡隱患早期識別和形變監測，但對于外部因素(如降雨)觸發的突發性滑坡，衛星雷達存在重訪周期較長的局限性，如Sentinel-1A為12 d，Sentinel-1A/B雙星為6 d。因此，隨著衛星數量增加，可通過多星組網構建衛星星座，縮短衛星雷達的重訪周期，提升突發性滑坡隱患的早期識別和形變監測能力。本研究在貴州水城縣滑坡隱患調查中所用技術、方法及實驗結論可為貴州省以及中國西南山區滑坡隱患調查與早期識別提供參考。

致謝：感謝英國利茲大學Andy Hooper教授提供StaMPS軟件，感謝中山大學蔣彌教授提供FaSHPS軟件，感謝歐空局提供Sentinel-1A數據，感謝美國航空航天局提供SRTM DEM數據。

圖5 各形變區豎直形變、水平形變、沿坡向合成形變量與方向Fig.5 Vertical and horizontal deformation, deformation and direction along the slope in each deformation area