森林下隐藏滑坡的测绘与早期预警：以香港大屿山为例

《Remote Sensing of Environment》：Mapping and early warning of hidden landslides under forests: A case in Lantau, Hong Kong

香港陡坡林地滑坡的InSAR检测方法优化与降雨触发机制研究

　　在热带和亚热带地区，极端降雨的增强导致了广泛的滑坡灾害，这在世界人口密度最高的香港尤为显著。香港的地形陡峭且多森林覆盖，这使得传统的合成孔径雷达干涉测量（InSAR）难以有效监测滑坡，因为山坡的滑动通常较小且隐藏在密集的树冠之下。本研究开发了一种新的检测框架，结合了两种方法：一是同质化幅度-相位细化（HARMIE），通过增强局部相位变化来保留微小位移；二是基于相位梯度的检测方法，将山坡响应与极端降雨联系起来。模拟和实际数据实验表明，HARMIE在保留局部相位细节和幅度方面优于传统方法。利用高分辨率的升轨和降轨Lutan-1（LT-1）SAR数据集（2023年7月至2024年10月），该框架成功地映射了由2023年10月极端降雨触发的广泛山坡滑动，识别率比基于幅度同质化的检测方法提高了27%，在捕捉狭窄的滑坡（如约10米宽）方面也表现出显著的改进。此外，研究还识别出了十处活跃的、隐藏在森林中的滑坡。

滑坡在地理上具有多样性，其成因通常涉及多种因素。例如，长期的前期降雨可能会引发坡地的季节性缓慢蠕动，而对于某些坡地，极端降雨可能会加速其不稳定性。这些现象表明，滑坡的发生不仅受到降雨的直接冲击，还可能受到前期降雨条件的影响。本研究首次利用L波段SAR技术对香港的森林覆盖的小规模滑坡进行InSAR映射，为森林覆盖的山地地形滑坡不稳定性的研究提供了新的视角，并推动了此类地区的滑坡预警系统的发展。

滑坡对香港的地理环境构成了持续的威胁。香港的地形陡峭，地质材料普遍易受侵蚀，同时受到强烈的热带季风降雨的影响。这种地理和气候条件使得滑坡的发生具有一定的频率和强度，尤其是在极端降雨事件中，这种现象更为显著。例如，2008年6月7日和2023年9月7日的两次极端降雨事件分别引发了362和214次滑坡，造成了人员伤亡、基础设施损毁和巨大的经济损失。因此，主动的滑坡地图绘制和对山坡运动的早期检测对于香港这一高密度建筑和人口的城市至关重要。

InSAR技术已被广泛应用于滑坡地图绘制，并在地质灾害预警中展现出其优越性，特别是在检测活跃的坡地不稳定性方面。然而，在植被覆盖的区域，由于显著的去相关效应，InSAR的相干信号往往受到干扰。通常通过同质化多视和相位链接来减少噪声，这将同质化像素聚集在一起以估计最佳的相位序列。这些方法有助于减少干涉条纹噪声，并有效提取地物参数。然而，传统的时序后处理方法常常受到相位解缠误差和显著的大气效应的阻碍。作为一种替代方案，相位梯度方法避免了相位解缠的需求，通过计算相邻像素之间的空间差异，可以抑制长波长的大气效应并增强变形边界定义。在本研究中，所提出的InSAR方法为揭示香港的滑坡灾害提供了有希望的解决方案。

然而，香港特有的滑坡特征和复杂的环境条件显著增加了InSAR在跟踪山坡运动中的应用难度。香港的山地通常被浅层的坡积物覆盖，这些材料由于松散和渗透性较强，容易在强烈或极端降雨下发生瞬时饱和，引发强烈的水力响应，并使山坡更容易发生浅层失稳。大多数记录的坡地滑动发生在这些薄层的坡积物上，经常表现为泥石流，这是该地区常见的滑坡模式。因此，香港的滑坡通常以浅层、通道式滑动为主要特征，其源宽通常小于20米，滑动深度通常小于3米。此外，由于香港的山地普遍被森林覆盖，植被导致的严重去相关使得滑动信号难以检测。无论由于密集树冠的视觉遮挡还是InSAR中的相干性丧失，小规模滑坡的特征，如微小的地面变形或滑坡阶地，通常都隐藏在观测中。这一挑战限制了InSAR技术充分发挥其技术优势的能力。

