CN102779346A - 基于改进c-v模型的sar图像变化检测方法 - Google Patents

Abstract

本发明公开一种基于改进C-V模型的SAR图像变化检测方法，主要解决现有技术检测结果精度低的问题。其实现过程为：对同一地域不同时间获取的两幅SAR图像，用对数比方法构造差异图；初始化水平集函数为符号距离函数，根据水平集函数值的正负将差异图区域分为两个区域；分别计算两个区域的灰度均值；根据水平集函数和两个区域的灰度均值，构建包含基于全局区域信息的能量函数和距离正则项的总能量函数；通过不断更新水平集函数的值获得总能量函数最小值，得到变化检测结果图。本发明用水平集函数表示轮廓曲线，通过更新水平集函数的值，使轮廓曲线朝着变化区域边界处演化，具有全局优化性和抗噪声性，提高了检测结果的精度，可用于灾害评估和环境检测。

Description

基于改进C-V模型的SAR图像变化检测方法

技术领域

本发明属于SAR图像变化检测领域，具体来说是一种基于改进C-V模型的SAR图像变化检测方法。该方法可应用于环境监测、军事侦察和灾情评估领域中，提高SAR图像变化检测结果的精确度。

背景技术

随着合成孔径雷达SAR技术的快速发展，SAR系统可以全天时、全天候获取图像数据，是较好的变化检测图像源。SAR图像变化检测是通过比较分析同一地域不同时刻获得的多时相SAR遥感影像，根据图像之间的差异来得到人们所需要的地物或目标随时间发生的变化信息。近年来逐渐出现了许多种多时相SAR图像变化检测方法，最常用的SAR图像变化检测技术主要由两步组成：首先对多时相SAR图像构造差异图，然后在差异图基础上提取变化区域。

对多时相SAR图像构造差异图是变化检测的一个关键问题，目前构造差异图的传统方法包括影像差值法、影像比值法、影像对数比值法和变化向量分析法等，其中影像对数比值法是影像比值法的扩充。SAR图像变化检测主要采用影像对数比值法来构造差异图，这种方法可以抑制将SAR图像中的乘性噪声转化成加性噪声，便于处理，提高在差异图基础上提取变化区域的准确率。

在差异图的基础上提取变化区域是变化检测的另一个关键问题。在SAR图像的变化检测方法的研究中，在差异图的基础上提取变化区域的传统方法是阈值方法。阈值方法是在差异图的基础上先确定阈值，再根据确定后的阈值将差异图中的区域分为变化类和未变化类两类，形成变化检测结果图。其中比较经典的一种方法是意大利G.Moser,S.B.Serpico等人提出的广义最小误差阈值方法即广义KI阈值方法，该方法假设变化区域与非变化区域的直方图统计分布符合广义高斯分布，通过最小化惩罚函数来自动确定阈值。该方法具有自适应性和实时性的优点，然而当差异图中变化类和非变化类的直方图有较大混叠时，变化类与非变化类的统计分布信息将不能被准确的模拟为广义高斯分布模型，用该阈值方法将会得到比较差的变化检测结果。

为了改进阈值方法，Yakoub Bazi等学者于2009年在会议文章“A VariationalLevel-set Method for Unsupervised Change Detection in Remote Sensing Images”将图像分割中的水平集方法应用于变化检测领域，提出了一种基于传统C-V模型的变化检测方法，通过演化一条初始闭合曲线将差异图分为变化类和未变化类两类区域。该方法用零水平集函数隐式的表示闭合轮廓曲线，避免了对闭合轮廓曲线演化过程的跟踪，将曲线演化转化成一个纯粹的求解偏微分方程数值解问题。初始闭合轮廓曲线可以选在图像中的任意位置，通过演化水平集函数来更新轮廓曲线的位置，经过数次演化后就可以正确的区分出目标和背景，检测出变化区域，具有全局优化特性。然而，基于传统C-V模型的变化检测方法，需要进行曲线演化过程中的初始曲线修正，即重新初始化为符号距离函数，造成了计算量增大，降低了变化检测速度，并且给检测结果带来许多伪变化信息。

发明内容

本发明的目的在于针对现有的SAR图像变化检测技术的不足，提出了一种基于改进C-V模型的SAR图像变化检测方法，以避免C-V模型变化检测方法中的重新初始化过程，提高SAR图像变化检测的效率和精度。

