西藏然乌地区地质灾害隐患点InSAR识别与监测

张本浩，魏云杰，杨成生，熊国华，董继红

（1. 长安大学地质工程与测绘学院, 陕西 西安 710054；2. 中国地质环境监测院（自然资源部地质灾害技术指导中心）, 北京 100081）

0 引言

我国地质构造运动频繁，地形结构复杂，地质灾害隐患多、分布广、高易发[1]。地质灾害往往具有强大的破坏力，严重威胁群众的生命和财产安全、阻碍当地经济发展[2]。如何准确地识别灾害隐患所在位置、可能发生的时间并提前预警，是当前地质灾害工作人员关注的首要问题[3]。因此，开展地质灾害早期识别与监测，确定灾害体空间分布及活动规律，对防灾减灾工作具有重要意义。

嘉黎断裂地处青藏高原东南缘，地势高陡，经过多次构造变形运动，形成了大量节理、断层等[4−5]。唐得胜等[6]通过对然乌湖到通麦镇典型冰湖溃决泥石流灾害链的成因分析，得出该区域为内外应力耦合作用形成地质灾害链的结论。郭佳宁等[7]通过总结并分析前人对藏东南地区滑坡地质灾害的研究得出林芝—然乌段共存在8个大型滑坡、6个大型崩塌，其中通麦—波密段有102个滑坡群[7]。赵聪等[8]利用遥感数据解译发现帕隆藏布江流域（波密—然乌段）冰川退缩，大量冻融松散物源进入沟道，导致流域内冰川泥石流大规模发育[8]。这类地质灾害一般规模较大、易堵河，对当地群众的生命财产安全构成严重威胁。如2017年8月3日，波密县倾多镇暴发泥石流；2018年10月，林芝东部色东普沟暴发大规模冰川泥石流；2020年6月，嘉黎县金翁错暴发冰碛湖溃决事件。目前对然乌地区的监测集中在波密县到然乌镇的冰川和泥石流地质灾害发育情况和构造部位的运动特征[4,9]，观测方法主要是光学遥感数据分析，大部分学者关注重点部位地质灾害发育情况[10−12]，对区域地质灾害形变特征和规律的分析研究较少。

然乌镇地区地质灾害易发，地形陡峭复杂、落差大、植被茂密，传统地面监测手段难以开展。合成孔径雷达干涉测量（Interferometric Synthetic Aperture Radar，InSAR）作为一种大范围地表形变监测技术，其穿透能力强、覆盖范围广、能快速获取区域形变信息，已经被广泛用在冰川[13]、滑坡[14]、泥石流[15]等的调查和监测，为地质灾害的防治工作提供了新的技术手段。因此文章在分析InSAR技术在该地区适应性的基础上，利用InSAR技术对然乌地区开展地质灾害隐患早期识别和形变监测研究，为该区域防灾减灾提供参考。

1 研究区概况及数据来源

1.1 研究区概况

研究区位于藏东南然乌地区（图1），属嘉黎断裂东部，川藏铁路康定—林芝段的必经路线。研究区内相对高差最大超过2 000 m。受季风环流影响，夏季湿凉、冬季干冷漫长多大风。由于该区域属于典型的高山峡谷地貌，且受嘉黎断裂带活动和气候的影响，地面切割强烈，地质灾害易发。研究区内灾害主要有泥石流、滑坡等。2016年田小平[16]通过野外调查结合文献资料，统计得到然乌—波密段泥石流共有118条，1994—2015年国道318线波密县至然乌镇段降雨泥石流和冰川泥石流总计64处。泥石流会引发堰塞湖、阻断交通、毁坏基础设施等[17]问题。

图1 研究区影像覆盖示意图Fig.1 Image coverage diagram of the study area

1.2 数据来源

为了利用InSAR技术对然乌地区开展地质灾害隐患识别和监测，文中收集了覆盖研究区从2018年1月至2020年11月的Sentinel-1A升轨SAR影像82景，从2018年1月至2020年10月的Sentinel-1A降轨SAR影像75景，表1列出使用的SAR数据基本参数。使用日本宇宙航空研究开发机构（JAXA）发布的AW3D30（https://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm） 全球数字地表模型数据消除InSAR干涉图中的地形相位[14]。同时为了分析SAR成像几何对地表观测的影响，AW3D30被用来计算SAR影像的叠掩和阴影区。为了改正基线误差，使用了Sentinel-1A的POD精密轨道卫星星历（Precise Orbit Ephemerides，POD）数据，数据下载网址（https://s1qc.asf.alaska.edu/aux_poeorb/），定位精度优于5 cm[17]。

