高山峡谷区滑坡灾害隐患InSAR早期识别——以雅砻江中段为例

戴可人 铁永波 许 强 冯 也 卓冠晨 史先琳

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

②(成都理工大学地球科学学院 成都 610059)

③(中国地质调查局成都地质调查中心 成都 610081)

1 引言

雅砻江地处青藏高原向四川盆地过渡地带，属青藏高原东部边缘侵蚀山原区强烈-极强烈山原亚区，属嵌入式河流，沟谷深切，沟壑纵横，相对高差达2000～3000 m，河谷为“V”形谷，属典型的高山峡谷地貌[1]。雅砻江沿岸岩层风化强烈，切割破碎，岸坡稳定性差，加之沿岸水电的频繁开发以及大中型工程的兴建，该流域所遭受的滑坡、崩塌及泥石流等地质灾害也越来越严重[1—4]，对当地人民生命财产及水利水电等基础设施安全带来严重威胁。在雅砻江流域沿岸高山峡谷区域开展滑坡灾害隐患早期识别，对该区域的地质灾害防范以及确保人民生命财产与水利水电等基础设施安全具有重要的意义。

由于需要监测面积较大，传统地质调查手段在滑坡灾害隐患排查工作中很难达到大范围覆盖。随着对地观测技术的不断进步，光学遥感解译是目前滑坡灾害隐患大范围早期识别主要手段之一，但其受云雾影响较大，且往往只能对较大地貌特征进行定性解译[5]。作为近三十年的新兴雷达遥感技术，合成孔径雷达干涉测量技术(Interferometric Synthetic Aperture Radar,InSAR)[6]具有覆盖范围广、监测精度高、全天时、全天候、空间分辨率高等特点，已被成功用于滑坡灾害监测(如文献[7—11])，近年来被逐步应用于灾害隐患早期识别并取得一些成功案例。张毅[12]基于Envisat ASAR数据，应用时间序列InSAR技术成功识别出白龙江流域133处活动斜坡；张路等人[13]利用InSAR技术成功识别出大渡河区域17处持续变形中的不稳定坡体；张亚迪等人[14]利用时序InSAR对芒康地区滑坡灾害隐患进行了探测；Dai等人[15]利用短基线DInSAR技术对汶川附近国道G317沿线不稳定滑坡体进行了有效识别；陆会燕等人[5]结合InSAR与光学遥感技术，成功探测出金沙江下游区域共计7处具有较显著形变的滑坡隐患。高山峡谷区域海拔落差大、地形复杂、植被茂密，给InSAR处理带来干涉失相干、大气延迟、几何畸变等挑战[16]。对于InSAR高山峡谷区域的滑坡识别效果与适用性还值得进一步分析与研究。

本文利用小基线集时间序列方法(Small BAseline Subset InSAR,SBAS-InSAR)，对雅砻江流域雅江县至木里县段高山峡谷区域进行了滑坡灾害隐患广域早期识别的研究，探究其潜在的滑坡隐患区域，并对识别结果进行进一步的分析讨论，分析雷达遥感在高山峡谷区滑坡探测中的效果与适用性，可为当地的防灾减灾提供有力的支持，并为高山峡谷区的滑坡灾害隐患早期识别提供思路与参考。

2 研究区域与数据—研究区域位置与地势特点

