基于改进邻域比和分类的SAR图像变化检测

王 平 王宜怀 刘长勇 彭 涛

1(武夷学院数学与计算机学院 福建 武夷山 354300)2(苏州大学计算机科学与技术学院 江苏 苏州 215006)3(认知计算与智能信息处理福建省高校重点实验室 福建 武夷山 354300)

0 引 言

合成孔径雷达(SAR)图像的变化检测，目标是生成一个变化图，表述两时相或多时相图在不同标定时间之间的变化情况，其广泛应用在环境监测、生态监测、城市发展研究、农林监测、自然灾害评估等领域[1]。SAR图像以微波成像，图像强度值与地面目标的后向散射成正比，易产生不易滤波处理的斑点噪声，进而对变化检测精度造成较大影响。如何在噪声抑制和变化细节保留之间进行均衡，提升变化检测精度，是处理过程中需要解决的重要问题。

从SAR图像变化检测技术流程来看，主要有图像预处理、变化信息发现、区域提取或变化类型确定等几个步骤。其中，信息发现和区域提取是变化检测的研究热点，且方法众多。从差异图(DI)角度进行变化检测处理是一种常用的途径，其办法是先构建两时相SAR图的DI，然后再对DI分析处理。

DI构建方法中，传统的方法有差值法、比值法(R)[2]、改进比值法(IR)[3]、对数比值法(LR)[4]、均值法等，但这些方法均没有利用像素点周围的邻域空间变化信息来减少斑点噪声的影响，进而造成精度较低。近年来提出的组合差异法的思路是将各种方法的优点结合起来，特别是寻求噪声抑制和图像变化细节保留之间的均衡[5]。MR方法在处理中引入像素点的邻域空间均值信息，能抑制斑点噪声，图像平滑，但未考虑区域异质性，使得变化边缘等异质性较大的地方检测不稳定[6]。Liu等[7]利用MR方法和差值法进行组合获取DI，具有较好的效果。而NR方法采用自适应斑点噪声滤波算法常用的局部邻域异质性测度来控制邻域空间信息对中心像素点变化检测的影响，其本质是将IR方法和空间邻域均值信息通过异质性测度来融合，但DI像素样本分布不利于DI分析[8]。

在DI分析中，常用方法是将DI的像素点分为变化类和不变化类。文献[5，7]针对DI采用k为2的K-means聚类，但对斑点噪声影响的鲁棒性较差，误检率较大。王建明等[9]采用基于自适应距离的FCM算法，将DI聚类出变化和不变化的区域类别，精度有所提高。文献[4]采用基于遗传算法改进的FCM算法对DI进行处理，识别变化类和不变化类，进而得到变化检测图。文献[10]在DI上基于多粒度级联森林和多尺度融化的方法进行处理，进而得到变化检测图。文献[11]设计了识别DI特征空间的层次级联FCM算法来得到变化和未变化的信息。但上述方法的样本分类依据仅仅来源于DI，未考虑DI中残留噪声的影响，仅仅是粗粒度的分类。文献[12]在利用FCM将DI像素样本聚类成三类的基础上，利用加速遗传算法进一步对不确定变化与否的像素样本进行分类进而获得变化检测图。文献[13]设计了层次FCM算法来分析DI得到样本标签，然后再将样本结合原始图像邻域块进行SVD网络分类，进而得到变化检测图，对DI的噪声残留有一定的抑制效果。崔斌等[14]在层次FCM聚类的基础上，利用卷积神经网络对中间类进一步进行分类处理，但该方法对中间类处理算法复杂，耗时较长。

基于此，本文依托IR方法的高效性和NR方法抗噪性强的特性，提出了新的邻域比方法INR，使得所构建的DI样本分布更利于分析。而在DI分析中，常用的直接聚类得到的检测结果忽略了DI中残留噪声的影响。本文利用层次聚类获得样本标签、输入时相图邻域块特征向量作为样本进行ELM网络分类训练而生成变化检测二值图，依托ELM网络简单有效的非线性映射能力，能解决聚类分析方法中有用信息丢失的问题，尤其是残留噪声影响的数据判别，实现对中间类细粒度处理，使得检测效果明显。

1 方法框架

