基于SBAS-InSAR技术的深切割高山峡谷区滑坡灾害早期识别

2022-04-29周定义左小清喜文飞

中国地质灾害与防治学报 2022年2期

周定义，左小清，喜文飞,2，肖波,3，毕瑞，范馨

（1.昆明理工大学国土资源工程学院，云南昆明 650093；2.云南师范大学地理学部，云南昆明650500；3.云南交通职业技术学院公路学院，云南昆明 650500；4.成都理工大学地球科学学院，四川成都 610059）

0 引言

中国是地质灾害高发的国家，各类地质灾害给人们的生命财产造成了巨大的损失[1]。滑坡作为一种主要的地质灾害，具有隐蔽性强、危害性大、突发性高等特点，广泛分布于中国山区和峡谷地带[2−4]。近年来，高山峡谷滑坡频频发生。因此，对高山峡谷地区的滑坡灾害开展早期识别，能够为防灾减灾事业及政府部门决策提供一种有效的手段。

地表形变是反应当前坡体稳定性及运动状态最直接的物理量，因此，监测地表形变可以为探测滑坡等地质灾害的隐患点提供重要信息[1]。传统的地表形变监测方法主要采用精密水准测量、全球定位系统 (GPS)等。但这些监测方法存在变形监测工作量大、费时、花费大、测点难以保存等缺陷，同时，传统的监测方法无法对大区域滑坡形变进行探测[5−6]。与传统地表形变相比，合成孔径雷达干涉测量 (Interferometric Synthetic Aperture Radar,InSAR) 技术是近几年发展起来的一种新型大地测量手段。具有覆盖范围广、穿透云层、全天候运作、精度高的特点，理论上可以获得非常精确的数字高程模型和毫米量级的地表形变信息，已被成功用于滑坡灾害监测[7−9]。目前，利用InSAR技术对滑坡早期识别的研究取得了一些成功范例：张路等[2]采用自主研发的相干散射体时序InSAR 方法，成功识别出了17 处持续变形中的不稳定坡体。Dai 等[10]利用InSAR 对潜在的滑坡进行早期识别，10 个潜在滑坡区域被早期发现，结果表明，利用InSAR 可以作为获取滑坡早期识别的有效手段。韩守富等[11]证明了InSAR技术在黄土高原地区地质灾害早期识别方面的适用性和准确性,可以应用于黄土高原地区地质灾害隐患识别预警。冯文凯等[12]利用SBAS- InSAR技术对金沙江流域沃达村滑坡进行地表形变监测，表明SBAS- InSAR技术在复杂山区地质灾害监测预警领域有较为广阔的应用前景，为类似老滑坡监测预警提供了新的思路与借鉴。戴可人等[13]利用时间序列InSAR技术对雅砻江流域雅江县—木里县段的高山峡谷区域进行了滑坡灾害隐患广域早期识别，成功探测到8 处隐患区域。以上方法能够有效识别滑坡灾害，但对于深切割高山峡谷区的滑坡早期识别，仅利用SAR 单轨道数据监测会致使SAR 成像几何畸变造成部分滑坡不能识别，只有通过升降轨数据结合的方式才能全面准确的对滑坡灾害进行早期识别。

文章利用SBAS-InSAR技术，采用升降轨数据结合的方式对深切割高山峡谷区的滑坡进行早期解译识别，并根据识别结果,对不同潜在滑坡类型进行了分析与讨论，为高山峡谷地区防灾减灾及滑坡灾害早期识别提供一种更为全面的方法。

1 SBAS-InSAR技术简介

小基线集（SBAS-InSAR）是在差分InSAR 基础上发展起来的一种新的时间序列分析方法，能够降低相位噪声和误差[14]。

假定在时间t1至ts内获取同一地区的S幅SAR 影像，然后根据干涉组合条件，在短基线距的条件下形成N幅干涉条纹图，且满足：

对tA和tB（tA

式中：r——斜距；

假定tk时刻和tk+1时刻不同干涉图间的形变速率为vk,k+1，则tA—tB间的累积形变可表示为式：

对N幅干涉条纹图进行三维时空相位解缠即可求出不同SAR 获取时间的形变速率。

2 研究方法

文中利用SBAS-InSAR技术对深切割高山峡谷区滑坡灾害早期识别，方法的重点有两个关键：①针对深切割高山峡谷区的滑坡如何识别；②如何保证识别结果的准确性。

2.1 升降轨数据结合的方式对深切割高山峡谷区的滑坡识别

对于深切割高山峡谷区的滑坡早期识别，由于地形条件等因素的影响，仅利用SAR 单轨道数据监测会致使SAR 成像几何畸变造成部分滑坡不能识别。为此，在技术上采用升降轨数据结合的方式，对深切割高山峡谷区滑坡灾害进行早期识别。弥补了SAR 单轨道数据识别存在的不全面、不准确等弊端。