为实现本发明的目的，包括如下具体步骤：

(1)对同一地域不同时间获取的两幅SAR图像X₁，X₂，用对数比方法构造差异图X；

(2)初始化水平集函数φ为符号距离函数形式，根据水平集函数值的正负，将差异图X的整个图像区域Ω分为内部区域Ω₁和外部区域Ω₂；

(3)设循环次数n的初始值为n=0,设第n次循环时的水平集函数φⁿ=φ,根据当前水平集函数φⁿ将差异图X分成内部区域

和外部区域

计算内部区域

的灰度均值c₁和外部区域的灰度均值c₂；

(4)结合步骤(2)和步骤(3)构建基于全局区域信息的能量函数E_G:

$E_G=\underset{{\mathrm{\Ω}}_1^n}{\∫}{|X(x,y)-c_1|}^{2}H\left(\mathrm{\φ}\right)\mathrm{dxdy}+\underset{{\mathrm{\Ω}}_2^n}{\∫}{|X(x,y)-c_2|}^2(1-H\left(\mathrm{\φ}\right))\mathrm{dxdy},$

$+\mathrm{\μ}\underset{\mathrm{\Ω}}{\∫}{\mathrm{\δ}}_0\left(\mathrm{\φ}\right)|\▿\mathrm{\φ}|\mathrm{dxdy}$

其中，X(x，y)表示差异图X中像素点的坐标，Ω表示整个差异图区域，c₁为内部区域

的灰度均值，c₂为外部区域

的灰度均值，表示对水平集函数求梯度，

为梯度范式，μ为调节长度约束项

的权重参数，H(φ)为Heaviside函数，Dirac函数δ₀(φ)为H(φ)的导数；

(5)构建可取代重新初始化过程的距离正则项：该距离正则项中的函数

表达式为：