假设有提前配准预处理的两时相SAR图像I1和I2，则本文所提方法的总体框架如图1所示。

图1 本文所提方法的总体框架

本文所提方法的具体步骤如下：

(1) 基于INR的DI构建。

(2) 对DI进行层次FCM聚类，并得到聚类差异图，该矩阵元素值为0、0.5和1。其中，0和1分别代表DI对应像素点属于不变化类和变化类，0.5代表可能噪声残留影响的中间类。

(3) 根据FCM聚类结果的变化类和不变化类像素，且以这些像素点为中心，提取两时相图的k×k图像块并组合在一起作为邻域块特征向量，依托ELM网络简单有效的非线性映射能力，构造ELM网络的训练样本集，训练ELM网络，并将两时相图的所有k×k邻域块特征向量作为测试样本集。

(4) 根据训练好的ELM网络，对测试样本集进行分类，可实现将FCM聚类结果值为0.5的中间类像素进一步确定为变化类和不变化类，最终得到与该两时相SAR图关联的变化检测二值图。

2 相关要点原理

2.1 差异图矩阵(DI)计算

SAR图像的配准预处理虽然实现了图像斑点噪声滤波和对齐等处理，但其残留的乘性斑点噪声仍然会严重影响变化检测的精度，因而，好的DI生成方法能提高变化检测的精度。假设In(i,j)表示SAR图像In中处于(i,j)位置上的像素强度值，i∈[1,M]，j∈[1,N]，其中，M、N分别对应时相SAR图像In的维度，则R方法可以表示为：

(1)

式中：DR表示采用R方法得到的DI；DR(i,j)表示DR中处于(i,j)位置的DI的像素强度值。经过分析，R方法只能检测诸如像素点强度值增大或减小变化类型中的一种，基于R方法改进得到IR方法则可以实现各种变化类型的检测。IR方法可以表示为：

(2)

则DIR(i,j)值越大，代表该像素点对应的原始输入两时相SAR图的像素未发生变化，反之，则发生了变化。实际应用表明，像素比值法对校正和辐射测量的误差都具有较好的鲁棒性，且融合SAR图像的邻域强度均值信息生成的DI可以减少乘性或加性斑点噪声的影响，这也是提出MR、NR方法的依据。MR方法直接根据IR方法来改进，其表示为：

(3)

式中：u1(x)、u2(x)分别表示SAR图像I1和I2的以(i,j)位置为中心邻域空间x(大小为r×r)的强度均值。DMR(i,j)值越大，代表该像素点对应的原始输入两时相SAR图的像素未发生变化，反之，则发生了变化。分析式(3)可知，MR方法不仅可以检测各种类型的变化，且在一定程度上能抑制斑点噪声的影响，能提高两时相SAR图像变化检测精度。NR方法正是结合了IR和MR方法优点，采用异质性测度值θ来作为权重比，其表示为：

(4)

式中：Ωx表示两幅SAR图像以(i,j)位置为中心的邻域空间x的集合；I1(x)、I2(x)定义为取以(i,j)位置为中心的邻域空间x块矩阵像素强度值；θ=σ(x)/μ(x)是邻域Ωx的异质性测度，σ(x)是邻域Ωx的标准差，μ(x)是邻域Ωx的均值。较小的θ对应一个同质区域，则式(4)的后一部分所占比重大，反之，较大的θ对应一个异质区域，则式(4)的前一部分所占比重大。

基于NR思想，本文在IR的基础上引入异质性测度θ，实现融合中心像素点强度和局部邻域均值信息的线性组合，形成本文提出的INR方法，其表示为：

DINR(i,j)=

(5)

式中：DINR表示采用INR方法得到的DI；DINR(i,j)表示DINR中位置处于(i,j)位置的DI的像素强度值；u1(x)、u2(x)定义如前所述；θ1和θ2分别是SAR图像I1和I2是邻域x的异质性测度；常数C解决分母可能趋近于0问题，常根据C=(kL)2取值，L为SAR图像灰度级别，如8位灰度强度则取值为255，k一般取值为0.003。