表1 SAR数据基本参数Table 1 Basic parameters of SAR data

2 技术方法

2.1 基于InSAR堆叠技术的滑坡探测

本次实验使用干涉图堆叠技术（Stacking-InSAR）进行大范围地质灾害隐患调查。Stacking-InSAR技术由SANDWELL等[18]在1998年提出，通过对多幅解缠干涉图进行加权平均。在对解缠干涉图堆叠处理之前，需要将参与计算的干涉对重新采样到相同投影坐标系下。与传统的差分技术（D-InSAR）相比[19]，该技术能够明显的减弱数字高程模型误差和大气误差的影响[14]。因此，Stacking-InSAR技术用于然乌地区地质灾害隐患点识别具有较好的技术优越性。

Stacking技术假设任意一幅干涉图是独立、大气误差是随机且相等的、形变速率是线性的。在此基础上，将多幅解缠相位图按干涉对时间间隔进行加权平均，就可以得到研究区的形变相位信息。根据误差传播定律，加权平均后大气相位误差与解缠相位图数量之间并不是倍数关系，而是平方根关系。从而增加了形变和大气相位的信噪比，提高了形变监测精度，其数学模型可以表示为：

式中：V——年均形变速率；

M——参与加权平均计算的干涉图数量；

ΔTi——第i幅干涉图的时间基线；

φi——第i个干涉图的解缠相位值。

由于原始Sentinel-1A影像覆盖范围较大（250 km×190 km），为提高数据处理效率，仅选择了与研究区相关的区域进行处理。数据处理过程中，时空基线阈值分别为60 d和±200 m，升降轨影像分别产生368和317个干涉对，通过挑选相干性高且无明显解缠误差干涉对，最终分别选择205个和189个高质量干涉对进行Stacking技术处理。

为了抑制数据处理中噪声影响，提高干涉图的相干性，数据处理过程中对SAR影像进行了多视处理，多视比为4∶1，多视处理后SAR影像分辨率约为15 m，并采用16×16和8×8窗口进行两次自适应滤波[20]。对残余轨道误差，采用二次多项式拟合方法进行去除。我们将InSAR识别出的疑似隐患点，叠加到谷歌光学影像中，通过与实际地形条件对比分析，圈定最终的地质灾害隐患点，具体实验流程如图2所示。

图2 实验流程示意图Fig.2 Schematic diagram of experimental flow

2.2 二维形变解算

由于SAR系统侧视成像特点，导致单一轨道的SAR数据在高山峡谷区进行地质隐患点探测中会存在漏判情况[21]，而采用升轨和降轨数据组合的观测方法，则可以增加对地观测的可见范围，有利于减少隐患点的漏判，且通过不同方向的形变特征能够更好地反映隐患点的形变特征。文中采用多维小基线子集(Multidimensional Small Baseline Subset，MSBAS)技术，联合升、降轨数据对研究实验区开展地质灾害隐患点的二维形变监测。MSBAS技术是由SAMSONOV等[22]在2013年提出，和传统小基线技术（Small BAseline Subset，SBAS）有所区别，该技术可以计算隐患点的二维(水平东西向、垂直向)形变速率和时间序列，称为MSBAS-2D技术[23]。

MSBAS-2D技术利用两个不同轨道在相同时间段覆盖相同区域的影像，常用升轨影像和降轨影像，获取该区域的二维形变时间序列。其数学模型可以表示为：

θ——方位角；

φ——航向角；

λ——正则化参数；

VE、VU——表示水平东西向速率和垂直向速率；

当升轨和降轨数据获取时间存在差异时，利用吉洪诺夫正则化矩阵B和 参数 λ相乘，得到正则化的时间序列结果，类似于利用低通滤波的方法去除高频信号的噪声。再应用奇异值分解（SVD）方法获得二维形变速率，在形变速率的基础上，使用数值积分法获得形变时间序列[23]。

3 一维形变监测结果与分析

3.1 研究区几何畸变分析

