基于地基雷达图像的无监督变化检测

黄平平 任慧芳 谭维贤 段盈宏 徐 伟 刘 方

(内蒙古工业大学信息工程学院 呼和浩特 010051)

(内蒙古自治区雷达技术与应用重点实验室 呼和浩特 010051)

1 引言

我国地质构造复杂、地形地貌起伏变化大，极易发生滑坡、崩塌、泥石流等地质灾害。此外，人类的生产、生活活动也很容易诱发滑坡灾害的发生，大量岩土边坡的失稳或滑动屡见不鲜，给工程建设和人民生命财产带来了严重的危害和损失。灾害发生区域的变化检测分析对安全生产有着重要意义。从遥感数据中获取灾害信息是一种十分重要的灾害研究手段，而快速有效的变化检测方法对于掌握灾害情况尤为重要。地基雷达技术具有全天时、全天候、大范围、远距离、非接触等特点，在地质灾害监测、预警方面具有巨大优势。目前，地基雷达已广泛应用于滑坡、崩塌等地质灾害的监测，利用干涉测量原理地基雷达可以监测到目标区域发生的微小形变量，根据累计形变量、形变速度、形变加速度等对地质灾害进行实时监测。然而，在长期施工建设监测中，部分区域的变化并不会引起滑坡等灾害，但受人为因素、地质因素、气象因素等影响，导致监测区域的雷达图像失相干较为严重，给长期定量化监测带来了较大的难度。因此，迫切需要在定量监测的基础上，进一步开展变化检测方面的应用，为长期全面了解监测区域的动态变化提供有效信息。针对上述问题，本文将研究一种基于地基雷达图像的变化检测方法。

一般来说，图像变化检测技术可以分为两大类：监督式和无监督式。监督变化检测方法，需要先验信息；而无监督变化检测方法直接自动比较多时相的遥感图像，不需要任何额外信息，降低人为误差的影响，符合实际应用中先验变化信息缺失的现实情况，因此得到了广泛地研究[1]。雷达图像中的无监督变化检测一般包括3个步骤：(1)图像预处理；(2)获取差分图像；(3)分析差异图像及后处理[2,3]。

近年来雷达图像变化检测的研究主要集中在机载和星载合成孔径雷达(Synthetic Aperture Radar,SAR)图像上，基于星载SAR图像中的无监督变化检测，龚茂国等人[4]提出了一种基于邻域的比值法来产生差异图，虽然该算法优于产生差异图的传统方法，如差值法和比值法，但是其检测结果受噪声干扰严重；为了提高变化检测结果的准确率，使用均值比图像和对数比图像的互补信息来生成差异图，取得了较好的检测结果[5]。Sumaiya等人[6]采用了一种经济且简单的对数均值阈值进行星载SAR图像变化检测。杨祥立等人[7]运用D-S证据理论(Dempster-Shafer evidence theory,D-S)融合高分辨率SAR影像的相干/非相干差异特征进行变化检测，采用美国旧金山地区港口附近的TerraSAR-X数据切片进行实验验证。此外，在图像变化检测中一个重要的步骤是图像的分割聚类，Elguebaly等人[8]提出了一种基于广义高斯混合模型(Generalized Gaussian Mixture models,GGM)的高效无监督图像分割和变化检测算法，论证了贝叶斯GGM在图像分割中的应用潜力。邵宁远等人[9]提出了一种面向变化检测的SAR图像超像素协同分割算法，解决了变化检测中存在的多时相图像边界和空间对应关系不一致的问题。Carincotte等人[10]提出了一种新的模糊隐马尔可夫链(Hidden Markov ChainS,HMCS)，从而解决了模糊变化检测的统计方法。而最著名的模糊聚类算法是模糊C均值聚类(Fuzzy C-Means,FCM)分割算法[11]，它主要的缺点是在成像应用中缺乏对空间环境信息的整合，导致其对噪声和其它人为因素较为敏感。而FCM的变形，如鲁棒FCM(Robust FCM,RFCM)[12]和空间条件FCM(the Spatial Conditional FCM,SCFCM)[13]在星载SAR图像聚类中都得到了较好的效果。