分析式(5)可知，DINR(i,j)值越大，则两时相SAR图间的相似度越大，对应的像素点也越可能没有变化，反之，DINR(i,j)值越小则越可能发生了变化。这可解释为：若有一个以像素位置(i,j)为中心的同质区域，此时θ取值小，表明此区域在前后时相图像中是未变化的，则式中邻域x的均值信息u(x)占比重较大，而可能受斑点噪声干扰的中心像素点强度值I(i,j)占比重小，使得比值法的min/max模式能很好地识别出该像素点及其邻域空间的不变性；反之，如果是异质区域，则θ取值大，在式(5)中，中心像素点强度值I(i,j)占比重大，此时min/max模式又更能细致反映各种变化类型的像素点变化情况，实现了变化细节保留。

另外一个问题是异质性测度值归一化，式(5)中θ1和θ2的理想值为[0,1]，但是当σ(x)>μ(x)时，可能出现该值大于1情况，而造成DINR(i,j)取值可能为负，无法理想地表达变化关系，可以使用下式对异质性测度值进行统一归一化：

(6)

2.2 层次FCM聚类

根据所提的INR方法进行计算得到的DI，通过FCM聚类方法[15]，将其元素分3类，即变化类Ωc、中间类Ωi和不变化类Ωu。其中，变化类Ωc和不变化类Ωu为变化概率和不变化概率高的像素点，而中间类Ωi呈现变化与否的不确定性，其可能是噪声残留的影响。进一步，由于单一的FCM聚类可能会造成中间类Ωi的比重过大，使得有用信息的丢失过大，影响最终变化检测的精度，因而采用层次FCM聚类的方法，具体方法的步骤如下：

输入：差异图矩阵DI

(7)

(8)

(5) 转至第(3)步，直到t=m后退出。

输出：层次FCM聚类结果得到表示为带有标签{Ωc,Ωi,Ωu}的DI像素集合。

经过层次FCM聚类后，我们得到了一张标有{Ωc,Ωi,Ωu}的变化图。属于Ωc的像素变化概率较大且令其值为0标签，属于Ωu的像素不变化概率较大且令其值为1标签，这两种像素可以构造ELM网络的训练样本。而属于Ωi的标签值为0.5的像素为不确定变化与否的中间类，其进一步的标签确定将由ELM网络根据空间邻域块特征向量分类而完成，这样就解决了可能的DI残留噪声影响。

2.3 ELM二分类

根据2.2节层次FCM聚类后的结果，且为了进一步解决Ωi的类别问题，将Ωc和Ωu作为两输入时相图的邻域块特征向量的样本标签构造分类器的训练样本集，训练ELM网络[16]，并将所有邻域块特征样本构造成ELM网络的测试样本集，最终可得到变化检测图。具体方法的步骤如下：

(1) 生成以Ωc和Ωu类对应像素为中心的两输入时相图的邻域块特征，并构造标签为Ωc和Ωu的训练样本集。其中，邻域块特征向量构造方法是以对应位置像素为中心，取大小为k×k的像素块构成，将两输入时相图的邻域块特征按如图2所示(图中k取3)形式转换并连接为行向量形式，作为ELM网络的样本特征向量，再加上中心像素点的类别标签，即构成ELM网络的训练样本集。同理，按同样的方法按序选取所有的像素为中心，取大小为k×k的像素块构造图像邻域块特征，即构成ELM网络的测试样本集。

图2 两输入时相图邻域特征向量的生成方法

(2) 基于训练样本集对ELM网络进行训练。

(3) 基于测试样本集，利用训练好的ELM进行二分类，最终得到变化检测图。

3 实 验

3.1 数据集描述和实验设置

实验数据集使用三组真实且已配准预处理的SAR图像数据，分别是Bern地区数据、Ottawa地区数据和Sardinia地区的两时相图及手动真值图，如图3-图5所示。

图3 Bern地区数据

图4 Ottawa地区数据

图5 Sardinia地区数据

为了验证本文算法的有效性，具体实验设置如下：

(1) 四种DI构建方法的对比：本文所提的INR方法与IR方法、MR方法和NR方法，通过其所生成的DI及直方图进行分析。

(2) 变化检测二值图对比：将四种方法生成的DI按本文设计的DI分析方法取得变化检测二值图，并进行性能的对比。

(3) 本文方法与其他典型的几种变化检测方法对比，以评估本文方法的有效性。

