CN109919910A - 基于差异图融合和改进水平集的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于差异图融合和改进水平集的SAR图像变化检测方法，首先分别对两时刻SAR图像进行对数比运算和均值比运算得到对数比值图和均值比值图，而后融合对数比值图和均值比值图获取融合差异图；利用KI阈值对融合差异图进行预分类，获得初始标记场、两类数据集的均值及方差；根据初始标记场、两类数据集的均值及方差，计算融合差异图上每个像素点属于某一类的隶属度，进而构建自适应核函数；依据该核函数和图像的梯度信息，结合全局能量项和局部能量项控制水平集函数的曲线演化；最后根据收敛时的曲线划分变化区域和未变化区域。本发明不仅保护了边缘等细节信息，而且有效抑制了相干斑噪声，从而提高了SAR图像变化检测的精度。

Description

基于差异图融合和改进水平集的SAR图像变化检测方法

技术领域

本发明隶属于图像处理领域，主要解决合成孔径雷达(Synthetic ApertureRadar,SAR)图像变化检测问题，具体提供一种基于差异图融合和改进水平集的SAR图像变化检测方法。

背景技术

SAR图像变化检测是通过对不同时期、同一区域的两幅SAR图像进行对比分析，根据两幅图像之间的差异来获取地物或者目标在一段时间内的变化信息。由于SAR具有大范围成像、不受天气影响、全天时、全天候侦察的优点，因此SAR图像变化检测技术被广泛应用于城市规划、国土监测、战场信息感知等领域，它的广阔应用场景使得越来越多的国内外学者致力于开发高精度的SAR图像变化检测方法。

目前SAR图像变化检测方法可分为两类：先分类后比较和先比较后分类。后者思路简单、检测结果精度高，从而成为当前主流的变化检测方法，其基本流程范式分成三步：预处理、生成差异图、分析差异图。其中生成差异图和分析差异图是影响最终检测结果的关键步骤，也是SAR图像变化检测的重点研究内容。

差异图生成研究较多的方法包括基于简单代数运算的算法、基于变化矢量分析的算法、基于差异图融合的算法等。学者Bazi等人提出的对数比(log ratio detector,LRD)算法在SAR图像变化检测中的应用更为广泛，对数比算法可以有效地解决SAR图像中存在的乘性噪声，然而对数运算容易模糊边缘区域的像素值；学者Gong等人提出了融合对数比值图和均值比值(mean ratio detector,MRD)图的SAR图像变化检测算法，该算法有效地结合了两种差异图的优点。

得到差异图后，需要对差异图进行分析来得到二值变化检测结果。阈值法是通过某种阈值选择方法找出划分2类的最优阈值，其中学者Bazi等人提出利用KI(Kittler&Illingworth)阈值对差异图进行分析，这类分析方法计算复杂度低、方法简单，但存在检测精度低等问题；水平集方法是构建一个三维水平集函数，然后求解函数值为零的点构成曲线，从而获得图像分类的结果，学者Chan等人提出了基于区域的CV模型，具有良好的边缘检测能力，然而由于差异图易受噪声的影响，容易在区域中或边界上产生杂乱的噪点。

发明内容

本发明提出了一种基于差异图融合和改进水平集的SAR图像变化检测方法，该方法通过融合对数比值图和均值比值图来提升最终差异图性能，并在CV模型中引入带有先验信息的局部能量项，从而提升SAR图像变化检测方法的检测精度以及抗噪性能。

为实现上述目的，本发明所采用的技术方案如下：

基于差异图融合和改进水平集的SAR图像变化检测方法，其特征在于，包括如下步骤：

(1)输入原始图像：

输入两幅配准及校正后的同一地区不同时间的SAR图像X₁和X₂；

(2)获取差异图：

2a)分别对输入的SAR图像使用对数比值法和均值比值法得到对应的对数比值图和均值比值图，计算过程如下：

I_log＝|X₁(p)-X₂(p)|

式中I_log、I_ave为获取的对数比值图和均值比值图，p表示坐标为(i,j)的像素点，X₁(p)表示X₁中坐标为(i,j)对应像素灰度值，X₂(p)表示X₂中坐标为(i,j)对应像素灰度值，X'₁(p)和X'₂(p)分别为X₁和X₂中坐标为(i,j)对应像素的局部均值；

2b)融合对数比值图和均值比值图，得到最终差异图I，其公式如下：

