基于双边加权核图割的SAR图像变化检测

丁 稀，董 张 玉*，杨 学 志

(1.合肥工业大学计算机与信息学院/工业安全与应急技术安徽省重点实验室，安徽 合肥 230601；2.合肥工业大学软件学院/智能互联网系统安徽省实验室，安徽 合肥 230601)

0 引言

SAR(Synthetic Aperture Radar)图像具有全天时、全天候监测特性，已成为遥感图像变化检测[1]、分割等应用中的主要数据源之一[2]，在灾害监测、农业、军事等领域具有极高的研究和应用价值[3]。SAR图像变化检测的核心步骤是对两时刻SAR图像生成的差异图进行分析，从而获取最终的变化二值图[4]。由于差异图分析可视为二分类问题，因此学者们多利用影像分割算法进行差异图分析，如Bazi等对差异图进行广义高斯模型建模，并基于KI阈值算法选取最优分割阈值对差异图进行二分类[5]。然而，阈值法易造成分割绝对化，对相干斑噪声严重的SAR图像难以保证分割精度。近年来，部分学者引入聚类法进行差异图分析，获取变化类和非变化类的两个聚类中心，如Gong等在模糊聚类算法中引入MRF方法，根据局部区域中心像素和邻域像素的标签关系修正能量函数，有效减少了FCM算法的耗时，并在一定程度上削弱了散斑噪声的影响[6]。此外，由于图割(graph cuts)模型具有良好的空间上下文利用能力和影像分割能力，逐渐成为SAR图像变化检测研究的热点。例如：Moser等将图割算法与MRF相结合对差异图进行分析，基于每个像素点与其周围像素点的标签类型极可能相同的假设，对差异图的分布模型进行参数估计，并设计新的MRF能量函数，最后利用图割对该能量函数进行优化，从而获取变化检测结果[7]；Salah等将核理论和图割算法相结合提出核图割(KGC)图像分割算法，在差异图上构建图结构，利用核函数将差异图转到特征空间，使传统图割的分段常数模型更契合SAR图像，最后设计能量函数完成图像分割，从而提高差异图的可分性[8]；Gong等进而提出基于局部拟合搜索模型的核图割变化检测方法，该模型逼近直方图，可为核图割提供优异的初始化结果，提升了图割算法精度[9]。

上述算法在变化检测任务中表现良好，但对空间上下文信息的挖掘程度不够，仅从像素级角度考虑邻接顶点的相似程度，并未从区域级角度考虑以像素为中心的区域块之间的关系，导致差异图分析过程中邻域信息丢失。鉴于此，本文提出双边加权核图割算法，即在均值比差异图上构建双边加权核图割模型，利用双边权重将邻接像素点的像素级信息和区域级信息结合，解决图像块邻域信息难以利用的问题，进而依托该双边加权核图割模型，提出新的图割能量函数，并利用能量最小化算法对该能量函数进行求解，获取变化二值图。

1 研究方法

本文提出的基于双边加权核图割的SAR图像变化检测算法流程(图1)为：1)将经过配准的同一地区不同时刻的两幅SAR图像X1、X2通过高斯滤波器进行降噪平滑预处理，以减少相干斑噪声的影响；2)使用均值比差异算子[10]生成差异图，利用邻域信息降低噪声对像素点的干扰，为后续差异图分析提供初始分类结果；3)根据得到的差异图构建双边加权核图割模型，图像各像素点为图割模型顶点，图像的像素级信息和区域级信息为图割模型的边权重；4)依托双边加权核图割模型设计相应的图割能量函数；5)依据网络流理论获取图割模型的最大流/最小割，从而完成SAR图像变化检测任务。

图1 基于双边加权核图割的SAR图像变化检测流程Fig.1 Flowchart of SAR image change detection based on bilateral weighted kernel graph cuts

1.1 双边加权核图割模型

图割算法本质是将图像中的每个像素与有限标签集合(变化类和非变化类[11]二元标签集合)进行一一配对。该算法将图像建模成一个带权无向图，图中顶点与图像的像素点或关键特征点一一对应，图的边对应像素间按传统四邻域、八邻域或K近邻原则联系起来的邻接边，而边上的权重则反映某顶点的连接顶点对其影响力的大小[11]。本文根据生成的差异图构建图割模型：定义P为差异图内所有像素点的集合，p为集合内某像素点，N为图像内所有八邻域像素点对(pn，pk)构成的集合，无向加权图G={V，E，ω}(其中V为图的顶点集，包含图像内的所有像素点以及源点S(变化类标签)和汇点T(非变化类标签)(式(1))；E为图的边集，用于描述顶点间的邻接关系，根据图割理论，边集可分为两部分：一是邻接像素点对N的n-links边，二是像素点与源点和汇点连接的t-links边(式(2))；ω为像素点间的相似度权重)。一般情况下，距离越近的像素彼此间的影响力越强，仅考虑单个像素的变化不具备物理意义。因此本文提出假设：空间最近的像素点对其中心像素点的影响最大，其影响不仅与自身强度相关，还与其所在的区域信息有关。依据以上假设，本文利用像素相似度权重ωI和区域相似度权重ωp计算ω(式(3))，其中，像素相似度权重ωI主要考虑两顶点间的强度差异，常利用指数函数计算(式(4))。

