基于MSBAS技術的金沙江上游色拉滑坡形變分析

熊國華，楊成生，朱賽楠，董繼紅，張 勤

（1.長安大學地質工程與測繪學院，陜西 西安 710054；2.中國地質環(huán)境監(jiān)測院（自然資源部地質災害技術指導中心），北京 100081）

0 引言

滑坡堵江事件是斜坡巖土體因崩塌、滑坡及其轉化為碎屑流等而造成江河堵塞和回水的現(xiàn)象[1]。我國滑坡堵江事件多發(fā)生在西南山高坡陡的峽谷地區(qū)，這些地區(qū)經歷強烈隆升運動后形成了地形起伏變化巨大的地形地貌。在地震、暴雨、人工開挖、河流強烈快速下切、庫區(qū)水位驟然升降等諸多因素的作用下，河谷兩側時常發(fā)生大規(guī)模堵江滑坡事件。相比于單一滑坡災害而言，堵江滑坡的危險性通過時空上的延拓而大幅增加[2?3]。近年來，國內外許多學者對堵江滑坡做了大量的研究。柴賀軍等[4]分析了堵江滑坡在時間和空間上的發(fā)育規(guī)律，對我國堵江滑坡的區(qū)域分布及時間分布規(guī)律進行了總結。許強等[5]運用衛(wèi)星遙感、無人機航拍、地面合成孔徑雷達監(jiān)測等技術手段對岷江支流四川新磨滑坡進行綜合分析，揭示了滑坡的運動過程和成因機制，并對滑坡周邊直接受主滑坡動力作用而產生的欠穩(wěn)定巖土體特征和危險性進行了分析評價。劉傳正等[6]采用多因素賦值統(tǒng)計研判對西藏林芝市米林縣雅魯藏布江左岸色東普溝上游發(fā)生的滑坡進行監(jiān)測，預測出其可能引發(fā)崩滑-碎屑流造成雅魯藏布江再次堵江的臨界條件。王群等[7]運用永久散射體干涉測量技術（Permanent Scatter InSAR，PS-InSAR）和偏移量跟蹤技術（Offset Tracking）對金沙江白格滑坡形變進行監(jiān)測，獲取白格滑坡在發(fā)生滑動前的運動特征等。朱賽楠等[8]運用差分合成孔徑雷達干涉測量技術和偏移量跟蹤技術（Offset Tracking）分析了金沙江上游的色拉滑坡在白格滑坡發(fā)生前后的形變特征，并結合地質結構、地層巖性及降雨等分析了滑坡的形成機理。由于堵江滑坡發(fā)生的地方隱蔽，多是發(fā)生在深切峽谷中，占比達到94.1%[9]，而傳統(tǒng)的監(jiān)測手段很難做到大范圍、全天候的監(jiān)測。近年發(fā)展起來的時序InSAR技術（MT-InSAR）彌補了部分常規(guī)滑坡識別與監(jiān)測手段的不足，具有對微小形變敏感的獨特優(yōu)勢，成為滑坡識別與位移測量的重要手段之一[10?14]。

野外調查顯示色拉滑坡近期活動明顯，且具有較高的堵江風險[15]。文章在收集覆蓋色拉滑坡區(qū)域2018年1月—2020年4月的Sentinel-1A/B衛(wèi)星升降、軌數(shù)據的基礎上，采用相位堆疊技術（Stacking InSAR）分析了滑坡體的形變特征，然后使用多維小基線集技術（Multidimensional Small Baseline Subset，MSBAS InSAR）獲取了滑坡的二維動態(tài)形變序列，成功捕獲了滑坡體形變加速的時間點，同時結合降雨數(shù)據開展了長時序形變監(jiān)測及形變原因研究，為金沙江流域堵江滑坡的監(jiān)測預警研究提供一定的參考。

1 研究區(qū)域及數(shù)據

1.1 研究區(qū)概況