I＝aI_log+bI_ave

式中I为融合差异图，a、b为融合的加权系数，a、b决定了两种差异类型在整体差异度所占的权值大小且满足a+b＝1；对于不同的原始图像，所得到的对数比值图和均值比值图表现力也不同，为了最大程度地融合差异图，需要通过实验来确定最佳的a、b值；

(3)KI阈值预分类：

利用KI阈值法推导出差异图I的最佳分类阈值T、两类数据集的均值和方差，通过阈值T对差异图中的每个像素标记为变化类w_c或者非变化类w_u，得到初始标记场W＝(w₁,w₂,…,w_n)，其中n为图像像素点数，并记w＝{w_c,w_u}；

(4)构建自适应核函数：

4a)假设数据集类别标记的先验概率符合Gibbs分布，局部概率为MLL模型，类别w的条件概率分布为高斯分布，根据步骤(3)中得到的初始标记场、类别w像素集的均值和方差，计算当前状态下每个像素点的标记场能量及特征场能量，其公式分别如下所示：

式中U(w_p)为像素点p的标记场能量，β是一个常数，用于控制MRF模型中背景信息的影响；N_p为以像素p为中心的局部邻域w_p和w_q表示相邻点像素标签，U(p|w)为像素点p的特征场能量，μ_w和σ_w分别为类别w像素集的均值和方差，I(p)是差异图I像素点p的灰度值；

4b)计算每个像素属于类别w的条件概率，确定每个像素点的隶属度，计算公式如下：

式中P(w|p)表示像素点p的隶属度，Z为局部窗口内的归一化因子；

4c)分别利用两类隶属度矩阵构建自适应核函数，例如当局部窗口大小是3×3时，以像素点p为中心的自适应核函数可以表示为：

式中K_i(p)(i＝u,c)表示由类别w_i的隶属度构建的核函数，N_p,q(q＝1,2,…,8)是以像素点p为中心的邻域,Z为局部窗口内的归一化因子；

(5)改进的水平集方法分类：

5a)使用步骤(4)中得到的核函数构建一个新的局部能量项，将局部隶属度和曲线梯度信息加入到CV水平集能量函数中，这可以对水平集曲线的演化进行约束，其水平集函数的梯度下降流方程为：

式中δ(φ)为Dirac函数，φ为水平集函数，α、β、ν和μ分别为各个能量项的权重系数，E_g和E_l分别为全局能量项和局部能量项。

其中：

E_g＝-|I-c₁|²+|I-c₂|²

E_l＝-K_c*|I-c₁|²+K_u*|I-c₂|²

β＝1-α

式中I表示融合差异图，c₁和c₂分别为差异图演化曲线C轮廓内外的灰度均值，K_u和K_c分别是未变化类和变化类的核函数，参数α、β是基于水平集函数的梯度信息自适应得出，对于靠近边界区域的像素，α的值将远小于β的值，以避免在边界区域中使用过多的局部信息；相反，如果像素具有相对低的梯度值，这意味着α的值大于β的值，则局部能量将在水平集演化中占主导地位；

5b)令K＝0，初始化水平集函数φ，定义一个初始曲线C，令φ在曲线内部为大于零的值，曲线外部为小于零的值，并计算此时曲线C内部和外部的灰度均值c₁和c₂；

5c)根据梯度下降流方程更新水平集函数，直至函数收敛，否则返回5b)；函数收敛时水平集函数φⁿ对应的零水平集即为两类区域的分界线，将变化区域的灰度值设为1，非变化区域的灰度值设为0，得到SAR图像变化检测结果图；

(6)输出变化检测结果。

本发明与现有技术相比具有如下优点：

1)本发明由于利用对数比值图和均值比值图融合得到最终差异图，有效地利用了两种差异图的互补信息，从而提高检测算法性能。

2)本发明由于在CV模型中加入带有先验信息的局部能量项，这些先验信息会引导能量函数至全局最优，使得变化检测的精度得到了提升。

实验结果表明，本发明与KI阈值法、基于马尔科夫随机场(Markov RandomField,MRF)的FCM(MRFFCM)算法和基于CV模型的水平集算法相比，本发明改善了边界等细节的保持能力，而且有效提高了检测精度和抗噪性能。

附图说明

图1为本发明方法的整体流程图；

图2为Ottawa地区两幅SAR图像及变化参考图；