在植被覆盖区域检测小规模滑坡信号时，当前的InSAR方法面临以下困难。传统的基于幅度同质性的方法不考虑相位变化，因此难以选择同质化像素以保留与微小山坡运动相关的局部相位。在标准的圆对称复高斯（CCG）分布假设下，分布式散射体（DS）的复协方差矩阵被视为表征单视复数（SLC）SAR向量的充分统计量，并用于推断最佳相位历史。实际上，幅度同质化方法通常从基于幅度同质性的像素邻域中推导出复协方差矩阵的极大似然估计（MLE）。然后，复协方差矩阵的估计值 [Math Processing Error]$\stackrel{?}{\Gamma }$ 被视为足够近似于相干矩阵，用于推断最佳相位序列。然而，对于小规模滑坡监测，当局部相位不一致时，多视干涉相位在 [Math Processing Error]$\stackrel{?}{\Sigma }$ 的参数中可能被低估。基于相位的方法虽然可以更好地保留局部相位异质性，但会导致所选像素包含分散的幅度样本，同样会导致 [Math Processing Error]$\stackrel{?}{\Gamma }$ 的低估。

为了评估上述低估现象，我们进行了蒙特卡洛模拟，考虑了三种情况：(1) 完全同质化集，其中幅度和相位样本是从独立同分布（i.i.d.）概率密度函数（p.d.f.）中抽取的，(2) 具有30%相位异质性，(3) 具有30%幅度异质性。具体模拟参数和结果分析见补充材料（第S1.2节）。如图3所示，完全同质化集的相干估计在样本数量增加时表现出轻微的正偏差（灰色阴影在黑色实线之上），这与之前的研究结果一致。相比之下，相位或幅度异质性的引入导致了更大的偏差和相干估计的潜在低估，随着样本数量的增加，这种偏差更为显著。这些发现强调了对同质化像素选择的严格要求，以减少相干估计偏差。

在小规模滑坡监测中，局部相位不一致的存在使得多视干涉相位的估计容易受到低估。基于幅度的方法选择幅度同质化的像素，但可能包括方向各异的SLC向量，这会降低平均向量的幅度，从而削弱 [Math Processing Error]$\stackrel{?}{\Gamma }$ 的幅度。基于相位的方法则没有对幅度同质性施加限制，这导致所选像素包含分散的幅度样本，同样导致 [Math Processing Error]$\stackrel{?}{\Gamma }$ 的低估。因此，必须采用更加严谨的同质化像素选择方法，以减少相干估计偏差。

本研究通过模拟实验进一步验证了这一点，我们模拟了两组32幅图像（每组700×700像素），分别包含无噪声的相位和幅度数据。相位模拟的图像4、14和30如图5所示（第一列）。我们为图像分配了一个高斯椭圆相位作为变形信号，其变化率约为1.3厘米/年，模拟为指数变化并具有0.5厘米的波动幅度。在图像中加入了规则的几何形状的DEM误差相位。为了增强相位变化，每组四个较大的相位矩形被一条相位不连续线分割成两个相邻块，具有不同的值。幅度图像则采用类似的布局，但在变形区域包含密集排列的对比矩形，以强调异质性，如图5所示。最终，基于固定相干矩阵生成了相同大小的SLC噪声图像堆栈，并将相位和幅度图像耦合在一起。完整的模拟设置见补充材料（第S1.4节）。