色拉滑坡位于西藏自治區(qū)昌都市貢覺縣沙東鄉(xiāng)金沙江右岸，總體地勢由西北向東南傾斜（圖1）。該滑坡距上游葉巴灘水電站23 km，距下游拉哇水電站約63 km。色拉滑坡前緣高程約2 649 m，后緣高程約3 342 m，相對高差達693 m，滑坡后緣至斜坡后緣陡緩交界處，左側以沖溝為界，右側以山脊基巖出露處為界，前緣至基巖出露處?；抡w坡向132°，總體呈陡坡地貌[8,15]。光學影像圖1（c）顯示色拉滑坡總體形態(tài)近似呈舌狀，上部較窄，下部寬，前緣已有明顯變形，中部有拉張裂縫，呈現(xiàn)灰褐色、灰白色。從光學影像上可以識別出滑坡壁，滑坡拉張裂縫、變形區(qū)等信息。

圖1 研究區(qū)數(shù)據覆蓋范圍示意圖Fig.1 Schematic image coverage of the study area

色拉滑坡在構造上位于歐亞大陸西南活動大陸邊緣中段、即松潘—甘孜陸緣活動帶和昌都陸塊，構造線由北西轉向近南北向，區(qū)內斷層和褶皺發(fā)育，近東西向的色協(xié)龍斷裂和近南北向的洛冷登—巴巴斷裂在滑坡東北方向交匯[8]?；聟^(qū)出露的地層主要為第四系滑坡堆積層，第四系殘坡積層及二疊系崗托巖組片麻巖。根據地質資料顯示，滑坡區(qū)主要為第四系滑坡堆積層，巖性主要為碎石塊土，碎塊石母巖為片麻巖?；马敳亢蛢蓚葏^(qū)域為第四系殘坡積層，巖性主要為粉質黏土夾碎石。在滑坡前緣和右側山脊處出露部位的巖性主要為二疊系崗托巖組片麻巖，區(qū)域內巖體糜棱巖化和蝕變作用嚴重。

1.2 SAR數(shù)據

文章使用的SAR影像數(shù)據來自歐空局（European Space Agency）2014年4月3日發(fā)射的對地觀測衛(wèi)星Sentinel-1A/B(哨兵數(shù)據)，數(shù)據覆蓋區(qū)域如圖1（a）所示。本文所使用的數(shù)據時間跨度為2018年1月—2020年4月，分別來自于升軌Track 99和降軌Track 33，數(shù)據詳細信息見表1。采用了美國地質調查局(USGS)提供的空間分辨率為30 m的SRTM DEM數(shù)據，作為外部數(shù)據來消除地形相位的影響。同時通過POD精密定軌星歷數(shù)據(POD Precise Orbit Ephemerides)來糾正干涉組合中的基線誤差。

表1 研究區(qū)SAR數(shù)據集主要參數(shù)Table 1 Main parameters of the SAR data sets used in this study

2 研究方法

2.1 相位堆疊InSAR技術

相位堆疊InSAR技術（Stacking InSAR）是對多幅差分干涉圖進行加權平均，獲取平均形變速率的一種方法[16?18]。其基本假設為：在每一幅干涉圖中，地表形變近似看作線性變化，殘余大氣誤差是隨機變化。在進行相位堆疊之前，通過對差分干涉圖進行挑選，以相干性和解纏結果為標準，挑選相干性且解纏結果較好的干涉對。在進行相位堆疊（Stacking InSAR）后，疊加的形變相位為對應時間段的累積形變相位，而疊加的大氣相位并沒有通過干涉圖數(shù)量的增加而不斷累加，而是與干涉圖的數(shù)量成平方根倍增長。由此使得形變相位在疊加圖上占的比例較大，而大氣誤差所占貢獻較小，達到抑制大氣噪聲和DEM誤差，提高形變信息對大氣相位的相對精度的作用[19?20]。

Stacking InSAR技術可以在數(shù)據較少情況下，對多幅干涉圖進行疊加，獲取累積形變量。其中，每幅干涉圖的相位變化速率（Vi）的標準差、干涉相位和單個干涉組合的時間基線（?Ti）之間的關系為：

把時間基線作為權重因子，求取相位形變速率（ph_rate）：

式中：wi——第i幅干涉圖的權，wi=

phi——單個干涉圖的解纏相位值。

從而可以獲取年平均形變速率，其標準差為[21?23]:

相位的標準差為：