$p\left(\right|\▿\mathrm{\φ}\left|\right)=\left\{\begin{array}{cc}\frac{1}{{\left(2\mathrm{\π}\right)}^2}(1-\cos \left(2\mathrm{\π}\right|\▿\mathrm{\φ}\left|\right)),& |\▿\mathrm{\φ}|\≤1\\ \frac{1}{2}{\left(\right|\▿\mathrm{\φ}|-1)}^2& |\▿\mathrm{\φ}|\≥1\end{array}\right.;$

(6)根据步骤(4)中构造的全局能量函数E_G和步骤(5)中的距离正则项D，构建总能量函数E：

$E=E_G+D$

$=\underset{{\mathrm{\Ω}}_1^n}{\∫}{|X(x,y)-c_1|}^{2}H\left(\mathrm{\φ}\right)\mathrm{dxdy}+\underset{{\mathrm{\Ω}}_2^n}{\∫}{|X(x,y)-c_2|}^2(1-H\left(\mathrm{\φ}\right))\mathrm{dxdy},$

$+\mathrm{\μ}\underset{\mathrm{\Ω}}{\∫}{\mathrm{\δ}}_0\left(\mathrm{\φ}\right)|\▿\mathrm{\φ}|\mathrm{dxdy}+{\mathrm{\μ}}_0{\∫}_{\mathrm{\Ω}}p\left(\right|\▿\mathrm{\φ}\left|\right)\mathrm{dxdy}$

其中μ为调节长度约束项的权重参数，μ₀为距离正则项D的权重参数。

(7)通过最小化总能量函数E更新水平集函数φ，获得新的水平集函数φⁿ⁺¹，用该水平集函数φⁿ⁺¹对差异图X进行分割；

(8)判断当前的水平集函数φⁿ⁺¹能否使总能量函数E达到最小值，如果能，停止迭代,此时的水平集函数φⁿ⁺¹对应的零水平集就是变化区域的边界，将变化区域的像素灰度值赋值1，非变化区域的像素灰度值赋值0，得到变化检测结果图，否则,返回步骤(3)，用当前的水平集函数φⁿ⁺¹替代步骤(3)中的φⁿ，继续循环迭代。

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

1、本发明采用的基于全局区域信息的水平集方法，由于该水平集方法的区域拟合能量项的驱动力非常强，曲线演化速度非常快，因而提高了SAR图像变化检测的检测速度，增强了抗噪声能力，使变化检测结果达到了亚像素级别，提高了SAR图像变化检测的检测精度。

2、本发明由于采用距离正则项来维持水平集函数演化过程的平稳性，确保水平集函数演化时不偏离符号距离函数，避免了传统C-V模型的重新初始化过程，消除了重新初始化过程给检测结果带来的伪变化信息，提高了变化检测结果的精确度。

3、本发明由于采用了水平集函数隐式表示图像区域中的闭合轮廓曲线，避免了对闭合轮廓曲线演化过程的跟踪，同时由于将基于差异图的全局区域信息的轮廓曲线演化过程转化成一个纯粹的求解偏微分方程的数值解问题，具有全局优化性，即使差异图中的变化区域的边缘呈模糊或离散状，仍然可以获得理想的变化检测结果。

附图说明

图1是本发明方法的实现流程图；

图2是本发明仿真使用的Ottawa地区两幅SAR图像；

图3是用对数比方法对Ottawa地区两幅SAR图像构造的差异图；

图4是用本发明方法和现有的广义KI阈值方法、基于C-V模型的方法对Ottawa地区SAR图像进行变化检测的仿真结果图；

图5是本发明仿真使用的黄河入海口地区的两幅SAR图像；

图6是用现有对数比方法对黄河入海口地区两幅SAR图像构造的差异图；

图7是用本发明方法和现有的广义KI阈值方法、基于C-V模型的方法对黄河入海口地区SAR图像进行变化检测的仿真结果图。

具体实施方式

以下参照附图，对本发明的实现方案及优点进行详细描述。

实施例1,

参照图1，本实例的实现步骤如下：

步骤1，对Ottawa地区不同时间获取的如图2所示的两幅SAR图像X₁，X₂，用对数比方法构造如图3所示的差异图X，其中图2(a)表示SAR图像X₁，图2图(b)表示SAR图像X₂。

将如图2(a)所示的图像X₁中位于i行j列的像素点(i,j)的灰度值X₁(i，j)和相对应的图2(b)所示的图像X₂中位于i行j列的像素点的灰度值X₂(i，j)，通过对数比运算

得到差异图X中位于i行j列的像素点(i,j)的灰度值X(i，j)，若X(i，j)为0，则表示图像X₁中的该像素点没有随时间的变化而发生变化，否则，认为图像X₁中的该像素点发生了变化；对图像X₁和图像X₂中每个位于i行j列的像素点的灰度值从左到右，从上到下都进行对数比运算，得到如图3所示的差异图X。

步骤2，将水平集函数φ初始化为符号距离函数形式，根据水平集函数值的正负，将差异图X的整个图像区域Ω分为内部区域Ω₁和外部区域Ω₂。

在差异图X上做一个矩形轮廓曲线，初始化水平集函数φ＝±d。当φ>0时，表示矩形轮廓曲线的内部区域Ω₁，当φ<0时，表示矩形轮廓曲线的外部区域Ω₂，d表示差异图X中的像素点到矩形轮廓曲线的欧式距离，当φ=0时，表示差异图中的像素点在矩形轮廓曲线上。

步骤3，设循环次数n的初始值为n=0,设第n次循环时的水平集函数φⁿ=φ,根据当前水平集函数φⁿ将差异图X分成内部区域

和外部区域

通过统计方法计算内部区域

的灰度均值c₁和外部区域

的灰度均值c₂：

$c_1=\frac{\underset{{\mathrm{\&Omega;}}_1^n}{\&Integral;}X(x,y)\mathrm{dxdy}}{\underset{{\mathrm{\&Omega;}}_1^n}{\&Integral;}\mathrm{dxdy}},$ $c_2=\frac{\underset{{\mathrm{\&Omega;}}_2^n}{\&Integral;}X(x,y)\mathrm{dxdy}}{\underset{{\mathrm{\&Omega;}}_2^n}{\&Integral;}\mathrm{dxdy}},$

其中，X(x，y)为差异图X中像素点的坐标。

步骤4，结合步骤2和步骤3构建基于全局区域信息的能量函数E_G为：

$E_G=\underset{{\mathrm{\&Omega;}}_1}{\&Integral;}{|X(x,y)-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}+\underset{{\mathrm{\&Omega;}}_2}{\&Integral;}{|X(x,y)-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy},$

$+\mathrm{\&mu;}\underset{\mathrm{\&Omega;}}{\&Integral;}{\mathrm{\&delta;}}_0\left(\mathrm{\&phi;}\right)|\&dtri;\mathrm{\&phi;}|\mathrm{dxdy}$

其中，

表示对水平集函数求梯度，

为梯度范式，μ为调节长度约束项

的权重参数，H(φ)为Heaviside函数，δ₀(φ)为Dirac函数。

步骤5，构建可取代重新初始化过程的距离正则项

该距离正则项中的函数

表达式为：

$p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)=\left\{\begin{array}{cc}\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}(1-\cos \left(2\mathrm{\&pi;}\right|\&dtri;\mathrm{\&phi;}\left|\right)),& |\&dtri;\mathrm{\&phi;}|\&le;1\\ \frac{1}{2}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^2& |\&dtri;\mathrm{\&phi;}|\&GreaterEqual;1\end{array}\right..$

