多尺度子空间融合谱聚类的SAR图像变化检测

张建龙，杨亚东

(西京大学 信息工程学院，西安 710123)

0 引言

图像变化检测是研究从同一地理位置、不同时间获取的2幅或者多幅图像之间发生变化的一种技术[1]。合成孔径雷达(synthetic sperture radar，SAR)是一种有源微波成像传感器[2]。SAR图像具有分辨率高、能全天候工作的特点，因此SAR图像变化检测技术被广泛应用于灾害评估、土地利用、城市规划以及军事侦察等重要领域。

SAR图像变化检测方法一般可分为2类：有监督方法和无监督方法[3]。有监督方法是基于样本的图像分类方法，分为训练和测试2个阶段。其中，训练阶段主要利用大量具有监督信息的样本完成对分类器模型的参数优化；测试阶段，要将待分类样本输入分类器以获得分类结果。由于监督信息的约束，因此相比无监督方法而言，有监督方法能够获得更高检测精度。然而，在遥感图像变化检测领域，要获取大量监督信息的样本，需海量专业人士的手工标定工作，实际中难以实现，因而该类方法的研究相对较少。无监督方法则无需监督信息，根据输入图像建立数学模型，采用聚类或者其他迭代优化方式获得变化检测结果。无监督方法不受监督信息限制，因此在实际应用中获得了更多的关注。目前，无监督变化检测方法主要有以下几种：①直接比较法。通常将2幅不同时相的SAR图像直接做差值或比值，再通过选取阈值来对像素点是否发生变化做出判断[4]。②基于图像变换的方法。慕彩红等[5]通过小波融合的方法构造差异图，并且采用PCA方法提取差异图像的特征。Li等[6]采用Gabor Wavelet方法提取输入图像的特征，并且生成差异图像。③基于图像建模的方法。根据图像的变化区域和非变化区域的不同统计信息，建立适用于变化检测的数学模型。窦方正等[7]将变化检测问题转化为图像像素二分类模型，提取图像多尺度特征，并利用深度置信网络得到分类结果。Zhou等[8]利用条件随机场(conditional random field，CRF)模型，将变化检测问题定义为一个标签问题，用以区分不同图像中的变化类和不变类。④基于图像分割的方法。张羽君[9]采用多层动态排序统计区域合并方法获取差异图像的超像素分割结果，以像素块为基本操作单元完成后续变化检测流程。Gong等[10]采用超像素分割进行变化特征提取，以实现对变化区域和非变化区域的分类。

上述SAR图像变化检测方法主要有以下不足：①遥感图像中地物尺寸、外形等物理信息混杂多变，语义提取困难，单一尺较度检测易造成误检，增大虚警率；②差异图能够将时相图像的变化信息统一映射到差异空间，不同机理的差异图能够捕获不同的变化信息，单一差异图易造成信息丢失和误匹配，导致检测性能下降。

针对这些不足，本文提出了一种基于形态多尺度子空间融合谱聚类的SAR图像变化检测方法。首先，生成多差异图丰富变化信息；其次，引入形态学扩展尺度子空间，分离不同地物的几何特征映射；最后，利用随机采样谱聚类的子空间融合算法对子空间分配权重，子空间特征信息分布被优化，从而获得最终的变化检测结果。

1 本文研究方法

本文研究方法的原理框图如图1所示。首先，利用输入2个时相的图像生成对数比值差异图和均值比值差异图；其次，采用MP(morphological profile)操作对双通道差异图像进行多尺度扩展，寻找SAR图像中不同尺寸地物的结构不变性，保留图像中变化目标的更多几何细节；最后，采用随机采样谱聚类子空间融合算法(AANSC)，实现各差异图像特征描述矩阵的最优化融合，生成变化检测结果。

图1中，X1、X2是经过配准后的2个时相的图像，分辨率均为I×J。

图1 本文研究方法的原理框图

1.1 生成差异图像

在变化检测中，差异图对最终变化检测结果影响很大，因此差异图的获取是一个关键步骤。对于SAR图像，为了消除相干斑噪声的影响，通常采用对数比值法[11]和均值比值法[12]来获得2个时相图像的差异图XLD和XMD。2种方法的算子分别如式(1)、式(2)所示。

(1)

(2)

式中：μ1和μ2分别表示图像X1和X2的局部均值。

对数比值法可以将SAR图像中的乘性噪声转换成加性噪声，可更加准确地反映背景区域(非变化区域)；均值比值法在对应像素灰度值相比之前加入了邻域信息，使得生成的差异图像比较接近实际情况，同时也降低了斑点噪声的影响。

图2(a)和图2(b)分别为Ottawa数据的对数比差异图和均值比差异图。

图2 Ottawa数据的2种模态差异图

1.2 基于MP的差异图像多尺度扩展