图3为Bern地区两幅SAR图像及变化参考图；

图4为本发明应用于Ottawa数据集的变化检测结果图；

图5为本发明应用于Bern数据集的变化检测结果图。

具体实施方式

按照图1所示，基于差异图和改进水平集的SAR图像变化检测方法，包括以下步骤：

步骤1.输入原始图像：输入两幅配准及校正后的同一地区不同时间的SAR图像X₁和X₂；

步骤2.获取差异图：使用对数比值法和均值比值法得到对数比值图和均值比值图，而后设置融合权重系数a＝0.7、b＝0.3，得到融合差异图I；

步骤3.KI阈值预分类：对融合差异图进行KI阈值分类得到初始标记场W以及两类数据集的均值和方差；

步骤4.构建自适应核函数：根据初始标记场W计算当前状态下每个像素点的标记场能量U(w)，通过两类数据集的均值μ_w和方差σ_w(w＝w_u,w_c)计算当前状态下的每个像素点的特征场能量U(I|w)，而后利用标记场能量U(w)和特征场能量U(w|I)得到当前状态下每个像素点为类别w的条件概率P(w|I)，P(w|I)也即每个像素属于类别w的隶属度。在两类隶属度对的基础上，分别构建3×3大小窗口的核函数K_i；

步骤5.改进的水平集方法分类：对经典CV模型进行改进，在计算全局能量的同时，考虑利用核函数K_i卷积全局能量项作为带有先验信息的局部能量项，两者共同约束能量函数。确定初始演化曲线，如图4(a)和图5(a)，根据梯度下降流方程更新水平集函数，直至函数收敛。函数收敛时水平集函数φⁿ对应的零水平集即为两类区域的分界线，将变化区域的像素点灰度值设为1，非变化区域的像素点灰度值设为0，得到SAR图像变化检测结果图；

步骤6.输出变化检测结果。

通过以下两组真实SAR图像实验进一步说明本发明的有效性。

1.实验设置：

实验环境：MATLAB R2016b，Intel(R)Core(TM)i3-2130CPU3.40GHz，Window10专业版。

实验数据：图2(a)和图2(b)为Radarsat-1获取的Ottawa地区水灾前后的SAR图像，图像大小为290×350，分辨率为10m，图2(c)为真实变化参考图；图3(a)和图3(b)为ERS-2获取的Bern城市的SAR图像，图像大小为301×301,分辨率为30m，图3(c)为真实变化参考图。

实验参数：融合的加权系数a＝0.7、b＝0.3，标记场能量计算窗口N_p和核函数窗口N_p，q均为3×3，梯度下降流方程中ν＝1，μ＝0.01×255×255，Δt＝0.1。

实验对比方法：KI方法、MRFFCM和基于CV模型的水平集方法。

实验结果评价指标：漏检数(False Negatives，FN)、误检数(False Positives，FP)、正确检测率(Percentage Correct Classification，PCC)以及Kappa系数(kappacoefficient，KC)。

2.实验结果与分析：

实验1对Ottawa地区水灾前后的SAR图像分别利用KI方法、MRFFCM方法、CV方法和本发明方法进行变化检测。变化检测结果图如图4所示，为了展示水平集函数的曲线演变过程，在本发明方法进行变化检测时分别选取了迭代次数n＝0,60,150,200时的曲线演化图，如图4(a)～(d)所示，图4(e)～(f)所示分别为KI方法、MRFFCM方法、CV方法和本发明方法的变化检测结果图。

实验2分别利用上述的四种变化检测方法对Bern城市的两幅SAR图像进行变化检测。实验结果如图5所示，其中5(a)～(d)是迭代次数n＝0,60,150,200时的曲线演化图，图5(e)～(h)分别为KI方法、MRFFCM方法、CV方法和本发明方法的变化检测结果图。两组实验的定量分析见表1所列。

表1 Ottawa及Bern地区SAR图像变化检测结果

从图4、图5可以看出，本发明方法相比于其他变化检测方法能获得更好视觉效果的变化检测结果图，其结果图中不存在明显噪点，边缘保持效果良好；另一方面，由表1可知，本发明方法的检测正确率以及KC均高于其他图像变化检测方法，并且本发明方法能够平衡错检数和漏检数，使得KC明显提升，获得更好的图像变化检测结果。

Claims (1)

1.基于差异图融合和改进水平集的SAR图像变化检测方法，其特征在于，包括如下步骤：