传统的C-V模型水平集方法中，为了保证轮廓曲线演化稳定，需要对水平集函数φ进行重新初始化处理，否则演化过程会产生剧烈的震荡，导致图像分割效果不理想。但重新初始化过程会造成最后演化得到的轮廓曲线偏移边界处的副作用，且计算量巨大。为了避免重新初始化的问题，同时又保证轮廓曲线演化过程的稳定，本发明引入距离正则项替代重新初始化过程。该项能够使水平集函数φ自动保持为近似的符号距离函数，有效控制水平集函数φ与符号距离函数的偏差。

距离正则项是李春明等学者于2010年在期刊文章“Distance Regularized Level SetEvolution and Its Application to Image Segmentation”中提出的，并将其应用于基于图像边缘信息的测地线活动轮廓模型(GAC)中，替代该模型的曲线演化过程中的重新初始化过程。本发明将距离正则项添加到传统C-V模型的能量函数中，避免重新初始化过程，从而提高变化检测结果的精确度。

步骤6，根据步骤4中构造的全局能量函数E_G和步骤5中的距离正则项D，构建总能量函数E：

$E=E_G+D$

$=\underset{{\mathrm{\&Omega;}}_1^n}{\&Integral;}{|X(x,y)-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}+\underset{{\mathrm{\&Omega;}}_2^n}{\&Integral;}{|X(x,y)-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy},$

$+\mathrm{\&mu;}\underset{\mathrm{\&Omega;}}{\&Integral;}{\mathrm{\&delta;}}_0\left(\mathrm{\&phi;}\right)|\&dtri;\mathrm{\&phi;}|\mathrm{dxdy}+{\mathrm{\&mu;}}_0{\&Integral;}_{\mathrm{\&Omega;}}p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)\mathrm{dxdy}$

其中μ为调节长度约束项的权重参数，μ₀为距离正则项D的权重参数。

步骤7，根据步骤6构造的总能量函数E演化水平集函数φ，获得新的水平集函数φⁿ⁺¹，用该水平集函数φⁿ⁺¹对差异图X进行分割。

(7.1)对总能量函数E应用变分法

得到水平集函数φ的梯度下降流方程

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}={\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)[\mathrm{\&mu;}\&times;\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)-{(X(x,y)-c_1)}^2+{(X(x,y)-c_2)}^2]+{\mathrm{\&mu;}}_0\&times;\mathrm{div}(d_p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)\&dtri;\mathrm{\&phi;}),$

其中，

为总能量函数E对φ的偏导数，div(·)代表散度算子，δ_ε(φ)是Dirac函数δ₀(φ)的正则化形式，

为中间变量，表达式为：

