InSAR滑坡监测研究进展

朱建军，胡 俊，李志伟，孙 倩，郑万基

1. 中南大学地球科学与信息物理学院，湖南 长沙 410083； 2. 有色金属成矿预测与地质环境监测教育部重点实验室(中南大学)，湖南 长沙 410083； 3. 湖南师范大学地理科学学院，湖南 长沙 410081

降雨、融雪、地震及人类工程活动，伴随滑坡本身的重力作用，会使滑坡体的剪应力超过其抗剪强度，打破其本身的应力平衡状态，致使滑坡变形加速，从而引发一系列地质环境问题和灾难(如山体崩塌、泥石流、堰塞湖和建构筑物损毁)。开展滑坡变形监测对于了解滑坡致灾因素，评估灾害影响，分析滑坡变形机理具有至关重要的作用。传统手段一般采取人工实地调查的方式，通过布设GNSS监测点或地下测斜仪来实现滑坡形变监测。这种方法虽然监测精度高、时间采样率快，但是只能监测滑坡在某些点位上的形变情况，空间分辨率很低，不能反映滑坡整体的形变[1]。在一些人力难以到达的高陡山区往往会出现漏判的情况，致使一些“高位远程”滑坡未能被有效地发现和监测，从而引发意想不到的灾难。此外，监测仪器的布设往往是在地表产生了显著变形后才开展，因此监测结果无法追溯滑坡的历史变形情况，制约了滑坡参数的反演和机理解译工作。

合成孔径雷达干涉测量技术(interferometric synthetic aperture radar,InSAR)是近50年来逐步发展起来的一种空间大地测量手段，该技术通过主动发射微波信号获取地表相位，通过一系列的信号处理后得到地面高程或者地表形变信息。InSAR技术具备空间分辨率高、不受云雾影响、全天时监测和覆盖面积广等特点，已被广泛应用于包括地震、火山、地表沉降和冰川运动等领域的监测工作中[2]。得益于该技术空间覆盖广的特点，在InSAR技术发展的早期，就有学者利用该技术进行滑坡监测工作，并取得了良好效果，其监测精度也在试验中得到了很好的验证[3]。然而，受限于早期SAR卫星发展不成熟、SAR数据少、成像模式单一及轨道质量不理想等原因，利用InSAR技术开展滑坡监测仅在少数区域得到了应用。

近年来，随着新一代SAR卫星的入轨运行，InSAR技术在滑坡监测领域呈现出了井喷式的发展，已逐步成为了滑坡监测领域的常用手段之一。宽幅成像模式为实现广域滑坡探测奠定了基础，多星协同观测带来的超短重访周期可以提供更多滑坡变化的细节，理想的轨道控制确保了全球陆域上都可以有丰富的观测数据。基于这一系列新特征，SAR工作者们在进行滑坡精细监测[4-5]、滑坡参数反演[6-7]、滑坡预测预报[8-9]和滑坡机理解译[7,10]等方面都做出了令人振奋的成果。

然而，受限于SAR成像几何和微波电磁特性的影响，即使利用新一代SAR卫星开展滑坡监测，地形起伏导致的几何畸变、密集植被覆盖产生的失相关噪声、水汽变化引起的大气延迟误差等因素依然制约着InSAR技术在滑坡监测领域的应用。针对这一系列问题，国内外学者在近年来开展了大量的攻关工作。为系统地梳理该领域的研究进展，本文总结了当前InSAR滑坡监测技术在数据选择、精度提升和维度扩展等方面的主要制约因素，并归纳目前主流的解决方案，进而探讨目前InSAR滑坡监测亟待攻克的技术难点及未来发展方向，以期为读者提供一个了解该领域最新研究进展的途径。

1 InSAR滑坡监测的数据选择

滑坡往往发生在地形起伏较大、植被覆盖茂密的山区，由于SAR卫星侧视成像几何的固有限制及植被体的高动态特性，InSAR技术在山区滑坡形变监测中更容易受到几何畸变和时间失相干的制约，导致观测失败或结果不可靠[11]。因此，在开展InSAR滑坡形变监测前，选择合适的SAR数据至关重要。针对滑坡区地形和地面覆盖类型，需要考虑不同的因素，对其进行可行性分析，选择最合适的数据，从而充分发挥InSAR技术在滑坡形变监测中的最大性能。例如，当滑坡面向东面时，选用升轨、小入射角的卫星监测更有利；当植被覆盖茂密时，长波SAR数据比短波SAR数据具有更强的穿透性，能够更好地抵御时间失相干。鉴于此，本节将详细阐述地形和植被对InSAR滑坡形变监测结果的影响，为数据选择提供参考依据。

1.1 几何畸变的影响：成像几何的选择

几何畸变的形成是由于SAR卫星固有的侧视成像模式所导致的，根据卫星和地面的相对几何关系，可将其划分为透视收缩、叠掩、阴影和分辨率增强4种类别[11]。当滑坡面向卫星时，如果入射角大于坡度角，会导致投影到SAR影像上的斜坡长度被压缩，形成透视收缩(图1(a)中AB)，且入射角越小，透视收缩越严重；如果入射角小于坡度角，会导致投影到SAR影像上的斜坡顶部和底部颠倒，产生主动叠掩(图1(b)中EF)。此外，受到主动叠掩的影响，相邻区域也会产生被动叠掩现象，其中，靠近卫星的区域称为近被动叠掩(图1(b)中DE)，远离卫星的区域称为远被动叠掩(图1(b)中FG)。由于叠掩导致更多的散射回波相互混叠在一起，叠掩区域会在SAR图像上表现得比周围更亮。当滑坡背向卫星时，如果入射角小于坡度角的余角，会导致投影到SAR影像上的斜坡长度被拉长，导致分辨率增强(图1(a)中BC)，且入射角越小，分辨率增加越明显；如果入射角大于坡度角的余角，会导致卫星无法照射到背坡区，因而产生主动阴影(图1(c)中HI)；另一方面，由于传感器和地表之间存在障碍而无法获取到有效信号的区域称为被动阴影(图1(c)中IJ)，阴影在SAR图像上表现为黑色。在InSAR滑坡形变监测中，分辨率增强区能够获得最好的结果；透视收缩区由于分辨率降低，可能会导致形变监测结果质量下降；而叠掩和阴影区域由于信号混叠或缺失，导致干涉相位不连续，会造成相位解缠误差[12](表1)。由此可见，准确识别几何畸变区域的范围，能够有效辅助InSAR滑坡形变监测的数据选择，提高形变监测结果的质量。