雅砻江是金沙江最大支流，为横断山区北南向的主要河系之一[1]。本文选取雅砻江流域雅江县-木里藏族自治县段作为研究区域(如图1所示)。该区域位于雅砻江中游段，四川省甘孜藏族自治州、凉山彝族自治州境内，处于中国最高一级阶梯向第2级阶梯云贵高原和四川盆地过渡地带，属横断山系北段川西高山高原区，亦属于青藏高原的一部分。该段位于甘孜阿坝褶皱带，出露地层以三叠系浅变质砂板岩为主，有少量花岗岩及石灰岩分布。分布巨厚的中上三叠系浅变质岩系，砂岩，板岩构成北西-南东向紧密褶皱，褶皱轴部及断层带中有少量二叠系灰岩分布，并有零星燕山期花岗岩出露，且出露下古生界至上古生界碳酸盐岩类、浅变质岩及玄武岩等[1—4]。研究区域海拔3000～5000 m，峰峦重叠、沟壑纵横、河谷深切，为典型高山峡谷地貌，为崩塌、滑坡等斜坡地质灾害的发生孕育了条件。大部分山区沟道巨大的地形高差，使处于高处的松散碎屑物质拥有巨大的势能，陡急的山坡和沟床为坡面和沟床松散堆积物能量的释放和势能转化为动能提供有利条件，有利于滑坡、泥石流等灾害的形成，导致该区域内滑坡、崩塌、泥石流等地质灾害频发[17—20]。

图1 研究区域Fig.1 Study area

本研究共计使用了50景Sentinel-1卫星数据(每一个时间点需要两景数据上下拼接)，时间范围覆盖2017年9月9日—2018年10月10日共计396天。Sentinel-1是欧洲航天局(European Space Agency,ESA)发射的由两颗卫星组成地球观测卫星星座，载有C波段合成孔径雷达，可提供白天、夜晚及各种天气下的连续影像。表1列出了Sentinel-1卫星主要参数。

表1 Sentinel-1卫星SAR影像数据主要参数Tab.1 Main parameters of Sentinel-1 SAR datasets

Sentinel-1影像数据集的时空基线图如图2所示，时间基线阈值为48 d，空间基线阈值小于150 m。为了消除与减弱由轨道误差引起的相位误差，在数据处理过程中还使用了欧洲空间局精密轨道数据对轨道信息进行纠正。同时采用日本宇航局ALOS WORLD 3D 30 m空间分辨率的数字高程模型(Digital Elevation Model,DEM)，以消除或减弱地形相位的影响。

图2 SAR数据集时空基线图Fig.2 Spatial and temporal baselines of SAR datasets

3 时间序列InSAR方法

传统的差分干涉测量技术(Differential InSAR,DInSAR)易受到时间失相干和空间失相干的影响，不能准确获取形变信号。为降低时空失相干的影响，Berardino等人[21]提出了短基线集，即SBASInSAR技术，其核心思路是将同一个地区多幅SAR影像配准得到短基线对，采用奇异值分解法(Singular Value Decomposition,SVD)将多个短基线集联合求解，得到目标区域的整个时间位移序列以及地表形变平均速率。具体原理及流程如下：

假设获取了覆盖同一区域N+1幅SAR影像，按时间序列排序如式(1)

按照一定规则，选取某一影像作为主影像进行配准，得到M个干涉对，M满足式(2)

通过设定时空阈值，得到相干性较好的M对干涉对，再通过精密轨道文件及高精度数字高程模型去除轨道误差、平地效应及地形相位的影响后得到M幅差分干涉图以及解缠后的形变相位。

假设t0时刻作为起始时间，将该时刻研究区域的位移视为0，则第i幅干涉图中(1≤i≤m)的某个像素相对于起始点的相位可以表示为式(3)

其中，φitopo表示地形相位，φiatm大气延迟造成的相位，φinoise相干噪声引起的相位。分别可表示为

λ为波长；θ为雷达入射角；Δh为外部DEM误差；r为雷达到目标物体的斜距；d(t2)与d(t1)分别表示t2与t1时刻像元相对于起始时间t0在雷达视线方向的形变累积量。假设干涉对影像时间间隔内地表形变满足线性变化，将对相位时间序列的求解转变为相位变化速率的求解，可得式(5)

则Δφdef可表示为：Δφdef=B，其中，矩阵B为系数矩阵，每行一一对应干涉像对。矩阵的元素中主影像系数为1，辅影像系数为—1，其他系数为0。

因此，在不考虑大气相位与噪声的情况下，联合式(1)—式(5)可得式(6)

图3 SBAS-InSAR时序分析流程图Fig.3 Flowchart of SBAS-InSAR time series analysis

