基于相干性的InSAR时间序列方法追溯四川雅安地区汉源滑坡灾前形变

徐晓雪,季灵运,3*,张文婷,刘传金

(1. 中国地震局第二监测中心，陕西 西安 710054; 2. 国家遥感中心 地质灾害部，北京 100036;3. 防灾科技学院 地球科学学院，河北 三河 065201)

0 引 言

合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)技术由于全天时、全天候、覆盖范围广、监测精度高和空间分辨率高的优点，在山区监测活跃滑坡方面显示出巨大的潜力。InSAR技术可以提供滑坡位置、边界和运动特征的关键信息，为滑坡的灾前时空形变追溯提供了新的技术途径。自Achache等首次使用差分合成孔径雷达干涉测量(Different InSAR,D-InSAR)技术对法国阿尔卑斯地区La Clapiere 滑坡进行监测以来，为了提高形变测量的精度和可靠性，在随后的研究中诸如永久散射体InSAR(PSI)方法、InSAR小基线子集(Small Baseline Subset,SBAS)方法和混合时序InSAR方法等多种时间序列InSAR方法被提出来。其中，SBAS方法具有获取微小形变信息和长时间序列缓慢地表变形的优势，能够在一定程度上克服时间失相干和空间失相干的不利影响，在复杂山区滑坡形变监测方面具有较为广阔的应用前景。很多学者采用SBAS方法进行滑坡的早期识别与形变监测，如Manunta等使用SBAS方法对意大利中部翁布里亚地区进行大范围滑坡识别与监测；Zhao等使用SBAS方法对美国西部森林覆盖区域的滑坡进行调查，获得了滑坡的时间序列形变特征；Liu等采用SBAS方法对美国加利福尼亚州三只熊滑坡的滑前形变进行监测，清晰地厘定了森林山区缓慢移动滑坡的运动特征。随着SAR卫星的相继发射和大量SAR数据的存储，基于雷达遥感数据的SBAS方法在未来将会发挥更加重要的作用。

滑坡是中国西南山区最主要的地质灾害之一，由于地形复杂多变且滑坡破坏区域较大，往往难以防范，滑坡的高精度监测对于了解其运动学特征和减少滑坡灾害具有重要意义。2020年8月21日，四川省雅安市汉源县富泉镇中海村6组发生山体滑坡(简称“汉源滑坡”)，滑坡体总方量达到500×10m，属特大型滑坡。此次滑坡造成群众房屋受损及部分人员伤亡，省道435交通中断，给人民生命财产安全带来威胁。追溯滑坡体灾前形变状态，对后续开展针对性的工程地质调查及滑坡诱因分析具有重要作用。因此，本文基于升轨、降轨哨兵一号(Sentinel-1)卫星影像，求解滑坡体的InSAR时间序列，研究四川雅安地区汉源滑坡的灾前演变情况。由于汉源滑坡发生在大渡河流域高山河谷地区，坡体表面作物覆盖，极大地影响了InSAR数据处理中干涉图的相干性，所以采用基于相干性的小基线集时间序列算法，以适应低相干分析，更好地获得形变场，从而厘定其时空变化规律，为灾后稳定性评估和未来滑坡灾害的监测预警提供思路。

1 研究区概况

汉源县位于四川省西南部，隶属雅安市管辖(图1中矩形框代表Sentinel-1影像覆盖范围)，为四川盆地与青藏高原之间的过渡地带，省道306经县境西北顺着大渡河穿过，与成昆线相接。从地形来看，汉源县位于大渡河流域高山河谷地区，河道下切，山高坡陡，是滑坡易发区。此次滑坡所在的区域属典型中低山斜坡地貌，为台阶状地形，斜坡前缘较陡。坡面上既有陡坡又有缓坡，滑坡体整体表面坡度约25°，最陡地方的坡度达40°，为利于滑坡失稳的地形地貌条件。滑坡体整体呈长条形舌状，宽约280 m，斜长约880 m。

红色五角星及红色圆圈为滑坡位置；蓝色方框为SAR数据覆盖范围图1 四川雅安地区汉源滑坡位置和SAR数据覆盖Fig.1 Location of Hanyuan Landslide in Ya’an Area of Sichuan and Coverage of the SAR Datasets

