基于SBAS-InSAR技术的木场古滑坡变形特征分析

张钟远，徐世光，邓明国，曾营，4，李超，5，晁江琴

(1.昆明理工大学国土资源工程学院，昆明 650000；2.哈尔滨工业大学重庆研究院，重庆 400020；3.云南地矿工程勘察集团有限公司，昆明 650000；4.西南交通大学土木工程学院，成都 610000；5.青岛市勘察测绘研究院，青岛 266000；6.昭通学院地理科学与旅游学院，昭通 657000)

滑坡灾害严重威胁着人们的生命财产安全，因此非常有必要对滑坡灾害开展监测工作。合成孔径雷达(synthetic aperture radar，SAR)技术可随时主动成像进行观测，具有穿透性强、不受云雾等自然现象影响的特点，且有数据精确度高的优点，近年来众多学者将合成孔径雷达干涉测量(interferometric SAR，InSAR)技术应用于滑坡监测中，取得了不错的成效[1-3]。

主要的星载SAR系统有日本的对地观测卫星(ALOS-2)，欧洲航空航天局的哨兵1号(Sentinel-1A)和国产高分3号(GF-3)等[4]。通过星载SAR系统获取研究区的影像数据，采用InSAR技术获取研究区地表变形速率和时间序列结果，识别潜在滑坡隐患，监测古滑坡变形特征，从地表形变角度有效监测滑坡近期变形状态，从而分析预测滑坡灾变趋势，可为地质灾害预警提供数据支持[5-7]。对于已发生的滑坡灾害研究，可通过时序InSAR技术追溯滑坡的变形过程，重现滑坡不同演化阶段的地表位移，实现一些重大滑坡地质灾害的“复盘”，结合数值模拟等其他研究手段进而揭示滑坡演化机理[8-10]。近年来，中国学者提出“天-空-地”一体化地质灾害监测技术，且成效显著，使得InSAR技术得到大力推广[11-13]。常用的InSAR技术有合成孔径雷达差分干涉测量(differential InSAR，D-InSAR)，永久散射体雷达干涉测量技术(permanent scatterer InSAR，PS-InSAR)和短基线子集干涉测量技术(small baseline subset InSAR，SBAS-InSAR)，D-InSAR与PS-InSAR在植被覆盖区容易造成失相干，而SBAS-InSAR技术用于监测地表形变，可以克服时间与空间失相干现象，在植被覆盖区具有一定的优势。

基于此，获取了2018年6月19日—2020年12月17日期间覆盖木场古滑坡的58景哨兵1号(Sentinel-1)降轨数据，采用SBAS-InSAR技术对木场古滑坡进行形变监测，结合木场古滑坡的地质特征分析滑坡时序变形趋势，促进InSAR技术在滑坡中的应用，以期为滑坡防灾减灾工作提供科学依据。

1 木场古滑坡基本概况

1.1 形态特征

古滑坡位于镇康县木场镇木场村，分布于木场镇南片河北侧的山坡，目前古滑坡整体处于稳定状态。但是，近年来受当地居民修筑房屋、乡村公路开挖边坡的影响，在古滑坡体范围内发育了次级滑坡H1、H2和H3(图1)。古滑坡总体呈扇状分布，分布高程1 760～1 845 m，轴向长为150 m，均宽约350 m，坡面面积约为60 746 m2，滑坡物质组成主要为含碎石粉质黏土，平均厚11.46 m，钻孔揭露最大厚度为17.20 m，古滑坡全貌如图1所示。该滑坡体积约7×105m3，为牵引式土质中型滑坡。古滑坡后缘根据最后一级错落台坎1、台坎2确定(图2)，错落坎后基岩中风化灰岩，属稳定岩层，不易产生变形，滑坡中后部发育两级斜坡滑移后的错落坎，陡坎3发育在古滑坡后缘，分布高程约为1 839 m，连续发育长45 m，走向NE45°；陡坎4发育在古滑坡体中前部，分布高程约为1 789 m，连续发育长55 m，走向NE55°。古滑坡两侧滑动变形迹象不明显，边界根据后缘错落台坎的结束位置及地形条件作综合分析和确定，从地形上看，古滑坡右侧：错落坎与地势低洼处相连，将低洼处定为古滑坡的右侧边缘较为合理；古滑坡左侧：大致以左侧的岩土界线为界；滑坡前缘特征不明显，受早期人类工程活动影响，造成坡体临空。探槽见有古滑坡变形滑动裂隙(图3)，其现已被泥土完全充填。该裂隙面张开5～10 cm，倾向坡外，根据裂缝延伸趋势，古滑坡此前变形失稳应在坡体前缘坡脚处剪出，因此推测滑坡前缘位于乡村公路的北部斜坡坡脚。

