基于InSAR技术的川西高山峡谷区地质灾害早期识别研究——以小金川河流域为例

唐尧， 王立娟， 廖军， 邓琮

(1.四川安信科创科技有限公司，四川 成都 610045; 2.四川省安全科学技术研究院，四川 成都 610045； 3.四川省地震与地质灾害应急技术保障中心，四川 成都 610045； 4.四川交通职业技术学院，四川 成都 611130)

0 引言

我国地质灾害多发频发，灾害事件总量偏大。川西地区地处第一地势阶梯与第二地势阶梯的过渡带，区内安宁河—鲜水河—龙门山等构造活动活跃，地震频发，断块的差异化隆升以及河流强烈的下切侵蚀作用使得高山峡谷地貌发育，同时，大量的震裂山体与高陡河谷岸坡临空面的侧向卸荷作用使得该区域地质灾害异常发育，受灾害影响严重。频发的地质灾害事故导致了大量的人员伤亡与财产损失，对社会民生造成了严重影响，极大制约了当地社会经济的发展[1]。前期的灾害隐患识别是开展灾中应急救援处置、防灾减损的有力途径[2]，但受川西高山峡谷地区恶劣的交通与地形条件限制[3]，以及灾害隐患信息化建设的滞后，地质灾害的预防与应急救援工作面临极大的挑战，因此亟需开展该类地区地质灾害早期识别关键技术的研究[2,4]。

目前，地质灾害早期识别主要采用野外人工排查、光学可见光遥感、无人机低空航摄等技术，缺乏对灾害分布、形变特征与时间演化特征过程的分析[5]。新型遥感对地观测技术的出现，为实现地质灾害隐患的早期识别提供了重要的技术支持，其中合成孔径雷达干涉测量(interferometric synthetic aperture radar，InSAR)作为一种新型的空间对地观测技术，具有全天时、全天候、大面积观测地表精确形变信息的能力，近年来已广泛应用于我国西南山地、西北高原及三峡库区等区域的地质灾害监测识别研究中[1,2,5]。本文采用InSAR技术，利用多时相合成孔径雷达数据，对川西高山峡谷区开展地表多时相、长时序形变监测，旨在实现小金川河流域地质灾害早期隐患的识别分析。

1 研究区概况及数据源

1.1 研究区概况

研究区位于四川省西部重要水系小金川河流域(图1)，属于大渡河上游重要支流，总面积约1 171 km2，区内地形属典型的高山峡谷地貌(强震频发、河谷深切、水系发育)，行政上隶属于甘孜州丹巴县与阿坝州小金县，包括墨尔多镇、半扇门镇、太平桥乡、宅垄镇、新桥乡、沙龙乡、美兴镇、崇德乡及八角镇等乡镇。流域最高海拔5 080 m，最低海拔1 847 m，相对高差约3 200 m，形成25 °～70 °的坡度变化。域内有重要国道G350分别通往丹巴县城与小金县城，经G248、S210通往国道G317与G318。

区内地质环境脆弱，滑坡、崩塌、泥石流等地质灾害频发(图1)，当地自然资源部门提供的资料表明，区内分布有滑坡、崩塌、泥石流等地质灾害共400余处，主要沿小金川河沿岸及次级支流展布。随着近年来人类活动的日益增强，受小金川河沿线梯级水电开发及相关基建项目施工的影响，流域内地质灾害发育有加剧的趋势，多处灾害复活变形或形成灾害链，如2020年的“6·17”梅龙沟泥石流-堰塞湖灾害链、阿娘寨复活滑坡群及“7·6”城隍庙沟泥石流等影响较大的典型地质灾害，给当地居民生命财产造成了重大损失。

图1 研究区地理位置(左)及流域地质灾害分布(右)Fig.1 Geographic location of the study area(left) and geological hazards in the watershed(right)

1.2 数据源