星载SAR适用于大范围对地成像，其成像范围可达数千平方公里，基于星载SAR图像变化检测适用于监测区域内发生的大面积、特征明显的变化。地基雷达通常针对特定的目标区域成像，多数情况下发生变化的特征不够明显，将星载SAR图像变化检测技术应用到地基雷达图像变化检测中，会造成差异图群分性比较差和聚类分割结果存在较多噪声点等缺陷。本文提出了一种无监督地基雷达图像变化检测方法，本方法通过对相干系数图和均值对数比值图进行非下采样轮廓波变换(NonsubSampled Contourlet Transform,NSCT)和局部能量法得到合成差异图，得到的合成差异图具有较好的群分性。然后采用主成分分析法(Principal Component Analysis,PCA)提取特征向量，增强抗噪性。最后使用改进的模糊C均值聚类(FCM)对特征向量进行聚类得到最终的变化检测结果。文章安排如下：本文第2部分提出了一种地基雷达图像变化检测方法；第3部分利用地基雷达实测数据进行结果分析并验证该方法的有效性；第4部分对全文做出总结。

2 地基雷达图像变化检测算法

本文所提基于地基雷达图像的无监督变化检测方法主要分为4个步骤：

(1) 初始差异图获取。先取两幅已配准的地基雷达图像X1与X2，通过相干系数算子以及均值对数比算子生成相干系数差异图和均值对数比差异图；

(2) NSCT变换。对得到的两幅初始差异图，进行NSCT变换并融合低频系数，然后使用融合的低频系数分别与初始差异图的高频系数进行NSCT反变换得到两幅融合结果图，对这两幅融合结果图进行像素能量合成得到最终的合成差异图；

(3) PCA提取特征向量。利用PCA算法提取出合成差异图的像素特征；

(4) 改进的FCM聚类。最后通过一种加相似性惩法项的模糊C-均值聚类算法(FCM)聚类，获取变化检测结果。

本文的算法流程如图1所示。

图1 地基雷达图像无监督变化检测流程图Fig.1 Flow chart for unsupervised change detection based on ground-based radar image

2.1 初始差异图获取

由于地基雷达图像自身存在乘性斑点噪声，对数比值法能够将乘性斑点噪声转化为加性噪声，便于有效去除斑点噪声。此外，对数运算作用于两幅地基雷达图中能够对像素的灰度值的比值做非线性拉伸，可以在一定程度上增强变化和非变化区域的对比。为了充分考虑像素点周围的特征，采用均值对数比获取的差异图为Xd，计算公式为

其中，m1(i,j)与m2(i,j)分别为图像X1和X2在像素点(i,j)处的3×3邻域均值。

雷达图像之间的相干系数用来衡量影像之间的相似程度，它提供了检测区域重要的变化信息。采用相干系数图，可以得到变化区域非常明显和背景信息较为平滑的差异图[14,15]。图像X1和X2的相干系数定义如式(2)所示

式中，*表示共轭复数，X1和X2为经过配准滤波的地基雷达图像，相干系数Xm(i,j)的值在0～1之间，0表示失相干，1表示完全相干。

采用对数比值图和相干系数图相结合的方法，综合利用各自的优势，能够得到群分性较好的融合差异图，为后续的图像聚类分割提供重要依据。

2.2 NSCT变换

NSCT是基于非下采样金字塔滤波器组(Non Subsamlped Pyramid Filter Banks,NSPFB)和非下采样定向滤波器组(Non Subsampled Directional Filter Banks,NSDFB)进行的，每个子带图像与原始图像在形状和大小上都是完全一样的，因此，NSCT具有多分辨率、局部化、方向性、各向异性和位移不变性等特点[16]。该方法能较好地捕捉图像的几何细节，保持目标信息和边缘信息[17]。

基于NSCT的图像融合算法通常将图像进行NSCT分解，针对不同频率域的特点选择不同的融合规则。对数比值图Xd和相干系数图Xm的融合结果图为Fd，具体实现步骤如下：

(1) NSCT变换。轮廓波变换使用“9-7”型拉普拉斯金字塔滤波器，使用“pkva”作为方向组滤波器，每个金字塔级方向滤波器组分解层数为4。经NSCT变换得到差异图Xd和Xm对应的低频系数Ld和Lm，高频系数Hd和Hm，高、低频系数反映了高、低频图像信息。

(2) 加权平均法。由于图像的低频信息通常反映图像的结构信息[18]，因此在融合过程中，为了凸显变化区域，针对低频系数采取加权平均的方式，对差异图Xd和Xm经NSCT变换得到的低频系数Ld和Lm进行融合，得到融合后的低频系数Lf，计算过程为[19]

其中，α(α ≥0)为加权平均的融合系数，本文实验中α=0.3。

(3) NSCT反变换。用融合的低频系数Lf分别与差异图Xd和Xm的高频系数进行NSCT反变换，得到融合结果图Fd和Fm[19]。

因为地基雷达图像的相干系数图的背景相对平坦，可以采用局部能量法抑制背景信息[5]，得到最终的合成差异图Fh为

