0.   引言

滑坡堵江事件是斜坡岩土体因崩塌、滑坡及其转化为碎屑流等而造成江河堵塞和回水的现象^[1]。我国滑坡堵江事件多发生在西南山高坡陡的峡谷地区，这些地区经历强烈隆升运动后形成了地形起伏变化巨大的地形地貌。在地震、暴雨、人工开挖、河流强烈快速下切、库区水位骤然升降等诸多因素的作用下，河谷两侧时常发生大规模堵江滑坡事件。相比于单一滑坡灾害而言，堵江滑坡的危险性通过时空上的延拓而大幅增加^[2-3]。近年来，国内外许多学者对堵江滑坡做了大量的研究。柴贺军等^[4]分析了堵江滑坡在时间和空间上的发育规律，对我国堵江滑坡的区域分布及时间分布规律进行了总结。许强等^[5]运用卫星遥感、无人机航拍、地面合成孔径雷达监测等技术手段对岷江支流四川新磨滑坡进行综合分析，揭示了滑坡的运动过程和成因机制，并对滑坡周边直接受主滑坡动力作用而产生的欠稳定岩土体特征和危险性进行了分析评价。刘传正等^[6]采用多因素赋值统计研判对西藏林芝市米林县雅鲁藏布江左岸色东普沟上游发生的滑坡进行监测，预测出其可能引发崩滑-碎屑流造成雅鲁藏布江再次堵江的临界条件。王群等^[7]运用永久散射体干涉测量技术（Permanent Scatter InSAR，PS-InSAR）和偏移量跟踪技术（Offset Tracking）对金沙江白格滑坡形变进行监测，获取白格滑坡在发生滑动前的运动特征等。朱赛楠等^[8]运用差分合成孔径雷达干涉测量技术和偏移量跟踪技术（Offset Tracking）分析了金沙江上游的色拉滑坡在白格滑坡发生前后的形变特征，并结合地质结构、地层岩性及降雨等分析了滑坡的形成机理。由于堵江滑坡发生的地方隐蔽，多是发生在深切峡谷中，占比达到94.1%^[9]，而传统的监测手段很难做到大范围、全天候的监测。近年发展起来的时序InSAR技术（MT-InSAR）弥补了部分常规滑坡识别与监测手段的不足，具有对微小形变敏感的独特优势，成为滑坡识别与位移测量的重要手段之一^[10-14]。

野外调查显示色拉滑坡近期活动明显，且具有较高的堵江风险^[15]。文章在收集覆盖色拉滑坡区域2018年1月—2020年4月的Sentinel-1A/B卫星升降、轨数据的基础上，采用相位堆叠技术（Stacking InSAR）分析了滑坡体的形变特征，然后使用多维小基线集技术（Multidimensional Small Baseline Subset，MSBAS InSAR）获取了滑坡的二维动态形变序列，成功捕获了滑坡体形变加速的时间点，同时结合降雨数据开展了长时序形变监测及形变原因研究，为金沙江流域堵江滑坡的监测预警研究提供一定的参考。

1.   研究区域及数据

1.1   研究区概况

色拉滑坡位于西藏自治区昌都市贡觉县沙东乡金沙江右岸，总体地势由西北向东南倾斜（图1）。该滑坡距上游叶巴滩水电站23 km，距下游拉哇水电站约63 km。色拉滑坡前缘高程约2649 m，后缘高程约3342 m，相对高差达693 m，滑坡后缘至斜坡后缘陡缓交界处，左侧以冲沟为界，右侧以山脊基岩出露处为界，前缘至基岩出露处。滑坡整体坡向132°，总体呈陡坡地貌^[8,15]。光学影像图1（c）显示色拉滑坡总体形态近似呈舌状，上部较窄，下部宽，前缘已有明显变形，中部有拉张裂缝，呈现灰褐色、灰白色。从光学影像上可以识别出滑坡壁，滑坡拉张裂缝、变形区等信息。