图1 木场古滑坡全貌

图2 古滑坡后缘错落坎

图3 探槽揭露古滑坡裂缝

1.2 古滑坡堆积体物质结构特征

古滑坡可分为滑坡体、滑动面(带)和滑床，现分述如下。

1.2.1 滑坡体

根据调查、槽探及钻探揭露可知，滑坡土层组成主要由粉质黏土组成，局部含有块石土，粉质黏土偶见碎石：稍湿～湿，松散～稍密。碎石成份主要为强风化泥质灰岩、白云岩，碎石粒经一般0.5～3 cm，最大块石粒经可达4 cm，多呈次棱角状～棱角状，碎石含量约占总量的5%，平均厚11.46 m，钻孔揭露最大厚度达17.20 m，该层在滑坡体上均有分布，如图4、图5所示。根据地质测绘及钻探显示，滑坡体与下覆基岩呈犬齿交错接触。

图4 粉质黏土

图5 块石土

1.2.2 滑动面(带)

古滑坡主要是中层土质滑坡，由于下覆永德组灰岩呈一峰一谷交替分布，古滑坡滑带在峰位置沿基覆界面滑动，在谷位置滑面从土体中剪出。大部分钻孔中均揭露有滑带土(图6)，滑带土为含碎石黏土，受地下水侵染呈灰褐色，湿，呈软塑状，碎石含量约35%，碎石呈角砾状，磨圆差，滑带土与上下土层或基岩物性有明显区别。受地下水侵染呈灰褐色，湿，呈软塑状，碎石含量约35%，碎石呈角砾状，磨圆差。

图6 钻孔揭露的滑带土

1.2.3 滑床

通过工程地质测绘及钻探资料显示永德组灰岩较为完整，工程地质条件较好。根据工程地质测绘及钻探显示，滑坡区下覆基岩呈峰-谷峰出现，呈犬齿交替状横向分布的滑床有利于古滑体的稳定。

1.3 次级滑坡H1、H2、H3概况

1.3.1 H1

H1滑坡，位于木场古滑坡体南西部，2012年当地居民对斜坡开挖拟修筑房屋，形成坡长约155 m。最大高差达17 m的陡坎，坡面角度45°～65°(图7)。据调查，2014年7月以前，滑坡H1所在边坡未见有明显的变形迹象，2014年8月，经访问，该边坡在强降水期间，坡体的粉质黏土在饱水的状态下受其自身重力的影响而发生垮塌，特别是2014年8月夜间强降雨，边坡出现不同程度的垮塌。在强降雨期间，该边坡经常发生垮塌，垮塌长度1～2 m，垮落后，当地居民又进行了及时清理，加剧了该边坡的失稳，后又进行了边坡加固。在2020年8月19日凌晨，在强降雨条件下加固失效，发生大面积滑坡，此次滑坡事件造成15幢民房不同程度损毁，如图8所示。

图7 滑坡H1发生前的边坡(拍摄于2014年)

图8 滑坡H1发生后图片(拍摄于2020年)

1.3.2 H2

H2滑坡呈舌状展布，分布高程1 758～1 972 m，轴向长为55 m，宽约50 m，坡面面积约3 156 m2。滑坡物质组成主要为含碎石粉质黏土，一般厚5～9 m，最大厚度达11 m，滑坡体平均厚度约7.5 m，该滑坡为一浅表层地质牵引式滑坡；体积约2.37×104m3，按其规模属于小型滑坡，滑坡全貌如图9所示。滑坡体总体滑向125°。