2.2 多維短基線集技術

多維小基線集（Multidimensional Small Baseline Subset, MSBAS）技術是由SAMSONOV等[24]提出的多維形變時序分析方法。多維小基線集-二維（MSBAS-2D）技術是基于不同的入射角與方位角，選擇覆蓋公共區(qū)域的影像，依據其幾何關系，對視線向形變進行分解，從而獲得東西向和垂直向形變特征的方法。該方法考慮到由于SAR衛(wèi)星飛行的幾何原因，導致其對于南北向的形變不敏感，因此MSBAS-2D忽略了滑坡南北向的形變。通常，監(jiān)測獲取的視線向形變只是真實形變的分量，面對復雜的地形條件，無法準確反映斜坡面的真實形變情況。MSBAS-2D技術有效地解決了單一影像只能獲取一維視線向形變的弊端，為滑坡不同方向的形變特征提取提供了可能。

MSBAS-2D的關鍵步驟有：①輸入數(shù)據準備：一組或多組升軌和降軌數(shù)據形成的差分解纏相位圖，將升軌與降軌的差分解纏相位進行地理編碼，選擇高相干的解纏相位通過重采樣至升降軌影像的公共區(qū)域；②創(chuàng)建時間矩陣，并將失相干的像素去除掉；③通過奇異值分解（SVD）對每一個像素進行求解，獲取每一個像素的二維形變速度分量，然后通過對變形速率進行數(shù)值積分重構形變時間序列。其利用多軌道數(shù)據集獲取二維形變時間序列的反演矩陣如下：

式中：θ 、φ——分別為方位角和入射角；

A——獲取的SAR影像的時間間隔；

λ——正則化參數(shù)；

L——零階、一階、二階差分算子；

VE、VU——東西向、垂直向的地表形變速量；

其中，系數(shù)矩陣用C來表示，則C=從而可以求得東西向的形變速度分量（VE）和垂直向的形變速度分量（VU） ：

文章首先分別對覆蓋色拉滑坡的升、降軌Sentinel-1影像進行了處理，獲得各自視線向形變相位；隨后選取干涉質量較好的升軌和降軌結果進行地理編碼，并重采樣到一個公共網格；最后，基于MSBAS-2D方法計算獲取了公共區(qū)域二維時間序列變形結果。多維小基線技術（MSBAS-2D）流程如圖2所示。

圖2 MSBAS-2D處理流程（虛線框為可選步驟）Fig.2 MSBAS-2D processing flow (dotted lines are optional steps)

3 形變結果分析

3.1 單軌形變速率監(jiān)測結果與分析

本文分別對升、降軌的哨兵數(shù)據進行Stacking InSAR處理，獲取研究區(qū)域的視線向年平均形變速率（圖3）。圖3中正值表示形變朝向傳感器方向，負值表示形變背離傳感器方向。結合光學遙感圖像可以看出，在色拉滑坡所處的區(qū)域，植被覆蓋密度不大，所以干涉圖的相干性較高，干涉效果好。

圖3 色拉滑坡2018年1月—2020年4月視線向形變速率圖Fig.3 Line of sight deformation rate map of Sela landslide from January 2018 to April 2020

結合升、降軌視線向年平均形變速率圖可以看出，由于升、降軌衛(wèi)星的飛行姿態(tài)不同導致InSAR所獲得的形變特征不一致。從升軌形變速率圖3（a）可以看出，色拉滑坡在近兩年的時間里，滑坡前緣形變速率較快，中后部有兩處活動性滑坡，最大年形變速率可達?100 mm/a以上。降軌形變速率圖顯示，色拉滑坡近兩年的主要形變區(qū)也是集中在滑坡前緣，最大年形變速率同樣達100 mm/a以上，但其形變區(qū)范圍明顯小于升軌衛(wèi)星監(jiān)測的結果。由圖3可見，基于單軌SAR影像的InSAR監(jiān)測結果雖然可以實現(xiàn)對活動性滑坡的有效識別，但由于受到視線向對滑坡形變觀測視角的影響，所觀測結果并不能對滑坡體的真實變形進行很好地展示。

3.2 二維形變監(jiān)測結果與分析