由于SAR系统的侧视成像模式，在较小范围内，可将微波看作为平面波且入射角视为固定。若地面不平坦则表示地表存在起伏，导致几何畸变现象的产生（顶底倒置、阴影、透视收缩）（图3）[24]。当隐患点位于叠掩区时（图3b），会造成滑坡体形变特征被压缩且点数量较少；当隐患点位于阴影区域时（图3c），雷达卫星不能接收到该区域的任何后向散射信息，则导致隐患点的漏判。因此，同一区域采用不同视角的卫星（如升降轨联合观测），则可以提高隐患点的识别范围。我们根据升降轨数据的成像参数和地形数据计算了研究区的叠掩阴影分布如（图4）。从图4中可以看出，升轨受到叠掩和阴影几何畸变影响，且在整个影像上均有分布，其面积占比分别为22.72%和0.86%；降轨受到叠掩和阴影几何畸变影响，且空间分布集中在西部，面积占比分别为7.45%和3.77%。升轨和降轨数据几何畸变空间分布和占比面积均存在差异，其原因主要是卫星的航向角和入射角的不同，导致几何畸变的空间分布和占比面积存在差异。对位于几何畸变区域中的隐患点，将会在InSAR识别结果中进行综合判断以减少误判。

图3 SAR成像几何畸变示意图（a:顶底倒置；b:透视收缩；c:阴影）Fig.3 Geometric distortion diagram of SAR imaging (a: top bottom inversion; b: fluoroscopic contraction; c: shadow)

图4 然乌地区升降轨叠掩阴影分布图Fig.4 Overlay shadow distribution of lifting rails in Ranwu area

3.2 大范围识别结果

相比较D-InSAR技术，Stacking技术能够明显的削弱数字高程模型误差和大气误差的影响[25]，研究中利用Stacking技术获取了然乌地区沿卫星视线方向年均形变速率结果。结合然乌地区地形和地质灾害的活动特性，选择然乌地区人工建筑物作为稳定区，计算标准差为1.96 mm/a，最大形变速率为8.3 mm/a。因此选择形变量级±10 mm/a的区域作为稳定区，当形变速率超过这个阈值时，则将其判定为地质灾害疑似隐患点。然后，结合光学影像和地形条件，对识别出的地质灾害疑似隐患点进行对比分析，最终共圈定潜在隐患点67处，隐患点分布如图5所示。升轨和降轨影像分别识别出隐患点41处和45处，其中共同识别19处。灾害类型主要为滑坡和冰川，部分典型隐患点InSAR结果如图6所示。从识别结果的空间分布看：识别隐患点主要集中在318国道附近，且在然乌附近分布较多的隐患点；识别隐患点在空间分布上有所差异，升轨影像主要集中在然乌到吉达乡，降轨影像隐患点主要集中在然乌到玉普乡。分析原因是由于升降轨影像入射角差异和叠掩阴影空间位置分布不同导致的。

图5 然乌地区地质灾害隐患点InSAR判识结果分布图Fig.5 Distribution of InSAR identification results of geological hazard hidden spots in Ranwu area

图6 典型隐患点InSAR结果及遥感影像（图a和图b分别对应图5中1号和2号点）Fig.6 InSAR results and remote sensing images of typical hidden trouble points (Fig.a and b correspond to no.1 and no.2 in Fig.5 respectively)

3.3 典型隐患点时空形变特征

为了更好地分析识别的隐患点在时间和空间上的形变特征，选择形变特征明显、对附近居民点或交通路线存在隐患的两个隐患点，分别是日曲滑坡（1号点）、莫格曲冰川（2号点）（图5），提取得到了这两个隐患点的时间序列结果。日曲滑坡位于古绕村日曲北岸陡峭山体上，坡体西侧为318国道。根据该滑坡体在2018年1月至2020年11月的时间形变序列和整体的年均形变速率结果（图7），可以看到该滑坡体处于等速蠕滑阶段，结合其剖线图（图8）发现，在距离山顶0.56 km处滑坡的形变速率达到最大值−94.8 mm/a，随着坡度增大，形变速率增大且两年多时间累计形变量达到278 mm。莫格曲冰川位于莫格曲东岸，通过光学影像发现冰川上部存在两个冰湖，根据2018年1月至2020年11月时间形变序列和整体的年均形变速率结果，可以看出相比日曲滑坡，莫格曲冰川体的形变速率更大，结合其剖线图（图7）发现该冰川形变主要集中于冰川中部，最大年均形变速率处距离山顶0.64 km，且在雷达视线方向达到107 mm/a，两年多时间累计形变量达472 mm。根据这两个隐患点的剖线结果，可以看出变形趋势均是先增大后减小，呈漏斗状，最后形变速率均趋近稳定。

图7 典型隐患点时间序列Fig.7 Time series of typical hidden trouble points

图8 典型隐患点剖线图Fig.8 Section diagram of typical hidden danger points (Riqu landslide on the left and Mogequ Glacier on the right)

4 二维形变监测结果与分析

4.1 二维时序形变监测结果