本文所提算法的参数设置为邻域空间x均值矩阵 大小r取值为3，层次FCM聚类的第一轮聚类数为3，第二轮聚类数m为7，聚类预设参数σ1取值为1.12，σ2取值为1.10，邻域块特征矩阵大小k取值为5，ELM分类器输入层节点数为50，隐含层节点数为10，输出层节点数为2。数据对比采用较为广泛使用的定量分析方法，假设图像总像素个数为N，则定量数据包括漏检像素数FN、误检像素数FP、总错误数OE、正确率PCC和KAPPA系数。

OE=FN+FP

(9)

PCC=(N-OE)/N

(10)

KAPPA=(Pr(a)-Pr(e))/(1-Pr(e))

(11)

式中：TP为确定变化的像素且在变化检测中被正确检测为变化的像素数；TN为确定未变化的像素且在变化检测中被正确检测为未变化的像素总数；Pr(e)=((TP+FN)×(TP+FP)+(TN+FN)×(TN+FP))/N2，Pr(a)=(TN+TP)/N。KAPPA系数值越趋近于1，表明变化检测的效果越好，根据其一致性等级划分，一般为当其取值为0.81～1时，得到的变化检测结果跟手动真值图更接近且趋于一致。

3.2 实验结果及分析

(1) 四种DI构建方法的对比：本文INR方法和IR方法、MR方法和NR方法对两组数据处理得到的DI及其直方图分别如图6-图8所示。

(a) IR方法及其直方图

(b) MR方法及其直方图

(c) NR方法及其直方图

(d) 本文INR方法及其直方图图6 Bern地区数据DI对比

(a) IR方法及其直方图

(b) MR方法及其直方图

(c) NR方法及其直方图

(d) 本文INR方法及其直方图图7 Ottawa地区数据DI对比

(a) IR方法及其直方图

(b) MR方法及其直方图

(c) NR方法及其直方图

(d) 本文INR方法及其直方图图8 Sardinia地区数据DI对比

差异图DI反映两时相SAR图像之间的变化情况，尤其需突出变化区域和轮廓，且还应有高的噪声抑制鲁棒性。在其直方图中，将会呈现两个波峰的数据分布，其中未发生变化像素数的波峰应高、集中且更靠近零灰度值附近，而发生变化像素数对应的波峰应该更靠近1灰度值附近。从得到的DI及其直方图可以看出，如图6(a)、图7(a)和图8(a)所示，由于IR方法按单像素比的min/max求得DI，未考虑邻域空间信息和区域异质性，对图像残留斑点噪声的抑制不够，细微的强度值变化也会在DI中反映出来，所得DI的直方图毛刺较多，且直方图分布也表明样本数据峰值分布不合理，尤其是噪声残留较多的图6(a)，没有明显地体现变化像素和未变化像素的区分，势必影响后续的DI分析。对于MR方法，其采用邻域均值比的min/max求得DI，能较好地抑制斑点噪声的影响，如图6(b)、图7(b)和图8(b)所示，得到的DI平滑度好，但忽略了区域异质性，对变化轮廓细节信息的丢失较为严重，因而边缘表现不好，所对应的直方图显示其数据峰值分布也不太好。对于NR方法，其结合了邻域空间信息和区域异质性测度，融合了IR方法和MR方法的优点，去噪和细节保留也较好，但从其所生成的DI直方图可以看出，如图6(c)、图7(c)和图8(c)所示，波峰数据的分布仍然存在有不太好的现象，对进一步的DI分析有影响。而本文所提的INR方法中引入异质性测度来加权考虑中心像素灰度值和邻域灰度均值信息，并以IR方法的min/max方式求得DI，从图6(d)、图7(d)和图8(d)可以看出，图像平滑性增强的基础上仍然突出保留了变化轮廓细节，其直方图显示该法数据峰值分布最好，毛刺少，这将更利于后续的DI分析处理求得变化检测图。

(2) 变化检测图对比：将上面的DI分别按本文的DI分析方法得到的变化检测二值图，如图9-图11所示，定量分析的数据对比如表1-表3所示。

图9 Bern地区变化检测二值图对比

图10 Ottawa地区变化检测二值图对比

图11 Sardinia地区变化检测二值图对比