4 滑坡灾害隐患早期识别结果

4.1 总体识别结果

通过短基线集SBAS-InSAR时间序列方法对雅砻江中段雅江县-木里县流域进行了形变监测。考虑到数据处理过程误差的累积与处理的效率，将雅江县-木里县流域分为5段，分别提取了各段的雅砻江沿岸坡体形变特征信息，获得的雷达视线方向上的年均形变速率如图4所示，其中红色负值代表的是目标地物远离卫星方向位移，蓝色正值代表的是目标地物靠近卫星方向位移。5段的详细划分如图4(a)所示，具体划分为雅江县-恶古乡(A段)，恶古乡-牙衣河乡(B段)，牙衣河乡-三岩龙乡(C段)，麦地龙乡-田镇村(D段)，卡拉乡-三桷垭乡(E段)。本文主要针对雅砻江沿岸对河流有危险的坡体进行In-SAR结果解译，早期识别地质灾害隐患，并进行相干性、几何畸变等复核检验，共识别出8处正在发生蠕变的不稳定坡体，其中B段4处，即鲁日、日阿、日衣、木恩；D段共中铺子村、麻撒村2处；E段共阳山村、独家村2处。基于SAR成像参数与坡度坡向数据，我们对8处坡体进行了几何畸变分析与验证(依据5.2节、图8与图9)，这8处坡体均处于非几何畸变区，形变监测结果可靠。