V=P∪{S，T}

(1)

E=N∪{{p，S}，{p，T}}

(2)

ω=ωI·ωp

(3)

ωI(pn，pk)=exp{-|log[I(pn)]-log[I(pk)]|}

(4)

式中：I(pn)表示像素点pn的灰度值；pk表示像素点pn的八邻接像素点。

利用图像块的结构信息可有效降低图像噪声，增加算法的鲁棒性，已被应用于图像去噪和分割等领域[12]。针对SAR图像变化检测过程中未充分利用图像块的区域级空间上下文信息问题，本文受非局部均值去噪算法[13]启发，提出新的区域相似度公式，即采用区域局部均值的归一化比率衡量两个成对图像块的强度相似度，以减少SAR图像相干斑噪声的影响；然而均值信息作为图像块的特征信息有一定局限性，当两个图像块均值相差不大时，均匀图像块和包含边界的图像块之间的相似性很难通过均值归一化比率衡量。为解决图像块边缘信息利用困难的问题，本文将图像块的边缘强度映射(Edge Strength Maps,ESM)矩阵[14]作为边缘细节信息引入区域相似度特征度量中，以增强边界信息精准分割能力。最终得到的区域相似度权重ωp的计算公式为：

(5)

(6)

式中：ESD(pn,pk)为像素点pn和pk的边缘相似度距离，本文利用高斯伽马窗算法[14]获取两时刻SAR图像的边缘信息，得到ESM矩阵，该矩阵各元素与像素点一一对应，值越大，对应像素点位于边缘的可能性越大，反之该像素点处于均匀区域的可能性越大；h为边缘相似度距离参数；d为区域块的半径；μ(pk)、μ(pn)分别为以像素点pk、pn为中心的区域块的平均值；H(pn)、H(pk)分别为图像ESM矩阵中以像素点pn、pk为中心的图像块，当H(pn)、H(pk)包含不同的边缘信息时，其ESM矩阵差异越大，边缘相似性距离值越大，当二者具有相同的边缘或均处于均匀区域时，其边缘相似性距离为0。

1.2 图割能量函数

在双边加权核图割模型构建完成后，两时刻SAR图像上的八邻域像素点可被连接，其对应的相似度信息可通过邻接边获取。依托以上信息，本文提出双边加权核图割能量函数F(L)，计算公式为：

F(L)=D(L)+αR(L)

(7)

(8)

(9)

(10)

式中：D(L)为图割能量函数的一元区域项，主要描述像素点属于某类标签L={Lc,Ln}的概率：Lc为变化类标签，Ln为非变化类标签，本文采用经典的分段常数图割模型[11]，假设μl为分段常数模型参数,即源点、汇点的实际值，在给定观测数据I(pn)下，使用参数μl与像素点pn的灰度值差异进行计算；R(L)为图割能量函数的二元边界项，用于衡量邻接像素点间的相似度，并对非邻接的像素点对进行不连续惩罚；α为常系数，用于平衡一元区域项和二元边界项的比重，一般取0<α<1；δ(·)为狄拉克函数。

由于SAR图像具有非线性特性，分段常数模型难以直接表示变化类与非变化类的分布，因此采用核函数将图像转换到特征空间，以此提高数据的可分性。根据Mercer定理[15]，核函数本质是数据在高维空间中的内积，因此，可以利用核函数计算能量函数的一元项。假设φ(·)是从数据观测空间到更高维(可能是无限维)特征/映射空间的非线性映射，同时引入高斯径向基核函数(RBF)，即：K(i，j)=exp(-‖i-j‖/2σ2)，则在核诱导空间中能量函数计算公式为：

(11)

‖φ(I(pn))-φ(μ)‖2=K(I(pn),I(pn))+
K(μ,μ)-2K(I(pn),μ)

(12)

该能量函数的一元项由核函数隐式映射到特征空间，二元去斑邻域项综合衡量像素相似度和区域相似度，并在区域级信息引入边缘相似度特征。至此，双边加权核图割能量函数建立完毕，其最佳切割方法可依据α-β移动交换标签优化算法[11]求得。

图2展示了图割模型结构和算法处理过程：不同颜色的圆圈表示不同灰度的像素点，点S和点T分别表示变化类标签和非变化类标签，每一顶点均与两终端顶点S和T连接，形成t-links边，表示为虚箭头，该边上的权重反映各像素点属于两类标签的概率；邻域内的顶点相连，形成n-links边，表示为实线，边上的权重反映像素点间的相似度。当图割模型建立后，为完成图像分割任务，需寻找一条割线(Cut，由虚线表示)将图中顶点分隔为两部分，而割线所包含的边称为割集；当割集所在边的权重之和最小时，即为图像的最终分割结果，该割线被称为最小割[16]。

图2 图割算法示意Fig.2 Schematic diagram of graph cuts

2 实验结果与分析

