基于SBAS-InSAR技术的金沙江流域沃达村巨型老滑坡形变分析*

冯文凯 顿佳伟 易小宇 张国强

(地质灾害防治与地质环境保护国家重点实验室(成都理工大学)，成都 610059，中国)

0 引 言

随着人类活动深度和广度的日益加深和拓展，我国西部地区大型突发地质灾害日益频发。这些灾害表现出隐蔽性强、危害性高，具有显著的突发性，严重影响了社会经济发展和人民安居乐业，给地质灾害防治与地质环境保护带来了巨大挑战。

2017年6月24日清晨6时许，四川省茂县叠溪镇新磨村新村组后山约450×104￣￣m3山体突发顺层高位滑动，造成80余人失踪或遇难，引起国内外广泛关注(许强等， 2017)。灾害发生后，意大利Tre Altamira (2017)通过InSAR(Interferometric Synthetic Aperture Radar，合成孔径雷达干涉)技术对2014年以来45景Sentinel-1雷达卫星影像进行了分析，证实新磨村滑坡的显著形变主要发生在滑前的数月之内。2018年10月11日凌晨4时许，金沙江流域西藏昌都市江达县波罗乡白格村发生特大库岸岩质滑坡，总方量约3165×104￣￣m3，滑坡体迅速阻断金沙江，回水造成多条道路中断，川藏滇等地紧急转移安置2.1万余人(冯文凯等， 2019)。许强等(2018)通过分析白格滑坡灾前雷达卫星数据，指出该滑坡灾害发生前一年，滑源区斜坡最大位移量约25 m。这些案例充分说明InSAR技术对处于高位、交通困难、人迹罕至的隐蔽性地质灾害的识别和监测研究具有显著优势。

InSAR技术是近20 a来兴起的一种新型对地观测技术，具有全天候工作、覆盖范围广、探测精度高、非接触和综合成本低等优点，已被广泛地应用于滑坡(李凌婧等，2017)、地面沉降(张静等，2018)、地裂缝(赵超英等，2011)等地质灾害的形变监测。InSAR技术中的小基线(Small Baseline Subsets， SBAS)技术具有获取微小形变信息和长时间序列缓慢地表变形的优势(莫玉娟等， 2018)，在复杂山区地质灾害预警和监测领域有较为广阔的应用前景。在基于SBAS-InSAR技术对滑坡早期识别与形变监测的研究领域中，很多学者取得了相应的实质性成果。余睿(2014)采用SBAS-InSAR技术对甘肃省西和县慢速型滑坡进行全局和局部的监测，揭示了影响西和县滑坡体移动的主导因素； 刘筱怡等(2017)以鲜河水断裂带内蠕滑型滑坡为研究对象，采用SBAS-InSAR技术获取蠕滑型滑坡时间序列形变特征，分析其发展变化规律，对鲜水河断裂带以及类似构造活动地区的地质灾害防治起到了重要作用； 聂兵其等(2018)基于SBAS-InSAR技术对丹巴县城进行滑坡形变探测及隐患识别，验证了InSAR监测结果与现场调查的一致性。可以说，SBAS-InSAR技术为地质灾害隐患点的普查和分析研究提供了重要支撑。在充分利用星载雷达对地观测技术的前提下，开展针对性的工程地质调查和复核工作，一方面有利于正确认识InSAR解译数据的规律性，另一方面也有利于灾害预警研判的准确性，从而大大提升我国科学防灾减灾的技术实力。

沃达村巨型老滑坡是区域地质灾害调查隐患数据库的隐患点，位于江达县岩比乡沃达村、金沙江右岸，距离白格滑坡直线距离42 km。本文在搜集沃达村老滑坡区域升轨哨兵一号(Sentinel-1A)2017年3月30日至2019年9月28日的39景数据的基础上，采用SBAS-InSAR技术分析了老滑坡堆积体的形变特征，并进一步针对性地开展了现场工程地质调查和复核工作，在此基础上分析了SBAS-InSRAR技术在复杂山区地质灾害预警和监测领域的适用性，为类似老滑坡监测预警提供了新的思路与借鉴。

1 沃达村老滑坡概况

沃达村老滑坡位于西藏自治区昌都市江达县岩比乡沃达村，坐标东经E93°26′06″、北纬N31°26′06″(图 1)。前期区域地质灾害调查表明，该地段地质历史时期曾发生滑坡灾害，其滑坡堆积体在降水、地震和水位波动情况下易发生复活，应对其重点监测研究。