从图5结果可以发现，然乌镇附近安目错北岸的升轨、降轨数据均存在明显的冰川群形变信息，其中包括吉穷隆冰川、迫弄冰川、瓦巴村冰川（图5中3、4、5号点）等多处典型隐患点。因此对该区域采用MSBAS-2D技术处理，反演该区域隐患点的二维形变结果(图9），其中水平东西向正值和负值分别代表向东和向西运动；垂直方向正值和负值分别代表向上和向下运动。由于该区域没有实地监测数据，为了评价MSBAS-2D结果的精度，选取升轨、降轨数据的解缠参考区域，通过计算得到标准差在2.2 mm/a内，小于隐患点判识的阈值（±10 mm/a），证明了反演结果的可靠性。从图9(a)可以发现，三个冰川(吉穷隆冰川、迫弄冰川、瓦巴村冰川)在水平东西向主要表现为正值，最大能达到90 mm/a，表明这些冰川在水平方向均向东运动，结合区域实际地形可以发现这三个冰川的朝向均为偏向东，即沿坡体方向向下运动；从图9(b)可以发现，这三个冰川在垂直向均表现为负值，最大超过−90 mm/a，即沿近似坡向的向下的主滑方向运动。

针对然乌区域3个典型冰川，我们选择了5个特征点(图9中A、B、C、D、E)，并提取了他们的形变时间序列结果(图10)。结果显示：吉穷隆冰川在水平方向为正值，近3年时间内水平东西方向最大累计形变量达到了202 mm，且在该冰川前缘垂直方向同样表现为加速变形趋势；迫弄冰川是3个典型冰川中面积最大的冰川，通过光学影像计算该冰川面积约为2.26×106m2，冰川形态为狭长条带状，两年多的时间内水平东西方向最大累计形变量为283 mm，且该冰川前缘同样表现为水平方向加速位移趋势；瓦巴村冰川运动方向近似平行东西方向，从形变时间序列结果上发现该冰川在东西和垂直方向变形状态仍然表现为加速，其中水平和垂直方向最大累计形变量均超过150 mm，且都呈现加速趋势。因此，建议后续加强对这三个典型冰川的形变监测，关注该冰川后期形变趋势。

4.2 形变外界驱动因素分析

瓦巴村冰川位于瓦巴村北方的山坡上，面积约63×104m2，在升降轨年均形变速率结果中均发现该滑坡形变特征明显。文中选取瓦巴村冰川来分析其形变特征，为此，在该坡体上选取了一处特征点（图9中E点），获得了其二维时间序列结果（图10）。从结果上可以发现，无论水平方向还是垂直方向，其累积形变量均超过150 mm，且2020年6月以后表现为加速趋势。

图9 重点地区冰川群二维年均形变速率图Fig.9 Time series of horizontal east-west and vertical deformation in key areas

图10 重点区域水平东西向和垂直向形变时间序列Fig.10 Time series of horizontal ew and vertical deformations in key areas

为了探讨降雨和温度对瓦巴村冰川形变速率的影响，收集了与InSAR覆盖时间段一致的降雨量和温度数据。经对比分析，降雨量和温度变化对冰川形变存在一定的影响(图11)。在每年雨季期间，冰川的形变量会产生加速趋势，其他月份冰川较稳定，且随着温度的升高，冰川变形也表现出部分加速趋势。

图11 瓦巴村冰川二维时间序列月降雨量和温度关系Fig.11 Relationship between monthly rainfall and temperature in twodimensional time series of Waba glacier

5 结论

文章通过采用Stacking技术和MSBAS技术对然乌镇附近区域地质灾害进行识别和形变特征的研究。

（1）在影像覆盖范围内共识别出67处灾害隐患点，其中19处疑似隐患点为升降轨影像共同识别。识别出的滑坡体主要集中于然乌镇318国道和201省道两侧，且最大滑坡速率可达−168 mm/a，威胁两侧道路的安全。利用MSBAS技术获取了然乌镇附近安目错北岸冰川群二维形变时间序列，分析表明冰川位移沿坡体方向向下运动，垂向位移速率超过90 mm/a，由于严重的位移失相干，冰川体物源区的位移未能获取。选取其中的瓦巴村冰川进行分析表明，其中水平和垂直方向最大累计形变量均超过150 mm，且都呈现加速趋势。需加强对该冰川后续形变监测，同时研究表明降雨和温度对冰川位移速率有一定的影响。

（2）文中仅利用了InSAR方法对然乌区域地质灾害隐患点识别与监测技术方法的适用性进行了研究分析，缺乏野外实地调查资料。通过野外调查对InSAR结果进行验证，并采取不同角度分析然乌区域地质灾害隐患点的形变特征是下一步将要做的事。