其中，Ed(i,j)和Em(i,j)分别表示差异图Fd和Fm的局部能量，Li,j表示以像素(i,j)为中心，大小为3×3的邻域，Fd(q)和Fm(q)表示邻域内第q个像素值[10]。

2.3 PCA提取特征向量

主成分分析(PCA)，也称Karhunen-Leve变换(简称K-L变换)，是一种线性变换。它是一种均方误差最小的最佳正交变换，是在统计特征基础上的多维线性变换。本文采用PCA方法对融合差异图进行特征提取，由于采用分块的思想，使得该方法具有一定的抗噪能力[20]。

将合成差异图分为n个大小为3×3的相邻的小块，小块之间互不重叠。每个小块可以看作一个矩阵，然后将每个小块矩阵转换为列向量Pt，其中t=1,2,…,n，并计算全部列向量的均值向量

将所有列向量重构为32×n的矩阵S，计算S的协方差矩阵C

其中，T 表示矩阵转置。协方差矩阵C是一个大小为32×32的矩阵。对协方差矩阵进行特征值分解，求出特征值和特征向量，按特征值从大到小的顺序排列，选出对应的特征向量，这些特征向量形成一组正交基，取特征向量的前r个列向量构成矩阵A,A是一个大小为32×r的矩阵。取合成差异图Fh中每个像素所在的3×3邻域小块，将每个小块转化为列向量Yη,η=1,2,…,k,k为合成差异图Fh中像素的总个数。将所有的列向量映射到矩阵A上，即

其中，Yη为每个像素代表的特征向量，AT为A的转置矩阵，这样每个像素都用一个R维的特征向量来表示，[Y1Y2…Yk]构成特征向量空间Q，Q即为一个r×k矩阵。

2.4 改进的FCM聚类

2.4.1 FCM算法

FCM是一种无监督的模糊聚类方法，传统FCM[11]目标函数的形式为

其中，uηi是像素点Yη在i类的隶属度矩阵，c为分类数，m为模糊权重(一般为2),η为图像的像素点个数，Ci为聚类中心。

uηi表示Yη隶属于第i类的概率，介于0和1之间。因此，约束必须满足

考虑到式(10)中的约束，通过拉格朗日数乘算子计算当式(9)中的目标函数达到最小值时，uηi为

聚类中心Ci为

当用于图像聚类时，FCM没有考虑到相邻图像像素之间的空间关系[21]。采用变形的FCM为鲁棒模糊C-均值(RFCM)[12]的目标函数采取了如式(13)所示形式

其中，Nη是像素点Yη的空间邻域集合，Pi是除i类外的所有类的集合，β是控制附加空间惩罚项效果的正偏差。

增加的空间惩罚项鼓励像素点在同一类中拥有有较高的聚类资格。最小化式(13)得到隶属度函数为

Ci的计算过程为式(12)，因为惩罚项不依赖于它。

FCM的另一个变型也在其目标函数中包含了空间关系，它是在目标函数中引入的空间条件FCM(也称SCFCM)[13]。具体来说，它的目标函数采用如式(15)的形式

最小化式(15)得到隶属度函数为

Ci的计算过程为式(12)。

2.4.2 改进的FCM算法

为了弥补FCM和其变形RFCM和SCFCM分割率低以及对于地基雷达图像噪声敏感等问题。本文对FCM的目标函数进行改进，在原始的FCM的目标函数中加入相似性惩罚项。目的是使图像中相类似的像素点以更大概率聚类到相同的类，提高聚类效果并抑制噪声点。改进后的聚类算法的目标函数的表达式为

其中，Yη为输入每个像素的特征向量，分类数c=2(变化类和非变化类)，m模糊权重(一般为2),β控制惩罚项的积极因子，本文中β=0.2615。Nη为Yη及其周围邻域的一个集合，l为此集合中的一个元素。

改进的目标函数在保留传统FCM目标函数的基础上，添加了相似性惩罚项。需要强调的是，在式(17)的邻域中，邻域是以相似性而不是以空间定义的，这样可以找到每个输入图像的最相似邻域，以鼓励相似的像素以更大的概率聚类到同一类中。通过拉格朗日数乘算子计算当式(17)中的目标函数达到最小值时，uηi为

由于相似度惩罚项不包含Ci，聚类中心采用式(12)中给出的FCM算法的相同形式。改进的FCM聚类算法计算过程如表1所示。

通过比较每个像素的隶属度uηi的大小，uη1和uη2分别代表第η个像素对应两类隶属度，采用如式(19)所示的判决条件对uη1和uη2进行二值化处理

R为最终变化检测结果图。

3 实验结果与分析