图1 几何畸变示意Fig.1 Illustrations of geometric distortion

表1 不同几何条件下InSAR观测性能

近年来，国内外大量学者对几何畸变区域的有效识别和提取进行了更为深入的研究，主要有两种不同的策略。第1种策略直接通过传感器和地表的相对几何关系来精确确定几何畸变的范围，具有代表性的方法是叠掩阴影图法(LSM)[13]、斜距指数法(R-Index)[14]和相邻梯度法(P-NG)[15]，后两种方法是在前者的基础上发展而来的。LSM方法能够通过卫星位置、方位角和地面高程等参数来模拟几何畸变主动区和被动区，但该方法计算较为复杂，效率较低[13]。R-Index法是通过定义雷达斜距和地距的比值来识别几何畸变区域[14]。随后经过其他学者的不断改进，R-Index法已经发展得较为完善，能够较准确地识别几何畸变范围[12]。P-NG方法是基于几何关系和相邻像素的梯度来确定几何畸变的范围[15]。R-Index和P-NG两种方法效果相当，不同之处在于R-Index法效率更高，但在结果中将叠掩和阴影区合并，没有对两者进行明确的区分，也没有区分主动区和被动区；而P-NG方法将主动区和被动区准确区分开来，同时还考虑了几何畸变区域的多重混叠，但缺点是逐像素运算效率较低。第2种策略基于SAR信号的统计特性。文献[16]基于叠掩区域信号满足高斯分布的假设来识别叠掩区域。文献[17]利用核密度函数与恒虚警率相结合的方法来确定几何畸变范围。基于统计的方法虽然能够取得一定效果，却无法避免阈值的设置，因此基于几何的方法更为准确可靠。此外，除了以上两种策略，还有GIS模型法[18]、光谱技术[19]等方法。目前，准确识别几何畸变范围的方法已经相对成熟，能够有效避免SAR数据的错误选择。

克服几何畸变影响的方法主要有两种。一种是有效识别几何畸变区的范围，并对叠掩和阴影区域进行掩模，以降低这些区域对相位解缠的影响，然而这种方法必然会导致监测范围的降低[20]。另一种方法是利用升降轨数据叠掩和阴影区的互补性来扩大监测范围[21-22]。图2为甘肃省舟曲的Sentinel-1卫星升降轨的几何畸变结果。图2(a)为升轨数据的几何畸变结果，图2(b)为降轨数据的几何畸变结果。由图2可以看出，升轨和降轨的几何畸变区域相互补充，融合两者能够有效扩大监测范围。然而，仍有部分区域难以同时获取升降轨数据，该方法依旧存在很大的挑战。一般而言，分辨率增强区的质量最好，透视收缩区的分辨率降低，可能会导致监测结果不可靠，如何准确评估透视收缩区与分辨率增强区的信号水平及其可用性，尚缺少相关研究。此外，随着极化SAR技术和立体SAR技术的发展，目前已有部分研究试图通过统计信号分离[23]、SAR三维重建[24]等技术来恢复叠掩区的观测信号，但尚未发展成熟，因此如何准确恢复叠掩区的混叠信号仍然是当前研究的难点。

图2 P-NG方法获取的甘肃舟曲地区的几何畸变结果Fig.2 Illustrations of geometric distortions based on P-NG method in Zhouqu, Gansu province

1.2 植被覆盖的影响：数据波段的选择

滑坡灾害频发区往往被植被所覆盖，植被冠层引起的微波信号衰减、植被高动态变化引起的时间失相干都会导致InSAR技术无法准确捕捉到植被覆盖层下的滑坡形变信号，极易出现监测盲区。由植被引起时间失相干的原因有很多，包括植被的季节性生长变化、自然风所导致的植被姿态扰动、降雨引起的介电常数变化、人为活动或自然灾害所引起的土地变迁(例如伐木、林火)等[25]。此外，研究证明，植被区的失相干程度随着树木的高度和密度的增加而增加[26]。因此，如何实现植被区InSAR滑坡形变监测是当前亟待解决的难题。

在研究中，学者们发现，长波SAR数据具有更强的穿透能力及抵御时间去相干的能力，可以获取植被层下更加丰富的信息，有利于相干性的准确估计。文献[27]证明了在高植被覆盖下滑坡形变监测中，L波段比C波段具有更强的穿透性，能够更好地抵御失相干的影响，获得更高质量的形变结果。文献[28]在对美国加州植被区进行失相干分析时发现，L波段的ALOS数据在森林区能够在超过两年的时间内保持较高的相干性，而C波段的ERS数据则不到6个月就会产生严重的失相干。文献[29]分析了L、C和X 3个不同波段数据在草地环境下的失相干速率，证实了L波段数据能够获得更高的相干性，且高重访周期和高分辨率条件也能够提高相干性。文献[30]分别使用C波段的RADARSAT、ENVISAT数据和L波段的ALOS数据对荷兰芬丹地区的农田进行了监测，结果表明穿透性更强的L波段数据可以更好地克服时间去相干的影响，相比C波段数据更适宜用于植被区地表形变监测。文献[31]使用L波段数据，采用相干散射体测量的方法识别了四川丹巴县滑坡隐患点，证明L波段数据具有更好的相干性，其在我国西部山区的滑坡普查和监测中更具优势。图3为X波段的TerraSAR卫星、C波段的Sentinel卫星和L波段的ALOS-2卫星在三峡植被覆盖的滑坡区差分干涉结果，所用数据见表2。由图3可以看出，随着波长的增加，观测点的密度增加，相位质量也显著提高。由此可见，在开展植被区InSAR滑坡形变监测时，选用长波SAR数据能够获得更好的结果。

表2 图3中所用数据参数

图3 三峡滑坡差分干涉相位Fig.3 Interferograms of different sensors in the Three Gorges areas

目前，除了在轨运行的L波段日本ALOS-2卫星、阿根廷SAOCOM卫星和我国陆地探测一号 (LT-1) 卫星，未来会有L波段的德国TanDEM-L卫星、日本ALOS-4卫星及欧空局P波段的BIOMASS卫星等长波雷达卫星发射，将为开展植被区长波InSAR滑坡形变监测提供充足的数据支撑。