2.2 引入高分辨率光学影像等作为辅助识别保证准确性

为保证识别结果的准确性，仅凭形变监测结果无法区分是否为潜在的滑坡导致形变，为此，引入高分辨率光学影像对形变区域进行辅助识别，结合形变范围、高程、坡度、植被覆盖和坡体是否具有滑坡特征等进行识别，避免过度依赖形变结果导致的误判等问题。同时高分辨率光学影像能起到验证识别结果的作用，早期的滑坡一般会有一定的滑动痕迹，但不代表所有早期滑坡都有滑动痕迹。综上，引入高分辨率光学影像等作为辅助识别能有效保证识别结果的准确性。

3 滑坡灾害早期识别试验

3.1 研究区概况和数据来源

东川小江流域为世界典型暴雨泥石流区，被称为“泥石流的天然博物馆”，该地区泥石流的发生往往与滑坡密切相关，是属于暴雨型滑坡泥石流，即先滑坡，再经暴雨冲刷后形成泥石流[15]。文中选取小江沿线两侧高山峡谷作为研究区（图1），以小江为界，河谷凹陷，形成“V”字型，属于典型的深切割高山峡谷，东侧为牯牛寨山，最高峰海拔4 017.3 m；西部为拱王山，最高峰海拔4 344.1 m。该区域地势陡峭，其独特的地形和地质构造，致使局部区域暴雨多、土质松软、水土流失严重，导致该区域内地质灾害频发。

图1 研究区位置Fig.1 Location of study area

研究数据是从欧空局（European Space Agency,ESA）下载的40 景C 波段Sentinel-1A 升降轨影像，时间跨度为2018年4月25日—2019年4月20日，极化方式为VV，成像模式为IW。数据参数如表1所示。为了提高影像轨道精度，引入了POD 精密定轨星历数据，使用JAXA（日本宇宙航空研究开发机构）提供的ALOS WORLD 3D 30 m 空间分辨率的数字高程模型（Digital Elevation Model,DEM），用于去除地形相位影响。拼接后的DEM 如图2所示。

图2 研究区DEMFig.2 Digital elevation model of study area

表1 Sentinenl-1A 数据参数Table 1 Sentinenl-1A data parameters

3.2 地表形变监测试验

选取时间跨度2018年4月25日—2019年4月20日的Sentinel-1A 斜距单视复数（Single Look Complex，SLC）影像，其中升降轨数据各20 景。利用SARscape软件进行处理，升轨和降轨数据选取日期为2018-08-11和2018-08-13 的影像作为超级主影像，通过设置临界基线和时间基线，生成103 和100 对像对，设置多视数为1∶4 可以较好地抑制斑点噪声，采用Minimum Cost Flow 解缠方法和Goldstein 滤波方法做干涉工作流，最终生成干涉图，调整删除不理想的数据，在研究区生成较为理想的部分升降轨干涉图如图3所示。从图3 可以看出，针对山区地区，运用升降轨数据处理后得到的干涉图相干性较为理想。

图3 部分升降轨干涉图Fig.3 Part of the lifting rail interference diagram

经过轨道精炼和重去平，估算和去除残余的恒定相位和解缠后还存在的相位坡道，进行第一次反演、第二次反演，最后对序列信息进行地理编码获得研究区域2018年4月25日—2019年4月20日雷达视线方向（LOS）的形变速率图。如图4所示，正值表示形变朝着卫星的方向，负值表示形变沿着远离卫星的方向。图4（a）为采用升轨数据获取的地表形变速率图，图4(b)表示采用降轨数据获取的地表形变速率图。从图4 可以看出，采用升轨获取的形变主要分布在小江东侧，最大形变速率为−162.074 mm/a，降轨获取的形变主要分布在小江西侧西北方向，最大形变速率为−120.425 mm/a。致使不同轨道得到的形变结果不同，是由于升轨数据飞行方向大致从南到北，雷达视线（LOS）方向位于右侧，能够很好地将峡谷两侧由西向东的地表形变监测出来，相反，降轨数据的飞行方向与之相反，能够将峡谷两侧由东向西的地表形变监测出来。为此，利用不同轨道数据可以互补，使得监测结果更为准确全面，能够避免单一轨道带来的几何畸变等问题。