图 1 滑坡地理位置图Fig. 1 Geographic location map of the landslide

图 2 综合光学遥感解译图Fig. 2 Integrated optical remote sensing interpretation map

从综合光学遥感解译图(图 2)分析，老滑坡体纵长约1.9 km，横宽约1.7 km，主滑方向为NE30°； 老滑坡堆积体局部复活，复活部分纵向长约0.88 km，横宽约1.08 km，面积约0.925 km2。滑坡体平面形态为半圆状，前缘呈由缓变陡的陡崖形态，缓部与江面高差最大约240 m，前端临空面陡倾，坡脚即为金沙江河谷，高程约2950 m，中后部相对平缓。据前期调查显示，该老滑坡上部覆盖层为第四系滑坡堆积物，含碎石、砾石黏性土，厚约20～30 m，下伏基岩为变质页岩，岩层产状约为240°∠80°。

图 3 数据覆盖范围图Fig. 3 Data coverage chart

2 SAR数据获取及处理方法

2.1 数据选取及覆盖概况

本文使用的单视复数SAR影像数据来自欧空局(European Space Agency)2014年4月3日发射的对地观测卫星Sentinel 1A，数据覆盖区域如图 3所示。波段为C波段，所用影像为合成孔径雷达99轨道数据，时间跨度从2017年3月30日至2019年9月28日， 24 d一景，共39景，轨道方向为升轨，成像中心入射角为36.2°，成像模式为IW宽幅模式，地面分辨率为5 m×20 m，幅宽为250 km×250 km，极化方式为VV，方位向和距离向采样间距分别为2.329 m、13.956 m(表 1)。本研究采用了日本JAXA提供的30 m分辨率ALOS World 3D(AW3D30)DEM作为外部数据来消除地形相位的影响。

表 1 沃达村滑坡Sentinel 1A数据列表Table 1 Data list of sentinel 1 in Dawo Village landslide

2.2 技术原理与形变速率获取

2.2.1 SBAS-InSAR技术原理

为减少时间和空间去相干，Berardino et al.(2002)在Casu(2000)研究的基础上，提出了短基线差分干涉法。该方法的基本原理是利用具有较短时-空基线的影像对产生干涉图提高相干性，再进行解缠干涉相位。依据相干像元相位和观测时间的关系，采用奇异值分解方法将多个小基线集合数据联合起来求解(许军强等， 2019)，有效地解决了各数据集之间空间基线过长造成的时间不连续问题，提高了监测的时间分辨率，从而求得影像序列间地表形变速率的最小范数最小二乘解(何秀凤等， 2012)。

其基本原理：假设有N+1景单视复数影像，其成像时间为(t0，t1，….，tn)，设某一差分解缠图为j，某像元(x，r)对应的解缠相位可表示为式(1)：

δφj(x，r)=φ(tB，x，r)-φ(tA，x，r)

(1)

(2)

由于在处理过程中产生M幅干涉图，根据式(2)可得M个方程，以矩阵形式表示为式(3)：

δφ(x，r)=Aφ(x，r)

(3)

其中，A为M×N系数矩阵；φ(x，r)为(x，r)点在N个时刻对应的未知形变相位构成的矩阵，当M≥N时，通过最小二乘法得出式(4)：

φ(x，r)=(ATA)-1ATδφ(x，r)

(4)

当M

图 4 影像时间基线Fig. 4 Image time baseline

图 5 影像空间基线Fig. 5 Image spatial baseline

2.2.2 视向形变速率获取

本文使用ENVI中SArscape模块，用SBAS算法对Sentinel 1A数据进行干涉处理，共生成109个干涉对，时间基线为75 d，临界基线百分比为2%，各像对时间基线和空间基线的连接方式如图 4和图 5所示，图中绿色的点代表 39景影像的成像时间，黄色点代表作为参考主影像的成像时间。对合成的干涉图做第1次反演、第2次反演、地理编码以及栅矢转换，生成地表沿雷达视向的平均形变速率(mm·a-1)，最终得到时间序列上的形变结果。具体处理流程如图 6所示。

图 6 SBAS处理流程Fig. 6 SBAS-InSAR processing flow