为了能够更准确地选择合适的数据，最直接有效的方法是定量衡量植被参数与InSAR相干性之间的关系。目前，部分学者对此开展了初步的探索。文献[32]发现植被的高度和密度是影响相干性的重要因素。一般情况下，InSAR相干性随植被高度和密度的增加而减小。此外，不同类型的植被覆盖下，InSAR的相干性也不同。随后，越来越多的研究发现相干性与植被指数之间存在着显著的关系[33-36]，这为InSAR失相关的定量预估提供了契机。然而大多数研究局限于定性分析，无法得到相干性和植被参数之间的定量关系。针对此问题，文献[37]利用ALOS和Landsat 5数据建立了归一化植被指数NDVI与InSAR时间失相干之间的三阶段模型，该模型能够在SAR数据获取前预先定量评估不同植被覆盖度下InSAR数据的表现性能，从而指导选择合适的SAR数据，该研究为植被区滑坡InSAR形变监测的数据选择提供了一种新的思路。综上所述，不同的研究区域环境千差万别，基于外部数据或经验模型事先评价不同数据的相干性，可为数据的选择提供参考依据。

2 InSAR滑坡监测的精度提升

星载SAR卫星运行高度一般在距地500 km以上，在利用InSAR技术开展滑坡监测的过程中，微波将穿过大气层和植被覆盖，二者也制约了InSAR监测滑坡的精度。具体而言，在利用InSAR开展城市变形监测时，由于城市具备较多的建构筑物，其每个像素中一般都有一个永久散射体(permanent scatters,PS)，但是在滑坡监测中，一般难以获取足够多的PS点，以至于InSAR监测点的密度较低。因此，往往需要采用分布式散射体(distributed scatters,DS)作为观测对象，在对DS点的观测相位进行优化后，方可获取密度较高的滑坡监测结果。另一方面，滑坡一般分布于地形陡峭的山区环境下，其大气环境较为复杂，在对应的InSAR观测相位中往往也会存在严重的大气误差，从而降低InSAR的监测精度。据此，本节将分别从InSAR观测相位优化和大气噪声抑制两方面对目前InSAR滑坡监测精度改善的相关工作进行介绍。

2.1 观测相位优化

滑坡体属于自然目标，通常位于非城市地区，地表以稀疏的植被覆盖为特征并且表现出典型的低相干性，使得传统时序InSAR探测到的测量点的分布密度低，造成监测困难。DS表示一个包含所有相近散射幅度值的基本散射体且没有主导散射体的散射单元，对应像素的散射值等于所有基本散射体后向散射的矢量和[38]。DS一般都分布在低相干区域内，例如城市郊区或者山区。如果SAR影像分辨率约为10 m，那么大多数的自然地物散射点(森林、农田、裸露土壤、岩石表面)都可以被认为是DS；如果SAR影像分辨率达到几米，那么DS大多出现在具有低矮植被覆盖的干旱地区。综上所述，DS的特征与滑坡区域的特征具有良好的吻合性，如果能充分利用DS目标来对滑坡进行监测，便能够突破传统时序InSAR方法所受到的限制。一般而言，基于分布式目标的地表形变监测涉及两个关键性步骤：空间维DS目标同质像元的自适应提取与时间维相位一致性约束条件下的DS目标最优相位估计。

2.1.1 空间维DS目标同质像元的自适应提取

DS-InSAR将二阶统计量样本协方差矩阵(或相干矩阵)作为观测值进行输入，因此样本协方差矩阵(或相干矩阵)的正确估计将直接影响到最后对形变相位参数的解算精度。通常相干矩阵估计需要提供较大的窗口范围，但是其中会有许多与中心像元统计特性不同或非同质的像元，因此需要对同质像元进行自适应地识别和提取。当前的同质像元识别方案大多基于假设检验的基本理论，大致可分为非参数检验与参数检验两类。第1类非参数检验方法，最早由文献[39]提出，将双样本KS(Kolmogorov-Smirnov)方法应用在同质点选取上，通过对时序振幅样本的经验累积分布函数进行比较分析，从而确定其是否为同质像元。同时该文还提出了双样本AD(Anderson-Darling)方法进行同质点选取，AD方法在KS方法的基础上改进了对于尾部差异的探测。KS方法在剔除高振幅的离群值之后，其快速和便捷的优势便得以体现，面对大范围、长时序SAR数据的同质点选取时具有优势。AD方法和KS方法相比具有更好的表现[40]，其对较大振幅值的灵敏度更高，结果更稳定，但是耗时更长。第2类的非参数检验方法，应用较多的是广义似然比检验(generalized likelihood ratio test，GLR)[41]，其基于像元时序振幅信息服从瑞利分布的假设，通过对两个样本的方差统计量进行显著性分析，可以有效地用于同质像元的识别，但是当分布的形状参数产生变化时，其结果会受到很大影响[40]。文献[42]提出了快速同质像元识别(FaSHPS)算法，在基于振幅样本均值服从高斯分布假设条件下，将传统假设检验问题从显著性差异判断转换为置信区间估计的方法。只要置信区间确定，就可以通过简单的逻辑运算实现样本之间的相似性判断，很大程度上提高了同质像元识别的计算效率。通常情况下，非参数方法适用性更好，可广泛应用于各种样本分布类型，但可能不够精确。而对于参数方法，当其基本假设成立时，结果会优于非参数的方法，并且在进行小样本检验时结果也更好。但缺点是稳健性不足，容易受到样本分布差异的影响。目前而言，没有通用的自适应同质点选取方法，实际应用中需要兼顾研究区特征、样本数及算法运行效率等因素来选取对应的方法。

2.1.2 时间维相位一致性约束条件下的DS目标最优相位估计

假定存在3景已完成配准处理的SAR影像M0、Sm、Sn，其中M0作为主影像，其余两景作为从影像，将它们两两之间进行干涉可以得到三景干涉图。在单视条件下，会满足如下关系

φmn=W{φom-φon}

(1)

式(1)即为相位三角条件，也称为相位测量一致性，这是时序InSAR技术的最基本假设[43]。空间自适应同质滤波处理往往会引起干涉相位不一致[44]，因此，针对DS多视干涉相位而言，有必要采取一定的相位优化算法在满足相位一致性的条件下构建出一组单主影像优化相位值，以尽可能地减弱分布式目标去相干现象的影响。