(1)输入原始图像：

输入两幅配准及校正后的同一地区不同时间的SAR图像X₁和X₂；

(2)获取差异图：

2a)分别对输入的SAR图像使用对数比值法和均值比值法得到对应的对数比值图和均值比值图，计算过程如下：

I_log＝|X₁(p)-X₂(p)|

式中I_log、I_ave为获取的对数比值图和均值比值图，p表示坐标为(i,j)的像素点，X₁(p)表示X₁中坐标为(i,j)对应像素灰度值，X₂(p)表示X₂中坐标为(i,j)对应像素灰度值，X'₁(p)和X'₂(p)分别为X₁和X₂中坐标为(i,j)对应像素的局部均值；

2b)融合对数比值图和均值比值图，得到最终差异图I，其公式如下：

I＝aI_log+bI_ave

式中I为融合差异图，a、b为加权系数，a、b决定了两种差异类型在整体差异度所占的权值大小且满足a+b＝1；对于不同的原始图像，所得到的对数比值图和均值比值图表现力也不同，为了最大程度地融合差异图，需要通过实验来确定最佳的a、b值；

(3)KI阈值预分类：

利用KI阈值法推导出差异图I的最佳分类阈值T、两类数据集的均值和方差，通过阈值T对差异图中的每个像素标记为变化类w_c或者非变化类w_u，得到初始标记场W＝(w₁,w₂,…,w_n)，其中n为图像像素点数，并记w＝{w_c,w_u}；

(4)构建自适应核函数：

4a)假设数据集中类别标记的先验概率符合Gibbs分布，局部概率为MLL模型，类别w的条件概率分布为高斯分布，根据步骤(3)中得到的初始标记场、类别w像素集的均值及方差，计算当前状态下每个像素点的标记场能量和特征场能量，其公式分别如下所示：

式中U(w_p)为像素点p的标记场能量，β是一个常数，用于控制MRF模型中背景信息的影响；N_p为以像素p为中心的局部邻域w_p和w_q表示相邻点像素标签，U(p|w)为像素点p的特征场能量，μ_w和σ_w分别为类别w像素集的均值和方差，I(p)是差异图I像素点p的灰度值；

4b)计算每个像素属于类别w的条件概率，确定每个像素点的隶属度，计算公式如下：

式中P(w|p)表示像素点p的隶属度，Z为局部窗口内的归一化因子；

4c)分别利用两类隶属度矩阵构建自适应核函数，例如当局部窗口大小是3×3时，以像素点p为中心的自适应核函数可以表示为：

式中K_i(p)(i＝u,c)表示由类别w_i的隶属度构建的核函数，N_p,q(q＝1,2,…,8)是以像素点p为中心的邻域；

(5)改进的水平集方法分类：

5a)使用步骤(4)中得到的核函数构建一个新的局部能量项，将局部隶属度和曲线梯度信息加入到CV水平集能量函数中，这可以对水平集曲线的演化进行约束，其水平集函数的梯度下降流方程为：

式中δ(φ)为Dirac函数，φ为水平集函数，α、β、ν和μ分别为各个能量项的权重系数，e₁和e₂分别为全局能量项和局部能量项；

其中：

e₁＝-|I-c₁|²+|I-c₂|²

e₂＝-K_c*|I-c₁|²+K_u*|I-c₂|²

式中I表示融合差异图，c₁和c₂分别为差异图演化曲线C轮廓内外的灰度均值，K_u和K_c分别是未变化类和变化类的核函数，参数α、β是基于水平集函数的梯度信息自适应得出，对于靠近边界区域的像素，α的值将远小于β的值，以避免在边界区域中使用过多的局部信息；相反，如果像素具有相对低的梯度值，这意味着α的值大于β的值，则局部能量将在水平集演化中占主导地位；

5b)令K＝0，初始化水平集函数φ，定义一个初始曲线C，令φ在曲线内部为大于零的值，曲线外部为小于零的值，并计算此时曲线C内部和外部的灰度均值c₁和c₂；

5c)根据梯度下降流方程更新水平集函数，直至函数收敛，否则返回5b)；函数收敛时水平集函数φⁿ对应的零水平集即为两类区域的分界线，将变化区域的灰度值设为1，非变化区域的灰度值设为0，得到SAR图像变化检测结果图；

(6)输出变化检测结果。