图  1  研究区数据覆盖范围示意图

Figure  1.  Schematic image coverage of the study area

下载: 全尺寸图片 幻灯片

色拉滑坡在构造上位于欧亚大陆西南活动大陆边缘中段、即松潘—甘孜陆缘活动带和昌都陆块，构造线由北西转向近南北向，区内断层和褶皱发育，近东西向的色协龙断裂和近南北向的洛冷登—巴巴断裂在滑坡东北方向交汇^[8]。滑坡区出露的地层主要为第四系滑坡堆积层，第四系残坡积层及二叠系岗托岩组片麻岩。根据地质资料显示，滑坡区主要为第四系滑坡堆积层，岩性主要为碎石块土，碎块石母岩为片麻岩。滑坡顶部和两侧区域为第四系残坡积层，岩性主要为粉质黏土夹碎石。在滑坡前缘和右侧山脊处出露部位的岩性主要为二叠系岗托岩组片麻岩，区域内岩体糜棱岩化和蚀变作用严重。

1.2   SAR数据

文章使用的SAR影像数据来自欧空局（European Space Agency）2014年4月3日发射的对地观测卫星Sentinel-1A/B(哨兵数据)，数据覆盖区域如图1（a）所示。本文所使用的数据时间跨度为2018年1月—2020年4月，分别来自于升轨Track 99和降轨Track 33，数据详细信息见表1。采用了美国地质调查局(USGS)提供的空间分辨率为30 m的SRTM DEM数据，作为外部数据来消除地形相位的影响。同时通过POD精密定轨星历数据(POD Precise Orbit Ephemerides)来纠正干涉组合中的基线误差。

表  1  研究区SAR数据集主要参数

Table  1.  Main parameters of the SAR data sets used in this study

轨道方向	升轨	降轨
轨道号	99	33
入射角	33.8°	39.3°
方位角	−10.4°	−170.0°
成像模式	IW宽幅模式	IW宽幅模式
极化方式	VV	VV
影像数量	66	64
影像时间间隔ΔT	12	12
影像时间范围	2018年1月—2020年4月	2018年1月—2020年4月

下载: 导出CSV 

| 显示表格

2.   研究方法

2.1   相位堆叠InSAR技术

相位堆叠InSAR技术（Stacking InSAR）是对多幅差分干涉图进行加权平均，获取平均形变速率的一种方法^[16-18]。其基本假设为：在每一幅干涉图中，地表形变近似看作线性变化，残余大气误差是随机变化。在进行相位堆叠之前，通过对差分干涉图进行挑选，以相干性和解缠结果为标准，挑选相干性且解缠结果较好的干涉对。在进行相位堆叠（Stacking InSAR）后，叠加的形变相位为对应时间段的累积形变相位，而叠加的大气相位并没有通过干涉图数量的增加而不断累加，而是与干涉图的数量成平方根倍增长。由此使得形变相位在叠加图上占的比例较大，而大气误差所占贡献较小，达到抑制大气噪声和DEM误差，提高形变信息对大气相位的相对精度的作用^[19-20]。

Stacking InSAR技术可以在数据较少情况下，对多幅干涉图进行叠加，获取累积形变量。其中，每幅干涉图的相位变化速率（$V_i$）的标准差、干涉相位$\left({ph}_i\right)$和单个干涉组合的时间基线（$\mathrm{\Delta }{T}_i$）之间的关系为：

把时间基线作为权重因子，求取相位形变速率（$ph\mathrm{\_}rate$）：

式中：$w_i$——第i幅干涉图的权，$w_i={\mathrm{\Delta }T}_i^{-1}$；

$p{h}_i$——单个干涉图的解缠相位值。

从而可以获取年平均形变速率，其标准差为^[21-23]:

相位的标准差为：

2.2   多维短基线集技术