对DS的历史相位参数的估计方法大致可分为3类，第1类为极大似然估计法(maximum likelihood estimaton,MLE)，最早由文献[45—46]提出，是一种渐进的无偏估计方法。MLE基于DS去相干噪声为高斯分布的假设，利用时序上所有的干涉相位信息，构造极大似然函数最大化准则以求取DS历史相位的最优值。但估计相干矩阵时直接采用了窗口内的所有像元，并未考虑同质像元的空间分布。文献[47]在此基础上引入同质点自适应选取技术，将其发展为SqueeSAR方法。此后，文献[48]发展了加权的MLE估计，使得估计结果更为稳健，且能够抵御去相关影响。MLE的估计精度高，但是由于解算时涉及高度非线性函数的迭代求最优解问题，运行效率受到了限制。第2类主流的DS历史相位参数估计方法是特征值分解(eigenvalue decomposition,EVD)方法，其灵感来源于极化SAR[49]。EVD直接对协方差矩阵进行特征分解，并选择最大特征值对应的特征向量作为DS历史相位估计值。文献[50]通过将EVD输入量由协方差矩阵替换为相干矩阵，从而减弱了振幅不平衡对于结果的影响。EVD方法能在一定程度上对不同类型的形变信号进行分离，并且计算效率高，但是同时也损失了更多的信息，降低精度。第3类方法是采用整数最小二乘(integer least squares,ILS)来估计DS缠绕相位最优解[51]，它采用Fisher信息指数进行定权，相比于MLE和EVD，ILS的估计结果精度最高，并且能利用协方差矩阵对优化结果进行质量控制。但缺点是计算的代价过于高昂，难以大范围推广应用。对于大数据量的运算效率问题，文献[44]提出了EMI方法，通过结合MLE和EVD方法各自的优点，在提高计算效率的同时改善结果精度。文献[52]进一步将EMI与序贯处理的思想相结合，将相干矩阵在时序上进行分块，使得海量SAR数据的DS历史相位的实时优化成为了可能。事实上，如果从模型的角度出发，上述算法的数学模型可以导出为相同形式[53]，而本质区别则在于其各自的协方差矩阵的获得方式和对权值确定方式的不同[54]。

此外，对于滑坡体监测而言，由于现有外部数字高程模型的精度有限，地形残差总是存在于差分干涉图中，地形误差有可能被误认为地质灾害诱发的变形信号。在时序InSAR技术中，地形残差是基线相关的分量，通常在时序InSAR中建模为未知参数与形变同时估计。在参数模型合适和长基线的条件下，时序InSAR技术能够准确地估计出相干目标上的形变和地形残差[55]。另一种替代方法是通过干涉图组合形成基线极短的伪干涉图来削弱地形残差的影响。如果在时序处理中只涉及基线极短的干涉图，由于与DEM相关的相位相当小，地形残差可以安全地忽略，此时形变将成为时序InSAR模型中唯一的未知参数[56]。

本文以青海省拉西瓦水电站果卜岸坡滑坡为例，利用44景2019年8月—2021年1月的Sentinel-1数据对比分析了观测相位优化的效果。试验采用单主影像构网，具体影像参数见表3，分别采用传统的PS-InSAR方法和DS-InSAR的优化方法对试验区进行形变解算，其中DS-InSAR方法采用极大似然估计(MLE)优化后相位值。无论在形变区域还是非形变区域上，DS-InSAR方法的监测点数量都大于PS-InSAR方法获取的监测点数量，证明了DS-InSAR在滑坡监测领域具有无可比拟的优势(图4)。

表3 图4中所用SAR影像的参数

图4 利用PS-InSAR方法和DS-InSAR方法得到的果卜滑坡区域的形变速率Fig.4 Deformation rate of the Guobu landslide derived by PS-InSAR method and DS-InSAR method

2.2 大气噪声抑制

在不同时间获取SAR影像时，电磁波传播路径上的温度、湿度、大气压等会有差异，从而使干涉图上出现大气相位屏。大气相位在空间上相关性较强，通常呈现出随高程分布或片状的形态；而在时间上的相关性较低，没有显著规律可循。大气相位的存在严重制约了InSAR的高精度地表形变监测，目前去除大气相位的方法可分为外部数据辅助的方法和利用大气相位时空特性的方法。

利用外部数据的方法一般会利用如GNSS、地面气象数据、MODIS、MERIS等数据，进行InSAR大气噪声的估计和抑制[57]。文献[58]提出利用参考点的气象观测值及对流层延迟模型来校正大气误差。但是，气象站的分布非常稀疏，有时难以获取足够多的数据，并且对流层模型建立不精确，以上都会带来额外的误差。GNSS气象学的发展为利用GNSS评估和校准InSAR测量中的大气影响提供了机会。文献[59]将Taylor的“冻结流”假设应用于InSAR大气校正，试验结果表明，加入SAR采集前后的GNSS观测值，可以进一步提高InSAR的精度。然而，GNSS站的空间分辨率通常远低于InSAR数据，这对应用GPS观测校正InSAR测量造成了限制。相较于GNSS数据，MODIS数据具有更高分辨率，对于InSAR的大气效应建模和校正是非常有用的。文献[60]首次使用InSAR和MODIS数据校正埃特纳山和洛杉矶上空的大气延迟。此外，文献[61]提出了一种将MODIS和GNSS数据集成用于InSAR大气校正的方法，在洛杉矶地区的试验表明，该方法显著地降低了InSAR相位中的大气信号，校正后InSAR测量中的地球物理信号更加突出。尽管MODIS具有较高的分辨率，但其测量结果对云层的存在非常敏感，这大大限制了MODIS在多云地区的应用。ENVISAT卫星上的MERIS分辨率高达300 m，精度高于MODIS。因此，MERIS为精确模拟大气对于ASAR数据造成的影响提供了机会[62]。以洛杉矶地区为例，文献[63]研究表明MERIS水汽数据可以显著降低SAR干涉图中的大气影响。

尽管外部数据可以有效提高InSAR形变测量精度，但其与SAR数据的空间分辨率的不一致性、测量时间的不同步性都使得该方法无法推广到所有场景中。因此不依赖于外部数据的方法，如时空滤波法、大气与高程建模法、独立成分分析(ICA)信号分离法相继提出，仅通过SAR数据本身就能估计大气相位。

