高山峽谷區滑坡災害隱患InSAR早期識別
——以雅礱江中段為例

戴可人 鐵永波 許 強 馮 也 卓冠晨 史先琳

①(地質災害防治與地質環境保護國家重點實驗室(成都理工大學) 成都 610059)

②(成都理工大學地球科學學院 成都 610059)

③(中國地質調查局成都地質調查中心 成都 610081)

1 引言

雅礱江地處青藏高原向四川盆地過渡地帶，屬青藏高原東部邊緣侵蝕山原區強烈-極強烈山原亞區，屬嵌入式河流，溝谷深切，溝壑縱橫，相對高差達2000～3000 m，河谷為“V”形谷，屬典型的高山峽谷地貌[1]。雅礱江沿岸巖層風化強烈，切割破碎，岸坡穩定性差，加之沿岸水電的頻繁開發以及大中型工程的興建，該流域所遭受的滑坡、崩塌及泥石流等地質災害也越來越嚴重[1—4]，對當地人民生命財產及水利水電等基礎設施安全帶來嚴重威脅。在雅礱江流域沿岸高山峽谷區域開展滑坡災害隱患早期識別，對該區域的地質災害防范以及確保人民生命財產與水利水電等基礎設施安全具有重要的意義。

由于需要監測面積較大，傳統地質調查手段在滑坡災害隱患排查工作中很難達到大范圍覆蓋。隨著對地觀測技術的不斷進步，光學遙感解譯是目前滑坡災害隱患大范圍早期識別主要手段之一，但其受云霧影響較大，且往往只能對較大地貌特征進行定性解譯[5]。作為近三十年的新興雷達遙感技術，合成孔徑雷達干涉測量技術(Interferometric Synthetic Aperture Radar,InSAR)[6]具有覆蓋范圍廣、監測精度高、全天時、全天候、空間分辨率高等特點，已被成功用于滑坡災害監測(如文獻[7—11])，近年來被逐步應用于災害隱患早期識別并取得一些成功案例。張毅[12]基于Envisat ASAR數據，應用時間序列InSAR技術成功識別出白龍江流域133處活動斜坡；張路等人[13]利用InSAR技術成功識別出大渡河區域17處持續變形中的不穩定坡體；張亞迪等人[14]利用時序InSAR對芒康地區滑坡災害隱患進行了探測；Dai等人[15]利用短基線DInSAR技術對汶川附近國道G317沿線不穩定滑坡體進行了有效識別；陸會燕等人[5]結合InSAR與光學遙感技術，成功探測出金沙江下游區域共計7處具有較顯著形變的滑坡隱患。高山峽谷區域海拔落差大、地形復雜、植被茂密，給InSAR處理帶來干涉失相干、大氣延遲、幾何畸變等挑戰[16]。對于InSAR高山峽谷區域的滑坡識別效果與適用性還值得進一步分析與研究。

本文利用小基線集時間序列方法(Small BAseline Subset InSAR,SBAS-InSAR)，對雅礱江流域雅江縣至木里縣段高山峽谷區域進行了滑坡災害隱患廣域早期識別的研究，探究其潛在的滑坡隱患區域，并對識別結果進行進一步的分析討論，分析雷達遙感在高山峽谷區滑坡探測中的效果與適用性，可為當地的防災減災提供有力的支持，并為高山峽谷區的滑坡災害隱患早期識別提供思路與參考。

2 研究區域與數據—研究區域位置與地勢特點