图9 滑坡H2全貌

目前滑坡H2在天然状态下处于基本稳定状态，在暴雨或地震等外力作用下易沿滑面产生滑动，在坡体、后缘及两侧发育有大量张裂缝，裂面粗糙，无擦痕，在坡体发育的横张裂缝大致呈圆弧形，弧尖背向滑动方向，延伸较长，而纵张裂缝比较顺直，成梭形，延伸不远，后缘和两翼裂缝连续分布并有明显位移。

滑坡H2裂缝L5位于滑坡体左侧，为边界剪切裂隙，裂缝宽为0.5～1 m，倾角72°，深度可见2 m，延伸长度17 m，走向300°，裂隙面偶见碎石，粒径一般为0.5～1.5 cm，最大为3 cm，主要为粉质黏土；L6分布于滑坡体后缘，变形模式属滑移-拉裂。错落高度0.5～2.0 m，延伸长度24 m，走向40°，错落面偶见碎石，粒径一般为0.5～1.5 cm，最大为3 cm，主要为粉质黏土；L7位于滑坡体右侧，为边界剪切裂隙，裂缝宽为0.05～1.5 m，倾角近乎直立，深度可见1.8 m，延伸长度为18 m，走向340°，裂隙面偶见碎石，粒径一般为0.5～1.5 cm，最大为3 cm，主要为粉质黏土；L8位于滑坡体中上部右侧，为滑移-压致拉裂，裂缝宽为0.05～0.25 m，倾角近乎直立，深度可见0.5 m，延伸长度8 m，走向60°，裂隙面偶见碎石，粒径一般为0.5～1.5 cm，最大为3 cm，主要为粉质黏土；L9位于滑坡体中部，为滑移-压致拉裂，裂缝宽为0.05～0.30 m，倾角近乎直立，深度可见0.5 m，延伸长度为13 m，走向45°。裂隙面偶见碎石，粒径一般为0.5～1.5 cm，最大为3 cm，主要为粉质黏土；L10位于滑坡体中后部，为滑移-压致拉裂，裂缝宽为0.10～0.25 m，倾角近乎直立，深度可见1.0 m，延伸长度10 m，走向50°，裂隙面偶见碎石，粒径一般为0.5～1.5 cm，最大为3 cm，主要为粉质黏土；L11位于滑坡体前部，错落高度为0.30～1.20 m，倾角近乎直立，深度可见1.0 m，延伸长度为10 m，走向70°。裂隙面偶见碎石，粒径一般为0.5～1.5 cm，最大为3 cm，主要为粉质黏土；L12位于滑坡前缘的乡村公路内侧的滑坡体前缘，为鼓胀裂缝，受滑坡推力的影响，公路护脚墙部分被摧毁，公路坡下即为滑坡剪出口，如图10所示。

图10 滑坡H2裂缝

1.3.3 H3

滑坡H3边位于木场古滑坡体北东部，临近木场村村委会。滑坡方向106°，由于当地居民对斜坡开挖拟修筑房屋，形成两级平台，平台与上下地形形成陡坎，第一陡坎高程1 769～1 774 m，长约45 m，第二陡坎1 774～1 982 m，长约35 m，两处陡坎临空面角度55°～75°，近乎直立。前缘出现鼓胀裂缝，后缘发育拉张裂缝(图11)，为滑移-压致拉裂，裂缝宽为0.10～0.25 m，倾角近乎直立，深度可见1.0 m，延伸长度10 m，走向50°。裂隙面偶见碎石，粒径一般为0.5～1.5 cm，最大为3 cm，主要为粉质黏土。据实地调查，该滑坡在强降雨期间，坡体的粉质黏土在饱水的状态下受其自身重力的影响而发生局部垮塌，每次垮落长度为0.5～2 m，垮落后，当地居民又进行了及时清理，为滑坡提供了临空面，及时清理反而促进了滑坡H1的发生。