我国西南某省发生地震之后山体严重垮塌，形成堰塞体，为确保堰塞体整治工程安全、顺利进行，在现场部署地基雷达LSA对堰塞体治理施工过程进行监测。经现场调研，地基雷达LSA布置在距堰塞体300 m正对面，雷达的监测范围完全可以有效覆盖堰塞体。表2为地基雷达LSA系统参数，图2所示为监测现场光学图。监测期间，施工地出现持续性强降水，导致边坡发生滑坡。如图3(a)为滑坡前雷达图，数据获取时间为2019年7月15日，图3(b)为滑坡后雷达图，数据获取时间为2019年7月19日。截取滑坡前后数据切片对进行实验，如图3(c)和图3(d)所示。

表1 改进的FCM算法计算步骤Tab.1 Improved FCM algorithm calculation steps

图4 初始差异图Fig.4 Initial difference images

采用均值对数比值算子和相干系数算子对图3(c)和图3(d)的实验数据进行处理，得到均值对数比值图和相干系数图，如图4所示。将相干系数图和对数比值图经过NSCT处理融合(方法1)产生的合成差异图Fh与通过邻域均值比和对数比值图经过NSCT处理融合(方法2)产生的差异图进行对比，如图5(a)和图5(b)所示，其中colorbar代表图像像素值域，方法2为星载SAR图像变化检测生成差异图最常采用的方法[19]。将两种方法产生的差异图显示在3维图中，3维图如图5(c)和图5(d)所示，X,Y轴代表图像像素的位置信息，Z轴代表差异度值越大发生变化的可能性越大，可以明显发现通过相干系数图和对数比值图经过NSCT处理融合的差异图群分性最好，这就代表了本文产生差异图的方法更容易进行聚类分割，进而区分变化区域与非变化区域。

对两时相的地基雷达图像进行变化检测，本文方法产生的合成差异图采用改进的FCM算法与传统的FCM算法的变形(包括RFCM和SCFCM)进行聚类结果对比，由于采用FCM算法进行聚类，无法得出聚类结果，在本文中不作分析。其结果如图6所示，图6(a)为RFCM聚类结果，图6(b)为SCFCM聚类结果，图6(c)为采用本文改进的FCM聚类算法得到的变化检测结果，即式(19)的变化检测结果R，其中蓝色的像素值为0，表示变化区域；黄色的像素值为1，表示未变化区域。

图6(a)中可以明显地看到，RFCM算法聚类结果含有较多噪声点，无法聚类出结果。图6(b)中SCFCM算法去除较多噪声点，但依然有少量噪声存在，如图6(b)中红色矩形区域。而图6(c)中本文算法很好地滤除了周围的噪声点，得到的变化检测结果更加准确。通过实地考察，证实实验数据对应区域发生明显的滑坡，图7所示为监测现场光学图，图7(b)中黄色勾选区域发生了明显的滑坡。此外，图6(c)和图6(b)相比，椭圆红框中检测到的变化区域存在较大差异，本文方法检测到该部分并没有发生明显的变化区域，在滑坡现场实地考察过程中发现滑坡前后此处为地质条件较为稳定的岩石，并没有随着滑坡的发生有所变化，只有两块岩石中间的部分小区域发生浮土下滑，如图7(b)中最小的黄色框区域，由此可确定检测到的变化结果与实地考察结果相符。地基雷达图像变化检测结果与地形3维数据融合结果如图8所示，根据融合图，可以准确地定位变化区域，如图8中红色框区域，通过融合图，进一步验证了本文方法的有效性，在实际应用中，可以为安全生产提供有效监测数据。

图5 差异图对比Fig.5 Difference images comparison

图6 聚类结果对比Fig.6 Clustering result comparison

图7 监测现场光学图Fig.7 Optical image of monitoring site

4 结束语

图8 变化检测结果与地形3维数据融合图Fig.8 Effect of fusion of change detection results and UAV tilt photography

对于地基雷达图像采用均值比以及对数比都无法得到群分性较好的差异图，本文提出使用相干系数与均值对数比得到群分性较好的差异图，使得地基雷达图像在变化检测过程中更容易进行聚类分割。此外，本文对传统FCM聚类算法进行改进，使图像中相类似的像素点以更大概率聚类到相同的类，以提高聚类效果并抑制噪声点。通过将所提出的方法应用于地基雷达实测数据，并将改进的FCM与FCM的另外两种变形RFCM和SCFCM聚类算法进行聚类结果对比，结果表明，本文所提方法具有更好的聚类效果，并且变化检测结果在保留变化区域的同时噪声点明显减少。未来的工作中，拟进一步探索地基雷达图像的无监督变化检测方法，以得到更为精确的变化检测结果，为安全生产、排险、救援提供数据参考。