欢迎访问中国地质灾害防治与生态修复协会

当前位置：首页> 专题 > 技术推广 > 地灾监测

基于数字形态学的InSAR形变主动滑坡自动识别研究：以中国白鹤滩水库为例

来源：智慧灾害调控 发布时间 2024-12-12

摘要

光学遥感与实地调查无法满足主动滑坡探测对精度和时效性的要求。干涉合成孔径雷达（InSAR）技术因其具有大范围探测和对地表形变高灵敏度的优势，近年来已成为监测滑坡的主流方法。然而，如何快速且准确地从InSAR观测结果中提取滑坡边界，仍是灾害防控与流域管理的关键问题。本研究首先建立了一种基于InSAR观测结果的主动滑坡自动识别方法，以快速提取形变斜坡。在识别过程中，以InSAR形变图为目标，利用图像梯度边缘信息确定最佳阈值以提取变形像素。其次，采用对比度受限自适应直方图均衡化（CLAHE）算法增强弱变形图像的对比度。最后，应用形态学规则优化分割结果，确保其接近自然滑坡的边界。为验证该方法的有效性，选取滑坡频发的金沙江流域下游白鹤滩库区作为测试区。在升轨与降轨的观测条件下，分别识别出336处与590处滑坡。通过无人机（UAV）和实地调查验证，该方法在无需样本训练的情况下，识别精度达76%，单一滑坡的最大交并比（IOU）提升了0.3。研究结果表明，该自动识别方法能够在大空间尺度及复杂地形条件下快速识别危险的主动滑坡。

方法

研究区域

白鹤滩水电站位于金沙江下游（如图1所示），是中国第二大水电站，总装机容量为1600万千瓦，坝高289m，库区正常蓄水位825m，库容206亿立方米，水位变化一般为60m。库区位于青藏高原和川西高原东缘向云贵高原和四川盆地的过渡地带。地质构造特征主要为褶皱和断层，地形地貌以海拔4000米以上的高原和山地为主，呈现典型的V型峡谷地貌。库区坡体海拔900米以上区域的坡度一般为30°至50°，年平均气温为21.7℃，年均降水量为715.9毫米，年降水天数约为100天。在汛期前水库排空以及汛期后蓄水过程中，水位每年长期且剧烈波动，对库区岸坡（特别是古滑坡堆积体）的稳定性极为不利，使其成为滑坡活动的多期活跃区。

图1 研究区域及现场照片

研究方案

本研究构建了一种基于形态学规则的滑坡自动识别方案，以InSAR形变值为目标对象（如图2所示）。在InSAR观测阶段生成研究区域的形变值，随后在预处理与初始分割阶段，首先计算图像梯度幅值的均值和标准差，以确定最佳阈值。利用该阈值，计算机能够快速去除正常形变，仅保留异常形变区域的像素点。为增强形变区域的显著性，本研究引入对比度受限方法，并通过直方图分布特性调整形变像素，使目标像素被“均匀”分配到概率密度分布中，从而拉伸图像的动态范围，突出地表弱形变区域。

此外，在利用形态学优化变形斜坡边界时，采用膨胀操作填充识别区域内的空值并平滑其边界线；随后通过腐蚀操作消除特征区域外的噪声信息，保留有效形变信息。该方法不仅能够增加特征区域的信息丰富度，还能增强形变区域位置的显著性。因此，通过图像的地理编码和参数校核，可以输出变形斜坡的地理信息数据。

图2 InSAR滑坡形态识别流程图

数据和InSAR处理

本研究采用了Sentinel-1A卫星在升轨和降轨条件下的观测影像，时间跨度为2021年1月至2022年3月，数据参数如（表1）所示。观测范围为宽幅干涉模式(IW)，极化模式为VV。升轨与降轨影像的中心视角为39.64°，相邻影像的最小时间间隔为12天，最大时间间隔为24天。选用Stacking-InSAR技术作为InSAR计算方法，以获取地表平均形变速率。首先，对时间基线和空间基线较短的SAR数据进行D-InSAR计算，以获得高质量的干涉相位。其次，去除层状大气效应，并对湍流大气进行Gurkin空间滤波。最后通过叠加多场景干涉结果，有效消除低频形变，突出高频主动滑坡的空间形态特征。在得到的形变图中，每个像素值代表地表的相对形变量，采用“hls”色带赋值，形变量的值范围表示为[−π, π]。在此过程中，本研究使用日本宇宙航空研究开发机构（JAXA）提供的空间分辨率为12.5米的数字高程模型（DEM）作为辅助数据，用于InSAR观测中的地形误差消除和地理编码。

表1 研究区域SAR数据的基本参数

实验

形变像素值提取

在InSAR结果的计算过程中，不同的数据和处理方法会导致检测目标体的分布存在差异，但在形变值图中，变形斜坡与背景之间的灰度值差异并不明显。由于观测结果是基于雷达反射信号干涉处理以获取地表形变特征，在InSAR观测结果图像中，变形斜坡的像素值与周围环境存在显著差异，因此可以利用图像梯度边缘信息来确定最佳阈值。在本部分中，本文构建了梯度边缘方法（Gradient Edge Method，简称GA），并采用一阶微分等效算子（Sobel算子）计算输入图像在水平方向和垂直方向的梯度矩阵。随后，使用高斯滤波器对梯度矩阵中的孤立像素进行滤波处理。接着，计算梯度矩阵的均值和标准差，并最终选择两者中的较大值作为最佳阈值m。应用Sobel算子后，原始图像中x方向的梯度为水平方向梯度，y方向的梯度为垂直方向梯度，其计算公式如下：