汉源县区域岩土体可分为松散岩、碎屑岩、碳酸盐岩、变质岩和岩浆岩5种岩类，易滑岩组的存在在降雨诱发下易于发生区域滑坡事件。汉源滑坡发生后，多部门开展了现场成因调查。从地质结构来看，滑坡体表层的土壤为耕作土和砂性土，结构十分松散；土层下面的岩石是风化很强烈的砂泥岩，容易失稳变形。汉源县由于所处的独特地理位置，具有年降雨丰沛、降雨集中和连续多日强降雨的特点。从2020年8月10日起，汉源县经历了两轮强降雨天气，岩土体含水率接近饱和，斜坡整体稳定性降低，容易发生地质灾害。降雨通过表层土壤下渗到岩石，慢慢浸润岩土体，超过稳定临界值之后就发生了滑坡。强降雨可能是此次滑坡发生的诱发因素。

2 数据来源与分析方法

2.1 研究区SAR数据与干涉处理

本文所采用的SAR数据为Sentinel-1卫星影像。Sentinel-1卫星影像不仅可以实时免费下载，而且覆盖范围广，重访周期短，轨道精度高，可以进行全天候、全天时的高分辨率监测。由于SAR采用侧视成像方式，而且西南山区地形复杂多变，致使覆盖汉源滑坡区域的降轨Sentinel-1 SAR影像在滑坡区域形成阴影，雷达接收不到后向散射信息，严重影响滑坡监测的准确度。因此，本文只针对覆盖汉源滑坡区域的升轨Sentinal-1 SAR影像进行讨论，影像覆盖范围如图1中矩形框所示。本研究所采用的数据分别为59景、56景两个升轨的Sentinel-1影像数据，不同轨道SAR数据的时间跨度统一为滑坡发生之前的2018年10月至2020年8月，时间间隔为12 d。使用美国航空航天局(NASA)发布的30 m空间分辨率SRTM1 DEM地形数据模拟并消除地形对干涉相位的贡献，与SRTM3数据相比具有更高的空间分辨率。

数据处理使用开源InSAR处理软件GMTSAR进行。由于此次滑坡的空间范围较小，在雷达坐标距离向和方位向不进行降采样，多视比设置为1∶1；经过地理编码后，单个像素的空间分辨率约为25 m。采用较小的多视比得到的差分干涉图相干性较低，存在较多斑点噪声，因此，对原始差分干涉图进行滤波，以减少噪声，提高相干性。对于轨道误差改正，通过建立二次多项式模型来消除。由于滑坡位于复杂地貌山区，茂密的植被使SAR影像难以保持较好的相干性，也很难在植被茂密的滑坡上提供有用信号，干涉对在很大程度上受到时间失相干的影响，特别是临滑前夏季影像生成的干涉对失相干现象明显，导致了长时间跨度干涉图的时间失相干。在去除低相干干涉对后，最终选择P128轨道的54个干涉对和P26轨道的41个干涉对进行后续InSAR时间序列分析。这些干涉对的基线通常小于100 m，SAR数据的具体获取日期与时空基线关系如图2所示。

20181216是指2018年12月16日，其他依此类推图2 SAR数据时空基线图Fig.2 Spatial and Temporal Baselines of SAR Datasets

2.2 InSAR时间序列分析

SBAS方法采用较小的时空基线，在提取滑坡沿雷达视线方向(Line of Sight,LOS)的时间序列形变方面较为常用。SBAS方法的主要思路是基于一定的时空基线阈值，选取多个主影像组成自由组合的干涉对，若干涉构网形成一个唯一子集并充分连接，可采用最小二乘法求解时间序列信息。当存在多个子集时，依据各相干像元相位与观测时间的关系，利用奇异值分解(Singular Value Decomposition,SVD)方法将多个小基线集合数据联合起来求解，有效地解决了各数据集之间的时间不连续问题，从而求得时间序列形变、平均形变速率结果。然而，由于研究区植被覆盖茂密，传统SBAS方法容易受到失相干的影响，低相干像元通常使用相干性阈值进行掩膜，降低了时间序列结果的点位密度。