雅礱江是金沙江最大支流，為橫斷山區北南向的主要河系之一[1]。本文選取雅礱江流域雅江縣-木里藏族自治縣段作為研究區域(如圖1所示)。該區域位于雅礱江中游段，四川省甘孜藏族自治州、涼山彝族自治州境內，處于中國最高一級階梯向第2級階梯云貴高原和四川盆地過渡地帶，屬橫斷山系北段川西高山高原區，亦屬于青藏高原的一部分。該段位于甘孜阿壩褶皺帶，出露地層以三疊系淺變質砂板巖為主，有少量花崗巖及石灰巖分布。分布巨厚的中上三疊系淺變質巖系，砂巖，板巖構成北西-南東向緊密褶皺，褶皺軸部及斷層帶中有少量二疊系灰巖分布，并有零星燕山期花崗巖出露，且出露下古生界至上古生界碳酸鹽巖類、淺變質巖及玄武巖等[1—4]。研究區域海拔3000～5000 m，峰巒重疊、溝壑縱橫、河谷深切，為典型高山峽谷地貌，為崩塌、滑坡等斜坡地質災害的發生孕育了條件。大部分山區溝道巨大的地形高差，使處于高處的松散碎屑物質擁有巨大的勢能，陡急的山坡和溝床為坡面和溝床松散堆積物能量的釋放和勢能轉化為動能提供有利條件，有利于滑坡、泥石流等災害的形成，導致該區域內滑坡、崩塌、泥石流等地質災害頻發[17—20]。

圖1 研究區域Fig.1 Study area

本研究共計使用了50景Sentinel-1衛星數據(每一個時間點需要兩景數據上下拼接)，時間范圍覆蓋2017年9月9日—2018年10月10日共計396天。Sentinel-1是歐洲航天局(European Space Agency,ESA)發射的由兩顆衛星組成地球觀測衛星星座，載有C波段合成孔徑雷達，可提供白天、夜晚及各種天氣下的連續影像。表1列出了Sentinel-1衛星主要參數。

表1 Sentinel-1衛星SAR影像數據主要參數Tab.1 Main parameters of Sentinel-1 SAR datasets

Sentinel-1影像數據集的時空基線圖如圖2所示，時間基線閾值為48 d，空間基線閾值小于150 m。為了消除與減弱由軌道誤差引起的相位誤差，在數據處理過程中還使用了歐洲空間局精密軌道數據對軌道信息進行糾正。同時采用日本宇航局ALOS WORLD 3D 30 m空間分辨率的數字高程模型(Digital Elevation Model,DEM)，以消除或減弱地形相位的影響。

圖2 SAR數據集時空基線圖Fig.2 Spatial and temporal baselines of SAR datasets

3 時間序列InSAR方法

傳統的差分干涉測量技術(Differential InSAR,DInSAR)易受到時間失相干和空間失相干的影響，不能準確獲取形變信號。為降低時空失相干的影響，Berardino等人[21]提出了短基線集，即SBASInSAR技術，其核心思路是將同一個地區多幅SAR影像配準得到短基線對，采用奇異值分解法(Singular Value Decomposition,SVD)將多個短基線集聯合求解，得到目標區域的整個時間位移序列以及地表形變平均速率。具體原理及流程如下：

假設獲取了覆蓋同一區域N+1幅SAR影像，按時間序列排序如式(1)

按照一定規則，選取某一影像作為主影像進行配準，得到M個干涉對，M滿足式(2)

通過設定時空閾值，得到相干性較好的M對干涉對，再通過精密軌道文件及高精度數字高程模型去除軌道誤差、平地效應及地形相位的影響后得到M幅差分干涉圖以及解纏后的形變相位。

假設t0時刻作為起始時間，將該時刻研究區域的位移視為0，則第i幅干涉圖中(1≤i≤m)的某個像素相對于起始點的相位可以表示為式(3)

其中，φitopo表示地形相位，φiatm大氣延遲造成的相位，φinoise相干噪聲引起的相位。分別可表示為

λ為波長；θ為雷達入射角；Δh為外部DEM誤差；r為雷達到目標物體的斜距；d(t2)與d(t1)分別表示t2與t1時刻像元相對于起始時間t0在雷達視線方向的形變累積量。假設干涉對影像時間間隔內地表形變滿足線性變化，將對相位時間序列的求解轉變為相位變化速率的求解，可得式(5)

則Δφdef可表示為：Δφdef=B，其中，矩陣B為系數矩陣，每行一一對應干涉像對。矩陣的元素中主影像系數為1，輔影像系數為—1，其他系數為0。

因此，在不考慮大氣相位與噪聲的情況下，聯合式(1)—式(5)可得式(6)