由于差异图像中不同尺寸地物的几何结构信息都混杂在一起，难以分离，若仅采用单一尺度的工具进行检测，则难以覆盖所有的尺度变化信息，造成检测精度下降。考虑到这一问题，本文采用MP工具对差异图像进行多尺度扩展，采用不同尺寸的结构元算子SE(structure element)对差异图像进行形态学滤波，以寻找SAR图像中不同尺寸地物的几何结构不变性，达到分离变化信息，保留更多图像变化细节，提高图像解译度的目的，从而提高了变化检测的精度。

MP是一种提取和重建图像中目标物体几何结构的有效工具。它被定义为采用不同尺寸的结构元SE对图像进行一系列开操作和闭操作的运算。

SE的特定形状可以被定义为菱形(diamond)、矩形(rectangle)或圆盘形(disk)，用于和中心像素的相邻区域相互作用[12]。在无法获得变化区域的任何形状先验信息的情况下，选择一个等向性的SE是较为合理的。如圆盘形SE，其中心像素到各个相邻像素的距离相同[13]；相反，如果可以通过调查获取变化区域的形状信息，那么就可以根据地物的几何特征选择相匹配的SE算子的形状。如在市区，矩形SE可以更好地匹配建筑物的形状。由于本文算法使用的SAR图像数据集无法获取变化区域的先验信息，因此，选择圆盘形SE对差异图像进行多尺度扩展具有较好的鲁棒性[14]。

由于图像中经常存在大小不同的物体，MP中的开操作能够抑制图像中比SE尺寸更小区域的几何特征，而闭操作能够保留大于SE尺寸区域的几何特征，因此通过使用一系列不同尺寸的SE对图像进行多尺度扩展，可以探索图像的不同空间变化信息，从而获得不同结构的最佳响应。

对于灰度图像f，开操作OR和闭操作CR的定义分别如式(3)、式(4)所示[15]。

(3)

(4)

式中：i是SE的尺寸大小；δi(·)和εi(·)分别是膨胀和腐蚀操作。通过设置一系列不同尺寸的SE，MP可以用式(5)定义。

(5)

式中：K是SE尺寸的最大值(即i的最大值)。

(6)

1.3 随机采样谱聚类子空间融合算法

聚类是一种无监督的机器学习方法。聚类根据未知标签样本的数据集内部的数据特征，将数据集划分成一组类内相似度高、类间相似度低、互不相交的子集。聚类方法广泛应用于变化检测过程中，处理差异图像以生成最终的变化检测结果。常用的聚类方法有：K-means、FCM(fuzzy c-means)、谱聚类(spectral clustering，SC)等。

其中，SC是一种从图论角度出发的聚类方法。它的主要思想是将每个数据样本看成空间中的点，这些点用边连接起来形成图。聚类的目的就是对图进行切分，找到“最佳”的切分方式：子图内的边权重和尽可能高，子图间的边权重和尽可能低[16]。谱聚类具有对样本空间的形状不敏感并且能够收敛到全局最优解的优点。

传统的SC算法仅采用单一的相似度矩阵。然而，在许多实际应用中，可能存在多个潜在有用的特征，因此可能存在多个相似度矩阵。为了扩大谱聚类的应用范围，需要通过特征选择或特征融合，将不同的相似度矩阵聚合到一个相似度矩阵中[17]。文献[17]提出了AASC (affinity aggregation spectral clustering)算法，该算法旨在寻求多个相似度矩阵的最佳组合，然后利用多重相似度来提高算法的聚类效果。由于多尺度扩展产生多个子空间，AASC算法的思想可以扩展到遥感变化检测中，这也是本文算法的重要动机。

相似度矩阵Wk(k=1，…,m)的计算公式如式(7)所示。

(7)

式中：wij表示di和dj之间的距离，计算公式如式(8)所示。

(8)

(9)

(10)

式中：Dk-Wk是Wk的拉普拉斯矩阵；αk=fTDkf。

AASC算法的求解分3个步骤。

①初始化相似度矩阵权重系数vk=1/m，m为待融合差异图的数目。

③运用K-means对fi聚类，得到最终的变化检测结果。

传统谱聚类算法的空间复杂度通常为O(n3)。本文相似度矩阵W的维度为l2×IJ×m，直接进行计算会消耗很大内存。因此，谱聚类算法无法直接应用于高分辨率图像和视频问题中。

Nyström方法[18]通过抽取部分样本点，用小样本的特征向量推导出大样本的特征向量的近似值，当数据量很大时，它可以应用于使用谱聚类算法的场景。对于给定数量的采样点，其复杂度与图像分辨率成线性关系。因此，可将Nyström方法引入AASC算法框架中，构建随机采样谱聚类的子空间融合算法，解决该算法无法处理大规模数据的问题。

2)随机采样子空间融合算法介绍。在随机采样子空间融合算法中，随机从差异图像的N个像素点所对应的特征空间中采样n个样本，那么，融合后的相似度矩阵W可表示为式(11)的形式。

(11)

式中：A∈IRn×n为采样得到的n个样本的相似度矩阵；B∈IR(N-n)×n表示n个采样点和剩余N-n个样本点对应的相似度矩阵；C∈IR(N-n)×(N-n)是剩余N-n个样本点对应的相似度矩阵。