图 7 雷达视线方向和坡度方向几何示意图 (Cascini et al., 2010)Fig. 7 Geometric sketches of radar line- of -sight direction and slope direction(Cascini et al.， 2010)

2.2.3 二维形变转换

由于滑坡多沿斜坡面进行滑动，雷达视线方向的形变信息无法准确反映斜坡面的真实形变情况(戴可人等， 2019)，考虑雷达视线方向、坡度方向和垂直沉降方向的几何关系(图 7)，假设运动沿单位矢量指定的方向发生，采用式(5)、式(6)(Cascini et al.， 2010)和式(8)(Colesanti et al.， 2006)将视线方向的形变速率转化为坡度方向和垂直方向的形变速率。

(5)

VSlope=VLos/cosβ

(6)

式(6)中，cosβ可用式(7)(Cascini et al.， 2010)计算；

(7)

(8)

3 形变结果分析

结合遥感影像分析，该滑坡中部和后缘部位表面植被覆盖密度不大，相干点较多，干涉效果较好，有较明显的速率积累，滑坡部分地方植被覆盖茂密，造成影像的失相干，速率点较为稀少。

形变速率为负值表示形变点向远离卫星传感器方向运动，正值表示形变点向靠近卫星传感器方向运动。由InSAR视向形变分析结果可见，图 8中形变值均为负值，表明滑坡整体向着金沙江方向向下滑动，符合滑坡体的运动规律。从图 8中可以看出该滑坡形变较强部位位于滑坡中部和中前缘(图 8紫色虚线圈内)，其中紫色虚线圈内的红色区域形变量最大，形变速率≥80 mm·a-1。

图 8 2017年3月～2019年9月形变速率图(LOS方向)Fig. 8 Deformation rate map from March 2017 to September 2019(LOS Direction)

沃达村滑坡沿坡向年平均速率(VSlope)和垂直方向年平均速率(Vu)分别如图 9、图 10所示，解译出的强变形区与LOS方向大致一致，但在形变速率量值上存在较大差异，沿坡向形变速率最大。

图 9 2017年3月～2019年9月形变速率图(Slope方向)Fig. 9 Deformation rate map from March 2017 to September 2019(Slope Direction)

图 10 2017年3月～2019年9月形变速率图(Vertical方向)Fig. 10 Deformation rate map from March 2017 to September 2019(Vertical Direction)

4 实地工程地质调查与复核

为复核解译数据，开展了实地工程地质调查与复核工作。经实地踏勘，滑坡体前缘地形陡峭； 中前部坡度较大，树木生长茂盛； 中部为缓坡平台，主要为耕地，树木零星生长； 后部为古滑坡滑源区(图 11a)。受冲沟下切侵蚀及前缘陡倾临空面综合作用，滑坡前缘发育多处向冲沟、临空面方向的崩滑变形体(B01～B04)，崩滑变形体表部裂缝横张、错坎发育，裂缝张开2～10 cm，最大错距1～2 m，前端陡崖可见明显基覆界面(图 11b)。滑坡复活区中部发育多处小型滑动变形体(H01～H03)，坡体发育横切梯级陡坎，高度几十厘米到几米不等，延伸几十米到数百米不等，陡坎上部多耕作，有灌木丛生，表明陡坎形成时间较长(图 11c)； 可见明显马刀树，结合树龄推测稳定时间至少10 a以上(图 11d)。

图 11 沃达村老滑坡照片Fig. 11 Photographs of the old landslide in Woda Village

4.1 滑坡整体变形与形变速率对比分析

图 12 InSAR时序形变与降雨响应Fig. 12 InSAR time series deformation and rainfall response

复活区主要位于老滑坡体中下部，强烈变形区(Ⅰ现场)后缘位于老滑坡后缘陡缓交界处，两侧边界依地形延伸至老滑坡边界(图 8)，其内裂缝、下错台坎、局部滑塌现象明显。现场的变形迹象与SBAS-InSAR方法获取的形变速率图(图 8)特征有着较高的一致性。在滑坡视向形变速率图(图 8)形变量较大的红色区域取一点A作出该点历史累积位移形变曲线，与搜集的2017年1月～2019年8月月累计降雨曲线对比分析(图 12)：随着2017-05～2017-06、2018-05～2018-08、2019-06～2019-07 3个不同时间段内降雨量急剧增加，A点形变量明显加速； 降雨量较小的时间段内，A点形变量增加变缓，形变趋势与降雨具有较好的响应规律。据现场调查，滑坡所处金沙江水位波动最高高程为2930～2950 m，而滑坡局部复活区位置位于3220～3580 m，远离金沙江水位波动影响区。说明该滑坡中前部复活区主要受降雨影响，为老滑坡降雨复活型滑坡。