$d_p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)=\frac{p^{\&prime;}\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)}{|\&dtri;\mathrm{\&phi;}|}=\left\{\begin{array}{cc}\frac{1}{2\mathrm{\&pi;}|\&dtri;\mathrm{\&phi;}|}\sin \left(2\mathrm{\&pi;}\right|\&dtri;\mathrm{\&phi;}\left|\right),& |\&dtri;\mathrm{\&phi;}|\&le;1\\ 1-\frac{1}{|\&dtri;\mathrm{\&phi;}|}& |\&dtri;\mathrm{\&phi;}|\&GreaterEqual;1\end{array}\right.,$

其中，表示对水平集函数求梯度，

为梯度范式，

表示对函数

求导；

(7.2)根据梯度下降流方程

的离散化公式

计算新的水平集函数其中φⁿ⁺¹表示第n+1次循环迭代后的水平集函数，φⁿ表示第n次循环迭代后的水平集函数，Δt为时间步长；

(7.3)用水平集函数φⁿ⁺¹对差异图X进行分割，由φⁿ⁺¹的正负值得到新的内部区域

和新的外部区域

即φⁿ⁺¹<0,表示分割轮廓曲线的新的内部区域

φⁿ⁺¹>0,表示分割轮廓曲线的新的外部区域φⁿ⁺¹=0，表示此时的分割轮廓曲线。

步骤8，判断当前的水平集函数φⁿ⁺¹能否使总能量函数E达到最小值，如果能，停止迭代,此时的水平集函数φⁿ⁺¹对应的零水平集就是变化区域的边界，将变化区域的像素灰度值赋值1，非变化区域的像素灰度值赋值0，得到如图4(c)所示的变化检测结果图，否则,返回步骤(3)，用当前的水平集函数φⁿ⁺¹替代步骤(3)中的φⁿ，继续循环迭代。

实施例2

参照图1，本实例的实现步骤如下：

步骤一，对黄河入海口地区在不同时间获取的如图5所示的两幅SAR图像X＇₁，X'₂，用对数比方法构造如图6所示的差异图X'，其中图5(a)表示SAR图像X＇₁，图5图(b)表示SAR图像X'₂，具体实施步骤如实施例1中的步骤1。

步骤二，初始化水平集函数φ为符号距离函数形式，根据水平集函数值的正负，将差异图X'的整个图像区域Ω分为内部区域Ω₁和外部区域Ω₂，具体实施步骤如实施例1中的步骤2。

步骤三，设循环次数n的初始值为n=0,设第n次循环时的水平集函数φⁿ=φ,根据当前水平集函数φⁿ将差异图X分成内部区域

和外部区域

计算内部区域

的灰度均值c₁和外部区域

的灰度均值c₂，具体计算方法如实施例1中的步骤3。

步骤四，结合步骤二和步骤三构建基于全局区域信息的能量函数E_G:

$E_G=\underset{{\mathrm{\&Omega;}}_1^n}{\&Integral;}{|X(x,y)-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}+\underset{{\mathrm{\&Omega;}}_2^n}{\&Integral;}{|X(x,y)-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy},$

$+\mathrm{\&mu;}\underset{\mathrm{\&Omega;}}{\&Integral;}{\mathrm{\&delta;}}_0\left(\mathrm{\&phi;}\right)|\&dtri;\mathrm{\&phi;}|\mathrm{dxdy}$

其中，X(x，y)表示差异图X中像素点的坐标，Ω表示整个差异图区域，c₁为内部区域

的灰度均值，c₂为外部区域

的灰度均值，

表示对水平集函数求梯度，

为梯度范式，μ为调节长度约束项

的权重参数，H(φ)为Heaviside函数，Dirac函数δ₀(φ)为H(φ)的导数。

步骤五，构建可取代重新初始化过程的距离正则项：

该距离正则项中的函数表达式为：

$p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)=\left\{\begin{array}{cc}\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}(1-\cos \left(2\mathrm{\&pi;}\right|\&dtri;\mathrm{\&phi;}\left|\right)),& |\&dtri;\mathrm{\&phi;}|\&le;1\\ \frac{1}{2}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^2& |\&dtri;\mathrm{\&phi;}|\&GreaterEqual;1\end{array}\right..$

步骤六，根据步骤四中构造的全局能量函数E_G和步骤五中的距离正则项D，构建总能量函数E：

$E=E_G+D$

$=\underset{{\mathrm{\&Omega;}}_1^n}{\&Integral;}{|X(x,y)-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}+\underset{{\mathrm{\&Omega;}}_2^n}{\&Integral;}{|X(x,y)-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy},$

$+\mathrm{\&mu;}\underset{\mathrm{\&Omega;}}{\&Integral;}{\mathrm{\&delta;}}_0\left(\mathrm{\&phi;}\right)|\&dtri;\mathrm{\&phi;}|\mathrm{dxdy}+{\mathrm{\&mu;}}_0{\&Integral;}_{\mathrm{\&Omega;}}p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)\mathrm{dxdy}$

其中μ为调节长度约束项

的权重参数，μ₀为距离正则项D的权重参数。