InSAR分析采用的遥感数据主要为欧空局(European space agency, ESA)的Sentinel-1A影像数据及精密轨道数据(precise orbit ephemerides, POD)，包括26景C波段数据，IW模式单视复数影像(single look complex，SLC)，时相为2018-11—2019-12。轨道方向为降轨，轨道号为135，方位向分辨率9.32 m，距离向分辨率13.97 m，重访周期12 d，飞行方位角6.7°，视线方位角为-83.3°，视线入射角为39.5°，极化方式VV。差分干涉处理选用30 m空间分辨率的ASTER GDEM作为参考DEM(digital elevation model)，该数据由美国航天局(NASA)主导测量，其标称高程精度约为10～25 m。

影像数据需做必要的预处理，本次Sentinel-1A雷达影像数据导入与裁剪采用PIE-SAR 6.0软件完成。PIE是航天宏图自主开发的一款专业遥感图像处理软件，提供了面向多源、多载荷的遥感图像处理、辅助解译及信息提取功能，是一套高度自动化、简单易用的遥感工程化应用平台。

2 短基线集时序InSAR方法

2.1 时序InSAR方法原理

短基线集时序干涉测量方法(SBAS-InSAR)，最早于2002年由意大利学者Berardino等人提出，之后经多位学者加以改进完善，其主要利用多期影像对地表形变进行反演，对区域地表形变的时间演化有很好的分析能力。该方法通过时间基线和空间基线的阈值，将获取的SAR影像进行自由配对组合，生成多个组合，去除地形相位，生成一系列差分干涉图； 对差分干涉图的像元基于相干性筛选出高相干点，对这些高相干像元进行相位解缠，利用地面控制点进行轨道精炼，去除相干性较低和解缠错误的干涉像对，通过奇异值分解法和最小二乘法解得影像获取时刻的累计形变； 最后，利用时域的高通滤波和空间域的低通滤波去除大气影响，得到最终的形变结果。其方法原理为： 获取同一区域按照时间顺序t0,t1,…,tn排列的(n+1)幅SAR影像，选取其中一幅影像作为主影像，将其他SAR影像配准到这幅影像上，生成j幅多视差分干涉图。

该方法具有不受天气制约、可大面积监测等优势[1,6-7]。

2.2 时序SBAS-InSAR数据处理

时序InSAR方法在山区地质灾害早期识别中有较广泛的应用，可获得沿雷达视线方向(line of sight，LOS)的地表形变信息，而对于山区来说，LOS向的形变速率不足以真实反映斜坡的真正形变情况[6-8]。本节采用SBAS-InSAR的方法完成小金川河流域数据分析，数据处理基本流程包括数据(SAR与DEM)预处理、差分干涉计算、时间空间域形变估算及形变量计算等步骤，如图2所示。

图2 SBAS-InSAR数据处理基本流程Fig.2 Basic flow chart of SBAS-InSAR data processing

2.2.1 数据预处理

数据预处理包括计算所有Sentinel-1A影像对之间的时间和空间基线，经配准裁剪后，利用流域内DEM数据完成影像精配准，选择满足给定阈值的相对组合生产差分干涉图集[9-10]。本节利用欧空局发布的对应精密轨道数据(POD)和30 m分辨率的SRTM DEM进行影像轨道误差和地形相位校正。为了提高时间、空间相干性，保证干涉图质量，将最大空间基线阈值设为5%，最大时间基线阈值设为90 d，进行3D解缠。超级主影像为2018-11-01，共生成了230个干涉像对，数据配对情况如图3所示。

图3 数据集像对配对图与时空基线分布

2.2.2 差分干涉计算

差分干涉计算包括平地与地形相位去除、差分干涉图滤波、相关系数计算机相位解缠等内容。相位去除依据空间基线参数和地球椭球体参数计算平地相位，之后利用配准后的DEM计算地形相位，从干涉相位中去除平地和地形相位，获得差分干涉相位，逐像元计算差分干涉图[8-10]。对生成的120个干涉像对进行干涉处理，主要包括了相干性生成、去平、滤波和相位解缠等步骤。为了得到方位向和距离向一致的地面分辨率，将多视的视数比设置为1∶4； 为了削弱斑点噪声的影响，采用可以提高干涉条纹清晰度的Goldstein滤波进行噪声去除； 由于小金川河流域内地形起伏大，森林覆盖度较高，选用可以减小植被等因素的Delaunay MCF方法进行相位解缠。