在确定最佳阈值后，将源图像中的每个像素与该最佳阈值进行比较。如果像素i的值大于阈值 m，则该像素被归类为形变区域，并取其绝对值后赋值为d；否则，像素i被归类为背景区域，其赋值按以下公式计算。

在后续的数据处理阶段，本文采用Min-Max方法将原始数据映射到区间 [0,255]。该方法不会改变数据分布，同时能够为所有数据属性赋予相同的权重，从而使属性间的比较与聚合更加容易，并改善结果的收敛条件。为验证GA方法的分割效果，本文设计了对比实验。对比方法包括多模态直方图方法（HI）、一维交叉熵方法（EN）、OTSU方法以及高斯自适应方法（AG）。

弱信号增强

常用的图像增强算法通过调整图像的整体灰度级，容易导致灰度级的压缩，无法有效提升局部对比度。特别是对于灰度值范围不均的InSAR观测结果，局部弱形变无法被有效突出，反而引入了噪声，导致观测结果的失真。在本节中，基于现有灰度图像增强技术，本文增加了自适应调整过程，并引入对比度受限方法。在图像增强过程中，通过计算局部区域的直方图分布构建映射函数，并对亮度值进行重新分配以改善图像质量。超过阈值的部分被“均匀”分布到概率密度曲线中，从而限制累计直方图的过度增加。该过程拉伸了图像的动态变化范围，增强了图像对比度，并突出了地表弱形变区域。

此外，阈值的计算公式为：

需要使用计算得到的阈值对每个子块的灰度直方图进行截断，并将被截断的多余像素数平均分配到其他灰度级中。分配到每个灰度级的平均像素数计算公式为：

在对每个子块重新分配像素后进行直方图均衡化后，还需要通过双线性插值法计算每个子区域的像素值。在计算过程中，以每个子块的中心点作为参考点。例如，需要计算的点的像素值由其四个相邻的参考点确定。在计算过程中，首先在x方向进行线性插值操作，分别得到像素值：

随后在y方向再次进行插值操作，以获得点的像素值：

双线性插值仅适用于包含四个参考点的中心点，四个顶点位置的像素值通过映射函数转换获得，而边缘处剩余像素的计算需要通过相邻的两个参考点进行插值。在完成插值后，图像的块效应被消除，并可通过高斯双边滤波进一步处理，以提高信噪比。本文设计了对比实验验证CLAHE算法的有效性，对比方法包括直方图均衡化（HE）、拉普拉斯方法（LAP）、对数变换方法（LOG）以及幂律变换方法（PLT）。

数字形态学

在提取变形像素并进行图像增强后，构成变形斜坡的像素呈离散分布，且二值图像中仍存在噪声、孔洞等现象。基于数字形态学，通过定义选择结构元素，对待处理图像进行比较与匹配，从而提取目标图像的形状结构。随后，通过形态学算子组合计算，去除背景干扰信息，保持目标几何特性，并消除无关形状，从而简化原始图像结构，确保图像的空间结构特性及目标几何特性。因此，在形态学方法中引入膨胀、腐蚀、开运算和闭运算等步骤，可以解决斜坡边缘锯齿状和内部空值现象。

灰度腐蚀操作会使图像变暗，其灰度图像腐蚀操作定义如下：

当在其上移动时，将所包含的像素值中的最小值赋予的原点。这一操作能够去除二值图像特征区域中的噪声，填补空值，并使特征区域中的线条更加平滑。

经过灰度膨胀操作后，图像会变得更加明亮。灰度图像的膨胀操作定义如下：

当结构元素在图像上移动时，将图像上被覆盖的像素值中的最大值赋予的原点位置。该操作可以去除特征区域外的背景信息，仅保留有效特征。

灰度开运算首先对图像进行腐蚀操作，然后在腐蚀操作之后进行膨胀操作。灰度开运算的定义如下：

该步骤可以去除图像中的孤立像素、单个像素之间的连接以及边缘粗糙区域，从而平滑图像轮廓。

灰度闭运算与灰度开运算互为对偶，即先对图像进行膨胀操作，然后在膨胀操作后对图像进行腐蚀操作，灰度闭运算的定义如下：

该步骤能够平滑图像的轮廓，但与开运算不同的是，它可以填补目标区域中的空值并增强像素的空间相关性。

结论

从InSAR观测结果中快速且准确地提取滑坡边界是灾害识别中的迫切需求，而视觉解译难以满足灾害监测的时效性要求。利用无监督识别方法从InSAR观测中提取潜在滑坡，需要评估三个方面：最佳形变阈值的选择、变形斜坡内像素的不均匀分布，以及滑坡边界的优化。本研究提出了一种基于InSAR形变值图的新型自动识别方案。首先，基于梯度边缘信息构建阈值判定函数，用于提取地表变形像素。其次，采用对比度受限的图像增强算法拉伸图像的动态范围，突出地表弱形变区域。最后，进一步构建基于数字形态学的优化方法，用于去除干扰信息，保留遥感图像的空间结构和目标几何结构。

对比实验表明，以上步骤在变形像素提取、弱信号增强和边界优化方面均有改进。与监督学习算法相比，该方法不需要样本，直接使用形变值。对金沙江下游白鹤滩库区的验证结果表明，在升轨观测条件下和降轨观测条件下，该方法识别的滑坡变形总数分别达到336处和590处，平均识别精度为76%，单一滑坡的最大IOU提升了0.3。结果表明，本研究构建的方法能够通过突出地表弱形变区域并改善变形斜坡边界的识别能力，有效识别主动滑坡。该方法适用于大规模、短周期滑坡的识别，对于早期预警系统具有重要意义。

下一篇：地质灾害前沿|多周期振动下颗粒材料的相变：多次地震对引发滑坡的影响上一篇：纳克斯三维植被护坡方法