利用InSAR技术在植被覆盖地区进行地形测量非常困难，这主要是因为电磁场和（或）散射体的物理特性随时间的变化而造成的，植被覆盖度高地区，失相干严重，导致监测地面形变信息的能力较差，对形变监测精度的影响较大[16−17]。为考虑监测形变结果的有效性，引入归一化植被指数（Normalized Difference Vegetation Index，NDVI）对研究区植被覆盖数据进行分析。如图5所示，NDVI 结果被限定在[−1,1]，负值表示地面覆盖为云、水、雪等，对可见光高反射，0 表示有岩石或裸土等，正值表示有植被覆盖，且随覆盖度增大而增大。对比图4 和图5 可以发现，升降轨获取的部分形变区处于植被覆盖较高的地区，这正是由于植被造成失相干严重，致使这些区域监测结果不准确，在分析和识别滑坡时，需把其剔除，避免导致滑坡识别的误判。

图4 升降轨获得的地表形变速率图Fig.4 Surface deformation rate map obtained by lifting rail

图5 研究区NDVIFig.5 NDVI in the study area

3.3 滑坡风险区的总体识别

深切割高山峡谷地区同样存在非滑坡原因的沉降变形，为保证获取的形变信息为滑坡灾害，仅凭形变监测结果无法区分是否为潜在的滑坡导致形变，为此，引入高分辨率光学影像对形变区域进行辅助识别，结合形变范围、高程、坡度、植被覆盖和坡体是否具有滑坡特征等进行识别，避免过度依赖形变结果导致的误判等问题。由于不同轨道获取的形变数据能够对研究区不同方向的滑坡进行早期识别，分别对不同轨道获取的数据结合光学影像解译，升轨数据解译结果如图6所示，黑色区域为滑坡形变区域，红色虚线区域为对应光学影像滑坡灾害隐患，黄色箭头代表滑坡方向,共解译出11 处潜在滑坡，分别用H01、H02、······、H11进行编号，其中H01、H02、H08 和H09 位于现存泥石流沟处，H01 位于尖山沟老尖山、土岩子一带，H02 位于蒋家沟一带，H08、H09 位于大白泥沟和小白泥沟一带，其他潜在滑坡有一定的滑坡痕迹。识别出来的滑坡区域平均最大形变速率超过−40 mm/a，最大形变速率−120 mm/a，位于H06 岩子脚村上方。根据滑坡是否直接威胁附近周围村落，对11 处潜在滑坡进行危险划分，共识别出4 处高风险滑坡，分别是H06、H09、H10 和H11，H06 坡体对岩子脚村构成威胁，H09 坡体对大村和鲁纳窝村构成威胁，H10 对小多红村和妥托村构成威胁，H11 对小凹子村、厂上村、大麦地、小村子等构成威胁，这些区域一旦发生滑坡，极有可能对当地群众造成生命及财产威胁。

图6 升轨数据滑坡灾害识别结果Fig.6 Landslide disaster identification results of lifting data

降轨数据解译结果如图7所示，共解译出7 处潜在滑坡，分别用H12、H13、······、H18进行编号，滑坡区域平均形变速率均超过−40 mm/a，最大形变区域位于H18 大白泥沟和小白泥沟连接处，其最大形变速率为−94.11 mm/a。根据滑坡是否直接威胁附近周围村落，对7 处潜在滑坡区进行危险划分，共识别出1 处滑坡高风险区，编号为H16，该坡体对姑庄村和新寨田村构成威胁。对比图6 和图7 发现，不同轨道识别出来的潜在滑坡分布不同，利用升轨获取的滑坡主要分布在小江右侧峡谷中部，降轨获取的滑坡主要分布在小江左侧西北方向，其中不同轨道识别到的潜在滑坡有相同位置，但形变区域不一致，例如，H01、H13 共处于尖山沟一带，H02、H15 共处于蒋家沟一带，H08、H18 共处于大小白泥沟一带。由此证明，利用升降轨结合的方式能够有效识别深切割高山峡谷地区不同坡度方向存在的潜在滑坡，避免了单一轨道存在的识别结果不准确，不全面等问题。

图7 降轨数据滑坡灾害识别结果Fig.7 Landslide disaster identification results of rail descent data

3.4 典型灾害滑坡分析

经过对比光学影像发现，研究区识别出来的潜在滑坡可定义为三种类型，分别是处于泥石流区域的潜在滑坡（H09）、有滑坡痕迹的潜在滑坡（H06）和无滑坡痕迹的潜在滑坡（H16）。为有效了解不同类型潜在滑坡的形变趋势，分别对三种典型的潜在滑坡进行分析。

3.4.1 H09 典型滑坡