多维小基线集（Multidimensional Small Baseline Subset, MSBAS）技术是由SAMSONOV等^[24]提出的多维形变时序分析方法。多维小基线集-二维（MSBAS-2D）技术是基于不同的入射角与方位角，选择覆盖公共区域的影像，依据其几何关系，对视线向形变进行分解，从而获得东西向和垂直向形变特征的方法。该方法考虑到由于SAR卫星飞行的几何原因，导致其对于南北向的形变不敏感，因此MSBAS-2D忽略了滑坡南北向的形变。通常，监测获取的视线向形变只是真实形变的分量，面对复杂的地形条件，无法准确反映斜坡面的真实形变情况。MSBAS-2D技术有效地解决了单一影像只能获取一维视线向形变的弊端，为滑坡不同方向的形变特征提取提供了可能。

MSBAS-2D的关键步骤有：①输入数据准备：一组或多组升轨和降轨数据形成的差分解缠相位图，将升轨与降轨的差分解缠相位进行地理编码，选择高相干的解缠相位通过重采样至升降轨影像的公共区域；②创建时间矩阵，并将失相干的像素去除掉；③通过奇异值分解（SVD）对每一个像素进行求解，获取每一个像素的二维形变速度分量，然后通过对变形速率进行数值积分重构形变时间序列。其利用多轨道数据集获取二维形变时间序列的反演矩阵如下：

式中：$\theta \text{、}\phi$——分别为方位角和入射角；

A——获取的SAR影像的时间间隔；

$\lambda$——正则化参数；

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

$V_{\mathrm{E}}\text{、}{V}_{\mathrm{U}}$——东西向、垂直向的地表形变速量；

$\hat{\varphi }$——获取的干涉位移。

其中，系数矩阵用$\text{boldsymbol}C$来表示，则$\text{boldsymbol}C$ = $\left(\begin{array}{c}\begin{array}{cc}-\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{s}\mathrm{i}\mathrm{n}\phi A& \mathrm{c}\mathrm{o}\mathrm{s}\phi A\end{array}\\ \lambda L\end{array}\right)$从而可以求得东西向的形变速度分量（$V_{\mathrm{E}}$）和垂直向的形变速度分量（$V_{\mathrm{U}}$）：

文章首先分别对覆盖色拉滑坡的升、降轨Sentinel-1影像进行了处理，获得各自视线向形变相位；随后选取干涉质量较好的升轨和降轨结果进行地理编码，并重采样到一个公共网格；最后，基于MSBAS-2D方法计算获取了公共区域二维时间序列变形结果。多维小基线技术（MSBAS-2D）流程如图2所示。

图  2  MSBAS-2D处理流程（虚线框为可选步骤）

Figure  2.  MSBAS-2D processing flow (dotted lines are optional steps)

下载: 全尺寸图片 幻灯片

3.   形变结果分析

3.1   单轨形变速率监测结果与分析

本文分别对升、降轨的哨兵数据进行Stacking InSAR处理，获取研究区域的视线向年平均形变速率（图3）。图3中正值表示形变朝向传感器方向，负值表示形变背离传感器方向。结合光学遥感图像可以看出，在色拉滑坡所处的区域，植被覆盖密度不大，所以干涉图的相干性较高，干涉效果好。

图  3  色拉滑坡2018年1月—2020年4月视线向形变速率图

Figure  3.  Line of sight deformation rate map of Sela landslide from January 2018 to April 2020

下载: 全尺寸图片 幻灯片

结合升、降轨视线向年平均形变速率图可以看出，由于升、降轨卫星的飞行姿态不同导致InSAR所获得的形变特征不一致。从升轨形变速率图3（a）可以看出，色拉滑坡在近两年的时间里，滑坡前缘形变速率较快，中后部有两处活动性滑坡，最大年形变速率可达−100 mm/a以上。降轨形变速率图显示，色拉滑坡近两年的主要形变区也是集中在滑坡前缘，最大年形变速率同样达100 mm/a以上，但其形变区范围明显小于升轨卫星监测的结果。由图3可见，基于单轨SAR影像的InSAR监测结果虽然可以实现对活动性滑坡的有效识别，但由于受到视线向对滑坡形变观测视角的影响，所观测结果并不能对滑坡体的真实变形进行很好地展示。