步骤七，通过最小化总能量函数E更新水平集函数φ，获得新的水平集函数φⁿ⁺¹，用该水平集函数φⁿ⁺¹对差异图X进行分割，具体实施步骤如实施例1中的步骤7。

步骤八，判断当前的水平集函数φⁿ⁺¹能否使总能量函数E达到最小值，如果能，停止迭代,此时的水平集函数φⁿ⁺¹对应的零水平集就是变化区域的边界，将变化区域的像素灰度值赋值1，非变化区域的像素灰度值赋值0，得到如图7(c)所示的变化检测结果图，否则,返回步骤(3)，用当前的水平集函数φⁿ⁺¹替代步骤(3)中的φⁿ，继续循环迭代。

本发明的效果可以通过以下仿真进一步说明：

1、仿真条件：

本发明的仿真是在主频1.87GHZ的Intel Pentium CPU P6000、内存2.00GB的硬件环境和MATLAB R2009a的软件环境下进行的。

2、仿真参数

对于具有参考图的实验仿真图，可进行定量的变化检测结果分析：

①漏检个数：统计实验结果图中发生变化区域的像素个数，与参考图中变化区域的像素个数进行对比，把参考图中发生变化但实验结果图中检测为未变化的像素个数，称为漏检数；

②误检个数：统计实验结果图中未发生变化区域的像素个数，与参考图中未发生变化区域的像素个数进行对比，把参考图中未发生变化但实验结果图中检测为变化的像素个数，称为误检个数；

③总错误个数：漏检个数和误检个数的和；

④正确率：

3、仿真实验内容与结果分析

先用对数比运算构造出差异图，然后用本发明方法与现有的广义KI阈值方法、基于传统C-V模型的方法分割差异图，得到变化检测结果图，并对照标准参考图，定量分析检测结果，通过两组真实SAR图像的仿真实验完成。

实验1.Ottawa地区SAR图像数据

这里将本发明方法应用在如图2所示的Ottawa地区SAR图像上，并与现有的广义KI阈值方法、基于传统C-V模型的方法进行对比。

如图2所示的Ottawa地区两幅SAR图像，大小均为290×350像素，其中图2(a)表示1997年7月洪水季节Ottawa地区的地貌信息；图2(b)1997年8月洪水过后Ottawa地区的地貌信息。用图2(a)和图2(b)所示的两幅SAR图像用对数比运算构造差异图，如图3所示。

分别用现有的广义KI阈值方法、基于传统C-V的模型方法和本发明方法对图3所示的差异图进行分割，得仿真结果图，并与仿真结果的标准参考图进行对比，如图4所示。其中图4(a)表示广义KI阈值方法的仿真结果图；图4(b)表示基于传统C-V模型的方法的仿真结果图；图4(c)表示本发明方法的仿真结果图;图4(d)表示仿真结果的标准参考图。

从图4的仿真结果图可以看出，基于传统C-V模型的方法和广义KI阈值方法的仿真结果图都出现了一些白色的小斑点，造成误检数的增加，而本发明方法的仿真结果图轮廓清晰，消除了基于传统C-V模型方法仿真结果图中的噪点，提高了变化检测结果的精确度。

将图4所示的三种方法仿真结果的性能指标进行对比，如表1所示：

表1Ottawa地区各种算法变化检测结果

从表1中可以看出，广义KI阈值方法的仿真结果总错误个数为5308,传统C-V模型的仿真结果总错误个数3752，本发明的仿真结果总错误个数只有1766，很大程度上减少了检测结果的总错误个数，提高了变化检测结果的正确率，达到98.26%。本发明方法的仿真结果误检个数774处于三种方法的中间水平，比较均衡，漏检个数最少，显示了本发明方法的优越性。

实验2.黄河入海口地区SAR图像数据

这里将本发明方法应用在如图5所示的黄河入海口地区SAR图像上，并与现有的广义KI阈值方法、基于传统C-V模型的方法进行对比。

如图5所示的黄河入海口地区两幅SAR图像，大小均为257×289像素，其中图5(a)表示2008年6月该地区的地貌信息；图5(b)2009年6月该地区的地貌信息。用图5(a)和图5(b)所示的两幅SAR图像用对数比运算构造差异图，如图6所示。