时空滤波方法是最常用的去除大气噪声的方法之一，该方法在时序InSAR技术(如SBAS-InSAR和PS-InSAR)中被广泛应用。该方法认为形变具有时空上的高相关性，大气在空间上有较高的相关性而时间上无相关性。因此，时间域高通滤波结合空间域低通滤波可以将大气与形变分离开。当大气在时间上没有任何规律且形变解算模型与真实形变情况相一致时，这种方法有较好的效果。然而，实际的大气中的季节性变化，会使大气估计不准确，降低形变解算精度[64]。文献[65]指出应用滤波方法时，应适当考虑季节性的大气延迟。此外，合适的滤波窗口对准确的大气估计至关重要，文献[66]指出可以通过大气与形变的先验信息辅助确定窗口大小。针对InSAR干涉相位中的分层大气相位，则可以通过大气与高程建模法来去除。目前已有多种模型被用于描述大气延迟与高程的关系，如线性模型、三次多项式模型、三次样条模型、幂律模型。文献[67]考虑到大气的空间异质性，以幂律函数来估计大气。文献[68]将幂律方法进一步发展，将分块思想融入其中。文献[69]指出了幂律方法的局限性，提出从高分辨率天气模型中估算大气延迟的方法。文献[70]提出一种线性模型，应用于协方差矩阵加权的相位差，缓解分层延迟与部分湍流延迟。文献[71]通过四叉树将影像分块，小块内部相邻点的时序差分相位建模后再拼接来估计整张影像的高程相关大气。这种方法适合于形变区域有一定地形起伏，大气与高程相关的场景。不过实际的大气比较复杂，通过简单的数学模型较难实现完美的大气估计。

除了上述方法外，文献[72]通过独立成分分析(independent component analysis,ICA)的方法来估计大气相位。ICA是一种统计方法，用于将观测到的多维随机向量转化为尽可能独立的分量，其优势在于不需要建立观测值与未知数之间的函数方程且无须外部数据就能够实现很好的信号分离[73]。近年来，ICA在InSAR领域得到了广泛的应用，包括时序形变信号的分解[74-75]、地形残差估计[76]、大气信号分离[72]等。本文将ICA应用到美国南加州洛杉矶地区获取的Sentinel-1数据大气延迟改正中。因毗邻太平洋空气中水汽含量较高，该区域具有明显的大气延迟信号。为了使用ICA方法分离出大气成分，首先，需要确定独立分量的个数，利用特征值分解方法得到InSAR时序相位的每个主成分的方差百分比；然后，给定一个阈值，当方差百分比大于该阈值时，将其保留，否则视为误差去除。本文选取0.2%作为该阈值，前5个主分量的方差百分比大于该阈值。因此，保留前5个主成分，最终确定独立成分数量为5。利用ICA从时序相位中分离出5个独立分量，经过时空特征分析，判断第3个分量代表大气信号。对第3个分量运用ICA逆运算重建每一个干涉图的大气相位，将重建得到的大气信号从每一个干涉图中减去得到校正之后的干涉图。经过ICA大气校正前后的干涉图如图5所示，结果表明，在经过ICA改正后，大气相位得到了显著的抑制。

3 InSAR滑坡监测的维度扩展

滑坡运动发生在三维空间中，而InSAR只能得到滑坡三维形变在InSAR视线向(line-of-sight,LOS)上的投影，无法反映滑坡真实的运动情况，制约了利用InSAR进行滑坡机理解译、灾害预警的能力。为了重建真实的滑坡三维形变信息，一般需要至少3个成像几何差异显著的观测结果[77-78]。然而，这样的观测一般仅在机载SAR观测平台上可以获取[6]，在常用的星载SAR平台上一般只能获取两个成像几何差异显著(即升轨和降轨)的InSAR观测结果。为了突破这一限制，国内外学者近年来开展了探索性研究工作，主要可以分为两类：一是在经典InSAR相位观测的基础上，引入像素偏移量追踪技术(pixel-based offset tracking，POT)或者多孔径InSAR技术(multiple aperture InSAR，MAI)，获取滑坡的方位向形变信息，增加第3个方向的观测量；二是根据滑坡的运动特性，通过引入滑坡先验信息的方式降低获取滑坡三维形变信息所必须的观测量。两类方法在获取滑坡三维形变信息时各有优势，在应用时需要根据滑坡的具体特征及数据情况来选择。下文将针对两种方法分别介绍其研究现状。

3.1 升降轨InSAR与POT/MAI融合

目前在轨的SAR卫星基本都是以极轨飞行的方式运行的，观测方式一般以右侧视为主，因此，在多数区域，一般只能获取升轨和降轨两个成像几何差异显著的观测值。为了利用有限的观测获取滑坡三维形变信息，有学者提出可以引入POT技术的方位向观测。该方法最早由文献[79]提出，并被用于同震形变场获取。POT方法充分利用了SAR信号不受云雾影响的特点，基于SAR强度影像通过亚像元级影像配准计算形变前后SAR影像上同名像点的变化量同时获取LOS向和方位向形变。

利用该思想，有学者开始尝试将其应用于滑坡监测中。文献[80]使用升降轨的高分辨率TerraSAR-X数据获取了法国南部La Valette滑坡的三维时序形变结果，与GNSS的结果进行比对后，在水平向上监测误差为12 cm，垂直向上则为8 cm左右。文献[81]基于高分辨率TerraSAR-X数据，利用POT技术获取了植被覆盖区域的角反射器三维形变结果。POT技术除了可以在SAR数据上应用，也有学者尝试结合光学卫星进行滑坡三维形变监测。文献[8]提出结合ALOS-2/PALSAR-2数据和高分二号光学数据进行三维形变监测，基于该结果分析了白格滑坡的稳定性。为了实现亚像元级影像配准，多数研究采用的是归一化互相关系数法，但在滑坡的边缘，容易因为中心像元与其他像元存在显著运动差异而导致边缘区域形变估计精度的降低。为了解决这一问题，文献[82]提出了一种基于自适应窗口的POT方法，用于改善滑坡边界的形变结果，并获取了Slumgullion滑坡的三维时序形变结果。