3.2   二维形变监测结果与分析

为了进一步掌握该滑坡的形变特征，利用MSBAS-2D技术对Sentinel-l升、降轨数据进行处理，计算并获取了色拉滑坡的二维形变速率及时间序列形变结果。图4是利用MSBAS-2D获取的二维形变速率图，其中，图4(a)是东西向形变速率，正值表示向东运动，负值表示向西运动。图4(b)是垂直向形变速率，正值表示抬升，负值表示沉降。由图4可以看出，滑坡主要形变集中在滑坡前缘和中后部分。观测P2和P3处滑坡体，发现两处滑坡均出现不同程度的相干点缺失现象，尤其是在P2点的下方出现大范围失相干点缺失，分析其原因可能是该区域滑坡体由于形变量级过大导致的失相干。

图  4  色拉滑坡2018年1月—2020年4月二维形变速率

Figure  4.  Two dimensional deformation rate of Sela landslide from January 2018 to April 2020

下载: 全尺寸图片 幻灯片

由于缺少实地的监测数据作为外部参考数据，为了评价本文监测结果的准确性与可靠信，选取稳定区域作为研究区域的形变参考区，参考区位置如图4中所示。计算了东西向与垂直向的形变参考区的平均值、标准差。结果如表2所示，其中参考区的东西向与垂直平均形变速率约小于1.5 mm/a，远小于滑坡判识的阈值，且标准差在2 mm/a左右，说明在该区域没有出现较大的跳变，表明误差干扰较小，证明了该结果的可靠性。

表  2  形变参考区内的形变速率标准差

Table  2.  Standard deviation of deformation rate in deformation reference area /(mm·a⁻¹)

SAR数据	平均值	标准差
东西向	−0.489	2.115
垂直向	1.499	1.574

下载: 导出CSV 

| 显示表格

为进一步了解该滑坡的变形趋势，结合形变速率图及光学影像形变特征，在滑坡体上选取了位于滑坡前缘、滑坡体中后部和滑坡后缘的四个形变较为明显的点作为滑坡特征点，各点分布位置如图4所示。分别提取各特征点的累计形变时间序列(图5)。根据特征点的二维累计形变时间序列结果来看，P1、P3点在2018年11月滑坡形变速率出现突变，发生了明显加速。P2点在2018年6月—2018年11月，形变速率较快，在2018年11月后，形变呈现缓和趋势，P4点一直较为稳定，形变量级较小。在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  色拉滑坡特征点二维累积形变时间序列

Figure  5.  Two dimensional cumulative deformation time series of characteristic points of Sela landslide

下载: 全尺寸图片 幻灯片

根据时间序列监测结果，滑坡主要发生在滑坡前缘和中后部，后缘变形较小。结合滑坡各部位形变及地势特征分析，滑坡前缘受到河水冲刷的作用，发生强烈变形，使得滑坡中后部受到滑坡体前缘块体牵引的作用而发生形变。因此，初步判断该滑坡类型属于牵引式滑坡。