圖3 SBAS-InSAR時序分析流程圖Fig.3 Flowchart of SBAS-InSAR time series analysis

4 滑坡災害隱患早期識別結果

4.1 總體識別結果

通過短基線集SBAS-InSAR時間序列方法對雅礱江中段雅江縣-木里縣流域進行了形變監測。考慮到數據處理過程誤差的累積與處理的效率，將雅江縣-木里縣流域分為5段，分別提取了各段的雅礱江沿岸坡體形變特征信息，獲得的雷達視線方向上的年均形變速率如圖4所示，其中紅色負值代表的是目標地物遠離衛星方向位移，藍色正值代表的是目標地物靠近衛星方向位移。5段的詳細劃分如圖4(a)所示，具體劃分為雅江縣-惡古鄉(A段)，惡古鄉-牙衣河鄉(B段)，牙衣河鄉-三巖龍鄉(C段)，麥地龍鄉-田鎮村(D段)，卡拉鄉-三桷埡鄉(E段)。本文主要針對雅礱江沿岸對河流有危險的坡體進行In-SAR結果解譯，早期識別地質災害隱患，并進行相干性、幾何畸變等復核檢驗，共識別出8處正在發生蠕變的不穩定坡體，其中B段4處，即魯日、日阿、日衣、木恩；D段共中鋪子村、麻撒村2處；E段共陽山村、獨家村2處。基于SAR成像參數與坡度坡向數據，我們對8處坡體進行了幾何畸變分析與驗證(依據5.2節、圖8與圖9)，這8處坡體均處于非幾何畸變區，形變監測結果可靠。

