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]。

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

对$ t_{\text{A}} $和$ t_{\text{B}} $（$ t_{\text{A}} < t_{\text{B}} $）时刻两幅影像生成的第${{i}}({{i}}= $$ \mathrm{1,2},\cdots ,N)$幅干涉图，在去除平地及地形相位影响后，第$ x $个像素的干涉相位可表示为^[14]：

式中：r——斜距；

${\varphi ^i}_{{\rm{def}}}\left( {x,r} \right)$——$t_{\text{A}} — t_{\text{B}}$对应的斜距向变形；

${\varphi ^i}_{{\rm{topo}}}\left( {x,r} \right)$——地形相位误差；

${\varphi ^i}_{{\rm{aps}}}\left( {x,r} \right)$——大气相位误差；

${\varphi ^i}_{{\rm{orb}}}\left( {x,r} \right)$——基线轨道引起的相位误差；

${\varphi ^i}_{{\rm{noise}}}\left( {x,r} \right)$——噪声误差。

假定t_k时刻和t_k+1时刻不同干涉图间的形变速率为${{{v}}_{{{k,k + 1}}}}$，则$t_{\text{A}}— t_{\text{B}}$间的累积形变可表示为式：

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

2.   研究方法

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

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

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

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

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

3.   滑坡灾害早期识别试验

3.1   研究区概况和数据来源

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

图  1  研究区位置

Figure  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所示。

表  1  Sentinenl-1A数据参数

Table  1.  Sentinenl-1A data parameters

轨道 方向	成像 模式	波段	波长	入射角 /（°）	极化 方式	距离 分辨率/m	方位 分辨率/m	重访 周期/d
升轨	IW	C	5.63	34.17	VV	5	20	12
降轨	IW	C	5.63	39.35	VV	5	20	12

下载: 导出CSV 

| 显示表格

图  2  研究区DEM

Figure  2.  Digital elevation model of study area

下载: 全尺寸图片 幻灯片

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  部分升降轨干涉图

Figure  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）方向位于右侧，能够很好地将峡谷两侧由西向东的地表形变监测出来，相反，降轨数据的飞行方向与之相反，能够将峡谷两侧由东向西的地表形变监测出来。为此，利用不同轨道数据可以互补，使得监测结果更为准确全面，能够避免单一轨道带来的几何畸变等问题。

图  4  升降轨获得的地表形变速率图

Figure  4.  Surface deformation rate map obtained by lifting rail

下载: 全尺寸图片 幻灯片

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

图  5  研究区NDVI

Figure  5.  NDVI in the study area

下载: 全尺寸图片 幻灯片

3.3   滑坡风险区的总体识别

深切割高山峡谷地区同样存在非滑坡原因的沉降变形，为保证获取的形变信息为滑坡灾害，仅凭形变监测结果无法区分是否为潜在的滑坡导致形变，为此，引入高分辨率光学影像对形变区域进行辅助识别，结合形变范围、高程、坡度、植被覆盖和坡体是否具有滑坡特征等进行识别，避免过度依赖形变结果导致的误判等问题。由于不同轨道获取的形变数据能够对研究区不同方向的滑坡进行早期识别，分别对不同轨道获取的数据结合光学影像解译，升轨数据解译结果如图6所示，黑色区域为滑坡形变区域，红色虚线区域为对应光学影像滑坡灾害隐患，黄色箭头代表滑坡方向,共解译出11处潜在滑坡，分别用${\rm{H}}01\text{、}{\rm{H}}02\text{、}\cdots\cdots\text{、}{\rm{H}}11$进行编号，其中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  升轨数据滑坡灾害识别结果

Figure  6.  Landslide disaster identification results of lifting data

下载: 全尺寸图片 幻灯片

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

图  7  降轨数据滑坡灾害识别结果

Figure  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潜在滑坡形变速率图

Figure  8.  Deformation rate diagram of H09 potential landslide

下载: 全尺寸图片 幻灯片

图  9  P1—P3时间序列曲线与降雨量

Figure  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。该潜在滑坡区最高形变点高程为3606 m，地势陡峭，从年形变量可以看出该区域存在很大安全隐患，一旦发生滑坡，岩子脚村村民将面临极大的生命及财产安全。

图  10  H06潜在滑坡形变速率图

Figure  10.  Deformation rate diagram of H06 potential landslide

下载: 全尺寸图片 幻灯片

图  11  P4—P5时间序列曲线与降雨量

Figure  11.  P4—P5 time series curve and rainfall

下载: 全尺寸图片 幻灯片

3.4.3   H16典型滑坡

H16属于典型的无滑坡痕迹的潜在滑坡，位于姑庄村和新寨田上方。图12为该潜在滑坡区形变速率与三维影像的叠置图，从图12中可以看到，位于潜在滑坡上方有明显的形变，最大形变速率为−80.141 mm/a。通过光学影像看到山体表面并无明显的滑动痕迹。该潜在滑坡区最高高程为2304 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。两个特征点形变曲线大致一直，从潜在滑坡上方到中部形变速率逐渐减小。该潜在滑坡一旦发生滑坡，将直接威胁姑庄村和新寨田两个村。

图  12  H16潜在滑坡形变速率图

Figure  12.  Deformation rate diagram of H16 potential landslide

下载: 全尺寸图片 幻灯片

图  13  P6—P7时间序列曲线与降雨量

Figure  13.  P6—P7 time series curve and rainfall

下载: 全尺寸图片 幻灯片

4.   结论

文中基于SBAS-InSAR技术，采用Sentinel-1升降轨数据结合互补的方式对东川小江沿线两侧深切割高山峡谷区滑坡灾害进行早期识别实验，得出以下结论：

（1）升轨数据识别出的滑坡主要分布在小江右侧，降轨反之，这表明相比单一轨道获取的识别结果，利用升降轨结合的方式能够更全面的监测和识别高山峡谷滑坡，避免单轨道对高山峡谷区滑坡进行早期识别存在SAR成像几何畸变造等问题。

（2）在高山峡谷地区植被覆盖度过高容易致使失相干，同时，有部分地区并不属于滑坡导致的形变，为此，引入光学影像，结合形变范围、高程、坡度、植被覆盖和坡体是否具有滑坡特征等进行识别，避免过度依赖形变结果导致的误判等问题。

（3）通过H09、H06、H16 等3个典型潜在滑坡灾害分析，可以看到滑坡的形成和降雨量有一定的关系。

（4）文中共识别出18处潜在滑坡区，其中H06、H09、H10、H11、H16等5处为滑坡高风险区，证明利用该方法可作为高山峡谷区滑坡灾害识别的有效手段，但目前还有一定的不足，比如在高山峡谷地区SAR数据量过少致使部分滑坡并未获取到形变量等问题。