2.2.3 时间空间域形变估算

对干涉图的差分干涉相位进行时间域的线性变形相位估计，去除大气、噪声等残余相位，获得点目标的时间序列变形相位。经过相邻点间的参数估计、线性变形相位与残余高程相位计算、残余相位低通滤波、奇异值分解处理(singular value decomposition，SVD)、大气相位和非线性变形相位计算及时间序列变形相位计算等步骤[8-10]，在解缠后的干涉图像上选取75个远离形变区、相对稳定的点作为地面控制点(ground control point，GCP)。考虑本次监测时间周期较短(仅1年)等因素，将精炼次数设置为3，对残余的恒定相位及解缠后还残留的相位坡道进行去除。

2.2.4 形变量计算

形变量计算主要根据反演结果分析实现，反演是SBAS(small baseline subset)处理的核心步骤。第一次反演估算形变速率和残余地形，同时也会做二次解缠用来对输入的干涉图进行优化。第二次反演是在第一步的基础上，利用高通滤波和低通滤波估算和去除大气相位，得到更加纯净的时间序列上的最终位移结果，通过地理编码，最终获取研究区形变速率分布图。

3 地质灾害识别结果分析

3.1 雷达可视性分析

由于SAR采用侧视成像方式，即雷达波束是以一定的入射角照射地表，小金川河流域内山高坡陡，当使用SBAS技术对其地表进行监测时，形变结果易受地形起伏的影响。为了克服影像因区内地形起伏、地表坡度等引起的几何畸变，需要结合雷达影像成像规律和地表特征对研究区的可视性进行划分，以对基于SBAS技术获取的地表形变监测数据进行合理筛除。SAR影像常常会发生以下3类几何畸变： 阴影、透视收缩和叠掩[7,9-10]。根据Sentinel-1A卫星几何参数以及研究区30 m分辨率所获取的地形参数，可得出本文所用降轨轨道卫星在研究区的可视性分区方式，并划分了相应的可视性区域。

3.2 地表形变监测结果

采用短基线集时序InSAR分析方法，对覆盖小金川河流域的Sentinel数据进行时序干涉处理，获取了2018-11—2019-12期间该地区的地表形变特征信息，即雷达视线方向LOS(line of sight)的年平均形变速率，结合研究区可视性分析结果，剔除了透视收缩、叠掩和阴影不可信区域的形变信息，保留了可视区域和低敏感区域中的形变点。在小金川河主干道两侧及平坦的地区发生的可能性也很低，因此本次剔除了研究区密度较高的林地、灌木林、河流以及坡度小于5°的不可靠区域的形变点，保留其他区域可靠形变点。本次从形变速率数据集中共识别出500万个测量位置点(measurement points，MPS)，平均分布密度约为1 492 MPS/km2。

研究区年平均形变速率信息如图4所示，图中红色负值代表地表目标位移沿远离卫星视线方向，蓝色正值表示位移靠近或朝向卫星视线方向[10]。在1年多时间里，LOS方向的年平均形变速率达到-51.12～75.28 mm/a，分析区内形变异常分布规律，共判译出4处形变异常区，见图4(a)中Q1、Q2、Q3及Q4。Q1表现为形变负异常，位于丹巴县县城东南部，呈葫芦状，起始部位于小金川河与大渡河流域交汇处，主体沿大渡河流域展布，长约7.9 km。Q3位于小金县宅垄镇与丹巴县半扇门镇交界处，长约10.8 km，主要沿小金川河展布，呈带状分布，异常值明显区别于周边区域，分布面积较大，为本文形变异常重点关注区。Q2与Q4分别位于丹巴县半扇门镇与小金县八角镇，异常区面积相对较小，形变异常特征相对于周边区域较为显著。

图4 研究区形变速率及地质灾害隐患探测识别

3.3 地质灾害识别结果分析