本文采用基于相干性的SBAS方法进行InSAR时间序列分析，将干涉图的相位相干性引入到SBAS时间序列处理算法中，不丢弃干涉图中存在的部分噪声数据，尽可能多地保留干涉图中的像素，通过相干性对观测相位进行加权以减弱失相干的影响，提高了时间序列产品的空间分辨率。在时间序列处理过程中，通过SNAPHU软件进行相位解缠，选择较低的相干性阈值0.16，以获得尽可能多的相位信息，这可能会带来相位解缠误差，但是在时间序列分析中失相干像元点会被降权处理。这种基于相干性的SBAS方法使用加权矩阵根据每个差分干涉图中每个像素的相干性对观测相位进行加权，对差分干涉图中的噪声不太敏感，降低了失相干相位的权重，从而可以得到相干信号的密集空间覆盖。每个像素的加权最小二乘反演可以表示为

(1)

=diag{,,,,}

(2)

2.3 坡向形变提取

采用InSAR技术只能得到三维地表变形沿雷达视线方向的投影，然而大多数滑坡运动通常发生在平行斜坡方向。在一定的假设条件下，雷达视线方向变形可以转移到由特定的坡度角和方位角所确定的斜坡方向。本文将雷达视线方向测量结果投影到滑坡的斜坡方向，获得滑坡沿坡向的运动情况，平行斜坡运动可以更好地识别滑动体，提高对滑坡实际运动的量化分析。

经过投影转换后，坡向矢量大小总是不小于雷达视线方向。从雷达视线方向到坡向的转换可以很容易地通过放大比例因子来实现。放大比例因子()表达式为

(3)

式中：=[sinsin-sincoscos]，为雷达视线方向的形变单位向量；=[-coscos-cos·sinsin]，为坡向的单位向量；和分别表示雷达飞行方位角和入射角；和分别为坡角和坡向。

3 灾前形变结果分析

3.1 雷达视线方向形变速率与时间序列

采用基于相干性的InSAR时间序列方法成功监测到四川雅安地区汉源滑坡的灾前运动特征，各独立轨道的雷达视线方向平均形变速率结果如图3所示。其中，正值表示滑坡靠近卫星运动，负值表示滑坡远离卫星运动。升轨P128和P26轨道在2018～2020年变形分布相似，但变形幅度略有不同(图3)。这主要是因为滑坡位于P128轨道的近距离范围，而位于P26轨道的远距离范围，相对于同一目标而言，其雷达视线方向存在9.9°的偏差。两个升轨的平均形变速率结果显示，滑坡存在明显的滑前形变信号。在滑坡发生前的近两年时间，P128轨道的最大平均形变速率可达到-42 mm·年，P26轨道达到-42.1 mm·年，非形变区大多数监测点的平均形变速率位于-10～10 mm·年，两个升轨的形变值为负值表明滑坡远离卫星运动。从两个轨道的平均形变速率结果可以看出，位移变化较大处主要发生在河流沿岸，这与现场识别的活动滑坡位置一致。

为了更好地分析此次滑坡的滑前时空运动特征，对2018～2020年两个升轨的形变时间序列进行研究。图4为P128和P26两个轨道在雷达视线方向上的时间序列变形。从图4可以看出，不同轨道滑坡运动的时空变化规律基本一致，滑坡在整个监测期间存在持续变形。总的来说，滑坡在2018年底至2019年4月移动迅速，而在之后移动缓慢，基本处于稳定状态。滑坡运动在空间上并不均匀，随着时间的积累，滑坡形变区域逐渐形成形变中心，滑坡中心的累积形变可达到97 mm，而位于周边区域的累积形变相对较小，最大累积形变仅为70 mm。

图3 雷达视线方向平均形变速率结果Fig.3 Images of Average LOS Deformation Rate

两个轨道的时间序列形变结果表明，汉源滑坡的灾前形变过程经历了形变快速积累和稳定两个阶段。由于缺少可用的夏季干涉对，时间序列变形分析未获取2019年7月至9月及临滑前的时间序列信息，灾前形变分析并没有捕捉到临滑前的形变加速信号。

20181221是指2018年12月21日，其他依此类推图4 雷达视线方向累积形变结果Fig.4 Images of LOS Cumulation Deformation

图5 雷达视线方向单点累积形变统计结果Fig.5 Statistical Results of LOS Cumulation Deformation for Single Point