在获得模拟SLC堆栈的同质化像素集 [Math Processing Error]${\Omega }_{\mathrm{AP}}$ 后，我们需要找到一个最佳的权重 [Math Processing Error]$D_0^{\left(0\right)}$，以估计 [Math Processing Error]$\stackrel{?}{\Sigma }$。假设最小权重为 [Math Processing Error]$D_{0,\min }$，并约束 [Math Processing Error]$D_0^{\left(0\right)}$ 小于或等于 [Math Processing Error]$D_{0,\min }$。引入一个比例因子 [Math Processing Error]$\lambda$ 来缩放权重 [Math Processing Error]$D_0^{\left(0\right)}$ = [Math Processing Error]$\lambda D_{0,\min }$。通过调整 [Math Processing Error]$\lambda$ 的值，我们可以评估相干矩阵偏差的Frobenius范数 [Math Processing Error]${\left(\stackrel{?}{\Gamma }?\stackrel{～}{\Gamma }\right)}_F$。如图5所示，随着 [Math Processing Error]$\lambda$ 接近1， [Math Processing Error]${\left(\stackrel{?}{\Gamma }?\stackrel{～}{\Gamma }\right)}_F$ 逐渐减少并趋于收敛。因此，我们采用 [Math Processing Error]$D_0^{\left(0\right)}$ = [Math Processing Error]$D_{0,\min }$ 用于HARMIE中 [Math Processing Error]$\stackrel{?}{\Gamma }$ 的估计。

为了评估HARMIE的性能，我们将四种方法——DespecKS、KUiper的相位同质性测试（KUTPH）、准同质相位选择（Q-HPS）和HARMIE——用于在15×15的正方形窗口内对同质化像素进行聚类。Q-HPS遵循HARMIE的第二阶段，但不包含幅度同质性约束。四种方法的典型相位结果如图5所示，对应的同质化像素数量如图5d所示。进一步评估在四个空间异质性子区域（图6a–d）中进行。

总体而言，四种方法的相位估计结果通常反映了原始相位模式（图5）。然而，在图6的放大区域中，DespecKS和KUTPH的结果中出现了明显的模糊伪影。具体而言，仅依赖幅度同质性进行像素选择会无意中将幅度模式的痕迹印入相位估计结果中（见图6a的第一行）。这也会导致矩形内相位不连续的过度平滑，产生对尖锐相位边界的模糊。相比之下，Q-HPS和HARMIE方法不仅有效抑制了干涉条纹噪声，还准确重建了干涉条纹并突出了尖锐的相位不连续区域，展示了它们在处理局部相位变化方面的优势。每种方法获得的同质化像素数量，如图5d所示，也反映了它们在管理相位和幅度异质性方面的能力差异。

为了更定量地比较性能，我们提取了相位/幅度不连续区域周围的缓冲区（7像素宽），评估了相干偏差 [Math Processing Error]${\left(\stackrel{?}{\Gamma }?\stackrel{～}{\Gamma }\right)}_F$ 和相位估计误差的L2范数 [Math Processing Error]${\left(\stackrel{?}{?}??\right)}_2$。这些指标之间的关系如图5e所示，揭示了不同的性能特征。DespecKS和KUTPH对应的散点图显示 [Math Processing Error]$\stackrel{?}{\Gamma }$ 和 [Math Processing Error]$\stackrel{?}{?}$ 的估计误差之间存在弱或不明显的依赖关系，即使 [Math Processing Error]$\stackrel{?}{\Gamma }$ 的偏差相对较小，相位估计的不确定性依然显著。相比之下，Q-HPS和HARMIE的结果显示了一种明确的正相关关系，随着 [Math Processing Error]${\left(\stackrel{?}{\Gamma }?\stackrel{～}{\Gamma }\right)}_F$ 的减少， [Math Processing Error]$\stackrel{?}{?}$ 的不确定性也随之降低（见图5e第一行的插图）。值得注意的是，HARMIE方法在 [Math Processing Error]$\stackrel{?}{\Gamma }$ 和 [Math Processing Error]$\stackrel{?}{?}$ 偏差之间表现出更强的正相关性，相较于Q-HPS。图S3展示了图5e中散点的密度分布。HARMIE方法得出的 [Math Processing Error]$\stackrel{?}{\Gamma }$ 和 [Math Processing Error]$\stackrel{?}{?}$ 的偏差在低偏差区域更为密集。对于Q-HPS，散点密度分裂成两个集群，这可能由于幅度异质性对相干估计的影响所致。这些发现表明，同时施加幅度和相位同质性约束可以增强相干矩阵估计与相位重建之间的理论依赖关系。