参阅Colesanti C及张路等的研究成果[1,10]，当山地斜坡区域LOS方向形变速率绝对值超过8～12 mm/a，即可认定为不稳定状态。在基于SBAS-InSAR处理结果的平均形变速率信息与形变负异常区判译成果基础上，通过形变点速率、聚集特征圈定小金川河流域内地质灾害隐患点，之后借助国产高分等可见光高分辨率光学影像数据，结合地质灾害光学影像纹理特征，对隐患点进行核对筛查，完成潜在地质灾害隐患点识别。

本次共判译识别11处不稳定坡体，疑为地质灾害隐患点，如图4(b)所示，各隐患点编号、地理位置、坡向、隐患区面积、MPS平均点密度及最大LOS形变量等信息见表1。

表1 InSAR分析地质灾害隐患探测识别结果Tab.1 InSAR analysis of geological hazard detection and identification results

(续表)

经与已有地质灾害台账、Google地球、近期高分遥感影像、无人机航片及现场野外反馈照片等资料对比验证，11处隐患点中的6处： P1(阿娘寨滑坡)、P2(擦巴沟滑坡)、P5(大鹰寨滑坡)、P6(下槽滑坡)、P7(旦巴里滑坡)及P11(水井湾滑坡)为已知地质灾害点，其余5处隐患点在本文研究前尚不为人知。判译的各处隐患点主要沿小金川河流域分布(图5)，行政上主要位于小金川县宅垄镇(6处)、新桥乡(1处)、美兴镇(1处)与丹巴县太平桥乡(2处)、半扇门镇(1处)等； 面积大于0.6 km2的有5处； 最大LOS方向形变量介于-18.92～-36.56 mm/a； 平均空间密度为3 204～5 387 MPS/km2。

隐患点P2位于太平桥乡三木扎村，面积最大(1.127 5 km2)，最大形变量达-27.98 mm/a，成灾可能性较高; 区内分布有连片居民聚集区，危害性较大。P5与P3隐患点分别位于宅垄镇千家村与太平桥乡一只碉村区域，LOS方向形变量均较大，分别达约-36.56 mm/a和-34.06 mm/a，发展成灾的可能性高，隐患点面积约0.615 8 km2和0.553 5 km2。P4、P5、P6、P7、P8及P9等隐患点上部有居民聚集区，面积介于0.5 ～ 0.76 km2之间，下方分布有小金川河，有堵江次生堰塞湖灾害的可能性，应作为重点加强监测。P10与P11分别位于小金县城上游的新桥乡团结村与下游的美兴镇卯梁村，面积相对较小，约0.267 1 km2和0.464 6 km2，均在后缘部分形变量最大，分别达-19.15 mm/a和-18.92 mm/a，威胁对象主要为团结村与卯梁村聚集民居、周边乡村道路及耕地等。P10下部为水卡子沟，为小金川河流域次级子流域，一旦成灾，势必形成泥石流灾害物源，如遇强降雨等极端天气，极易堵江断流，威胁距沟口1.5 km的小金县城，建议加强雨季监管值守与预警。P11位于小金县城上游约1.5 km处，P11下部为小金川河美兴段与国道G350春厂段，灾害一旦发生有堵江断路的隐患，造成汶川与宝兴县通往小金县城的生命线中断，次生堰塞湖更是威胁小金县城安全，建议加强汛期监测排查与隐患治理。

分析2018-11-10—2019-12-23地质灾害各隐患点判译结果可知: P2、P3、P4、P5、P6、P7、P8、P9、P10、P11号隐患点持续变形，在1a多的时间里，累积形变在雷达视线方向上分别达到了-31 mm、-29 mm、-36 mm、-44 mm、-24 mm、-27 mm、-40 mm、-33 mm、-26 mm、-39 mm; 隐患点P2在2019-06之前相对稳定，累积形变量约-15 mm，随后开始形变加速； 隐患点P5、P8、P11在2018-10—2019-12间是累积形变量最大的3个区域，其中P5与P8在2019-07呈现变形加速，分别由-14 mm、-16 mm快速跃至2019-12的-44 mm与-40 mm，加速现象显著； P11则在2019-08出现累积形变加速，由-12 mm快速累积至-39 mm，加速迹象明显。

3.4 典型隐患点时序监测与验证