H09 属于典型的泥石流区域潜在滑坡，位于小白泥沟，该区域为暴雨型泥石流多发区。利用获取的形变速率值与三维光学影像叠置分析，其形变速率如图8所示，存在A、B 和C 三个潜在滑坡区，呈V 字型分布，滑动方向由坡面向下滑动，在每个潜在滑坡上方能明显的观察到断裂滑动痕迹，最大形变速率为−103.013 mm/a，位于潜在滑坡A 面上方。从图8 中可以看到，B、C 面分别位于小白泥沟两侧，其潜在滑坡属于长期滑坡，潜在滑坡上端形变速率大，下端由于滑坡堆积物导致一定的抬升。但A 面潜在滑坡还未形成真正意义的滑坡，一旦A 面滑坡，极有可能威胁大村和鲁纳窝村民的生命及财产安全。为更有效全面地分析该种类型潜在滑坡的特点，选取A 面进行详细分析，可以看到A 面呈现三处不均匀形变，有明显的拉张裂缝且有滑坡痕迹，其中P1 位于整个潜在滑坡顶端，P2 位于潜在滑坡中部，P3 位于潜在滑坡右侧。三个点均有可能在后期的降雨等因素下发展为新的潜在滑坡区。图9 为3 处潜在滑坡特征点形变量与当月降雨量的关系图，可以看到特征点的形变量与降雨量有一定的相关关系。研究区总体位于小江断裂带沿线，该区域断裂带宽5～20 km，呈现弱剪切强挤压活动特征[18]，以挤压穹起隆升变形为主[19]，从形变时序图9 可以看出，三个特征点在2018年5月先抬升后沉降，抬升可能是由于断裂带挤压活动致使，后经过雨水冲刷导致后期慢慢滑动沉降。P1、P2 和P3 点的平均形变量分别为−10.571 mm、−22.564 mm、−19.516 mm，年形变量分别为−49.063 mm、−58.740 mm、−62.635 mm。该区域月均降雨量为102 mm。随着时间的推移，形变量逐渐增加，到一定临界值便形成真正意义上的滑坡。

图8 H09 潜在滑坡形变速率图Fig.8 Deformation rate diagram of H09 potential landslide

图9 P1—P3 时间序列曲线与降雨量Fig.9 Time series curve of P1—P3 and rainfall

3.4.2 H06 典型滑坡

H06 属于典型的有滑坡痕迹的潜在滑坡，位于牯牛山岩子脚上方。图10 为该潜在滑坡形变速率与三维影像的叠置图，从图10 中可以白色虚线内两处潜在滑坡有明显的滑动痕迹，滑动方向为坡顶沿着坡脚移动。最大形变速率位于P5 附近，达到−127.093 mm/a。选取特征点P4、P5 形变量与月降雨量构建关系图（图11）。随着时间推移，形变量不断增大。特征点的形变量与月降雨量有一定的相关关系，从形变时序图11 可以看出，两个特征点在2018年5月先抬升后沉降，该潜在滑坡位于小江断裂带沿线，可能是由于断裂带挤压活动致使其抬升。当月降雨量增多加剧潜在滑坡的形变，这是由于降雨致使表层松散土体随着雨水滚落流失。P4、P5点的平均累积形变量分别为−37.562 mm、−42.054 mm，年形变速率分别为−105.903 mm/a、−112.469 mm/a。该潜在滑坡区最高形变点高程为3 606 m，地势陡峭，从年形变量可以看出该区域存在很大安全隐患，一旦发生滑坡，岩子脚村村民将面临极大的生命及财产安全。

图10 H06 潜在滑坡形变速率图Fig.10 Deformation rate diagram of H06 potential landslide

图11 P4—P5 时间序列曲线与降雨量Fig.11 P4—P5 time series curve and rainfall

3.4.3 H16 典型滑坡

H16 属于典型的无滑坡痕迹的潜在滑坡，位于姑庄村和新寨田上方。图12 为该潜在滑坡区形变速率与三维影像的叠置图，从图12 中可以看到，位于潜在滑坡上方有明显的形变，最大形变速率为−80.141 mm/a。通过光学影像看到山体表面并无明显的滑动痕迹。该潜在滑坡区最高高程为2 304 m，坡度大，滑动方向大致沿着山沟滑动。选取特征点P6、P7 形变量与月降雨量构建关系图，如图13所示，形变与降雨有一定关系。从形变时序图13 可以看出，P6、P7 两个特征点总是先抬升后沉降再抬升这一反复过程，在2018年5月可能是由于断裂带挤压活动致使其抬升，其他月份抬升微小可能是由于雨水冲刷导致土地堆积抬升，再经雨水冲刷致使微小沉降。P6、P7 平均累积形变量分别为−29.219 mm、−21.281 mm，年形变速率分别为−63.788 mm/a、−53.110 mm/a。两个特征点形变曲线大致一直，从潜在滑坡上方到中部形变速率逐渐减小。该潜在滑坡一旦发生滑坡，将直接威胁姑庄村和新寨田两个村。