我们的模拟实验确认了相位和幅度空间异质性对最终相位估计的不利影响。在比较的方法中，Q-HPS在捕捉空间异质性相位特征方面表现出色，而HARMIE在估计 [Math Processing Error]$\stackrel{?}{\Gamma }$ 和 [Math Processing Error]$\stackrel{?}{?}$ 方面达到了最准确的结果，这是由于其更严格和可靠的同质化像素选择，尽管所选像素数量较少。此外，HARMIE在抑制干涉条纹和保留相位异质性方面表现出很强的可靠性和有效性。

鉴于LT-1 SAR SLC在Lantau Island上的高空间分辨率——通常包括约10,000个像素在距离和方位角方向上，分辨率为1.7米和1.5米——计算效率是一个关键考虑因素。因此，采用一种更高效的参数形式的GLRT方法来构建幅度同质性像素集 [Math Processing Error]${\Omega }_A$。在完成第3.1.2节中描述的步骤后，应用HARMIE方法的剩余步骤来生成从噪声LT-1干涉图中优化的相位产品。为了评估HARMIE在现实世界中的性能，我们还使用 [Math Processing Error]${\Omega }_A$ 像素集作为比较基准生成了一个替代优化结果。

图7a和e展示了使用GLRT方法（[Math Processing Error]${\Omega }_A$）和HARMIE方法（[Math Processing Error]${\Omega }_{\mathrm{AP}}$ 和 [Math Processing Error]$W$）优化的下降和上升轨道LT-1干涉图。图7b-d和f-h展示了HARMIE方法的优化干涉相位的放大视图，其中列I和II分别对应基于幅度同质性（[Math Processing Error]${\Omega }_A$）和HARMIE（[Math Processing Error]${\Omega }_{\mathrm{AP}}$ 和 [Math Processing Error]$W$）方法的结果。根据第2.3节所述，从两个光学图像中划定的滑坡边界被用作评估相位异质性恢复的参考，并叠加在地理编码的干涉图上。确认的滑坡范围由实线标记，而虚线则表示通过InSAR变形结果识别出的持续滑动山坡区域（详见第5节和图10）。

如图7所示，HARMIE优化的相位表现出更大的幅度和更清晰的背景对比，与实际滑坡范围的对齐度更高，与基于幅度同质性优化的空间平滑相位相比。尽管上升轨道的俯冲效应在图7f中明显，但HARMIE方法有效地保留了滑动山坡的局部相位变化。此外，HARMIE得出的相位更准确地捕捉了局部相位异质性，如图7c、g和h所示。

对图7中几个代表性滑坡的相位异质性恢复进行了定量评估。对于每个滑坡，计算了滑坡范围内的像素相位与周围非滑坡背景的平均相位之间的 [Math Processing Error]$L_1$ 范数距离 [Math Processing Error]$\omega$，以衡量异质相位的幅度。为了比较基于幅度和HARMIE方法的结果（[Math Processing Error]${\omega }_{\mathrm{amp}}$ 和 [Math Processing Error]${\omega }_{\mathrm{HARMIE}}$），每个滑坡的值被标准化为 [Math Processing Error]${\omega }_{\mathrm{amp}}$ + [Math Processing Error]${\omega }_{\mathrm{HARMIE}}$。如表1所示，[Math Processing Error]${\omega }_{\mathrm{HARMIE}}$ 在选定的滑坡中平均比 [Math Processing Error]${\omega }_{\mathrm{amp}}$ 高13%。需要注意的是，HARMIE的两阶段选择导致像素数量减少，这可能会保留去相关异常，并略微高估之前的定量评估，尽管总体的改进依然显著。在严重去相关的像素中，相位分布为均匀分布时，HARMIE可能无法进行像素选择，导致检索的相位主要由原始空间随机性主导。

总体而言，在相干条件下，这些发现表明HARMIE在保留细尺度相位变化和幅度方面具有优越的能力，从而在检测微小地面运动和完整滑坡范围方面优于基于幅度同质性的优化方法。