(12)

该问题使用2步迭代方法求解指示向量f和权重向量v，求解步骤与AASC相同。

综上所述，本文方法由以下4个步骤组成。

1)利用输入图像X1、X2，根据式(1)和式(2)，计算差异图XLD和XMD。

3)将各差异图的特征描述矩阵Fk作为随机采样子空间融合算法的输入，计算得到聚类结果，并且恢复为分辨率为I×J的检测结果CM。

4)将检测结果CM与地面真值GT进行比较，得到本文方法的评价指标。

2 实验结果与分析

2.1 实验环境

本文仿真实验开发平台配置信息为Microsoft Windows7操作系统，3.30 GHz CPU，8.00 GB内存，仿真软件为Matlab R2016a开发平台。

2.2 实验数据集

选择2组SAR图像数据集进行算法验证。第1组数据集为Berne地区的2个时相的SAR图像，分别如图3(a)和图3(b)所示，图像分辨率为301像素×301像素；图3(c)为地面真值GT。第2组数据集是Ottawa地区的2个时相的SAR图像，如图4(a)和图4(b)所示，图像分辨率为290像素×350像素；图4(c)为地面真值GT。

图4 Ottawa地区的SAR图像

图3 Berne地区的SAR图像

2.3 评价指标

本文选取的评价指标有：漏检数FN(false negative)、虚检数FP(false positive)、检测总错误数OE(overall error)(式(13))、正确分类概率PCC(percentage correct classification)(式(14))、度量检测结果与参考图结果一致性的Kappa系数(式(15))。

OE=FP+FN

(13)

(14)

(15)

2.4 实验结果及分析

为验证算法性能，本文算法与以下DWT2-FLICMC[19]、MRFFCM[20]2种算法进行比较。2组数据集的检测结果如图5、图6所示。

图5 Berne地区的不同算法的检测结果

图6 Ottawa地区的不同算法的检测结果

由图5可以看出，图5(d)中绿色方框区域和图5(b)中绿色方框区域相比，漏检的像素点数有所减少，即本文算法的检测结果相较于DWT2-FLICMC算法来说，漏检率降低；图5(d)中红色方框区域和图5(c)中红色方框区域相比，虚检的像素点数大大减少，即本文算法和MRFFCM算法相比，虚警率大大降低。图5(d)和图5(a)相比，能较完整地显示变化目标并且保留图像大部分细节信息。综合考虑漏检和误检情况，本文算法的整体检测效果较好。但是，也可以看出，对于Berne数据，采用本文算法的检测结果丢失了部分细节信息，所以导致漏检率较高。原因分析为：该数据集中包含较多细小的变化区域，细节信息较为丰富。如何在检测结果中保持细节信息不丢失本身就是一个难点问题。本文算法采用MP操作对差异图像进行多尺度扩展时，MP中的开操作能够抑制图像中比SE尺寸更小区域的几何特征。因此，该图像中的细节区域受到抑制，导致最终的检测结果细节丢失，漏警率增大。

由图6可以看出，图6(d)中绿色方框区域和图6(b)中绿色方框区域相比，能更加完整地显示出变化区域，本文算法的检测结果相较于DWT2-FLICMC算法来说，漏检率大大降低；图6(d)中红色方框区域和图6(c)中红色方框区域相比，虚检的像素点数有一定减少，即本文算法和MRFFCM算法相比，虚警率有所降低。图6(d)和图6(a)相比，能较完整地显示变化目标并且保留图像大部分细节信息。考虑整体检测情况，本文算法的检测效果优于其他2个算法。

表1总结了本文算法和2种比较算法的评估指标。本文算法具有最低的FP和OE，最高的PCC和Kappa系数，以及最佳的检测结果。这不仅是因为本文算法采用MP工具对差异图像进行了多尺度扩展，将低维的差异空间映射到高维的多尺度子空间，达到了分离变化信息的目的，同时在多尺度子空间中寻找不同尺寸地物的结构不变性，提高了图像解译度，从而提高了检测精度。而且，本文算法利用随机采样谱聚类的子空间融合算法，该方法通过寻求权重融合系数的最优化组合，更好地利用各差异图之间的优势信息，解决单一差异图像易造成信息丢失和误匹配的问题，降低了检测结果中虚警率，使得该算法的抗噪性能得到明显改善，因此检测结果优于其他2种算法。

表1 算法性能指标

3 结束语

本文提出了一种基于形态多尺度子空间融合谱聚类的SAR图像变化检测方法。实验结果表明，该方法具有良好的检测效果，检测精度进一步提高，抗噪性能较好。在未来的工作中，本文将对MP进一步探索，并考虑用不同形状的形态学算子对差异图进行多尺度扩展，以提取更加丰富的变化特征。同时，本文方法仅在SAR图像数据集上进行测试，接下来的工作可以将该方法推广至多光谱及高光谱图像中，以提高算法的适应性和普遍性。