为了进一步了解滑坡的形变原因，以变形量较大且靠近金沙江的P3点为研究对象，将该点的东西及垂向形变与该地区近两年的月降雨量进行对比分析(图6)。该地区的降雨主要集中在每年的6—9月，在2018年6—9月，滑坡累计形变量较小，东西向和垂向形变呈小幅度波动，该阶段滑坡处于缓慢变形阶段。在2018年11月份，滑坡形变速率明显加快，分析其加速原因，认为导致滑坡加速的主要原因是上游白格滑坡的两次堵江事件。在2018年10—11月，色拉滑坡上游白格发生两次堵江事件并形成堰塞湖。第二次堵江后，堰塞体于11月12日开始泄洪，至13日坝体上下游水位贯通，堰塞湖险情解除^[25]。对比出现加速的时间点，初步分析色拉滑坡形变速率出现明显加速主要受白格滑坡二次泄流影响，由于堰塞湖疏通导致河流水体快速下切，使得色拉滑坡前缘受泄流洪水侵蚀，崩滑破坏加速，牵扯滑坡中后部出现局部崩滑现象。

图  6  P3点InSAR时序累计形变与降雨

Figure  6.  InSAR deformation and rainfall sequence response

下载: 全尺寸图片 幻灯片

同时，提取色拉滑坡在空间和时间上的东西向与垂直向的形变特征，结合斜坡坡度，绘制A-A′和B-B′两条剖线的累计形变量曲线（图7）。剖线位置如图4所示。根据剖线的时间序列结果可以看出，滑坡体在2018年11月—2019年3月期间，累计形变量级发生了较大的变化，出现了明显的变形加速趋势，滑坡前缘的部分坡体因在短时间内形变速率极速增大而出现了失相干现象。截至2020年3月，色拉滑坡复活区东西向最大累计形变达到165 mm，垂直向−102 mm。结合东西向形变、垂直形变和斜坡坡度关系，可以看出色拉滑坡的滑体坡度较陡，下部滑体坡度为33°~36°，且变形强烈，上部滑体变形迹象不明显，初步判定目前该滑坡仍处于表层变形阶段。

图  7  2018年1月—2020年4月期间色拉滑坡东西向与垂直向变形体

Figure  7.  East West and vertical deformations of Sela landslide from January 2018 to April 2020

下载: 全尺寸图片 幻灯片

结合2020年9月拍摄的该滑坡现场图片（图8），可以看出色拉滑坡已出现明显的拉张裂缝，局部坡脚形成高陡临空面，高度约在30～100 m，滑坡区域表层植被覆盖稀疏，岩体结构松散，中部拉张裂缝拓宽，主要变形方向为东西走向，与InSAR获得的二维形变特征基本一致。通过对比图5的滑坡各点的累积形变量，认为滑坡目前破坏形式以前缘散落、牵引拉张裂缝变宽，中后部出现局部滑塌现象。该滑坡目前仍处于匀速变形阶段，一旦遭遇极端天气，如强降雨和地震等，发生滑动导致堵江的几率很大，一旦堵江将可能对周边水电站和村庄造成较大损失，需要持续监测。

图  8  色拉滑坡现场照片

Figure  8.  Scene photos of Sela landslide

下载: 全尺寸图片 幻灯片

4.   结论

本文利用欧空局 Sentinel-1A/B 数据对西藏贡觉县色拉滑坡采用 InSAR 技术进行了形变监测研究，获取了该滑坡东西向与垂直向的二维形变监测结果，得出结论如下：

（1）从该滑坡形变速率监测结果分析，在 2018 年 1 月—2020 年 4 月，滑坡前缘形变速率最大，尤其是坡脚出现局部形变失相干现象。从形变时间序列结果来看，2018 年 11 月份，滑坡形变速率出现明显加快，其原因可能与色拉滑坡上游白格滑坡堰塞湖的两次泄洪事件有关。

（2）由于该滑坡累积形变最大的两处活动性滑坡区域位于滑坡前缘，随着江水对滑坡前缘坡脚的冲刷，已出现明显的变形。且滑坡的地质条件脆弱，在强降雨的条件下，出现失稳的可能性较大。

（3）本文从InSAR的角度对该滑坡进行了监测，研究成果较好地展示了色拉滑坡在 2018 年1 月—2020 年 4 月的形变特征，捕获了其形变加速的时间点，为该滑坡灾害的预警提供重要参考。