在森林覆盖的区域中检测小规模滑坡信号时，当前的InSAR方法面临以下困难。传统的基于幅度同质性的方法没有考虑到相位变化，因此在选择同质化像素时可能包含相位或幅度不一致的像素，这会导致对 [Math Processing Error]$\stackrel{?}{\Gamma }$ 的偏差。为了更有效地检测滑坡，我们开发了一种新的方法，该方法通过结合幅度和相位属性来选择混合同质像素，并和谐地调整其贡献以更好地近似 [Math Processing Error]$\stackrel{?}{\Sigma }$ 的MLE，从而恢复最佳相位历史。这种方法可以减少显著的去相关，同时保留局部相位异质性，以确保微小的坡地移动信号不会被淹没。我们随后定量相关连续相位梯度与每个获取周期内的极端降雨强度，识别出在空间上聚集的坡地移动易发区域。最终，这种提出的框架使我们能够对广阔森林区域的滑坡易发区进行映射，主动定位隐藏但活跃的滑动坡地，并支持对未来极端降雨事件的灾害缓解规划。

本研究的主要贡献在于提出了一种新的滑坡检测框架，该框架结合了HARMIE方法和基于相位梯度的检测方法，从而在高分辨率的L-band SAR数据集上提高了对小规模、森林覆盖滑坡的识别能力。HARMIE方法通过联合幅度-相位分析优化了同质化像素选择，并和谐地调整了像素间的相位不平衡，从而提升了相位估计的准确性。而基于相位梯度的检测方法则通过计算相邻像素间的空间差异，避免了相位解缠的需求，并且能够抑制长波长的大气效应，增强变形边界的定义。通过将这些方法应用于Lantau Island的LT-1 SAR数据集，我们成功地识别了2023年10月极端降雨触发的广泛滑坡，并且相比基于幅度同质化的检测方法，识别率提高了27%。这一成果表明，HARMIE方法在检测小规模滑坡方面具有显著的优势，能够更准确地捕捉微小的地面变形信号。

在实际数据实验中，我们利用了LT-1 SAR的升轨和降轨数据集，时间跨度为2023年7月5日至2024年10月11日（降轨）和2023年8月2日至2024年10月23日（升轨）。这些数据集与2023年9月-10月的两个热带气旋（台风海葵和台风康森）相关，这些气旋对研究区域产生了显著影响。HKO记录的累积降雨量超过400毫米和300毫米，分别对应台风海葵和台风康森。在Lantau Island，台风海葵期间最高每小时降雨量达到46毫米（HKA），而台风康森期间达到65毫米（CCH）。相比台风康森，台风海葵对Lantau Island的影响更为有限。2024年5月4日的一场严重降雨事件在CCH记录了每小时84毫米的降雨量，这表明存在局部降雨异常。其包含在GERP计算中会降低区域代表性，但会增强对空间异质性的敏感性，这可能减少因局部降雨触发的滑坡被遗漏的风险。此外，Lantau Island上广泛分布的坡积物，使坡地在极端水文气象条件下更容易发生广泛的滑坡，这也支持了使用GERP作为合理区域指标的合理性。

在滑坡检测过程中，我们首先通过计算升轨和降轨的相位梯度，与GERP进行相关分析。通过这一过程，我们能够识别出滑坡点，并将它们聚类。HARMIE方法通过引入相位异质性权重，有效提升了滑坡点的识别率。在图8中，我们展示了滑坡检测过程和结果，包括基于升轨和降轨优化的LT-1干涉图的相位梯度与GERP的相关分析。多个典型的梯度区域被放大，以突出极端降雨事件（2023年10月9日）对滑坡的影响。我们选择了四个代表性滑坡区域，并在图8中展示了它们的梯度时间序列，分别用深色和浅色圆圈标记。为了提高检测的准确性，我们还利用了Lantau Island上的LiDAR数据生成高分辨率的数字表面模型（DSM），并将其用于去相关区域的掩膜处理。