本文从变形范围、最大形变速率、坡体平均坡度、高程范围、植被覆盖情况、威胁对象这几个方面，对这8处正在变形区域进行了详细分析。根据隐患点的形变量级，坡体是否具备滑坡地质地貌特征、形变速率是否有加速、突变现象、是否直接对雅砻江或附近村落直接造成威胁作为判断依据，对风险等级进行了定性评估。其中5处隐患点处于高风险等级，倘若发生地质灾害，极大可能会对当地群众人身安全造成威胁，属于重点观测点值得特别关注。这5处主要分布于B段恶古乡-牙衣河乡与D段麦地龙乡-田镇村区域。恶古乡-牙衣河乡有3处风险等级高的区域，其中鲁日坡体，对鲁日村及村中的雅江县八衣绒乡中心小学、下次呷、次呷和上次呷的安全造成严重的威胁，并且具有堵塞雅砻江的风险；日阿坡体对附近的日阿村、科鲁村、道路与雅砻江构成直接的威胁；日衣坡体对日衣村与雅砻江构成直接的威胁。麦地龙乡-田镇村有2处风险等级高的区域，其中中铺子村坡体对麦地龙乡、中(铺子村与雅砻江构成直接的威胁；麻撒村坡体对附近的村落与雅砻江构成直接的威胁。详细结果列表如表2所示。从整体识别结果来看，InSAR技术在植被覆盖较低的区域，可以精确有效地识别出高山峡谷区滑坡灾害隐患区域，所发现的几处隐患点均为中低植被覆盖区。就此研究区域而言，滑坡灾害隐患区域最大形变速率普遍超过50 mm/year，且呈现出高位特征，平均高程大于2000 m，平均坡度大于55%，均可对当地村落与雅砻江构成直接的威胁。

图4 InSAR年均平均速率监测结果Fig.4 Mean velocity map derived from InSAR

4.2 重点区域分析解译

鲁日村与中铺子村是此次识别中具有直接威胁对象的两处重点区域，其所在的坡体年均形变速率如图5、图6所示，图中圆点代表监测InSAR有效相干点，点颜色代表年平均形变速率，其中绿色代表较弱的年平均形变速率，亮色彩(黄色、红色)代表显著的远离卫星视线运动，具体量级如图中颜色条所示。

鲁日村所在坡体位于雅砻江东岸(如图5(a)所示)，坡体植被较少，整体相干性较好。其在地形地貌上呈现阶梯状特征(如图5(b)所示)，其中鲁日村位于海拔2986～3085 m，上方存在高位的坡体；下次呷、次呷和上次呷位于鲁日村坡体的下方。根据形变监测的结果表明，鲁日村的形变主要位于鲁日村以及鲁日村至下次呷之间的坡体(图5(a)、图5(c)中的蓝色虚线所示)。若发生滑坡灾害，潜在的滑坡方向如图5(b)中红色箭头所示，可能对鲁日村及村中的雅江县八衣绒乡中心小学、下次呷、次呷和上次呷的安全造成严重的威胁，并且滑坡松散堆积物质可以直接经下次呷冲入雅砻江，具有堵江的风险。

表2 雅砻江潜在灾害点早期识别结果列表Tab.2 Early identification results on potential disaster sites along Yalong River

图5 鲁日村滑坡灾害隐患识别结果Fig.5 Early identification results of potential landslide geohazards in Luri village

图5(c)展示了3处特征点的时间序列形变曲线图，其中点P1位于鲁日村，点P2位于下次呷上方坡体，点P3位于西侧坡体。时序结果表明，3处累积形变量有所差异，但其趋势相似。下次呷(P2点)的累积形变最大，达到58 mm；鲁日村(P1点)累积形变次之，达到51 mm；西侧坡体(P3点)累积形变最小，达到44 mm，可见坡体整体都存在明显的形变。3个特征点的形变速率自2017年9月—2018年2月与2018年8月—2018年10月较缓，2018年2月—2018年8月较快。通过中国气象数据中心查询甘孜州地面观测站发现，当地气温在3月—10月月气温温度大于0 ℃且为当地雨季，气候变化可能是造成形变加速的主要自然原因。

中铺子村所在坡体位于雅砻江西岸，坡顶高程3000 m，坡底高程1900 m，平均坡度约60.4%，属于高位陡峭坡体。如图6(a)所示，由于坡体中上部植被较少，相干性较好，因此坡体整体中上部区域相干点的数量相对较多，但在河底部位，由于坡度更大，达65%，相干性点的数量较稀疏。坡体中部出现一处核心形变区，其范围约1700×1000 m。最大视线向形变量达60 mm/year。该区域的遥感解译结果如图6(b)所示，在两处红色虚线范围内存在古滑坡体，坡度较大，达到约65%，是坡顶存在不稳定的直接表现之一。在坡体的北侧，更是存在巴多沟崩塌、滑坡体，间接证明该坡体长期存在不稳定性。

图6(c)展示了3个特征点的时间序列形变曲线图，其中点P1位于坡体上部，点P2位于坡体下部左侧，点P3位于坡体中部。时间序列结果表明，3处形变其趋势相似，3个特征点的形变速率自2017年9月—2018年10月，形变量持续增加，形变速率稳定，呈现线性形变。P1,P2点的累积形变最大，达到120 mm；P3点累积形变次之，达到55 mm，整体都存在明显的形变。人类活动主要集中在江对岸，所以植被的稀疏与雨季降雨量的增大影响是造成坡体有滑坡趋势的主要原因。

图6 中铺子村滑坡灾害隐患识别结果Fig.6 Early identification results of potential landslide geohazards in Zhongpuzi village

5 识别结果分析验证与雷达影像适用性分析

5.1 识别结果对比验证

根据蠕变滑坡位移曲线特征[9]，InSAR形变结果是潜在滑坡灾害点早期识别的重要特征之一。这些潜在滑坡灾害点部分是发育于古滑坡上的不稳定坡体，部分滑坡表面不具备滑坡继续发展的地貌特征，还需要对这些坡体进行进一步验证与分析。在得到这些灾害隐患点后，接下来需要对InSAR早期识别结果与滑坡隐患历史资料(含地质调查资料、文献、滑坡发生报道等)、光学影像遥感解译结果进行对比(图7)。

雅砻江流域各个监测段共探测到形变坡体8处依次由北向南编号为1～8号(表3)。其中共有4处相对历史记录为新发现(占识别结果50%)，标记为新的潜在地质灾害点，其中有7处具有明显滑坡地貌特征(占识别结果87.5%)。我们以前述变形范围、最大形变速率、坡体平均坡度、高程范围、植被覆盖情况、威胁对象等作为指标，对这8处正在变形区域进行了危险等级快速定性评估，其中3处中危险等级潜在灾害点，5处高危险等级潜在灾害点。高危险等级潜在灾害点具备滑坡地貌地形特征，形变速率较大，直接威胁雅砻江或相近村落，需要对这些识别出的灾害点进行进一步调查及监测。

5.2 雷达影像几何畸变分析

图7 潜在地质灾害点验证与对比Fig.7 Verification and comparison on potential landslide geohazards

由于星载SAR采用侧视成像的方式，雷达波束斜向照射地表时会导致雷达图像出现距离向透视收缩、阴影或倒置等几何畸变[16]，在高山峡谷区域几何畸变尤为严重，将会造成无效的监测盲区[12,13]。雷达卫星在面对不同朝向、坡向时所会形成的各种几何畸变情况如图8所示。可以看出当坡面朝向卫星时，若坡角小于入射角时，沿坡面向下的变形在LOS方向上会表现为靠近卫星(图8情况①)，且坡体会发生透视收缩。当坡角与入射角相等时，会发生完全透视收缩，坡体无法被测量(图8情况②)。如果坡角过大(超过卫星入射角)，则会发生顶底倒置的成像，沿坡面向下的变形在LOS方向上会表现为远离卫星(图8情况③)。当坡面背向卫星时，若坡角较大(大于入射角的余角)，此时坡面会位于阴影区域，无法被卫星照射产生回波信号，整个坡面无法被测量(图8情况④)。当坡角与入射角互余时，此时所获得地距分辨率等于卫星斜距分辨率(图8情况⑤)。若坡角较小(小于入射角的余角)，此时将获得最接近于斜距分辨率(slant range resolution)的地距分辨率(ground range resolution)[23]，此时的地距分辨率将高于卫星照射平地情况，因此该朝向及坡度的坡体是最适合SAR卫星进行观测的情况。在实际应用中，联合升降轨的SAR影像，可以获取不同成像几何下的监测目标的形变特征，从而在一定程度上补偿单一成像几何带来的观测盲区，有利于提升滑坡灾害隐患早期识别的有效探测率[13,16,24]。

表3 雅砻江流域InSAR早期识别验证与对比Tab.3 Verification and comparison on early identification results from InSAR along Yalong River

图8 卫星入射角(LOS)观测方向与沿坡向形变关系(修改自文献[22])Fig.8 Relationships between the Line Of Sight (LOS) and the downslope displacements for different slope orientations (adapted from Ref.[22])

从定量关系上来看，根据Sentinel-1卫星轨道参数信息与坡体坡度信息进行测算，Sentinel-1升轨数据的卫星飞行方向是自南向北，沿着方位向—12.6°飞行，雷达入射角为36.8°。由于SAR采用侧视成像，当坡体朝向卫星的时候，若坡度角处于0～36.8°时，会出现透视收缩(距离压缩)的几何畸变；当坡体背向卫星的时候，若坡度角的绝对值大于53.2°时，会出现阴影的几何畸变。Sentinel-1降轨数据的卫星飞行方向是自北向南，沿着方位向—167.4°飞行，雷达入射角为39.7°。由于SAR采用侧视成像，当坡体朝向卫星的时候，若坡度角处于0～39.7°时，会出现距离压缩的几何畸变；当坡体被向卫星的时候，若坡度角的绝对值大于50.3°时，会出现阴影的几何畸变。

图9为基于升降轨Sentinel-1影像获取的雅砻江升降轨几何畸变图，红色代表的是距离压缩的区域，蓝色代表的是阴影区域。在雅砻江流域中，升轨数据中的距离压缩区域主要位于东南朝向的坡体，且主要分布在雅砻江的东岸；降轨数据中的距离压缩区域主要位于东北朝向的坡体，且主要分布在雅砻江的西岸。阴影区域主要分布在较为陡峭的背向雷达信号侧的坡体。

图10(a)、图10(b)分别为雅砻江日衣村-木灰村段Sentinel-1升轨数据、降轨数据几何畸变图，图10(c)、图10(d)分别为该区域坡度、坡向图。在升轨数据图10(a)中，红色区域相对卫星的入射角处于0°～36.8°，主要分布在坡面②与坡面④，代表的是距离压缩区域，对应图10(d)中雅砻江东岸褐色的区域；蓝色区域的入射角大于53.2°，代表的是阴影区域，主要分布在坡面①与坡面③，对应图10(c)中红色的大坡度的区域。

在降轨数据对应的图10(b)中，红色区域相对卫星的入射角处于0°～38.7°，主要分布在坡面①与坡面③，代表的是距离压缩区域，对应图10(d)中雅砻江西岸绿色的区域；蓝色区域的入射角大于50.3°，代表的是阴影区域，主要分布在主要分布在坡面②与坡面④，对应图10(c)中红色大坡度的区域。

图9 雅砻江流域Sentinel-1数据几何畸变分布图Fig.9 Geometric distortion of Sentinel-1 datasets along Yalong River

图10 雅砻江日衣村-木灰村段Fig.10 Yalong River Riyi village-Muhui village section

可以明显地看出，在坡向较复杂的区域，易出现严重的距离压缩几何畸变。在坡度较大的区域，易出现阴影区域，无法获取有效的监测信号。就雅砻江日衣村-木灰村区域而言，升轨数据在雅砻江东岸存在严重的距离压缩几何畸变，主要分布在坡面②与坡面④；在雅砻江高坡度的陡坡主要分布在坡面①与坡面③)，可能出现严重的阴影区域；降轨数据在雅砻江西岸存在严重的距离压缩几何畸变，主要分布在坡面①与坡面③；在雅砻江高坡度的陡坡(主要分布在坡面②与坡面④)，可能出现严重的阴影区域。因此，在现有SAR卫星系统飞行轨道差异性有限的情况下，高山峡谷区域单轨道监测可能会由于SAR成像几何畸变造成部分不稳定坡体被“漏检”，只有通过升降轨数据结合的方式，才能较为全面覆盖监测区域，在一定程度上补偿单一成像几何带来的几何畸变问题，才能实现全面准确的滑坡灾害隐患早期识别。