表1 Bern地区变化检测二值图数据对比

表2 Ottawa地区变化检测二值图数据对比

表3 Sardinia地区变化检测二值图数据对比

由图9可以看出，(d)的变化检测图轮廓和细节保留最好，与真实的手动真值图最接近，(c)和(b)次之，(a)虽然轮廓和细节保留也不错，但图噪声干扰太多，说明IR方法噪声残留较多。从表1数据对比来看，IR方法误检像素数最多，是造成其正确率低的主要原因；MR方法的漏检数较高，是因为其采用邻域均值比来构建DI，增强了DI的平滑性，但变化轮廓的边缘信息丢失增多；NR方法的DI取得了不错效果，说明其引入的邻域均值信息和区域异质性测度是有意义的；而本文的INR方法得到变化检测图各项数据最好，其PCC达到99.67%，KAPPA系数达到0.866 9，且比NR方法高出1个百分点。

由图10可以看出，(d)变化细节保留最好，尤其是上半部分的变化细节更清晰，(c)次之，(a)和(b)的漏检情况较多。从表2数据可以看出，本文INR方法得到变化检测图各项数据最好，漏检数和误检数均较低，PCC达到97.03%，KAPPA系数达到0.879 6，且比NR方法高了3个多百分点。对比数据均表明本文所提的INR方法来生成DI更有利于DI分析。

根据Sardinia地区的手动真值图可以知道其变化细节要求较高，由图11可以看出，仍然取得了较好的变化细节保留，(c)和(b)次之。从表3数据可看出，虽然几种方法的误检像素数均较高，但本文INR方法得到变化检测图各项数据最好，误检数相对较低，PCC达到97.65%，KAPPA系数达到0.812 0，且比NR方法高了2个多百分点。对比数据均也表明本文所提的INR方法来生成DI更有利于DI分析。

(3) 本文方法与其他典型的几种变化检测方法的对比：这些典型的SAR变化检测方法有K-means方法、RFLICM方法、SIFT方法、S-PCANet方法和MRFFCM方法，且为实现对比公平性，这些方法的DI构建均用本文所提的方法，其他参数设置均按其文献阐述取值。三组SAR图像数据的结果对比分别如图12-图14所示，定量分析的数据对比分别如表4-表6所示。

图12 Bern地区六种检测方法的变化检测图对比

图13 Ottawa地区六种检测方法的变化检测图对比

图14 Sardinia地区六种检测方法的变化检测图对比

表4 Bern地区六种检测方法数据对比

表5 Ottawa地区六种检测方法数据对比

表6 Sardinia地区六种检测方法数据对比

从上述三组SAR图像数据经六种方法处理得到的变化检测二值图对比可以看出，在这些方法中，本文方法和S-PCANet方法表现最稳定，因为这两种方法均对DI中残留的噪声进行了抑制，所得到变化检测二值图的轮廓和变化细节均较为清晰，更趋近于手动真值图。从表4-表6的数据对比进一步看出，本文方法的PCC值均在97%以上，KAPPA值也均大于0.81，数值居于前两位，尤其是Sardinia地区检测细节要求高，本文方法也获得了不错的结果，根据其一致性等级划分，所得到的变化检测结果均是有效的，趋近于手动真值图，这表明本文设计的DI分析方法是有效的，比其他典型的几种变化检测方法更优。

4 结 语

本文针对两时相的SAR图像变化检测问题，提出了一种新的邻域比方法INR来构建DI并设计了一种新的DI组合分析方法。实验表明，通过该INR方法构建的DI能更清晰地反映两时相图的变化性质，能有效抑制斑点噪声的影响，在图像平滑性增强的基础上仍然能突出变化轮廓细节，且使得DI的像素样本分布更加合理，更利于DI分析。通过本文设计的DI组合分析方法，利用层次聚类获得样本标签、邻域块特征向量作为样本进行ELM网络分类，得到变化检测二值图。实验结果再次表明INR方法有效的同时，数据对比的PCC值和KAPPA系数值也说明了本文所设计的DI分析方法的有效性，其能对DI中残留的噪声进行抑制，变化细节保留较好，使得变化检测精度得以提升。下一步，将继续对DI中残留噪声抑制进行分析处理，考虑新的特征提取方法。