通过图8的分析，我们可以看到，GP1-GP4的相位梯度时间序列几乎与 [Math Processing Error]$S_{\mathrm{desc}}^{\mathrm{ext}}$ 或 [Math Processing Error]$S_{\mathrm{asc}}^{\mathrm{ext}}$ 保持同步。所有主要的梯度峰值都与2023年10月9日的极端降雨事件相吻合，该事件由台风康森带来，预期在Lantau Island的不稳定坡地上发生显著的下滑运动。通过将这些梯度与GERP进行相关分析，我们能够区分滑坡点与其他受干扰的点。通过密度聚类，我们获得了全面的滑坡点库存。图8B展示了通过整合升轨和降轨结果获得的滑坡分布，叠加在Lantau Island的西部（主要滑坡影响区域）的阴影图上，并与通过LiDAR DSM确定的坡地滑动边界进行对比。检测到的滑坡集群（洋红色正方形）与2023年11月TDOP确定的滑动边界高度一致，证实了检测的可靠性。这种空间匹配，加上2023年10月GERP峰值与对应相位梯度峰值的时间一致性，表明这些坡地滑动是由2023年10月的极端降雨触发的。

在滑坡检测结果的比较和评估中，我们还使用了基于幅度同质性的优化LT-1干涉图生成了额外的检测集，以进行比较。89个已确认的坡地滑动点作为参考点，用于评估HARMIE和幅度方法的检测率。检测到的滑坡和未检测到的滑坡，以及两者均检测到的滑坡，如图9所示。HARMIE方法成功检测了72个滑坡（80.9%的检测率），其中包括幅度方法未检测到的26个滑坡。相比之下，幅度方法仅检测了48个滑坡（53.9%的检测率），其中只有两个未被HARMIE检测到。对图9中结构化结果进行的McNemar检验确认，HARMIE在检测率上比幅度方法取得了统计学上的显著改进。

图9显示，虽然一些滑坡被幅度方法检测到，但检测到的点通常比HARMIE方法更稀疏，并且与实际的滑坡轮廓相比不够一致。值得注意的是，HARMIE方法能够检测到宽度不足10米的小规模滑坡。此外，现场滑坡库存数据（2023年10月的坡地滑动）被用于验证结果（图9中的蓝色圆圈标记了参考点）。在图9b、d、h和i中，这些参考点被用于评估检测分布。需要指出的是，现场检查通常沿着可到达的道路或徒步路线进行，这使得一些真正的偏远且陡峭的滑坡区域未被记录。某些在Lantau Island南部记录的滑坡在InSAR结果中未被检测到，因为它们位于升轨LT-1覆盖范围之外，且在降轨视图中不够敏感。总体而言，当与提出的滑坡检测方法相结合时，HARMIE方法展示了对局部坡地位移信号的更好保留，并且在检测率上比幅度方法提高了27%，这证实了所提出的滑坡识别框架的有效性。

通过研究，我们能够更深入地理解滑坡的发生机制及其与降雨的相互作用。例如，我们发现长期的前期降雨可能会引发坡地的季节性缓慢蠕动，而某些坡地则可能在极端降雨的冲击下加速失稳。这表明，滑坡的发生不仅受到降雨的直接触发，还可能受到前期降雨条件的影响。在Lantau Island的特定区域，我们观察到滑坡点的分布与极端降雨事件的时间一致，尤其是在台风康森带来的降雨高峰期间，坡地滑动的显著增加表明降雨是触发滑坡的关键因素。此外，某些滑坡点的分布还与长期的降雨活动相吻合，这表明前期降雨可能在滑坡发生前对坡地进行了预处理，使其在随后的极端降雨中更容易发生滑动。

我们还发现，滑坡的动态过程可能与地形特征和地质条件密切相关。例如，某些滑坡点的分布显示出与降雨事件的时空一致性，而其他滑坡点则表现出与长期降雨的关联性。这表明，滑坡的发生可能受到多种因素的影响，包括降雨的强度、持续时间以及地形的稳定性。因此，滑坡预警系统需要综合考虑这些因素，以更准确地预测滑坡的发生。

本研究的结论表明，HARMIE方法在检测小规模滑坡方面具有显著优势，特别是在植被覆盖的区域。通过结合幅度和相位同质化像素选择，HARMIE能够更准确地保留局部相位变化，从而提高对微小地面运动的识别能力。同时，基于相位梯度的检测方法能够将坡地响应与极端降雨联系起来，从而提供更全面的滑坡分析。这些方法的应用不仅提高了滑坡检测的准确性，还为滑坡预警系统的开发提供了新的思路。未来，随着InSAR技术的不断发展，我们有望在更多复杂的地质环境中应用这些方法，以更有效地识别和监测滑坡风险。