分别用现有的广义KI阈值方法、基于传统C-V的模型方法和本发明方法对图3所示的差异图进行分割，得仿真结果图，并与仿真结果的标准参考图进行对比，如图7所示。其中图7(a)表示广义KI阈值方法的仿真结果图；图7(b)表示基于传统C-V模型的方法的仿真结果图；图7(c)表示本发明方法的仿真结果图;图7(d)表示仿真结果的标准参考图。

从仿真结果图7可以看出，广义KI阈值方法的仿真结果图漏检现象比较严重，得到的检测结果较差，本发明方法的仿真结果图轮廓很清晰，而且消除了基于传统C-V模型的方法仿真结果图中的一些噪点，对变化区域的检测最为准确。

将图7所示的三种方法仿真结果的性能指标进行对比，如表2所示：

表2黄河入海口地区各种算法变化检测结果

从表2中可以看出，本发明方法的仿真结果总错误个数3306，相比其他对比方法是最少的，广义KI阈值方法的仿真结果总错误个数7282最高，漏检现象比较严重，而本发明方法的仿真结果只有2060个像素点的漏检。基于传统C-V模型的方法仿真结果出现了较高的误检现象，改进后的本发明方法一定程度上减少了误检像素点的个数。本发明方法的仿真结果正确率达到了95.55%，高于广义KI阈值方法的94.77%和基于传统C-V模型方法的94.84%，显示了本发明方法的优越性。

Claims (5)

1.一种基于改进C-V模型的SAR图像变化检测方法，包括如下步骤：

(1)对同一地域不同时间获取的两幅SAR图像X₁，X₂，用对数比方法构造差异图X；

(2)初始化水平集函数φ为符号距离函数形式，根据水平集函数值的正负，将差异图X的整个图像区域Ω分为内部区域Ω₁和外部区域Ω₂；

(3)设循环次数n的初始值为n=0,设第n次循环时的水平集函数φⁿ=φ,根据当前水平集函数φⁿ将差异图X分成内部区域

和外部区域计算内部区域的灰度均值c₁和外部区域

的灰度均值c₂；

(4)结合步骤(2)和步骤(3)构建基于全局区域信息的能量函数E_G:

$E_G=\underset{{\mathrm{\&Omega;}}_1^n}{\&Integral;}{|X(x,y)-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}+\underset{{\mathrm{\&Omega;}}_2^n}{\&Integral;}{|X(x,y)-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy},$

$+\mathrm{\&mu;}\underset{\mathrm{\&Omega;}}{\&Integral;}{\mathrm{\&delta;}}_0\left(\mathrm{\&phi;}\right)|\&dtri;\mathrm{\&phi;}|\mathrm{dxdy}$

其中，X(x，y)表示差异图X中像素点的坐标，Ω表示整个差异图区域，c₁为内部区域的灰度均值，c₂为外部区域的灰度均值，

表示对水平集函数求梯度，

为梯度范式，μ为调节长度约束项

的权重参数，H(φ)为Heaviside函数，Dirac函数δ₀(φ)为H(φ)的导数；

(5)构建可取代重新初始化过程的距离正则项：该距离正则项中的函数

表达式为：

$p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)=\left\{\begin{array}{cc}\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}(1-\cos \left(2\mathrm{\&pi;}\right|\&dtri;\mathrm{\&phi;}\left|\right)),& |\&dtri;\mathrm{\&phi;}|\&le;1\\ \frac{1}{2}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^2& |\&dtri;\mathrm{\&phi;}|\&GreaterEqual;1\end{array}\right.;$

(6)根据步骤(4)中构造的全局能量函数E_G和步骤(5)中的距离正则项D，构建总能量函数E：

$E=E_G+D$

$=\underset{{\mathrm{\&Omega;}}_1^n}{\&Integral;}{|X(x,y)-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}+\underset{{\mathrm{\&Omega;}}_2^n}{\&Integral;}{|X(x,y)-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy},$

$+\mathrm{\&mu;}\underset{\mathrm{\&Omega;}}{\&Integral;}{\mathrm{\&delta;}}_0\left(\mathrm{\&phi;}\right)|\&dtri;\mathrm{\&phi;}|\mathrm{dxdy}+{\mathrm{\&mu;}}_0{\&Integral;}_{\mathrm{\&Omega;}}p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)\mathrm{dxdy}$

其中μ为调节长度约束项

的权重参数，μ₀为距离正则项D的权重参数。