图 13 2017年3～12月形变速率及剖面形变图Fig. 13 Deformation rate and profile deformation map from March to December 2017

图 14 2018年形变速率及剖面形变图Fig. 14 Deformation rate and profile deformation map of 2018

图 15 2019年1～9月形变速率及剖面形变图Fig. 15 Deformation rate and profile deformation map from January to September 2019

4.2 滑坡局部变形与形变速率对比分析

现场复核阶段对前缘临空面发育的4处崩滑变形体(B01～B04)和滑坡中、下部3处滑动变形体(H01～H03)进行了细致调查。4处崩滑变形体均为表层发生的溜滑，规模较小，树木歪斜的范围位于冲沟沟口，表明为集中降雨引起的变形。3处滑动变形体中， 2处(H02和H03)因为前部临空面变形导致出现牵引式滑动变形， 1处(H01)因降雨和人类工程活动因素发生牵引式滑动变形。下面主要分析滑动变形体H01的情况。

H01滑动变形体位于坡体的中后部，沃达村公路上侧，滑动方向与坡向一致(图 16a)。据访，H01滑动变形体前缘为沃达村村道，村道挡土墙于2013年左右发生隆起破坏(图 16b)； 2018年中旬，H01滑动变形体发生了较为强烈的变形，导致挡墙上方较高处田地出现多处裂缝及下错台坎(图 16c)。调查发现，其下部紧邻的村民房屋完好无破坏，表明其仅仅是近年来表层岩土体在降雨等作用下发生的局部牵引式滑动变形。

为细致分析H01滑动变形体的滑动性质与雷达形变速率信号之间的相关性，沿H01滑动变形体主滑方向布置剖面线，并依次取1#～4#监测点，利用图 9绘制其3年内滑坡历史累积位移形变曲线(图 16d)。如图所示， 4#监测点历史累积形变量较大，增加的速率也较大； 1#～3#监测点历史累积形变量增加幅度基本一致。前缘的形变速率远大于中部和后部的形变速率，说明H01滑动变形体属于牵引式滑动变形。

根据对滑塌体后缘下错台坎发育处沿坡向以及垂向形变速率点的分析：沿坡向的形变速率(图 9、图 16a)颜色各异，形变速率差异性很大，形变速率140 mm·a-1、160 mm·a-1、180 mm·a-1均有分布，与裂缝形成特征吻合； 沿垂向形变速率分析(图 10)，在后缘位置黄色形变速率(-80～-60 mm·a-1)附近出现多级的橙色形变速率(-100～-80 mm·a-1)，说明该区域发生多级下错沉降现象。两种现象印证了滑坡强变形区历史出现多级下错台坎，表现出蠕动变形，与地质调查的现象基本相符。

5 结 论

(1)通过对视向、坡向、垂向方向的形变速率分区以及形变点之间的差异性分析，初步将沃达村滑坡分为强烈形变区(Ⅰ雷达)和均匀形变区(Ⅱ雷达)，与现场判定的强烈变形区(Ⅰ现场)和弱变形区(Ⅱ现场)基本吻合。

图 16 H01变形体SBAS-InSAR监测与野外现象对比Fig. 16 Comparison of SBAS-InSAR monitoring data of H01 landslide with field review

(3)基于坡向形变速率图、垂向形变速率图中形变速率值的分布以及量值的差异性特征，初步判定了局部滑塌的性质以及多级台坎发育的位置，与野外地质调查吻合性较好。

(4)通过SBAS-InSAR技术分析以及现场调查复核，发现老滑坡复活变形迹象与SBAS-InSAR技术解译成果有着较好的一致性，在复杂山区地质灾害预警和监测领域有较为广阔的应用前景，为类似老滑坡的监测预警提供新的思路与借鉴。

致 谢本文在完成过程中得到了中国地震局地壳应力研究所、四川省国土空间生态修复与地质灾害防治研究院、江达县气象局等单位的帮助。欧洲太空局为本研究提供的Sentinel-1A数据，日本JAXA提供的30 m分辨率“ALOS World 3D—30 m” DEM。在此一并表示感谢。