利用POT技术获得的形变结果，其精度主要取决于数据的分辨率[83]。POT技术通过亚像元级影像配准技术获取像素偏移信息，目前该技术可以实现的最小像元变化探测度为1/10～1/20个像素[84]。因此，为了获取微小的滑坡蠕变及早期加速变形的形变信息，往往需要分辨单元更小的高分辨率SAR数据的支撑。图6是结合TerraSAR-X(图6(a))、Sentinel-1(图6(b))和ALOS-2/PALSAR-2(图6(c))卫星观测结果获取的新铺滑坡三维形变速率结果。其中，方位向形变是利用TerraSAR卫星的高分辨率聚束模式的数据，通过POT处理获得的。相较于Sentinel-1和ALOS-2/PALSAR-2数据，由于聚束模式的TerraSAR-X数据分辨率较高(方位向和距离向约为1 m)，在进行POT处理时，其对形变的敏感性远大于其他数据，可探测的最小形变量级达到了5～10 cm，因此，可以得到滑坡的方位向形变结果。进一步结合升轨Sentinel-1卫星和降轨ALOS-2/PALSAR-2卫星的InSAR观测结果后，利用最小二乘方法获取了该滑坡的三维形变结果(图6(d)—图6(e))。然而，通过融合POT观测虽然可以克服InSAR观测南北向形变不敏感的缺陷，但是获取的南北向形变精度相较垂直向和东西向仍然较低。

多孔径InSAR(multiple aperture InSAR，MAI)是一种基于相位的测量模式，该技术一般通过对SAR回波信号进行处理，将全孔径的SAR影像分为前视和后视两个部分，利用前视和后视干涉图在方位向形变上的差异性对两个影像进行差分处理后即可得到方位向形变。相较于POT获取的方位向形变结果，MAI获得的方位向形变精度更高，在地震、火山和冰川监测方面已有不少应用。在滑坡监测方面，文献[85]利用MAI技术结合DInSAR获取了露天矿边坡的二维时序形变结果，但是受限于MAI技术对失相关噪声较为敏感的原因，在滑坡监测上，该技术的适用性略差于POT技术，目前成功的案例较少。

3.2 先验信息约束的升降轨InSAR融合

POT只能监测形变显著的滑坡，而MAI对相干性要求非常高，因此对于地表环境复杂的缓慢滑坡而言，POT/MAI均难成为常规的形变观测手段。为进一步降低三维形变监测对于数据的严苛要求，有学者提出可以根据实际应用场景，利用先验信息约束来获取三维形变的方法[86-87]。在滑坡监测领域，最常用的先验信息约束方式主要有两种。

第1种是直接忽略南北向形变的方法，该方法利用InSAR观测对南北向形变不敏感的特点，假设南北向没有形变，利用升轨和降轨的InSAR观测即可得到滑坡东西向和垂直向的形变结果。基于该思想，文献[86]利用升轨和降轨的TerraSAR-X聚束模式数据获取了黑方台滑坡的二维形变结果分析了黑方台滑坡的破坏模式。

第2种方法是根据滑坡运动特点所提出的地表平行位移假设方法(surface-parallel flow,SPF)。该方法认为，地表上自由流动的物质主要受重力作用的影响顺着地势运动，因此，可利用地形作为先验信息构建三维形变分量内部之间的函数模型，实现仅以升轨和降轨InSAR观测即可获取滑坡三维形变的目的。该方法最早是由文献[87]提出，用于解决在极地获取三维冰川流速的问题，在后续的研究中，该方法已成为冰川监测中的常用方法。而在滑坡监测方面，该方法在理论方面比忽略南北向形变的处理手段更为严密。文献[88]利用ALOS-1升轨和ENVISAT降轨数据实现了我国甘肃舟曲滑坡的三维形变监测。文献[89]利用ALOS-1/2升轨和Sentinel-1降轨的数据，获取了美国克雷森特湖滑坡的三维形变结果，并在此基础上进一步反演了该滑坡的厚度信息。图7是在SPF方法约束下，利用升轨和降轨的Sentinel-1数据获取的InSAR观测结果解算的金坪子滑坡的三维形变信息。除了可以获取三维形变速率，文献[90]在SPF约束的基础上引入Tikhonov正则化的思想，用于解决在解算滑坡三维动态时序结果时升降轨观测时间不一致的问题，并利用COSMO-SkyMed升降轨观测得到了刚果地区Funu滑坡的三维时序形变结果。

因此，在充分利用先验信息约束的情况下，以有限的数据获取滑坡三维形变结果。但是，SPF方法需要获取精确的外部DEM来计算地形坡度。常用的DEM数据(如SRTMDEM和ASTERGDEM)的精度和分辨率均不高[91]，且获取时间常早于滑坡事件的发生，因此，DEM中的误差将不可避免地会引入函数模型中，导致不可靠的三维形变结果。此外，无论是忽略南北向形变的方法或者SPF约束，均需要至少两个方向的观测。然而，滑坡一般发生在山区，受限于SAR侧视观测的影响，在多数研究区都存在叠掩和阴影的问题，致使升轨和降轨中仅有一个轨道有观测结果，制约了滑坡三维形变的获取工作。

4 InSAR滑坡监测的发展方向

4.1 基于长波时序极化InSAR的滑坡形变透视监测

滑坡多发生在植被茂密的山区，植被的高动态变化所引起的时间失相干导致时序InSAR技术往往对于高植被覆盖度下的滑坡隐患识别与形变监测“束手无策”。长波SAR信号具有较强的穿透性，能够在一定程度上抵御时间失相干的影响[28-29]。未来，随着TanDEM-L、ALOS-4、BIOMASS、NiSAR等长波SAR卫星的入轨运行，植被覆盖区InSAR滑坡监测的失相干问题预计将会得到显著改善。

此外，有学者提出在时序InSAR中考虑引入极化测量信息来改善形变监测效果[92]。与单极化数据相比，多极化数据包含更多的地物散射信息[92]。同时，极化InSAR技术对植被体的空间结构、介电常数等物理属性更加敏感，不仅具备区分同一分辨单元内不同散射机制的能力，而且可以有效地刻画植被复杂散射过程，从而能够帮助重建更准确的地表相位[93]。目前，已有部分学者开展了时序极化InSAR形变监测的初步研究。一方面，通过数值优化的方法寻找最优散射机制来降低失相干的影响，进而可以最大限度地提高相干点的数量和质量，在低相干区域具有一定的效果[92,94]；另一方面，利用极化分解技术分离不同的散射机制，有效地追溯散射来源，从而获得更加准确的地表形变监测结果[95]。然而，现有研究仍停留在探索阶段，尚未充分挖掘极化信息的潜能。例如，极化数值优化方法仅仅在数学上获得最优解，没有考虑不同地物的极化散射差异和时序变化特征。更重要的是，大多研究使用的是短波数据，长波数据的复杂极化散射特性未得到充分研究。