(7)通过最小化总能量函数E更新水平集函数φ，获得新的水平集函数φⁿ⁺¹，用该水平集函数φⁿ⁺¹对差异图X进行分割；

(8)判断当前的水平集函数φⁿ⁺¹能否使总能量函数E达到最小值，如果能，停止迭代,此时的水平集函数φⁿ⁺¹对应的零水平集就是变化区域的边界，将变化区域的像素灰度值赋值1，非变化区域的像素灰度值赋值0，得到变化检测结果图，否则,返回步骤(3)，用当前的水平集函数φⁿ⁺¹替代步骤(3)中的φⁿ，继续循环迭代。

2.根据权利要求1所述的基于改进C-V模型的SAR图像变化检测方法，其中步骤(2)所述的初始化水平集函数φ为符号距离函数形式，是在差异图X上做一个矩形轮廓曲线，初始化水平集函数φ＝±d，当φ>0时，表示矩形轮廓曲线的内部区域Ω₁，当φ<0时，表示矩形轮廓曲线的外部区域Ω₂，d表示差异图中的像素点到矩形轮廓曲线的欧式距离，φ=0，表示差异图中的像素点在矩形轮廓曲线上。

3.根据权利要求1所述的基于改进C-V模型的SAR图像变化检测方法，其中步骤(3)中所述的计算内部区域

的灰度均值c₁和外部区域

的灰度均值c₂，通过统计方法中计算SAR图像灰度均值计算公式进行：

$c_1=\frac{\underset{{\mathrm{\&Omega;}}_1^n}{\&Integral;}X(x,y)\mathrm{dxdy}}{\underset{{\mathrm{\&Omega;}}_1^n}{\&Integral;}\mathrm{dxdy}},$ $c_2=\frac{\underset{{\mathrm{\&Omega;}}_2^n}{\&Integral;}X(x,y)\mathrm{dxdy}}{\underset{{\mathrm{\&Omega;}}_2^n}{\&Integral;}\mathrm{dxdy}}$

其中，X(x，y)为差异图X中像素点的坐标。

4.根据权利要求1所述的基于改进C-V模型的SAR图像变化检测方法，其中步骤(7)中所述的通过最小化总能量函数E更新水平集函数φ，获得新的水平集函数φⁿ⁺¹，通过以下步骤进行：

(7a)对总能量函数E应用变分法

得到水平集函数φ的梯度下降流方程

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}:$

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}={\&PartialD;}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)[\mathrm{\&mu;}\&times;\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)-{(X(x,y)-c_1)}^2+{(X(x,y)-c_2)}^2]+{\mathrm{\&mu;}}_0\&times;\mathrm{div}(d_p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)\&dtri;\mathrm{\&phi;})$

其中，

为总能量函数E对φ的偏导数，div(·)代表散度算子，δ_ε(φ)是Dirac函数δ₀(φ)的正则化形式，

为中间变量，表达式为：

$d_p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)=\frac{p^{\&prime;}\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)}{|\&dtri;\mathrm{\&phi;}|}=\left\{\begin{array}{cc}\frac{1}{2\mathrm{\&pi;}|\&dtri;\mathrm{\&phi;}|}\sin \left(2\mathrm{\&pi;}\right|\&dtri;\mathrm{\&phi;}\left|\right),& |\&dtri;\mathrm{\&phi;}|\&le;1\\ 1-\frac{1}{|\&dtri;\mathrm{\&phi;}|}& |\&dtri;\mathrm{\&phi;}|\&GreaterEqual;1\end{array}\right.,$

其中，

表示对水平集函数求梯度，

为梯度范式，表示对函数

求导；

(7b)根据梯度下降流方程

的离散化公式

计算新的水平集函数

其中φⁿ⁺¹表示第n+1次循环迭代后的水平集函数，φⁿ表示第n次循环迭代后的水平集函数，Δt为时间步长。

5.根据权利要求1所述的基于改进C-V模型的SAR图像变化检测方法，其中步骤(7)所述的用水平集函数φⁿ⁺¹对差异图X进行分割，是由φⁿ⁺¹的正负值得到新的内部区域

和新的外部区域即φⁿ⁺¹<0,表示分割轮廓曲线的新的内部区域

φⁿ⁺¹>0,表示分割轮廓曲线的新的外部区域

φⁿ⁺¹=0，表示此时的分割轮廓曲线。

Images