本文從變形范圍、最大形變速率、坡體平均坡度、高程范圍、植被覆蓋情況、威脅對象這幾個方面，對這8處正在變形區域進行了詳細分析。根據隱患點的形變量級，坡體是否具備滑坡地質地貌特征、形變速率是否有加速、突變現象、是否直接對雅礱江或附近村落直接造成威脅作為判斷依據，對風險等級進行了定性評估。其中5處隱患點處于高風險等級，倘若發生地質災害，極大可能會對當地群眾人身安全造成威脅，屬于重點觀測點值得特別關注。這5處主要分布于B段惡古鄉-牙衣河鄉與D段麥地龍鄉-田鎮村區域。惡古鄉-牙衣河鄉有3處風險等級高的區域，其中魯日坡體，對魯日村及村中的雅江縣八衣絨鄉中心小學、下次呷、次呷和上次呷的安全造成嚴重的威脅，并且具有堵塞雅礱江的風險；日阿坡體對附近的日阿村、科魯村、道路與雅礱江構成直接的威脅；日衣坡體對日衣村與雅礱江構成直接的威脅。麥地龍鄉-田鎮村有2處風險等級高的區域，其中中鋪子村坡體對麥地龍鄉、中(鋪子村與雅礱江構成直接的威脅；麻撒村坡體對附近的村落與雅礱江構成直接的威脅。詳細結果列表如表2所示。從整體識別結果來看，InSAR技術在植被覆蓋較低的區域，可以精確有效地識別出高山峽谷區滑坡災害隱患區域，所發現的幾處隱患點均為中低植被覆蓋區。就此研究區域而言，滑坡災害隱患區域最大形變速率普遍超過50 mm/year，且呈現出高位特征，平均高程大于2000 m，平均坡度大于55%，均可對當地村落與雅礱江構成直接的威脅。

圖4 InSAR年均平均速率監測結果Fig.4 Mean velocity map derived from InSAR

4.2 重點區域分析解譯

魯日村與中鋪子村是此次識別中具有直接威脅對象的兩處重點區域，其所在的坡體年均形變速率如圖5、圖6所示，圖中圓點代表監測InSAR有效相干點，點顏色代表年平均形變速率，其中綠色代表較弱的年平均形變速率，亮色彩(黃色、紅色)代表顯著的遠離衛星視線運動，具體量級如圖中顏色條所示。

魯日村所在坡體位于雅礱江東岸(如圖5(a)所示)，坡體植被較少，整體相干性較好。其在地形地貌上呈現階梯狀特征(如圖5(b)所示)，其中魯日村位于海拔2986～3085 m，上方存在高位的坡體；下次呷、次呷和上次呷位于魯日村坡體的下方。根據形變監測的結果表明，魯日村的形變主要位于魯日村以及魯日村至下次呷之間的坡體(圖5(a)、圖5(c)中的藍色虛線所示)。若發生滑坡災害，潛在的滑坡方向如圖5(b)中紅色箭頭所示，可能對魯日村及村中的雅江縣八衣絨鄉中心小學、下次呷、次呷和上次呷的安全造成嚴重的威脅，并且滑坡松散堆積物質可以直接經下次呷沖入雅礱江，具有堵江的風險。

表2 雅礱江潛在災害點早期識別結果列表Tab.2 Early identification results on potential disaster sites along Yalong River

圖5 魯日村滑坡災害隱患識別結果Fig.5 Early identification results of potential landslide geohazards in Luri village

圖5(c)展示了3處特征點的時間序列形變曲線圖，其中點P1位于魯日村，點P2位于下次呷上方坡體，點P3位于西側坡體。時序結果表明，3處累積形變量有所差異，但其趨勢相似。下次呷(P2點)的累積形變最大，達到58 mm；魯日村(P1點)累積形變次之，達到51 mm；西側坡體(P3點)累積形變最小，達到44 mm，可見坡體整體都存在明顯的形變。3個特征點的形變速率自2017年9月—2018年2月與2018年8月—2018年10月較緩，2018年2月—2018年8月較快。通過中國氣象數據中心查詢甘孜州地面觀測站發現，當地氣溫在3月—10月月氣溫溫度大于0 ℃且為當地雨季，氣候變化可能是造成形變加速的主要自然原因。

中鋪子村所在坡體位于雅礱江西岸，坡頂高程3000 m，坡底高程1900 m，平均坡度約60.4%，屬于高位陡峭坡體。如圖6(a)所示，由于坡體中上部植被較少，相干性較好，因此坡體整體中上部區域相干點的數量相對較多，但在河底部位，由于坡度更大，達65%，相干性點的數量較稀疏。坡體中部出現一處核心形變區，其范圍約1700×1000 m。最大視線向形變量達60 mm/year。該區域的遙感解譯結果如圖6(b)所示，在兩處紅色虛線范圍內存在古滑坡體，坡度較大，達到約65%，是坡頂存在不穩定的直接表現之一。在坡體的北側，更是存在巴多溝崩塌、滑坡體，間接證明該坡體長期存在不穩定性。

圖6(c)展示了3個特征點的時間序列形變曲線圖，其中點P1位于坡體上部，點P2位于坡體下部左側，點P3位于坡體中部。時間序列結果表明，3處形變其趨勢相似，3個特征點的形變速率自2017年9月—2018年10月，形變量持續增加，形變速率穩定，呈現線性形變。P1,P2點的累積形變最大，達到120 mm；P3點累積形變次之，達到55 mm，整體都存在明顯的形變。人類活動主要集中在江對岸，所以植被的稀疏與雨季降雨量的增大影響是造成坡體有滑坡趨勢的主要原因。

圖6 中鋪子村滑坡災害隱患識別結果Fig.6 Early identification results of potential landslide geohazards in Zhongpuzi village

5 識別結果分析驗證與雷達影像適用性分析

5.1 識別結果對比驗證

根據蠕變滑坡位移曲線特征[9]，InSAR形變結果是潛在滑坡災害點早期識別的重要特征之一。這些潛在滑坡災害點部分是發育于古滑坡上的不穩定坡體，部分滑坡表面不具備滑坡繼續發展的地貌特征，還需要對這些坡體進行進一步驗證與分析。在得到這些災害隱患點后，接下來需要對InSAR早期識別結果與滑坡隱患歷史資料(含地質調查資料、文獻、滑坡發生報道等)、光學影像遙感解譯結果進行對比(圖7)。

雅礱江流域各個監測段共探測到形變坡體8處依次由北向南編號為1～8號(表3)。其中共有4處相對歷史記錄為新發現(占識別結果50%)，標記為新的潛在地質災害點，其中有7處具有明顯滑坡地貌特征(占識別結果87.5%)。我們以前述變形范圍、最大形變速率、坡體平均坡度、高程范圍、植被覆蓋情況、威脅對象等作為指標，對這8處正在變形區域進行了危險等級快速定性評估，其中3處中危險等級潛在災害點，5處高危險等級潛在災害點。高危險等級潛在災害點具備滑坡地貌地形特征，形變速率較大，直接威脅雅礱江或相近村落，需要對這些識別出的災害點進行進一步調查及監測。

5.2 雷達影像幾何畸變分析

圖7 潛在地質災害點驗證與對比Fig.7 Verification and comparison on potential landslide geohazards

由于星載SAR采用側視成像的方式，雷達波束斜向照射地表時會導致雷達圖像出現距離向透視收縮、陰影或倒置等幾何畸變[16]，在高山峽谷區域幾何畸變尤為嚴重，將會造成無效的監測盲區[12,13]。雷達衛星在面對不同朝向、坡向時所會形成的各種幾何畸變情況如圖8所示。可以看出當坡面朝向衛星時，若坡角小于入射角時，沿坡面向下的變形在LOS方向上會表現為靠近衛星(圖8情況①)，且坡體會發生透視收縮。當坡角與入射角相等時，會發生完全透視收縮，坡體無法被測量(圖8情況②)。如果坡角過大(超過衛星入射角)，則會發生頂底倒置的成像，沿坡面向下的變形在LOS方向上會表現為遠離衛星(圖8情況③)。當坡面背向衛星時，若坡角較大(大于入射角的余角)，此時坡面會位于陰影區域，無法被衛星照射產生回波信號，整個坡面無法被測量(圖8情況④)。當坡角與入射角互余時，此時所獲得地距分辨率等于衛星斜距分辨率(圖8情況⑤)。若坡角較小(小于入射角的余角)，此時將獲得最接近于斜距分辨率(slant range resolution)的地距分辨率(ground range resolution)[23]，此時的地距分辨率將高于衛星照射平地情況，因此該朝向及坡度的坡體是最適合SAR衛星進行觀測的情況。在實際應用中，聯合升降軌的SAR影像，可以獲取不同成像幾何下的監測目標的形變特征，從而在一定程度上補償單一成像幾何帶來的觀測盲區，有利于提升滑坡災害隱患早期識別的有效探測率[13,16,24]。

表3 雅礱江流域InSAR早期識別驗證與對比Tab.3 Verification and comparison on early identification results from InSAR along Yalong River

圖8 衛星入射角(LOS)觀測方向與沿坡向形變關系(修改自文獻[22])Fig.8 Relationships between the Line Of Sight (LOS) and the downslope displacements for different slope orientations (adapted from Ref.[22])

從定量關系上來看，根據Sentinel-1衛星軌道參數信息與坡體坡度信息進行測算，Sentinel-1升軌數據的衛星飛行方向是自南向北，沿著方位向—12.6°飛行，雷達入射角為36.8°。由于SAR采用側視成像，當坡體朝向衛星的時候，若坡度角處于0～36.8°時，會出現透視收縮(距離壓縮)的幾何畸變；當坡體背向衛星的時候，若坡度角的絕對值大于53.2°時，會出現陰影的幾何畸變。Sentinel-1降軌數據的衛星飛行方向是自北向南，沿著方位向—167.4°飛行，雷達入射角為39.7°。由于SAR采用側視成像，當坡體朝向衛星的時候，若坡度角處于0～39.7°時，會出現距離壓縮的幾何畸變；當坡體被向衛星的時候，若坡度角的絕對值大于50.3°時，會出現陰影的幾何畸變。

圖9為基于升降軌Sentinel-1影像獲取的雅礱江升降軌幾何畸變圖，紅色代表的是距離壓縮的區域，藍色代表的是陰影區域。在雅礱江流域中，升軌數據中的距離壓縮區域主要位于東南朝向的坡體，且主要分布在雅礱江的東岸；降軌數據中的距離壓縮區域主要位于東北朝向的坡體，且主要分布在雅礱江的西岸。陰影區域主要分布在較為陡峭的背向雷達信號側的坡體。

圖10(a)、圖10(b)分別為雅礱江日衣村-木灰村段Sentinel-1升軌數據、降軌數據幾何畸變圖，圖10(c)、圖10(d)分別為該區域坡度、坡向圖。在升軌數據圖10(a)中，紅色區域相對衛星的入射角處于0°～36.8°，主要分布在坡面②與坡面④，代表的是距離壓縮區域，對應圖10(d)中雅礱江東岸褐色的區域；藍色區域的入射角大于53.2°，代表的是陰影區域，主要分布在坡面①與坡面③，對應圖10(c)中紅色的大坡度的區域。

在降軌數據對應的圖10(b)中，紅色區域相對衛星的入射角處于0°～38.7°，主要分布在坡面①與坡面③，代表的是距離壓縮區域，對應圖10(d)中雅礱江西岸綠色的區域；藍色區域的入射角大于50.3°，代表的是陰影區域，主要分布在主要分布在坡面②與坡面④，對應圖10(c)中紅色大坡度的區域。

圖9 雅礱江流域Sentinel-1數據幾何畸變分布圖Fig.9 Geometric distortion of Sentinel-1 datasets along Yalong River

圖10 雅礱江日衣村-木灰村段Fig.10 Yalong River Riyi village-Muhui village section

可以明顯地看出，在坡向較復雜的區域，易出現嚴重的距離壓縮幾何畸變。在坡度較大的區域，易出現陰影區域，無法獲取有效的監測信號。就雅礱江日衣村-木灰村區域而言，升軌數據在雅礱江東岸存在嚴重的距離壓縮幾何畸變，主要分布在坡面②與坡面④；在雅礱江高坡度的陡坡主要分布在坡面①與坡面③)，可能出現嚴重的陰影區域；降軌數據在雅礱江西岸存在嚴重的距離壓縮幾何畸變，主要分布在坡面①與坡面③；在雅礱江高坡度的陡坡(主要分布在坡面②與坡面④)，可能出現嚴重的陰影區域。因此，在現有SAR衛星系統飛行軌道差異性有限的情況下，高山峽谷區域單軌道監測可能會由于SAR成像幾何畸變造成部分不穩定坡體被“漏檢”，只有通過升降軌數據結合的方式，才能較為全面覆蓋監測區域，在一定程度上補償單一成像幾何帶來的幾何畸變問題，才能實現全面準確的滑坡災害隱患早期識別。