為了進一步掌握該滑坡的形變特征，利用MSBAS-2D技術對Sentinel-l升、降軌數(shù)據進行處理，計算并獲取了色拉滑坡的二維形變速率及時間序列形變結果。圖4是利用MSBAS-2D獲取的二維形變速率圖，其中，圖4(a)是東西向形變速率，正值表示向東運動，負值表示向西運動。圖4(b)是垂直向形變速率，正值表示抬升，負值表示沉降。由圖4可以看出，滑坡主要形變集中在滑坡前緣和中后部分。觀測P2和P3處滑坡體，發(fā)現(xiàn)兩處滑坡均出現(xiàn)不同程度的相干點缺失現(xiàn)象，尤其是在P2點的下方出現(xiàn)大范圍失相干點缺失，分析其原因可能是該區(qū)域滑坡體由于形變量級過大導致的失相干。

圖4 色拉滑坡2018年1月—2020年4月二維形變速率Fig.4 Two dimensional deformation rate of Sela landslide from January 2018 to April 2020

由于缺少實地的監(jiān)測數(shù)據作為外部參考數(shù)據，為了評價本文監(jiān)測結果的準確性與可靠信，選取穩(wěn)定區(qū)域作為研究區(qū)域的形變參考區(qū)，參考區(qū)位置如圖4中所示。計算了東西向與垂直向的形變參考區(qū)的平均值、標準差。結果如表2所示，其中參考區(qū)的東西向與垂直平均形變速率約小于1.5 mm/a，遠小于滑坡判識的閾值，且標準差在2 mm/a左右，說明在該區(qū)域沒有出現(xiàn)較大的跳變，表明誤差干擾較小，證明了該結果的可靠性。

表2 形變參考區(qū)內的形變速率標準差Table 2 Standard deviation of deformation rate in deformation reference area /(mm·a?1)

為進一步了解該滑坡的變形趨勢，結合形變速率圖及光學影像形變特征，在滑坡體上選取了位于滑坡前緣、滑坡體中后部和滑坡后緣的四個形變較為明顯的點作為滑坡特征點，各點分布位置如圖4所示。分別提取各特征點的累計形變時間序列(圖5)。根據特征點的二維累計形變時間序列結果來看，P1、P3點在2018年11月滑坡形變速率出現(xiàn)突變，發(fā)生了明顯加速。P2點在2018年6月—2018年11月，形變速率較快，在2018年11月后，形變呈現(xiàn)緩和趨勢，P4點一直較為穩(wěn)定，形變量級較小。在2018年1月到2020年4月期間，滑坡中后部分的P1在東西方向上累計形變量可達104 mm，垂直方向累計形變量可達?64 mm；滑坡中后部分的P2在東西方向上累計形變量可達164 mm，垂直方向累計形變量可達?102 mm；滑坡前緣的P3在東西方向累計形變量可達165 mm，垂直方向累計形變量可達?79 mm；滑坡后緣的P4在東西方向累計形變量可達20 mm，垂直方向累計形變量可達?8 mm。

圖5 色拉滑坡特征點二維累積形變時間序列Fig.5 Two dimensional cumulative deformation time series of characteristic points of Sela landslide

根據時間序列監(jiān)測結果，滑坡主要發(fā)生在滑坡前緣和中后部，后緣變形較小。結合滑坡各部位形變及地勢特征分析，滑坡前緣受到河水沖刷的作用，發(fā)生強烈變形，使得滑坡中后部受到滑坡體前緣塊體牽引的作用而發(fā)生形變。因此，初步判斷該滑坡類型屬于牽引式滑坡。