因此，基于长波数据的时序极化InSAR技术，为穿透植被覆盖实现林下滑坡形变监测提供了可能。如何充分利用长波极化SAR数据开展InSAR滑坡形变透视监测，将是未来的研究热点。但是在未来的应用中，仍有一些问题值得思考，例如，如何准确表征植被区的失相干特征，如何刻画植被区长波极化非对称散射过程，如何考虑极化散射机理进行精确地表相位重建等问题。

4.2 基于单轨数据的InSAR滑坡三维形变监测

如前所述，目前获取滑坡三维形变信息至少需要两个不同轨道的SAR观测结果。然而，在复杂的山区环境下，由于SAR卫星侧视成像导致的几何畸变的存在，在单个滑坡体上往往只能获取升轨或者降轨单个有效InSAR观测，限制了InSAR滑坡三维形变的获取。因此，如何利用单轨InSAR观测值获取滑坡三维形变信息已经成为了目前InSAR滑坡监测的主要瓶颈之一。

在利用单轨InSAR观测获取滑坡多维形变信息方面，目前已有少量的研究开展。文献[85]利用MAI技术获取边坡的方位向形变后，通过忽略东西向形变的方式，与InSAR观测进行融合获取边坡的二维形变。然而，MAI观测很容易受到监测场景的制约，主要表现为该技术对形变的敏感性较低，无法获取形变量级较小的滑坡蠕滑信号，而过大的形变会导致MAI观测失相关。利用POT技术进行监测同样可以获取滑坡的二维形变信息[82]，但是POT技术往往只能在特大形变量级的场景下使用，难以用于开展多数蠕滑型滑坡的常态化三维形变监测工作。因此，如何突破目前的技术瓶颈，改善利用SAR数据获取方位向形变的精度和可靠性，将成为利用单轨SAR数据获取滑坡三维形变信息的解决方案之一。

前文已述利用先验信息进行约束，可以降低InSAR监测三维形变所需要的观测数据量。因此，如何针对滑坡监测场景，进一步充分挖掘其运动本质，将成为扩展InSAR滑坡一维监测到多维监测的关键。目前一些其他形变监测场景中的相关研究可为这一问题的解决提供思路。例如，在矿区形变监测中，可以考虑在水平或近水平煤层开采条件下，矿区地表水平形变与垂直形变之间存在线性关系，即可构建东西向、南北向与垂直向之间的理论约束方程，此时仅需InSAR监测的一维LOS向形变，即可恢复矿区的三维形变场[96]。在地下流体开采监测中，则可以引入弹性半空间模型构建三维形变之间的函数关系，同样可以实现基于单轨InSAR的三维形变监测[97]。但是相较于矿区或开采区地面沉降，滑坡三维形变之间的约束关系目前较为单一，尚不足以支撑滑坡三维形变的获取。利用先验信息约束获取滑坡三维形变信息，还需要对滑坡变形机制机理开展进一步的挖掘，才有可能实现利用单轨SAR数据获取滑坡三维形变。

4.3 面向滑坡监测的InSAR设备载荷研发

星载InSAR技术因其监测范围广、全天时、全天候等优势，已经在滑坡监测方面有许多应用。然而目前在轨卫星的重访周期较长，以最常用的哨兵系列卫星为例，在双星联合观测的情况下，其最短周期也要6 d，难以满足滑坡实时监测预警的需求；其次，高分辨率的星载数据价格昂贵，用其进行时序监测成本高昂。面对这些问题，科研人员考虑将传感器搭载在其他平台(如机载SAR、地基SAR)以补充现有的星载InSAR技术，以实现高时间分辨率、低成本的滑坡监测。

机载SAR观测方向灵活，不受固定飞行轨道的约束；观测时间灵活，可及时地对突发滑坡开展应急测绘；搭载低频传感器可获取密集植被区的滑坡形变。机载SAR目前已在少数滑坡灾害监测中得到了应用。文献[98]利用机载InSAR数据对美国科罗拉多州的Slumgullion滑坡进行监测，证明了InSAR反演滑坡内部变形的动力学特征的能力。文献[99]通过NASA/JPL UAV SAR数据和像素偏移量追踪方法，等解算出加利福尼亚北部海岸山脉的上百个缓慢山体滑坡的三维地表形变，并推断出滑坡体的厚度、体积、几何缩放和摩擦强度等参数。在更低成本和更机动灵活监测的需求下，无人机载微小型SAR应运而生，如MiniSAR、BYUSAR、MFUSAR等，可以实现几千米范围内的亚米级分辨率成像。不过该技术仍然处于研发阶段，尚未像星载SAR一样得到广泛的应用。随着新材料、新部件、新算法的不断涌现，无人机载微小型SAR朝向小体积、低功耗、高分辨率、多模式成像、多极化成像等方向不断发展。

地基SAR则将雷达安置在地面上，时间分辨率可以可达到分钟级，最远监测距离可以达到10 km，形变监测精度理论上可达到亚毫米级，并且可以根据实地情况布设在最适合监测滑坡的位置，从而最大限度地捕捉地表位移信号。现有的地基雷达主要分为两种，一种以线状扫描方式成像，如IBIS-L、LiSA、FastGB-SAR、MIMO-SAR；另一种以弧形方式扫描成像，如ArcSAR、IBIS-ArcSAR、Arc-FMCW-SAR、GPRI-Ⅱ。无论哪一种，距离向的分辨率都是不变的，可以保证在亚米级，方位向分辨率随距离增大而降低，在距离雷达1 km处的方位向分辨率可实现优于10 m。相较于机载SAR，目前地基SAR的设备研发方面较成熟，并在滑坡监测方面有广泛的应用[100-103]。由于安装灵活，雾气、雨天均能进行监测，地基雷达在小范围内应急滑坡监测中发挥了重要作用。然而，现有地基雷达仍然存在不少问题亟待解决。首先，工作波段主要以Ku等高频波段为主，对于植被覆盖等复杂地表难以实施有效监测；然后，尺寸仍然偏大，需要进一步优化使其更加便携；最后，在无电源的观测场景中，供电需要人为参与，难以实现自供电。因此，地基雷达未来将在高穿透、高便携、低功耗等方面不断发展。