6 結論

本文使用Sentinel-1雷達影像，通過時間序列InSAR技術對高山峽谷區域雅礱江流域雅江縣-木里縣段進行了滑坡災害隱患廣域早期識別，成功探測到8處隱患區域，最大雷達視線向年平均形變速率可達95 mm/year，這些區域若發生滑坡，對當地居民的生命財產安全存在嚴重的威脅，同時存在堵塞雅礱江的風險。同時，利用滑坡隱患歷史資料與遙感解譯對InSAR早期識別結果進行了進一步的驗證與對比分析，綜合變形范圍、最大形變速率、坡體平均坡度、高程范圍、植被覆蓋情況、威脅對象等情況對這八處滑坡進行了危險等級評價。

由于星載SAR采用側視成像的方式，不可避免的會產生幾何畸變，在高山峽谷區尤為嚴重。本文對不同坡度坡體幾何畸變情況進行了詳細分析，展示了Sentinel-1升降軌數據在雅礱江高山峽谷區的具體幾何畸變區域與適用性。在地形復雜的高山峽谷區域，只有通過采用不同軌道數據結合的方式補償單一成像幾何帶來的幾何畸變問題，才能有效避免“漏檢”，實現全面準確的滑坡災害隱患早期識別。本文以雅礱江為例，揭示了時序InSAR技術運用于滑坡災害隱患早期識別的有效性與優勢，可為當地的防災減災提供有力的支持，結合幾何畸變的詳細分析討論，為高山峽谷區的滑坡災害隱患廣域早期識別提供范例與參考。