图11 滑坡H3裂缝

2 基于SBAS-InSAR的地表形变分析技术

SBAS-InSAR是Berardino等[14]于2002年提出的监测地面微小形变的时间序列InSAR技术方法，能有效克服D-InSAR在数据处理过程中的时空失相干、大气效应等问题，相比永久散射点干涉测量技术(PS-InSAR)，在一些人工建筑非常少的区域具有更大密度的监测效果[15]。

通过时间序列t0，t1，…，tn获得n+1幅SAR数据影像，对超级主影像进行配准，将垂直基线小于阈值的SAR数据影像归为一组。每组进行差分干涉处理，形成p组影像，得到m幅差分干涉图，若n为奇数，则m有

(1)

以t0为初始时刻，任意时刻ti(i=1,2,…,n)，相对于t0的差分干涉图中像元(x，r)的差分相位φ(ti，x，r)是未知参数，其中x，r分别为像元的行坐标、列坐标。监测量为数据处理过程中获得的差分干涉相位Δφk(x，r)，k=1，2，…，m。

对第k(k=1，2，…，m)幅差分干涉相位图中像元(x，r)组成方程，可表示为

Δφk(x，r)=φ(tb，x，r)-φ(ta，x，r)

φtopo(x,r)+φorb+φres(x,r)

(2)

式(2)中：λ为雷达波长；ta和tb分别为第k幅差分干涉图对应的主、从影像的获取时刻；d(tb，x，r)、d(ta，x，r)分别为ta、tb像元(x，r)相对于t0的雷达视线方向(line of sight，LOS)地表形变；φtopo(x，r)为DEM数据引起的地形相位；φorb为不精确SAR影像数据导致的轨道相位；φres(x，r)为残余相位。

假设d(t0，x，r)=0，则相位时间序列为

(3)

假设SAR影像中某像元组成的向量为所求参数，可表示为

φ=[φ(t1)φ(t2)…φ(tn)]T

(4)

Δφk=[Δφ1Δφ2… Δφm]T

(5)

假设E、S分别为主影像序列和从影像序列，E={E1E2…Em}，S={S1S2…Sm}，且满足Ek>Sk,k=1，2，…，m。

则差分干涉图相位组成观测方程Δφk=φ(tEk)-φ(tSk)，k=1，2，…，m，可转换为

aφ=Δφk

(6)

式(6)中：a为m×n维矩阵。

假设Δφ1=φ4-φ2，Δφ2=φ3-φ1，Δφm=φn-φn-2，则矩阵a可表示为

(7)

当所有影像同一组即p=1时，m≥n，矩阵a的秩为n，对式(6)采用最小二乘法可求解出φ的估算值φ′，可表示为

φ′=(aTa)-1aTΔφk

(8)

计算SBAS-InSAR技术的木场古滑坡形变的基本流程如图12所示。

3 SAR数据

用于地质灾害领域分析的合成孔径雷达(SAR)卫星主要有日本的对地观测卫星(ALOS)、欧空局的Sentinel-1、德国的TerraSAR-X和中国开发的高分系列卫星数据。

选用2018年6月19日至2020年12月17日期间的Sentinel-1A降轨影像数据(图13)，共计58景，数字高程模型(digital elevation model，DEM)采用美国国家航空航天局(National Aeronautics and Space Administration，NASA)提供的30 m分辨率的航天飞机雷达地形测绘使命数据(shuttle radar topography mission DEM，SRTM DEM)，精密轨道校正数据采用欧空局提供的精密定轨星历(precise orbit ephemerides，POD)数据。SAR卫星影像数据信息如表1所示。

表1 木场古滑坡SAR卫星影像数据信息

图13 Sentinel-1A数据覆盖范围图

4 SBAS-InSAR数据处理