6 结论

本文使用Sentinel-1雷达影像，通过时间序列InSAR技术对高山峡谷区域雅砻江流域雅江县-木里县段进行了滑坡灾害隐患广域早期识别，成功探测到8处隐患区域，最大雷达视线向年平均形变速率可达95 mm/year，这些区域若发生滑坡，对当地居民的生命财产安全存在严重的威胁，同时存在堵塞雅砻江的风险。同时，利用滑坡隐患历史资料与遥感解译对InSAR早期识别结果进行了进一步的验证与对比分析，综合变形范围、最大形变速率、坡体平均坡度、高程范围、植被覆盖情况、威胁对象等情况对这八处滑坡进行了危险等级评价。

由于星载SAR采用侧视成像的方式，不可避免的会产生几何畸变，在高山峡谷区尤为严重。本文对不同坡度坡体几何畸变情况进行了详细分析，展示了Sentinel-1升降轨数据在雅砻江高山峡谷区的具体几何畸变区域与适用性。在地形复杂的高山峡谷区域，只有通过采用不同轨道数据结合的方式补偿单一成像几何带来的几何畸变问题，才能有效避免“漏检”，实现全面准确的滑坡灾害隐患早期识别。本文以雅砻江为例，揭示了时序InSAR技术运用于滑坡灾害隐患早期识别的有效性与优势，可为当地的防灾减灾提供有力的支持，结合几何畸变的详细分析讨论，为高山峡谷区的滑坡灾害隐患广域早期识别提供范例与参考。