为验证InSAR技术时序监测方法判识地质灾害成果，本节以流域内隐患点P1(阿娘寨滑坡)为例，采用SAR数据集进行滑坡形变的长时序监测研究(图6)，研究灾害体的历史形变特征，印证灾害体的形变加速过程，以评估利用InSAR技术开展地质灾害隐患早期识别的可靠性。

(a) 隐患点P1形变速率 (b) 隐患点P1时间序列形变折线

隐患点P1总体呈圈椅状或舌状，即为已知的阿娘寨(烂水湾)滑坡，位于半扇门镇阿娘寨村，其西北缘正对梅龙沟，下方为小金川河曲折蜿蜒处，下部有国道G350通过，判译隐患区面积约0.651 5 km2，坡向约300°，滑区LOS方向的年平均形变速率为-20.12～7.05 mm/a，最大形变量-20.12 mm/a。如图6所示，滑坡形变迹象信息区域主要分为4个区域(图6(a)中a、b、c、d)，自滑坡中上部向外呈环形扩展，其中a区与b区属形变量较大区域，位于滑坡中后部及中部，呈手柄状，形变异常在隐患点较为明显，形变速率分别达-15 mm/a与-8 mm/a以上，成灾致灾可能性大，危害性较大； c区与d区属滑坡形变一般区，呈环形、葫芦状，位于形变较大区域外围，总体上窄下宽，为滑坡中心形变挤压周边所致，形变速率分别为-8.31～4.84 mm/a与-4.84～-0.84 mm/a，疑为缓慢沉降所致。

在2018-11—2019-12这一年多时间里，滑坡各区域累积形变在LOS方向总体呈沉降迹象，最大值达-27 mm。其中a区在2019-06-02前后有一次明显的沉降迹象，由-12 mm跌至-20 mm，之后略有抬升，到9月前后沉降速率加大，趋势更为显著，在12月达到极值，推测在该时间前后滑体裂缝形成并加剧所致； b区与c区累积形变量变化趋势与a区域相似，均在9月前后表现为沉降加速，最大累积形变值分别为-25 mm与-19 mm； d区相对其他3个区域累积形变相对较小，6月之前形变速率较小、迹象不甚明显，且在局部时间呈现抬升趋势，推测该处为滑坡中心部位物源挤压周边所致，其在9月后开始出现形变加速，至2019-12累积形变量达极值，为-11 mm。

隐患点P1滑坡沿小金川河右岸展布，位于河道拐角处，为土岩混合滑坡群，受上游强降雨诱发泥石流灾害影响，滑坡在2020-06-17加速形变、失稳、下坠，形成多次滑坡造成叠加效果，后在堰塞湖自然溃坝泄漏流水冲击作用下，点状滑坡互相作用诱发滑坡群，形成灾害链效应[2]。结合灾后无人机机载LiDAR监测数据(图7)，滑坡下部形变以增量趋势为主，变形量介于3～6 cm之间，推测原因疑为滑体上部老滑坡滑移物下坠堆叠所致，部分区域形变为微弱递增，介于1～2 cm，推测或为上部滑落物坠落至下部的靠小金川河处，后被堰塞湖下泄流水冲刷带走； 滑坡上部整体呈现下滑趋势，其在2020年6月18—20日的变形量在-5～0 cm之间，形变量较大，该隐患的存在对灾后应急救援抢险构成较大安全隐患，建议在灾后恢复重建过程中，利用InSAR及光学遥感等手段对其加强跟踪监测，保障重建施工安全。

针对上述短基线集InSAR分析判译识别的地质灾害结果，主要通过已有地质灾害点对比、可见光遥感影像解译成果佐证、实地现场调查验证及差分干涉处理结果定性比较等技术方法完成。

图7 典型隐患点(P1)现场照片及形变监测

3.5 灾害识别结果影响分析

对研究区小金川河流域的灾害识别结果，证明了SBAS-InSAR技术在地质灾害早期识别中的优势及有效性，该技术成果在川西高山峡谷区具有大范围推广应用潜力，但其实际应用效果尚受限于多种影响，具体包括如下几部分。

