1. 引言
2022年9月5日，四川省甘孜州泸定县发生Ms6.8级地震，震源深度16 km，最高烈度达到Ⅸ度，震中位于泸定县磨西镇(北纬29.59˚，东经102.08˚)。地震发生后，强烈震动诱发了大量次生滑坡灾害[1]。这些地震滑坡不仅在短期内造成严重的人员伤亡和基础设施破坏，还可能在震后持续一段时间内发生蠕动、局部崩塌甚至再次滑坡，给灾后恢复与重建工作带来长期隐患。因此，开展震后滑坡的持续形变监测，不仅有助于及时掌握滑坡区域的稳定性变化趋势，也为灾害预警、风险评估与应急响应提供重要的科学依据[2]。
合成孔径雷达干涉测量技术(Interferometric Synthetic Aperture Radar, InSAR)凭借其毫米级形变敏感度、广域覆盖能力及全天候、全天时观测优势，已成为震后地质灾害监测的核心手段之一[3]。其中，小基线集时序InSAR (Small Baseline Subset InSAR, SBAS-InSAR)方法通过优化干涉对的组合策略，有效降低了时空失相干带来的影响，可实现滑坡体长时间序列形变场的稳健反演，是震后滑坡持续监测的重要技术路径[4]。然而，传统SBAS-InSAR处理流程存在数据量大、处理效率低[5]，以及解缠误差与噪声干扰显著等问题，难以满足震后应急监测对时效性和精度的双重要求。
近年来，云计算架构与开源科学计算工具的协同发展，为InSAR技术的快速部署与高效应用提供了技术突破。由美国阿拉斯加卫星中心(Alaska Satellite Facility, ASF)研发的HyP3 (Hybrid Pluggable Processing Pipeline) [6]平台，集成了标准化预处理模块与弹性计算资源，可实现Sentinel-1数据的自动化辐射校正与干涉对生成。结合由加州理工学院喷气推进实验室(Jet Propulsion Laboratory, JPL)开发的MintPy (Miami InSAR Time-series software in Python) [7]时序分析框架，能够构建从原始SAR数据到地表形变产品的全流程处理体系。此外，针对时序InSAR处理中关键误差源，如相位解缠误差、噪声影响和参考点选取问题，仍需进一步优化策略，以提升震后滑坡形变监测的精度与可靠性。
基于此，本文以2022年泸定Ms6.8级地震为研究背景，利用HyP3-MintPy联合处理流程，对震区内典型滑坡区域开展震后持续形变监测。通过优化参考点选取、水体掩膜处理及参考日期筛选等关键步骤，提升InSAR解缠质量与形变反演精度，揭示滑坡震后的时空演化特征，为震后地质灾害监测与风险管理提供技术支撑与科学依据。
2. 研究区与数据源
2.1. 研究区概况
研究区范围包括此次地震烈度为Ⅸ度的全部区域和烈度为Ⅷ度的部分区域(图1(a))，涉及泸定、石棉两个县(图1(b))。研究区内主要覆盖大渡河、燕子沟、海螺沟以及田湾河流域，整体地形呈西高东低，为典型的高山峡谷地貌。区内地质构造复杂，发育多条断裂，历史地震活动强烈，作为本次地震的发震断裂—鲜水河断裂贯穿整个研究区。
Figure 1. Overview of the study area
图1. 研究区概况
2.2. HyP3 SAR数据
HyP3 (Hyperspectral Processing Pipeline Prototype)是由阿拉斯加卫星设施(ASF, Alaska Satellite Facility)提供的一项基于云计算的在线SAR数据处理服务，旨在为用户快速、便捷地生成高质量的分析就绪SAR数据产品(图2) [6]。
Figure 2. HyP3 data products ((a) original interferogram, (b) unwrapped interferogram)
图2. HyP3数据产品((a) 原始干涉图，(b) 解缠后的干涉图)
本文设置其多视系数为10:2，时间范围为2022年9月至2023年6月的升、降轨HyP3在线SAR数据，用于震后滑坡的形变监测。数据基本信息如表1所示。
Table 1. System resulting data of standard experiment
表1. 标准试验系统结果数据
参数 |
HyP3 |
7轨道方向 |
升轨/降轨 |
轨道号 |
26/135 |
幅号 |
93 (88)/493 |
时间基线 |
36 days |
波段 |
C |
极化方式 |
VV + VH |
影像数量(景) |
22/21 |
干涉对数量 |
57/54 |
影像获取时间段 |
2022/9/1至2023/6/1 |
3. 研究方法
3.1. SBAS-InSAR技术
小基线集技术(Small Baseline Subset Interferometric Synthetic Aperture Radar, SBAS-InSAR) [8]是差分雷达干涉测量(DInSAR)的一种改进方法，用于高精度提取地表形变信息。尽管小基线集技术在地表形变监测方面具有重要优势，但在实际应用中仍然存在一定的局限性，主要体现在以下三个方面。
1) 数据需求大，计算复杂度高。SBAS技术依赖多时相SAR影像构建干涉网络，以确保时间序列分析的精度和稳定性。因此，需要处理大量数据，涉及干涉图生成和时间序列分析等复杂计算，计算资源消耗大，处理时间长，对硬件设备要求较高。2) 解缠误差影响。相位解缠是时序InSAR处理中保证形变监测精度的关键步骤。然而，受噪声、相干性退化或复杂地形影响，可能出现相位跳变、错误解缠或累积误差，进而影响最终的形变计算结果。3) 噪声影响。时间序列InSAR处理中，噪声主要来源于大气延迟、轨道误差和热噪声等因素。其中，参考日期的选择对噪声积累影响显著。若参考日期处于高噪声水平，可能放大系统误差，降低形变估计的稳定性。
3.2. 快速形变监测流程及优化
上述问题在震后快速形变监测中尤其突出，为了应对这些局限性，本研究采用HyP3-MintPy数据处理流程，结合HyP3的自动化、高效预处理能力和MintPy的时序分析功能，在数据量管理、相干性保持、误差校正及计算效率等方面提供了解决方案。完成了研究区的震后滑坡形变监测，其流程如图3所示。
Figure 3. Post-earthquake deformation rapid acquisition process
图3. 震后形变快速获取流程
在HyP3的处理流程中，首先根据时间基线与空间基线自动选择合适的SAR图像对，并设置相应的滤波参数以提高干涉图的质量。随后将数据转换为GAMMA格式，并进行重叠Burst区域的重新拼接与处理，以保证干涉对的空间连续性。接着利用Copernicus全球DEM进行干涉图生成和共配准操作，进一步执行相位解缠，并进行选代配准优化，以提升几何一致性与相干性。最后，完成地理编码和产品打包，生成标准化的干涉产品，便于后续MintPy软件进行时序形变分析。为了进一步提高形变监测的精度与可靠性，本研究对关键误差源和噪声源进行了优化，重点改进以下关键步骤。
1) 解缠误差校正：参考点与水体掩膜的应用参考点的选取是时间序列形变分析和相位解缠的关键，能为干涉图提供稳定的相位基准。解缠从参考点向外扩展，距离越远，误差越容易积累，可能导致解缠错误或形变偏差。
此外，部分水体像素可能因满足相干阈值被错误解缠，尤其是在沿海或大内陆水体区，严重影响结果准确性。为避免此类问题，应使用水体掩膜排除水体区域，防止误差传播。因此，优化参考点策略并结合水体掩膜，可有效提升解缠质量与形变监测的可靠性。本研究中，通过闭合相位(Closure Phase)方法检测干涉图中的解缠误差[7]。通过计算干涉三角形(由三个影像组合成的干涉对)的闭合相位，可以检查是否存在解缠错误，闭合相位定义如下：
(1)
为三角形闭合相位，
，
，
分别为三个影像组合的解缠干涉相位。理论上，如果干涉图中没有解缠错误，闭合相位 应为零。然而，由于解缠过程中可能引入误差，闭合相位可能会偏离零，这部分误差以整数形式表示(称为整数模糊)。
(2)
将
的值重新映射到[−π, π]区间。等号右侧为整数模糊的计算方法，表示闭合相位中偏离零的整数倍数。
(3)
一个无误差的干涉图网络应满足
。如果存在解缠错误，会导致某些三角形的
。因此，通过计算所有三角形中非零整数模糊的数量
，可以用来评估解缠误差。
2) 基于残差相位评估的最优参考日期选取。计算每次采集的残差相位时间序列的均方根(RMS)，用于表征InSAR时间序列的噪声水平。计算中值绝对偏差(MAD)值，如果SAR采集的RMS值大于预定义的截止值，则将其标记为噪声。RMS值最小的噪声SAR采集可以解释为大气湍流最小的日期，选择并将其作为参考日期[7]。
(4)
式中，
，是网络反演期间基于时间相干性选择的可靠像素的集合，
表示雷达中心波长，
表示在i时刻的冗余相位。使用中值绝对偏差(MAD)作为基准，设定一个阈值(通常为3∙MAD，对应99.7%的置信度)，用于判断噪声水平。
(5)
其中
为中值，
为所有冗余相位的中值。如果某次采集的RMS值超过预定义的阈值，则标记该采集为“噪声采集”，将噪声最小的SAR数据采集设置为参考日期，以减少系统误差的累积。
3.3. MintPy参数设置
为高效开展HyP3数据的时序InSAR分析，本文采用MintPy软件，并结合研究区特点对各关键处理步骤进行了参数配置与优化。在干涉图网络构建阶段，为剔除低质量干涉图对反演结果的干扰，设置最小空间相干性阈值为0.3，并启用基于最小生成树(MST)的网络优化策略，以保障网络的连通性与稀疏性。在对流层延迟校正方面，选用PyAPS模块并调用ERA5全球气象再分析数据对大气分层延迟进行改正，该方法在四川地区复杂水汽条件下能有效提升相位精度。
时序反演过程中，采用最小范数约束的加权最小二乘法，设置反演权重为均匀分布，并引入最小速度正则化以增强结果的稳定性。同时，设定时间相干性阈值为0.1，用于筛除不稳定像元。针对高山峡谷地形中DEM误差导致的系统性相位残差，启用了地形残差校正模块，并基于整体几何信息进行拟合，以避免像素级建模带来的过拟合问题。此外，还启用了残差均方根(RMS)分析及时间参考日期优化功能，自动选取误差最小的参考时间点，从而进一步提升数据质量与时间序列一致性。
3.4. 震后典型滑坡选取
为了更系统地分析不同类型滑坡的震后形变规律，并为灾害监测与风险评估提供科学依据，本研究筛选出H1和H2两个典型滑坡区进行重点分析。
H1典型滑坡区位于大渡河沿岸(图4)，震前地表整体稳定，但震后发生大面积滑坡，导致地貌显著变化。滑坡区西南侧紧邻什月坪村，对周边居民区构成潜在安全隐患。此外，该滑坡区域的地表失稳可能进一步影响大渡河河道的稳定性，增加沿岸冲刷和次生灾害的风险，需要持续监测其形变演化特征。
Figure 4. Location of H1 typical landslide area
图4. H1典型滑坡区地理位置
H2典型滑坡区主要分布于陡峭山坡，震后滑坡范围明显扩大，地表裸露岩土增多，反映出强烈的地表扰动和地质不稳定性。该滑坡区域坡体高差较大，且其北侧紧邻湾东村，滑坡的持续形变可能对该区域的居民安全和基础设施稳定性构成潜在威胁，需加强监测与防护措施(图5)。
Figure 5. Location of H2 typical landslide area
图5. H2典型滑坡区地理位置
4. 实验结果与分析
4.1. 解缠误差改正与最优日期选取结果
从图6可见，高相干性区域(如裸地或稳定建筑物)通常被选作参考点，以降低解缠误差。然而，单纯选择高相干性区域作为参考点并不能保证整个研究区的解缠误差均匀降低。如图所示，尽管在高相干性区域选择参考点能够有效减少该区域附近的解缠误差，但研究区其他低相干性区域的误差仍然较高(图6(b))。研究区内地震滑坡主要分布在大渡河沿岸及沿岸附近，因此，为了提高解缠精度并减少形变估计的不确定性，参考点应优先选取在沿岸稳定区域。
Figure 6. (a) Spatial coherence map of the study area, ((b) (c)) Tint maps with different reference points in the study area
图6. (a)研究区空间相干性图，((b) (c))研究区不同参考点下的Tint图
沿岸平稳区域具有较高的相干性且地形变化较小，有助于降低解缠误差，从而提高形变监测精度。图6(b)展示了参考点选取在河沿岸时的误差分布，可见该策略有效减少了解缠误差的空间不均匀性。
如图所示，对比图7中使用水体掩膜前(a)和使用水体掩膜后(b)的解缠误差分布，可以直观地看到水体掩膜对解缠精度的影响。掩膜前，由于部分水面像素仍满足相干阈值，导致水体区域出现错误的相位解缠，尤其在河流沿岸，误差累积明显，表现为大量红色区域。
掩膜后，水体像素被成功去除，避免了错误相位传播，从而减少了解缠误差。特别是在M和N区域(图7(a)和图7(b))，可见掩膜后河流沿岸的解缠错误显著减少，蓝色区域增多，表明解缠精度和形变估计的可靠性均得到提升。这一结果表明，在沿海或大面积水体覆盖区域，应用水体掩膜可有效减少解缠误差，提高形变监测精度。
实验中计算了升轨和降轨干涉图中的所有RMS，如图8所示，最终将升轨参考日期设置为20230318，降轨参考日期设置为20230217。
Figure 7. Comparison of phase unwrapping error distribution. (a) Without water mask application, (b) With water mask application
图7. 研究解缠误差分布对比。(a) 未使用水体掩膜，(b) 使用水体掩膜
Figure 8. RMS of residual phase per acquisition
图8. 升降轨道每次采集的残差相位均方根图
4.2. 典型滑坡区形变监测结果
为了进一步分析其震后形变特征，获取了从2022年9月至2023年6月期间覆盖上述两个典型滑坡区的升降轨HyP3在线SAR数据，其多视系数为10:2，时间基线为60天，空间基线则严格控制在200米以内。基于这些参数，构建74对升轨干涉像对和54对降轨干涉像对。图9、图10展示了所有干涉像对的干涉基线连接情况。
Figure 9. Interferogram network of HyP3 ascending orbit data
图9. 升轨HyP3数据干涉图网络
Figure 10. Interferogram network of HyP3 descending orbit data
图10. 降轨HyP3数据干涉图网络
经过MintPy程序处理，获取了H1滑坡区在2022年9月至2023年6月震后期间的形变速率，结果如图11(c)所示。图11(a)和图11(b)分别展示了H1滑坡区在地震前后的光学影像图。
地震前，H1滑坡区植被覆盖良好，地表稳定(图11(a))。震后，该区域发生大规模滑坡，坡面大面积裸露，特别是南侧靠近河岸区域，受侵蚀影响明显(图11(b))。滑坡范围内发育多条冲沟，其中左下角区域形变显著，需重点关注其后续演化。图11(c)显示了H1滑坡区的形变速率，最大视线向形变速率达−13.2 cm/年，主要集中在滑坡体南部，表明该区域滑坡活动剧烈，存在持续蠕动或失稳风险。其中，黄色(正值)表示地表抬升，主要分布于高海拔区域；紫色(负值)表示沉降，集中于坡脚区域，反映滑坡体的整体下滑趋势。为了验证形变速率的可靠性，实验中计算其标准差(图11(d))，H1滑坡区的形变速率标准差约为0.89 cm/年。
Figure 11. H1 typical landslide area (a) Pre-earthquake optical/SAR image (b) Co-seismic image (c) Post-earthquake deformation velocity (d) Standard deviation of deformation velocity
图11. H1典型滑坡区地震前(a)后(b)影像图、形变速率图(c)、形变速率标准差图(d)
图12显示，H1滑坡区在2022年9月至2023年5月期间形变量持续增加，呈现“初始稳定–加剧沉降–持续蠕动”的演化趋势。前期变形较小，156天后沉降范围扩大，252天趋于稳定但仍有较大变形，表明滑坡体可能进入长期蠕动阶段。最大形变集中于滑坡体中央，132天后形变量显著上升，为滑坡活动的关键时期，需持续关注其演化趋势。
Figure 12. H1 typical landslide area post-earthquake time-series cumulative deformation map
图12. H1典型滑坡区震后时序累计形变量图
特征点P1选自2023-05-17的累积形变量图(图13)，位于H1滑坡区的典型变形区域。其时间序列形变量显示，该点在2022年9月至2023年5月持续沉降，最大累积位移达−11.87 cm，整体呈稳定下沉趋势。初期形变量较小，随后形变速率逐步增大，并在后期趋于缓和。
Figure 13. Time-series cumulative displacement of characteristic point P1 in H1 typical landslide area
图13. H1典型滑坡区P1特征点时序累计形变量图
时间序列拟合结果表明，P1形变速率为−10.71 cm/年，表现为持续沉降，加速度为18.12 cm/年2，反映出形变呈非线性加速趋势。这表明P1可能仍处于震后蠕变或滑坡调整阶段，需进一步关注其长期稳定性。
H2典型滑坡区的结果如图14(c)所示。图14(a)和图14(b)分别展示了H2滑坡区在地震前后的光学影像图。从地震前影像(图14(a))可见，H2滑坡区整体坡度较陡，地表主要为植被覆盖，未见明显滑坡迹象。地震后(图14(b))，该区域发生了大规模滑坡，坡面出现大范围裸露岩土，尤其是在滑坡体中央及坡脚区域，表明滑坡破坏较为严重，且坡体高差较大，可能存在进一步失稳的风险。
Figure 14. H2 typical landslide area (a) Pre-earthquake optical/SAR image (b) Co-seismic image (c) Post-earthquake deformation velocity
图14. H2典型滑坡区地震前(a)后(b)影像图、形变速率图(c)
图14(c)为H2滑坡区的形变速率分布，展现了该区域在震后时序InSAR监测下的视线向形变速率。最大沉降速率达−16.7 cm/年，主要集中于滑坡体中部及下部，说明该区域地表活动最为剧烈，存在持续蠕动或进一步破坏的可能。整体来看，H2滑坡区在地震后经历了大规模失稳和形变发展，需要重点关注其后续演化趋势及潜在的灾害风险。H2滑坡区的形变速率标准差约为1.21 cm/年(图14(d))。
图15展示了H2典型滑坡区在2022年9月2日至2023年5月24日期间的时序累计形变量，直观地反映了该区域的震后形变演化过程。
Figure 15. H2 typical landslide area post-earthquake time-series cumulative deformation map
图15. H2典型滑坡区震后时序累计形变量图
从时间演变来看，滑坡区的形变量随时间逐步增加，整体呈现初始缓慢–加速沉降–形变趋稳的趋势。在前72天，滑坡区整体变形较小，仅部分区域出现轻微沉降；132天后，滑坡体沉降范围明显扩大，表明坡体仍处于调整阶段；至264天，形变量趋于稳定，但仍维持较大沉降幅度，可能进入长期蠕变阶段。从空间分布来看，形变量的最大区域集中在滑坡体中央，表明该区域的地表活动最为剧烈。图中P2点(黄色三角标记)代表特征监测点，其在264天的时序监测中表现出持续沉降趋势，需重点关注其后续演化及潜在灾害风险。
P2特征点位于H2滑坡区的典型变形区域，其时间序列形变量(图16)展现了该点在2022年9月至2023年5月期间的形变演化趋势。P2点在震后持续发生沉降，最大累积形变量达到约−12.29 cm，整体呈现稳定的负向形变趋势。形变量初期较小，随后形变速率逐渐增大，并在后期趋于平缓，表明滑坡活动可能经历了一个由加速变形到逐步稳定的过程。
从P2点的时间序列拟合曲线来看，其形变速率为−18.33 cm/年，显示出显著的沉降趋势。加速度为3.13 cm/年2，表明形变速率在震后有所变化，可能受到滑坡内部应力调整或外部环境因素的影响。这一特征表明，P2所在区域仍可能处于震后调整期，需进一步关注其长期演化趋势和潜在灾害风险。
Figure 16. Time-series cumulative displacement of characteristic point P2 in H2 typical landslide area
图16. H2典型滑坡区P2特征点时序累计形变量图
5. 结论
本文以2022年四川泸定Ms6.8地震为研究区，基于HyP3平台获取的在线SAR数据，结合MintPy时序InSAR分析工具，对震区典型滑坡区域开展了震后持续形变监测研究。主要结论如下：
1) 本研究构建了基于HyP3-MintPy的震后滑坡监测流程，在确保监测精度的同时显著提升了处理效率与自动化水平，可有效满足震后应急监测的时效性需求。
2) 通过闭合相位法识别解缠误差，结合水体掩膜与高相干区域优化参考点选取，降低了解缠误差的空间不均匀性。同时，基于残差相位RMS与MAD方法选取噪声较低的参考影像，提高了形变反演的稳定性与准确性。
3) H1与H2两个典型滑坡区在震后均表现出持续沉降和非线性形变特征，最大沉降速率分别达−13.2 cm/年和−16.7 cm/年，局部区域存在持续蠕动与失稳风险，需持续监测与关注。