联合升降轨PS-InSAR技术的滑坡二维形变监测
武志刚1，陈恒祎2，赵超英2，马正龙1，田 野2
（1. 国家电投集团青海黄河电力技术有限责任公司，西宁 810000；
2. 长安大学 地质工程与测绘学院，西安 710054）
摘要:
针对InSAR库岸滑坡监测易受时空去相干、大气误差以及单一成像视角等因素的影响，难以满足滑坡形变规律与失稳特征的分析需求，甚至造成滑坡变形的错判或误判等问题。该文联合相干点目标分析法与升降轨Sentinel-1影像开展了果卜滑坡高精度二维时序形变监测，并结合Sgolay滤波器对滑坡失稳特征进行了分析。结果显示，果卜滑坡二维时序形变特征与GPS监测结果较为一致，相关系数大于0.9，且2016年—2021年坡体形变规律呈线性持续增加，垂直向最大累计沉降量达2 m以上。此外，坡体在雨季出现加速变形，而坡底受库区储水作用影响，形变存在滞后现象。结果表明，该文的技术方法可有效地获取库区滑坡多维形变规律，可为库区边坡的监测预警与安全防控提供技术支撑。
0 引言
库区河流蓄水与水电大坝修建等人类工程活动使得流域岩土体结构与应力平衡遭到破坏，易诱发坡体开裂、滑坡，甚至引起灾难性的坝体溃决[1-3]。由于库区坡体的失稳过程常具有明显的形变特征，因此可利用高效的监测手段获取地表的形变信息，从而为研究库区坡体失稳因素以及预警预报提供技术支撑。
目前，合成孔径雷达干涉测量（interferometric synthetic aperture radar, InSAR）目标由于其全天候作业、覆盖面广等优势，已在库区滑坡监测中得到广泛应用[4-8]。针对InSAR相位监测存在的时空去相干问题，基于永久散射体（persistent scatterers, PS）[9-11]或分布式散射体点（distributed scatterers, DS）目标[12-13]的时序InSAR技术通过高质量点目标选取、小时空基线干涉对组合或基于最大似然估计的相位优化等方式，可捕获滑坡体的有效形变信号。常见的时序InSAR方法包括斯坦福永久散射体技术（Stanford method for persistent scatterers, StaMPS）[14]、临时相干点InSAR（temporarily coherent point InSAR，TCPInSAR）[15]以及相干点目标分析法（interferometric point target analysis, IPTA）[11]等。然而，由于库岸边坡散射强度相对较低，以上几种基于单主影像的时序InSAR方法仅能获取稀少的监测点位，难以有效捕获坡体的失稳过程。此外，这些方法主要通过估计线性形变，进而叠加非线性成分来获取地表变形。而库区边坡的变形常常受降雨、库水位等因素影响，存在着较大的非线性成分。因此，本文针对InSAR库岸边坡的监测条件与误差特性，结合多主影像思想优化IPTA方法的技术流程，基于短时空基线干涉对与高相干点目标，通过IPTA回归分析迭代去除误差项，再基于加权最小二乘反演获取滑坡体的形变特征。
此外，库区边坡地形起伏较大，根据SAR卫星侧视成像特点，仅依靠单一轨道沿视线向（line of sight, LOS）形变通常会造成监测结果的歧义或误判。同时，单一LOS方向的监测结果难以满足复杂地形条件下库区滑坡形变解译的需求。InSAR获取的形变是东西、南北以及垂直方向沿LOS向的投影矢量总和，因此直接解算三维形变需要联合至少3个不同轨道LOS向观测，以解决观测方程的秩亏问题[4,16]。然而，由于南北向形变在LOS向的投影系数较小，对于形变方向主要集中在东西与垂直方向上的坡体，南北向形变可忽略不计。因此，可通过联合两个轨道观测实现滑坡二维时序形变反演，以揭示滑坡地表变形规律。
目前,针对果卜滑坡变形监测的相关研究主要包括联合光学遥感影像与InSAR的灾害调查[5-6]，基于SAR强度以及点目标的偏移量跟踪算法研究[7]，联合InSAR、带宽分频以及水文模型计算了多维变形、岩体的水力扩散系数以及厚度[8]等。然而，一方面果卜滑坡的形变特征需要进一步的回溯与更新；另一方面，受限于SAR数据的重访周期与相干性，前人研究主要基于单一数据来进行大梯度变形获取[7]，缺少基于多个成像视角来解算坡体的多维形变，且坡体的形变规律与失稳特征也需要进一步的探究。因此，本文将利用InSAR技术开展果卜滑坡的形变监测，获取坡体的多维变形规律，从而为分析其失稳因素以及灾害预警提供技术支撑。
综上所述，本文以位于我国黄河中上游的果卜滑坡为例，开展联合IPTA与两轨道SAR影像的滑坡二维形变监测研究。基于改进IPTA方法实现库区滑坡单一LOS形变监测，联合两轨道Sentinel-1存档数据反演二维形变特征，并结合GNSS监测验证形变结果的精度；提取时序形变的趋势项，从而开展滑坡体形变特征、影响因子的分析。
1 技术方法
1.1相干点目标分析法（IPTA）
IPTA与传统InSAR干涉解算方法采用的相位模型一致，即相位解缠后的干涉相位为形变相位、地形相位、轨道以及大气误差等相位的总和。
IPTA的主要步骤可分为[11]点目标选取、差分干涉相位提取、误差迭代改正、线性形变估计以及时序解算等。该方法的基本思路是选取具有稳定后向散射特性的目标点后，通过一维或者二维回归分析的方法估计点目标的线性形变速率、轨道误差与数字高程模型(digital elevation model,DEM)误差等，再基于时空域滤波分离残差中的大气延迟与非线性形变，从而可以获取每个目标点的时序形变。
由于库岸边坡散射强度相对较低，基于单主影像的PS方法仅能获取稀少的监测点位，而IPTA也可基于传统小基线子集算法（small baseline subset，SBAS）多主影像的时序解算思路，基于短时空基线的干涉对组合，选取点目标后提取对应的差分干涉相位进行相干点目标分析，并迭代估计轨道误差、DEM误差、大气误差，优化点位从而获取形变速率，最后基于奇异值分解（singular value decomposition，SVD）获取时序形变。
1.2 二维形变分解
由于升降轨Sentinel-1数据的入射角与航向角等几何参数不同，因此在获取各个轨道沿LOS向的形变后，可以基于式（2）进行果卜滑坡的多维分解。果卜滑坡的坡向整体上为东西朝向，且InSAR相位监测南北向的投影系数远小于其他两个方向[8, 16]，因此,本文忽略了南北向形变的影响，直接分解东西向与垂直向的二维形变场。此外，由于两组SAR数据的重访作业时间不同，而二维时序解算需要将时间进行统一。根据研究区域的形变特征与技术流程整体上的简便性，本文选取了三次样条插值进行升降轨时序的统一，进而实现后续的二维形变解算。
2 研究区域概况与数据处理
2.1 研究区域概况
果卜滑坡位于我国黄河流域中上游，边坡上下游两侧分别为我国的龙羊峡水电站与拉西瓦水电站。库区平均水位约为2 452 m，自2009年大坝蓄水后，距拉西瓦水电站上游500~1 300 m处的果卜边坡发生失稳，在坡体后缘形成了一条约21 m的拉裂带。受风化、降雨等因素影响，果卜边坡仍处于活动状态，威胁着大坝的安全生产。
2.2 数据概况
本文选取了果卜滑坡2016年1月18日—2021年8月10日的升轨与降轨Sentinel-1影像，其基本参数如表1所示。对于单个轨道的Sentinel-1数据，时间采样间隔一般为12 d，极化方式为VV极化。考虑到果卜滑坡DEM数据的时效性，本文选取了30 m分辨率的AW3D30（ALOS World 3D - 30 m）数据。
表1 研究区域Sentinel-1数据参数列表
此外，由于升降轨SAR影像的成像几何不同，可以提供不同的观测向量，从而分解获取滑坡体的多维形变场。根据表1中SAR数据在各方向上的投影系数，单个轨道LOS向形变在南北向的投影系数约为0.1，远小于东西向与垂直向，即表现为对南北向的形变信号不敏感，而果卜滑坡整体上为东西朝向，因此，本文仅针对滑坡的二维形变开展研究。
2.3 数据处理
如图1所示，本文联合2016—2021年升降轨Sentinel-1数据，基于IPTA和多几何成像关系获取了果卜滑坡的垂直向与东西向形变速率与时间序列，主要步骤如下： ①获取覆盖研究区域的升降轨Sentinel-1影像，并基于增强谱分级算法、精密轨道数据以及外部DEM配准SAR影像；②分别设置时空基线阈值36 d和200 m，生成差分干涉图，并剔除去相干严重的干涉对，干涉对组合结果如图2所示；③根据幅度离差与相干性进行联合选点，再提取点目标对应的干涉相位，选点结果如图4所示；④基于IPTA进行干涉点目标分析，改正轨道、大气、DEM等误差项，获取初始形变速率；⑤迭代移除各项误差后，基于最小二乘反演获取LOS向形变速率与形变时间序列；⑥将升轨、降轨SAR数据获取的LOS向形变时间序列进行空间域与时间域统一，即搜索同名点，并在时间域进行插值，反演坡体二维形变；⑦基于GNSS数据监测结果进行精度验证。
3 结果分析
3.1 LOS向形变监测结果
本文首先基于IPTA算法，联合升轨与降轨Sentinel-1数据分别获取了果卜滑坡不同成像几何视角下的形变速率与形变时间序列。本文基于幅度离差与相干性进行初始选点，采用幅度离差阈值为0.65，相干系数阈值为0.25，选点结果如图3所示。其中，底图分别对应SAR坐标系下的升轨与降轨SAR强度影像，红色点位则为初选的点位。整体上，果卜库岸边坡的坡顶与坡体上初选的点位分布均匀，地表高散射强度区域均有一定密度的点位覆盖，可被有效用于后续时序InSAR解算。
经过数据解算，获取了如图4（a）和图4（b）所示的果卜滑坡沿LOS向升轨和降轨形变速率。结果表明，两个轨道的观测形变范围基本一致，形变速率的量级存在较大差异，推测是由于成像几何视角不同。从图4（a）可以看出，坡体与顶部台塬位置的形变存在明显差异，即在坡体上表现为正值，在台塬位置则为负值。而对于图4（b）的降轨监测结果，滑坡体的形变整体保持负值。根据SAR影像的入射角与航向角，可计算获取如表1所示的投影系数。显然，垂直向与东西向投影系数较为接近，而南北向分量远小于其他方向。
由于坡体的滑动方向朝向升轨轨道的视线方向，而远离降轨轨道的视线方向，当水平向形变的投影分量大于垂直向时，会出现如图4（a）、图4（b）所示的误判现象。因此，需结合两者的几何成像差异，开展二维形变分解研究，获取果卜滑坡的多维形变特征。
3.2 滑坡二维形变分解结果
通过联合果卜滑坡升降轨的时序形变，分解获取了如4(c）~图4（d）所示的东西向与垂直向形变。本文在二维形变解算时，除考虑时间域的插值处理外，在空间域仅对两者的同名点进行了二维分解，因此图4（c）~图4（d）中点位的密度要小于4（a）、图4（b）。
对比4（c）~图4（d）中东西向与垂直向形变场，在形变量级上，坡体的东西向形变要大于垂直向的形变分量。由于果卜滑坡的坡顶与坡体均有明显的变形，且坡体两侧地形延伸方向不同，其水平移动的特征也可能不同。因此，为了尽可能地分析滑坡体的形变特征，分别在坡顶、坡脚、坡体两侧共提取了4条剖线（A-A’、B-B’、C-C’、D-D’），如图4（d）中4条白色剖线所示。图5为4条剖线对应位置的二维形变速率和地面高程信息。图5（a）显示，台塬位置整体上东西向形变小于垂直向形变，而坡体前缘东西向形变反而大于垂直向。综上所述，果卜滑坡坡顶的台塬位置垂直向形变大于东西向，而坡体的垂直向形变小于东西向，因此升轨Sentinel-1数据获取的形变速率在滑坡后缘呈负值，而在坡体前缘呈正值。
3.3 基于GNSS的精度验证
为检验基于Sentinel-1升降轨数据获取的二维形变时间序列的精度，本文获取了果卜滑坡2020年7月—2021年7月的GPS监测结果进行对比分析，GPS点的位置如图4（c）所示。由于二维形变场的解算主要基于IPTA获取的LOS向形变，仅时间域进行了插值处理，空间上的点为升降轨数据的同名点，因此图4（c）~图4（d）中点位密度较为稀疏。此外，InSAR获取的有效形变信号主要来源于库区边坡表面的强散射体，在与GPS监测结果进行对比时，同一位置InSAR可能由于散射强度较低点位缺失或难以保证精度，因此需要在对应点位选取一定大小的窗口，将内部所有点位的InSAR监测结果的均值与GPS进行对比，保证对比结果的可靠性。本文在GPS点对应的位置设定了2个像元大小的矩形窗口，将窗口内监测点位的均值作为InSAR的值，从而实现GPS与InSAR的对比，如图6所示。同时，本文计算了GPS点与对应InSAR时间序列的互相关系数。对于点GPS1，在东西向与垂直向上的相关系数分别为0.998 9、0.998 2，而对于点GPS2点相关系数分别为0.999 2、0.998 4。显然，两者趋势较为吻合。由于提取InSAR形变时序时采取了均值处理，因此在形变量级上有轻微的偏差。结果表明，联合IPTA与升降轨Sentinel-1数据获取的果卜滑坡二维形变精度较为可靠。
3.4 果卜滑坡失稳诱发因素分析
为进一步分析果卜滑坡的形变规律，分别从果卜滑坡的坡顶、坡体、坡底以及滑动面边缘位置选取了如图4（c）所示的特征点P1~P4，提取了4个点对应位置的形变时间序列，如图7（a）、图7（c）、图7（e）、 图7（g）所示。整体上4个点位的形变趋势较为一致，呈线性形变特征，其中点P3处在2016—2021年坡体最大形变量大于2 000 mm。点P1、P2以及P4处，东西向形变均大于南北向形变。根据特征点的形变时间序列可看出，果卜滑坡仍处在近似匀速下滑的活动状态，因此需要对其失稳因素进行进一步的探究。
Sgolay滤波器(Savitzky-Golay smoothing filter)是一种在时域内基于局域多项式最小二乘法拟合的滤波方法，在滤除噪声的同时可以确保信号的形状、宽度不变，被广泛运用于数据平滑除噪[18]。因此，本文在获取果卜滑坡二维时序形变后，利用Sgolay滤波器获取其趋势项并予以分离，最终分离结果分别如图7（b）、图7（d）、图7（f）、 图7（h）所示。其中，图7（a）、图7（c）、图7（e）、 图7（g）中线条分别表示P1~P4的4个特征点时序的趋势项，图7（b）、图7（d）、图7 (f)、 图7（h）为时间序列去除趋势项后的余项。
图7（a）、图7（c）、图7（e）、 图7（g）中模拟的趋势项与观测值基本一致，而图7（b）、图7（d）、图7 （f）、 图7（h）中去趋势项则表现为明显的周期性变化特征。其中，P1~P3点位于坡体的顶部与坡体中部，3个点位的去趋势项在雨季表现为加速的下降趋势，且周期与降雨量较为相似，证明了降雨对坡体变形的影响。P4点位于坡体底部，该点的去趋势项的周期性则存在一定的滞后特征，这与库区在蓄水期水位逐渐上升以及在雨季库水位下降较为吻合。
降雨一方面增大了岩土体重量，使得滑坡体的下滑力增强，抗剪强度降低，且降雨易通过后缘裂隙入渗，进一步加速滑坡的运动。除降雨作用外，库岸滑坡底部更易受库区水位的影响，在库水位下降时即雨季来临前的排水期出现下沉趋势。综上所述，果卜库岸滑坡的坡体和坡顶主要受降雨作用，在雨季出现加速变形，而坡底形变则受库区储水作用影响而出现滞后特征。
4 结束语
本文针对InSAR库区滑坡监测的误差特性与应用特点，基于优化的IPTA算法，联合升降轨Sentinel-1数据，获取了果卜滑坡2016—2021年的二维形变时间序列，根据收集到的GPS监测数据验证了InSAR二维分解结果的可靠性，并结合Sgolay滤波器获取并分离了滑坡体二维形变时间序列的趋势项，且与降雨等资料对滑坡的失稳因素进行了分析。结果显示，基于多轨道的PS-InSAR算法有效地获取了果卜滑坡的二维变形，整体上呈线性形变趋势，坡顶与坡体在雨季出现加速变形，而坡底受库区储水作用，在时间上存在一定的滞后特征。本文所提方法与技术流程可为库区滑坡监测的工程化应用提供技术支撑，也为黄河中上游库区的安全生产提供可工程化技术推广的应用案例。
【作者简介武志刚（1971—），男，甘肃天水人，教授级高级工程师，博士，主要研究方向为水利水电工程。
转自：“测绘学术资讯”微信公众号
如有侵权，请联系本站删除！