SBAS-InSAR技术对木场古滑坡区域Sentinel-1降轨影像数据处理技术流程如图12所示，其主要为对超级主影像进行地理编码配准，并划定研究区进行裁剪，生成连接图后，进行干涉处理，利用DEM数据移除干涉对进行地形相位，对差分干涉对进行自适应滤波，提高干涉数据质量，干涉图进行相位解缠处理采用最小费用流法(minimum cost flow，MCF)，去除大气相位，解算目标点时序形变，对结果进行地理编码，得到时序形变数据。

DEM为数字高程模型；GCP为地面控制点(ground control point);SAR为合成孔径雷达

4.1 生成连接图

本次干涉工作流处理输入58景Sentinel-1A的SLC数据。选择时间2019年12月11日的SAR影像为超级主影像，在整个处理中，超级主影像作为参考影像，然后进行影像配对。为了避免空间完全失相关，此处时间基线阈值调为365 d，生成连接图，得到各像对的时间基线图和空间基线图(图14)。

图14 时间和空间基线图

4.2 干涉处理

计算距离向、方位向值分别为4和1，以此获得方位向和距离向一致的分辨率。相位解缠采用MCF法进行，考虑到此次干涩对的整体相干性较好，不用对不理想相对进行剔除，解缠相关系数阈值设置为0.35，低于此相干值的区域不参与解缠计算。

采用30 m分辨率的SRTM DEM数据对干涉图进行掩膜处理，木场古滑坡地表起伏较大，侧视成像原理的雷达系统获取的数据会存在叠掩和阴影等现象，掩膜处理去除叠掩和阴影区域。

4.3 轨道精炼和重去平

对残余的恒定相位和解缠后还存在的相位坡道进行估算和消除，然后进行轨道精炼和重去平，通过斜距投影为选择GCP提供一个参考依据，在斜距数据上选择GCP定位，选择一个去平和滤波后的干涉图，可以在判断地形和形变区域时提供参考，此次选择30个点，利用GCP是对所有数据对进行重去平。

4.4 SBAS反演

SBAS技术需要两次反演，第一次反演用于估计形变速率和残余地形，第二次解缠优化输入干涉图。若选择无位移、二次方以及三次方模型，需要密集的连接图和高相干性才能获得准确的结果，因此选择了更稳定的线性模型进行反演。第二次反演对形变速率进行自定义大气滤波，估算并去除大气相位，得到最终的位移结果。

4.5 形变结果与地理编码

对SBAS反演的结果进行地理编码(地理坐标的统一)，将地理编码结果投影到研究区域，生成研究区平均形变速率(单位：mm/a)，如图15所示，累积形变量(单位：mm)，如图16所示。

5 木场古滑坡InSAR形变监测结果分析

5.1 木场古滑坡SBAS-InSAR形变分析

通过实际调查结合卫星影像(图17)可知，木场古滑坡区内地表植被以低矮灌木丛为主，覆盖度中等。SBAS-InSAR具有计算的相干点较多，干涉结果较好的特点，使得计算形变速率效果较好。

图17 滑坡卫星影像

根据木场古滑坡2018年6月19日—2020年12月17日期间雷达视线方向的形变结果，如图15所示，可以看出，木场滑坡正在发生强烈变形，形变向靠近SAR传感器方向运动时为正值；反之，形变向远离SAR传感器方向运动时为负值，可以看出滑坡中下部颜色较深，表明中下部变形较为强烈，滑坡体整体为负值，表明木场古滑坡整体上沿坡向方向向下滑移。依据滑坡形变速率分布与现场调查特征，将木场古滑坡分H1、H2和H3次级滑坡变形体，监测结果表明滑坡H1、H2和H3变形最为严重，最大形变速率均超过-25 mm/a，最大形变速率达到-26.52 mm/a。H1滑坡位于木场古滑坡南西侧，H2滑坡位于木场古滑坡北东侧，H3滑坡位于，木场古滑坡北东侧，在2018年6月19日—2020年12月17日期间最大累计形变量分别达到了-56、-49、-55 mm，如图16所示。

图15 InSAR形变速率

D1～D6为监测点

5.2 木场古滑坡变形规律与趋势分析