為了進一步了解滑坡的形變原因，以變形量較大且靠近金沙江的P3點為研究對象，將該點的東西及垂向形變與該地區(qū)近兩年的月降雨量進行對比分析(圖6)。該地區(qū)的降雨主要集中在每年的6—9月，在2018年6—9月，滑坡累計形變量較小，東西向和垂向形變呈小幅度波動，該階段滑坡處于緩慢變形階段。在2018年11月份，滑坡形變速率明顯加快，分析其加速原因，認為導致滑坡加速的主要原因是上游白格滑坡的兩次堵江事件。在2018年10—11月，色拉滑坡上游白格發(fā)生兩次堵江事件并形成堰塞湖。第二次堵江后，堰塞體于11月12日開始泄洪，至13日壩體上下游水位貫通，堰塞湖險情解除[25]。對比出現(xiàn)加速的時間點，初步分析色拉滑坡形變速率出現(xiàn)明顯加速主要受白格滑坡二次泄流影響，由于堰塞湖疏通導致河流水體快速下切，使得色拉滑坡前緣受泄流洪水侵蝕，崩滑破壞加速，牽扯滑坡中后部出現(xiàn)局部崩滑現(xiàn)象。

圖6 P3點InSAR時序累計形變與降雨Fig.6 InSAR deformation and rainfall sequence response

同時，提取色拉滑坡在空間和時間上的東西向與垂直向的形變特征，結合斜坡坡度，繪制A-A′和B-B′兩條剖線的累計形變量曲線（圖7）。剖線位置如圖4所示。根據剖線的時間序列結果可以看出，滑坡體在2018年11月—2019年3月期間，累計形變量級發(fā)生了較大的變化，出現(xiàn)了明顯的變形加速趨勢，滑坡前緣的部分坡體因在短時間內形變速率極速增大而出現(xiàn)了失相干現(xiàn)象。截至2020年3月，色拉滑坡復活區(qū)東西向最大累計形變達到165 mm，垂直向?102 mm。結合東西向形變、垂直形變和斜坡坡度關系，可以看出色拉滑坡的滑體坡度較陡，下部滑體坡度為33°～36°，且變形強烈，上部滑體變形跡象不明顯，初步判定目前該滑坡仍處于表層變形階段。

圖7 2018年1月—2020年4月期間色拉滑坡東西向與垂直向變形體Fig.7 East West and vertical deformations of Sela landslide from January 2018 to April 2020

結合2020年9月拍攝的該滑坡現(xiàn)場圖片（圖8），可以看出色拉滑坡已出現(xiàn)明顯的拉張裂縫，局部坡腳形成高陡臨空面，高度約在30～100 m，滑坡區(qū)域表層植被覆蓋稀疏，巖體結構松散，中部拉張裂縫拓寬，主要變形方向為東西走向，與InSAR獲得的二維形變特征基本一致。通過對比圖5的滑坡各點的累積形變量，認為滑坡目前破壞形式以前緣散落、牽引拉張裂縫變寬，中后部出現(xiàn)局部滑塌現(xiàn)象。該滑坡目前仍處于勻速變形階段，一旦遭遇極端天氣，如強降雨和地震等，發(fā)生滑動導致堵江的幾率很大，一旦堵江將可能對周邊水電站和村莊造成較大損失，需要持續(xù)監(jiān)測。

圖8 色拉滑坡現(xiàn)場照片F(xiàn)ig.8 Scene photos of Sela landslide

4 結論

本文利用歐空局 Sentinel-1A/B 數(shù)據對西藏貢覺縣色拉滑坡采用 InSAR 技術進行了形變監(jiān)測研究，獲取了該滑坡東西向與垂直向的二維形變監(jiān)測結果，得出結論如下：

（1）從該滑坡形變速率監(jiān)測結果分析，在 2018年1月—2020年4月，滑坡前緣形變速率最大，尤其是坡腳出現(xiàn)局部形變失相干現(xiàn)象。從形變時間序列結果來看，2018年11月份，滑坡形變速率出現(xiàn)明顯加快，其原因可能與色拉滑坡上游白格滑坡堰塞湖的兩次泄洪事件有關。

（2）由于該滑坡累積形變最大的兩處活動性滑坡區(qū)域位于滑坡前緣，隨著江水對滑坡前緣坡腳的沖刷，已出現(xiàn)明顯的變形。且滑坡的地質條件脆弱，在強降雨的條件下，出現(xiàn)失穩(wěn)的可能性較大。

（3）本文從InSAR的角度對該滑坡進行了監(jiān)測，研究成果較好地展示了色拉滑坡在 2018年1月—2020年4月的形變特征，捕獲了其形變加速的時間點，為該滑坡災害的預警提供重要參考。