本文提取了汉源滑坡形变中心区域P点[图3(a)]的雷达视线方向累积形变时间序列，结果如图5所示。单点形变时间序列结果显示，两个升轨表现出明显的滑前形变信号。滑坡运动经历了两个时段的形变：Ⅰ阶段为2018年10月至2019年4月，形变加速积累，其中P128轨道累积形变量达97 mm，P26轨道累积形变量达86 mm；Ⅱ阶段为2019年4月以后，滑坡运动速率几乎为0，并在2020年8月21日发生滑坡。由于失相干现象，缺少夏季的干涉对，在2019年7月至2019年9月和临滑前(即2020年6月至8月)并没有获得有效形变信号。2019年4月以后滑坡运动缓慢的具体原因还需结合后续工程地质调查展开。

3.2 坡向平均形变提取与结果验证

为了分析汉源滑坡滑前沿坡向的运动特征，首先通过研究区DEM数据提取滑坡的坡度和坡向。虽然不同轨道数据集测得的雷达视线方向的形变量级不同，但是平行坡向的分量是相同的。因此，在提取坡度和坡向之后，结合式(3)分别计算了P128轨道和P26轨道的放大比例因子。为了控制坡向的偏差，剔除了绝对放大值大于一定阈值(本文取值为2)的孤立点，便于更好地分析滑坡沿坡向的运动规律。

图6展示了P128轨道和P26轨道经放大比例因子逐像素放大校正后得到的平行坡向的平均形变速率结果。从图6可以看出，滑坡发生前在斜坡方向上存在明显的灾前形变信号，滑坡方向的最大平均形变速率可以达到43.9 mm·年，显著大于雷达视线方向。滑坡体中心为滑前形变最剧烈的区域，整体呈现出明显的漏斗状特征。滑坡沿坡向的形变结果显示，形变中心区域长约250 m，宽约250 m，形变中心海拔高度1 017 m，该区域是在滑坡发生前发生显著形变的范围。

由于没有获取研究区地面监测数据，本文使用相同时间跨度的独立SAR数据集来验证InSAR监测结果的精度。通过统计2个轨道独立观测的SAR数据结果重叠点的差值，验证InSAR结果之间的一致性。统计结果显示，2组相邻轨道滑坡区域重叠点的坡向速率差值符合正态分布，速率差值的标准差为2.9 mm·年，表明InSAR监测结果是可靠的。

综合此次滑坡发生前InSAR时间序列结果，汉源滑坡存在明显的灾前形变积累。在滑坡发生地为高山河谷的地形和砂泥岩地质结构背景下，经过强降雨的加载作用，可能打破边坡的稳定阶段。足够的降水浸透滑坡体基底面，降低有效抗剪强度；当边坡剪切应力超过剪切强度时，就可能触发滑坡。结合滑坡发生前该区域经历的强降雨过程，推测在灾前形变积累情况下，强降雨可能是此次滑坡发生的诱发因素。

4 结 语

本文利用多轨道Sentinel-1卫星影像，采用基于相干性的InSAR时间序列方法，对四川雅安地区汉源滑坡的灾前形变情况(2018年10月至2020年8月)进行追溯。首先，获取了滑坡沿雷达视线方向的平均速率与时间序列，结果显示滑坡雷达视线方向的最大平均形变速率为42.1 mm·年；汉源滑坡从2018年底便发生了变形，至滑坡发生前最大累积形变量达97 mm，共经历了形变加速和稳定两个阶段，形变在空间分布不均匀。然后，提取了滑坡沿坡向的形变，结果显示在坡向上同样存在明显形变信号，最大坡向平均形变速率达43.9 mm·年。本文研究结果为此次汉源滑坡提供了第一手、相对完整的形变测量资料，所揭示的滑坡时空运动特征对于评估滑坡危害及进一步实地调查至关重要。上述基于InSAR技术的灾前形变追溯方法也可以用来识别其他滑坡可能发生点的前兆运动。

图6 汉源滑坡平行坡向的平均形变速率结果Fig.6 Images of Average Deformation Rate of Hanyuan Landslide Parallel to Aspect

值得注意的是，由于研究使用的Sentinel-1 SAR影像数据植被穿透能力较弱，在研究区植被茂密的情况下相干性会受影响。本文采用的基于相干性的InSAR时间序列方法虽然在一定程度上可以提高形变结果的点位密度，但由于夏季植被非常茂密，临滑前干涉对相干性很差，所以本文形变追溯结果并未捕捉到临滑前明显的加速变形信号。这也说明InSAR技术在滑坡监测预警中容易受到观测环境影响，可能存在时间采样率不足等缺点。

欧洲航天局提供了Sentinel-1卫星影像，在此表示感谢!