为研究时序的木场古滑坡形变特征，分别从古滑坡的次滑坡H1、H2、H3选取形变特征点，提取形变特征点的时序变形量，绘制累积形变量曲线，横轴为时间，纵轴为对应时间的累积形变量。由图18所示的变形特征点形变量曲线可知，2018年6月19日—2020年12月17日，滑坡H1呈缓慢下滑趋势，最大累积形变量达到了-56 mm；滑坡H2最大累积形变量达到了-49 mm；滑坡H3最大累积形变量达到了-55 mm；D4点位于木场古滑坡体上部，D5点位于古滑坡体中上部，D6点位于古滑坡体中下部，监测点D4、D5、D6累积形变量分别达到了-26、-37、-54 mm。选取的这6个特征点，可能在累积形变量上存在一定差别，但在整体都呈下滑趋势明显，时间维度的形变趋势基本保持一致，说明木场古滑坡体处于持续发育阶段，尤其是次级滑坡H1、H2、H3发育尤为强烈，形变量在3处次级滑坡区最大。

通过美空局获取的月降雨量数据，绘制木场古滑坡区月降雨量柱状图(图18)，可以看出，2018年6月19日—2020年12月17日，该区雨量较大的月份集中在5—9月，在2018年7月、2019年6月和2020年7月的降雨量达到峰值，分别为380、273、459 mm。比较监测点的累积形变量曲线可发现：形变量在6—9月的形变量增长较快，表明在这期间木场古滑坡的变形与降雨量呈正相关，降雨量越大，监测点的形变量越大，古滑坡的变形越明显。

图18 古滑坡变形特征点形变量曲线与降雨量关系

6 结论

获取2018年6月19日—2020年12月17日期间覆盖木场古滑坡的58景Sentinel-1降轨数据，基于SBAS-InSAR技术对木场古滑坡进行形变监测，得出如下结论。

(1)形变速率结果显示木场古滑坡中下部变形较为强烈，滑坡体整体为负值，表明木场古滑坡整体上沿坡向方向向下滑移。

(2)木场古滑坡分为3个次级滑坡变形体分别为H1、H2和H3，监测结果表明：古滑坡在H1、H2和H3区域变形最为严重，最大形变速率均超过-25 mm/a，最大形变速率达到-26.52 mm/a，最大累计形变量分别达到-56、-49、-55 mm。

(3)监测点变形趋势基本一致，总体呈下滑趋势，表明木场古滑坡体处于持续发展阶段，尤其是次级滑坡H1、H2、H3发育尤为强烈，形变量在3处次级滑坡区最大。

(4)在监测期间木场古滑坡区雨量较大的月份集中在5—9月，其中在2018年7月、2019年6月和2020年7月的降雨量达到峰值，分别为380、273、459 mm。通过比较降雨量与监测形变可发现：形变量在6—9月份的形变量增长较快，表明在这期间木场古滑坡的变形与降雨量呈正相关，降雨量越大，监测点的形变量越大，古滑坡的变形越明显。

(5)结合现场调查资料可知，H1滑坡在2020年8月19日，强降雨条件下加固失效，发生大面积滑坡；H2滑坡后缘及侧缘形成了明显错落坎，迹象明显，滑坡前缘的公路挡墙被挤压变形，造成部分挡墙损坏及形成不规则裂缝；H3前缘出现鼓胀裂隙，后缘发育拉伸裂隙。现场变形迹象验证了SBAS形变监测的可行性。

笔者认为今后的滑坡监测预警发展趋势，除了从地表的形变监测考虑，还应该结合岩土体本身的性质，如木场古滑坡的岩土体具有软化性，在强降雨后软化岩土体需要经历时间的过程，随着时间的变化，岩土体强度折减，当强度降到一个临界值时，斜坡稳定性系数小于1，就会导致滑坡的发生，进而研究位移监测时序变化与岩土体强度衰减时效性的关系，不同尺度揭示滑坡的灾变过程，探寻更加精准的滑坡监测预警方法。