4.4 基于深度学习的InSAR滑坡隐患点自动识别

滑坡识别技术是滑坡监测的前提，至关重要，只有准确识别出滑坡体的位置，才能为后续滑坡监测提供条件。滑坡识别是通过对滑坡区域的形态及特征进行研究，确定滑坡的发生区域及分布情况的图像分析技术。传统的滑坡识别方法以人工目视解译为主，极大依赖于专业人员的知识技能与熟练程度，人为主观因素占据较大比重。目前全球共有400多颗在轨地球遥感卫星，面对如此海量的遥感影像数据，若都采用人工解译的方式进行滑坡识别，工作量巨大，难以实施。近年来，基于遥感图像的机器学习方法逐步实现了滑坡隐患的智能化、自动化、快速识别，特别是深度学习，因其拥有强大的图像学习、分析和预测能力，逐渐成为滑坡隐患自动识别的重要工具之一。

SAR影像凭借其全天候、全天时、不受云雾影响等特点被部分学者用于滑坡自动识别的研究当中。文献[104]结合高分二号影像获取的训练样本，构建基于BP神经网络的全极化SAR数据滑坡自动识别模型，实现滑坡体的自动快速识别。文献[105]利用深度学习方法从极化SAR数据中自动识别了滑坡区域。目前已有的研究只能自动识别出已发滑坡，但是对于正在蠕动的滑坡效果欠佳。InSAR作为可探测地表微小形变信号的技术，是识别滑坡蠕动的重要手段，将深度学习和InSAR技术进行结合，将对提高滑坡隐患点自动识别效果起到重要的作用。深度学习的出现为遥感影像自动分类技术提供了新的理论支撑，但是目前基于深度学习与InSAR相结合的研究仍处于初步探索阶段，在将二者进行融合时，仍有一些问题值得思考。首先，对比深度学习与InSAR结合的已有研究来看，其中主要集中在火山的自动识别中[106-109]，不同区域的火山形变在干涉图中具有非常相似的特征结构，而不同的滑坡在干涉图中形变并无明显特征，可能导致深度学习模型无法较好进行学习识别规律。然后，滑坡作为一种潜伏周期长、复杂的地质灾害类型，在短时间内的干涉图中难以发现明显的形变信号，因此对于人工的圈定训练样本范围带来了不小的挑战。最后，深度学习中缺乏针对SAR影像的训练样本集，由此则会产生样本不均衡问题，这个问题不管是深度学习在火山上的研究还是在滑坡中的研究都普遍存在。因此，未来还需要在滑坡特征学习上进一步探索，才能显著提升基于深度学习的InSAR滑坡隐患点自动识别的准确度。

4.5 形变大数据驱动的InSAR滑坡灾害预报预警

观测和分析各种滑坡发生的前兆现象，开展预报预警工作，可以在一定程度上避免滑坡灾害对人民生命财产安全的侵蚀。传统方法，如逆速法线性逼近、对数速率和对数加速度技术等数值模拟方法受测量误差、环境噪声及建模假设等因素的干扰效果明显，表现为模型抗差性能低、迁移应用能力弱。随着星载SAR技术的不断发展，目前已有多种不同波段、不同模式、不同分辨率的SAR卫星在轨服务，长期的稳定运行已累计获取了全球每个角落的海量存储SAR数据。深度学习等大数据挖掘和分析手段的出现，为利用InSAR数据开展形变大数据驱动的滑坡灾害预报预警提供了机会。

目前，在结合大数据挖掘和传统形变数据服务于灾害预报预警方面，已经有学者开展了探索性的研究。文献[110]提出一种基于双移动平均线(DMA)和长短期记忆(LSTM)网络的数据驱动时序变形预测方法，结果表明该耦合模型预测三峡库区白水河滑坡形变能够完成更高精度的形变预测，可用于早期预警系统。文献[111]将深度学习模型应用于地震诱发的山体滑坡形变预测中，结果表明深度学习模型在准确性方面优于朴素贝叶斯、逻辑回归、支持向量机及随机森林4种机器学习算法。文献[112]利用LSTM模型对白水河滑坡进行了单因素和多因素预测，结果表明具有多因素数据辅助的滑坡预测性能优于单因素深度学习模型的预测性能。文献[113—114]联合自适应噪声的完全集成经验模态分解和嵌入注意力机制的LSTM模型对由GNSS获取的滑坡形变数据实现了预测任务，并达到了较高精度，表明了深度学习方法对于滑坡位移预测存在巨大潜力。这些研究为开展形变大数据驱动的InSAR滑坡灾害预报预警提供了有益参考，仍有一些难点需要考虑：①滑坡相较于地面沉降、建筑膨胀大坝变形而言，其形变趋势的影响因素繁多且机理复杂，导致参与预测数据的选择和其最佳组织方式仍需进一步探究；②突降暴雨、人为因素或地震等大型灾害都有可能引起滑坡的瞬间垮塌，但由于此等现象并不多见，能够用于训练的突发性滑坡数据库并不丰富，这是基于深度学习的滑坡预测的又一瓶颈；③目前针对其他应用的形变预测都是基于单点的时序预测，然而滑坡形变势必还需要考虑空间关系，因此滑坡预测向空间维联合时间维方向上的拓展是亟待突破的一大难点。

5 结 论

InSAR技术为突破传统滑坡监测手段中所面临的瓶颈问题提供了新的契机。近年来，星载SAR数据日趋丰富，极大地推进了InSAR技术在滑坡监测领域上的应用和研究。针对InSAR滑坡监测过程中所面临的失相关噪声严重、大气误差扰动和无法获取三维形变等问题，相关学者提出了一系列的解决思路，极大丰富了InSAR技术理论。在可预见的未来，又将有一批新的SAR卫星入轨运行。新的数据将带来更短的重返周期、更高的空间分辨率及更丰富的极化方式。这一系列新特征不仅有助于InSAR技术在开展滑坡变形监测时突破植被覆盖导致的监测精度低、几何畸变导致的监测维度单一等问题，而且有望促进InSAR技术在机载/地基InSAR设备载荷研发、滑坡隐患点自动识别、滑坡灾害预报预警等方面的应用，从而进一步提高我国防灾减灾的能力。