(1)影像几何畸变。由于卫星SAR数据采集系统的固有观测成像模式，导致其对高地形起伏区的影像容易出现几何畸变现象，形成盲区。通常在面向SAR卫星的地表坡体区域出现叠掩或收缩透视的可能性更大，多个雷达散射回波信号更易叠混，较难获取有效地形形变信息； 相比而言，背向卫星的地表坡体出现叠掩现象较少，更适用于获取观测坡体形变信息，稳定性较好。影像几何畸变的大小主要与卫星传感器和观测地表坡体的相对几何关系有关，特别是卫星入射角和坡体的坡向、坡角等关键参数[1]。可综合采用降轨、升轨观测模式以减少盲区，提升地表形变观测效率。既有雷达卫星(民用)观测波段主要集中于C(λ=5.6 cm)、X(λ=3.1 cm)与L(λ=23.5 cm)波段，波长越长代表穿透植被冠层更容易，时间相干性更强，地表形变观测灵敏度越低； 波长越短则地表形变观测灵敏度越高[10]。

(2)对流层大气传播延迟。雷达卫星信号在传播过程中受到大气对流层水分含量时空变化的影响，回波信号易出现传播延迟现象，形成雷达回波差分相位测量误差。当前星载SAR卫星的大气相位传播延迟的去除，已成为利用时序InSAR技术分析地表形变与灾害早期判识的重要瓶颈[1,10]。提取高山峡谷区地表形变信息的关键是对流层大气传播延迟的定量估算与去除。针对本文研究的小金川河流域等典型川西高山峡谷区，与地表高程起伏相关的对流层垂直分层延迟相位占据主导地位，且具有季节变化特征与时间相关性，使用传统的时空滤波等方法在估算与去除时难度较大。

(3)可探测最大形变梯度。时序InSAR技术通常基于雷达波相位观测，在观测地表形变的空间梯度与时间形变速率时面临上限阈值，即空间上临近观测点的地表形变差异值需小于等于波长的四分之一(λ/4)，才能保证有效的相位解缠成果； 而时间形变速率取决于雷达数据获取时间间隔(Δt)与空间形变梯度[1,10-11]。波长较长的雷达波与高空间分辨率可提升地表形变空间梯度与较大时间形变速率，较短时间获取周期也是提高形变时间形变速率的有效方法。但若地表形变量过大则较难获取到相干散射体，容易导致失相干现象的出现。

4 结论与讨论

(1)采用SBAS-InSAR方法完成了小金川河流域的数据分析，梳理总结了数据处理的基本流程，包括数据(SAR与DEM)预处理、差分干涉计算、时间空间域形变估算及形变量计算等步骤。

(2)采用SBAS-InSAR分析方法，对覆盖小金川河流域区域进行时序干涉处理，获取了2018-11—2019-12的年平均形变速率(LOS方向)为-51.12～75.28 mm/a。分析了区内形变异常的分布规律，共判译出4处形变异常区，其中： Q1形变异常区位于丹巴县县城东南部，沿大渡河流域展布，长约7.9 km； Q3形变异常区长约10.8 km，主要沿小金川河展布，呈带状分布，异常值明显区别于周边区域，分布面积较大； Q2与Q4形变异常区分别位于丹巴县半扇门镇与小金县八角镇，异常区面积相对较小，异常特征显著。

(3)基于平均形变速率信息与形变负异常区判译结果，结合地质灾害光学影像纹理特征，共判译识别11处不稳定坡体，疑为潜在地质灾害隐患点，其中6处隐患点为已知地质灾害点，其余5处隐患点在本文研究前尚不为人知。这些隐患点主要位于小金川县宅垄镇(6处)、新桥乡(1处)、美兴镇(1处)与丹巴县太平桥乡(2处)、半扇门镇(1处)。

(4)选择隐患点P1(阿娘寨滑坡)开展长时序监测研究，并通过灾后LiDAR监测数据进行验证，结果表明利用InSAR技术开展川西高山峡谷区地质灾害隐患早期识别具有较好的可靠性。监测结果与实际情况存有一定误差，分析其可能受几何畸变、对流层大气传播延迟及可探测最大形变梯度等多种因素影响，这将是下一步研究改进的方向。