2.1 实验数据与参数设置

实验选取Ottawa和Mexico fire两组真实SAR数据，且均已经过预处理并达到亚像素级的配准精度(表1)。算法参数中，高斯滤波参数如果过大，会加深滤波程度，导致检测结果出现过平滑现象，如果过小，检测结果零碎点过多，抗噪性减弱，故本文选取常见的高斯滤波参数0.5；邻域加权斑块尺寸过大会导致算法复杂度变高和出现过平滑现象，所以本文选择较小的斑块尺寸d=3，在提高算法运行效率的同时，能保证良好的平滑降噪性能；边缘相似度参数h在合理范围内使边缘距离权重由边缘相似性距离值决定，过大或过小都会使该权重失去意义，本文设置h=0.002。

表1 SAR图像数据的先验信息Table 1 Prior information of SAR image dataset

选取基于修正马尔科夫随机场的模糊C均值算法[6](MRFFCM)、基于核的图割算法(Kernel-induced Graph Cuts,KGC)、PCA-Kmeans算法[17]和FLICM聚类算法[18]作为对比实验组，采用变化检测常用的指标[4]对各算法精度进行评价,包括：1)错检率(FPR)：变化类型的像素数占实际非变化类型像素总数的比例；2)漏检率(FNR)：遗漏的变化类型像素数占实际变化类像素总数的比例；3)总体精度(OA)：检测正确像素数占图像像素总数的比例；4)Kappa系数：综合反映变化检测的二值结果图与实际参考图的趋近程度。FPR和FNR越小，OA和Kappa系数越接近1，表明检测结果精度越高、越接近参考图[19]。

2.2 实验对比分析

2.2.1 Ottawa数据集变化检测结果分析 第一组数据集呈现了Ottawa地区洪水前后的真实地物变化情况(图3a、图3b)，从真实变化参考图(图3c)可以看出，变化区域(图中白色区域)多为线状或带状，且块状变化区域内有很多碎点和不规则边缘，因此变化检测算法要具有良好的抗噪性和准确的边缘定位能力。由5种算法检测结果(图3d-图3h)对比可知，PCA-Kmeans算法检测结果良好，但在中上部区域漏检点较多；FLICM算法碎点较多，中上部区域也存在较多漏检点，对变化区域和非变化区域表征均不理想；MRFFCM算法在非变化区域(图中黑色区域)存在大量碎点，抗噪性能不佳；KGC算法的噪声抑制能力较好，但变化区域和非变化区域边界模糊，可分性较差；本文算法合理利用了像素级和区域级空间上下文信息，抗噪性能较好且变化区域边缘定位准确。

图3 Ottawa数据集变化检测结果Fig.3 Change detection results of the Ottawa flood dataset

由5种算法在Ottawa数据集上的定量检测结果(表2)可以看出，PCA-Kmeans算法的错检率较低，但漏检率较高，影响了检测精度；FLICM和MRFFCM算法的错检率较高，表明其抗噪能力不强；KGC算法错检率较低，但漏检率较高，二值变化图出现过平滑现象，导致最终的检测结果不理想；本文算法OA和Kappa系数均最高，验证了算法的有效性。

表2 不同算法Ottawa数据集变化检测结果比较Table 2 Comparison of change detection results of the Ottawa flood dataset for various algorithms

2.2.2 Mexico fire数据集变化检测结果分析 第二组数据集呈现了Mexico地区火灾前后的真实地物变化情况(图4a、图4b)，与前文分析类似，PCA-Kmeans、FLICM算法检测结果中变化区域与非变化区域边界表现良好，但在较小变化区域表现不理想，漏检率较高；MRFFCM和KGC算法背景区域表现良好，但变化区域边界较平滑，且丢失部分较小变化区域，致使漏检率较高；本文算法(图4d)确保漏检较少，同时对变化区域边缘定位准确。由表3可以看出，本文算法的漏检率最低，总体精度OA和Kappa系数最高，进而验证了本文算法在检测精度和边缘检测能力上的优异性。

表3 不同算法Mexico fire数据集变化检测结果比较Table 3 Comparison of change detection results of the Mexico fire dataset for various algorithms

图4 Mexico fire数据集变化检测结果Fig.4 Change detection results of the Mexico fire dataset

3 结论

本文针对SAR图像固有的相干斑噪声问题，提出了双边加权核图割算法，该算法特点为：1)图割结构上，从像素级角度和区域级角度挖掘图像的空间上下文信息，并将两者结合作为图割的双边权重二元项，提升了抗噪能力；2)能量函数上，在去斑邻域项中引入边缘相似度特征，通过边缘强度映射矩阵获取成对斑块的边缘相似距离，缓解图割模型边缘细节信息丢失严重问题。与其他4种算法对比结果表明，本文算法总体精度和Kappa系数更高，边缘细节信息更丰富，且具有良好的抗噪能力。本文提出的SAR图像变化检测算法，主要利用像素的强度特征和边缘特征，尚未研究其他特征信息，今后将进一步挖掘图像多种特征，并将其引入双边核图割模型中，探究多种特征对变化检测的影响。