CN103839256B - 基于小波分解的多尺度水平集的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于小波分解的多尺度水平集的SAR图像变化检测算法，其实现过程为：（1）首先对两幅配准后的同一地区不同时相的SAR图像采用对数比操作算子做出差异图；（2）通过平稳小波分解（SWT）对差异图进行多层小波分解，得到不同分辨率的图像；（3）对低分辨率的图像通过水平集算法进行初步分割，所得到的分割图像的轮廓作为更高一层分辨率图像水平集算法的初始化曲线；（4）按照步骤（3）逐层进行，直到对原高分辨率图像进行水平集分割得到最终的分割图像。本发明中多尺度的应用增加了对噪声的鲁棒性，显著地提高了检测效果及检测的正确率，加速了变化检测过程。

Description

基于小波分解的多尺度水平集的SAR图像变化检测方法

技术领域

本发明属于遥感图像处理领域，涉及遥感图像的变化检测，具体提供一种基于小波分解的多尺度水平集的SAR图像变化检测方法，可用于对SAR图像变化的检测。

背景技术

与可见光、红外传感器比较，合成孔径雷达（syntheticapertureradar,SAR）具有分辨率高，全天时、全天候工作的特点，是较好的变化检测图像源。遥感图像变化检测是比较分析同一地区不同时间获得的两幅遥感图像，根据所得到的图像之间的差异进行地物变化和目标变化信息的提取。随着遥感技术的发展,SAR图像变化检测的技术需求日益广泛，目前，全球环境变化加剧，城市急速发展，洪水、地震等自然灾害频频发生，这些都需要及时掌握相关动态信息，为相关决策部门提供支持，SAR图像的种种优点为快速响应提供了技术支持和应急保障。

SAR图像变化检测主要有两种思路：（1）先分类后比较，即先对两幅配准后图像分别进行分类，通过比较分类结果得出变化部分。但该方法工作量很大，并且存在精度方面的缺陷。该方法是一种比较直观的方法，其优点是可以回避所用多时相遥感图像数据应获取环境条件不同和传感器不同所带来的辐射归一化问题，并可直接获取变化的类型，数量，位置,但此方法很难获得不同时相图像的不变信息具有相同类别的分类结果,且会夸大变化程度。同时,由于分类累积误差问题降低了变化检测精度。（2）对差异图进行分类，即先对两幅配准 图像做差异图，再对所得的差异图进行分类比较。这种方法的研究空间比较大，思路简单明确，检测精度较高，是当前比较流行的方法。很多经典和前沿的方法都是基于这种思路进行的。

国内外学者们对变化检测方法进行了大量的研究,目前SAR图像变化检测经典的方法有:（1）基于简单代数运算的变化检测，经典的变化检测方法包括图像差值法、图像比值法、对数比值法；（2）基于图像变换的变换检测，经典的变化检测方法包括主成份分析、变化向量分析法、相关分析法图像变换法；（3）基于图像分类的变化检测方法。在经典的算法基础上，针对直接代数运算没有考虑到像元的空间上下文信息，误检和漏检现象严重，近年来出现了在此基础上的基于上下文空间关系的变化检测，主要从两个方面进行改进。一方面单尺度上利用差异图像的上下文关系进行进一步的优化得到变化图，如利用马尔可夫随机场对差异图进行分析，继而求得最优变化检测结果；另一方面利用差异图像的小波多尺度信息进行优化得到变化图。

目前的SAR变化检测算法大部分为中低分辨率、象元级、单极化的SAR影像变化检测。随着SAR技术的日益成熟，图像质量逐渐提高，分辨率逐渐增强，获取数据的能力和精度越来越高，图像获取越来越便捷，SAR变化检测的前景会更加广泛。

发明内容

本发明的目的是提供一种利用小波分解多尺度水平集进行SAR图像变化检测的算法，以便降低噪声的影响，从而提高变化检测精度，同时加速了变化检测的过程。

本发明的技术方案是，基于小波分解的多尺度水平集SAR图像变化检测方法， 其特征是：包括如下步骤：

步骤101：开始基于小波分解的多尺度水平集SAR图像变化检测方法；

步骤102：使用对数比方法构造两幅不同时间同一地区的图像I₁,I₂的差异影像图DI并对其进行平稳小波分解获得不同尺度的图像；

步骤103：在尺度水平为L的图像（即最低分辨率的图像）上运行水平集算法；

步骤104：由低分辨率图像到高分辨率图像上依次进行水平集演化分割，每一次的低分辨率图像的分割结果作为下一层高分辨率图像水平集演化的初始化曲线，直到对原始图像进行水平集演化；

步骤105：获得最终的分割结果；

步骤106：结束基于小波分解的多尺度水平集SAR图像变化检测方法。

所述的步骤102，包括如下步骤：

步骤201：开始平稳小波分解SAR图像的方法；

步骤202：分别取两时相SAR图像I₁,I₂在同一位置x上的像素I₁(x)和I₂(x)，通过对数比方法构造差异图DI(x)：

DI(x)＝log(I₁(x)/I₂(x))

步骤203：设置尺度水平L，这里取L=3；

步骤204：利用平稳小波分解对所得到的对数比差异图DI(x)进行L-1层分解，从而获得不同分辨率的图像；

步骤205：结束平稳小波分解SAR图像的方法。

所述的步骤103，包括如下步骤：

步骤301：开始运行水平集算法；

步骤302：使用Kittler阈值分割算法对差异图DI(x)进行初始分割，得到初始分割曲线φ⁰，令k=0，φ=φ^k；

步骤303：以曲线φ^k为初始化曲线进行水平集演化，最小化能量函数F^MS(I,C)，得到一个轮廓曲线C来分割图像，其中F^MS(I,C)定义如下：

$F^{\mathrm{MS}}(I,C)=\underset{\mathrm{\Ω}}{\∫}{|X-I|}^2\mathrm{dxdy}+\mathrm{\λ}\underset{\mathrm{\Ω}\backslash C}{\∫}{\left|\mathrm{\ΔI}\right|}^2\mathrm{dxdy}+\mathrm{\μ}\left|C\right|$

其中λ和μ为正的常数，|C|是曲线C的长度，Ω表示图像区域，I为去除边界的图像区域，分割曲线若闭合曲线C没有位于两个同质区域的边界时，能量函数并不能达到能量最小，只有当闭合曲线C位于两个同质区域的边界时，能量函数才能达到能量最小。但是这个优化问题属于非凸问题，所以考虑一种特殊情况，I为分段常数函数，上述式子化为：

$\begin{array}{c}{F}^{\mathrm{MS}}(C,c_1,c_2)={\∫}_{\mathrm{inside}\left(C\right)}{|X-c_1|}^2\mathrm{dxdy}+\\ {\∫}_{\mathrm{outside}\left(C\right)}{|X-c_2|}^2\mathrm{dxdy}+\mathrm{\μ}\left|C\right|\end{array}$

式中式子中c₁和c₂为曲线C内部和外部的灰度均值，前两项称为全局能量项，最后一项作为调整项。为了便于上式的求解，将轮廓曲线C用函数φ刻画如下：

$\left\{\begin{array}{c}C=\&PartialD;\mathrm{\&omega;}=\{(x,y)\&Element;\mathrm{\&Omega;}:\mathrm{\&phi;}(x,y)=0\},\\ \mathrm{inside}\left(C\right)=\mathrm{\&omega;}=\{(x,y)\&Element;\mathrm{\&Omega;}:\mathrm{\&phi;}(x,y)>0\},\\ \mathrm{outside}\left(C\right)=\mathrm{\&Omega;}\backslash \stackrel{\&OverBar;}{\mathrm{\&omega;}}=\{(x,y)\&Element;\mathrm{\&Omega;}:\mathrm{\&phi;}(x,y)<0\}\end{array}\right.$

能量函数式子可以表达为：

$\begin{array}{c}{F}^{\mathrm{MS}}(\mathrm{\&phi;},c_1,c_2)={\&Integral;}_{\mathrm{\&Omega;}}{|X-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}+\\ {\&Integral;}_{\mathrm{\&Omega;}}{|X-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}+\mathrm{\&mu;}{\&Integral;}_{\mathrm{\&Omega;}}|\&dtri;H\left(\mathrm{\&phi;}\right)|\mathrm{dxdy}\end{array}$

其中H为Heaviside函数，可近似表达为H_ε：

$H_{\mathrm{\&epsiv;}}=\frac{1}{2}(1+\frac{2}{\mathrm{\&pi;}}\arctan \left(\frac{z}{\mathrm{\&epsiv;}}\right))$

因此上述最小化问题的求解可通过下面的Euler-Lagrange方程:

$\frac{\mathrm{d\&phi;}}{\mathrm{dt}}=\&dtri;H\left(\mathrm{\&phi;}\right)[\mathrm{\&mu;}\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)-{(X-c_1)}^2+{(X-c_1)}^2]$

求解时，先将φ固定不变，求与c₁和c₂相关的最小化，得到c₁ ^k和c₂ ^k，再通过下式计算φ^k+1：

${\mathrm{\&phi;}}^{k+1}={\mathrm{\&phi;}}^k-\mathrm{\&Delta;t}\frac{{\&PartialD;F}^{\mathrm{MS}}}{\&PartialD;\mathrm{\&phi;}}({\mathrm{\&phi;}}^{k-1},c_1^k,c_2^k)$

其中Δt为时间步长；

步骤304：重新初始化φ使其为曲线的符号距离函数；

步骤305：达到收敛即停止迭代过程，否则转向步骤303；

步骤306：得到最终分割图；

步骤307：结束运行水平集算法。

本发明与现有的技术相比具有如下有益效果：

1、本发明是通过平稳小波分解实现的多尺度，平稳小波分解具有时移不变形，能够有效的避免下采样带来的混叠，从而使变化区域的检测更加准确，降低了漏检率；

2、低分辨率图像上噪声和错误边界会得到一定程度的减少，因此低分辨率图像上的水平集演化分割对噪声更加鲁棒，避免了闭合曲线收敛到局部最小或错误边界上；

3、对于给定类型的输入图像，参数μ值只需要在允许的操作范围内调整，由于通过对数比获得的差异图进行了归一化处理，实验结果表明对参数μ并不敏感。

附图说明

图1是基于小波分解的多尺度水平集SAR图像变化检测主流程图；

图2是平稳小波分解的流程图；

图3是水平集算法的流程图；

图4是第一组实验仿真图，图4（a）和图4（b）的拍摄时间分别为1999.04和1999.05，大小均为301×301，图4（c）为参考图；

图5是对应图4的本发明方法（MSLSK）从低分辨率到高分辨率图像上水平集的分割结果；

图6是对应图4的YakoubBazi提出的多分辨率水平集方法（MLSK，参见BaziY.MelganiF.Al-SharaiHamedD,“UnsupervisedChangeDetectioninMultispectralRemotelySensedImageryWithLevelSet Methods,”IEEETransactionsonGeoscienceandRemoteSensing.2010,48(8):3178-3187）从低分辨率到高分辨率图像上水平集的分割结果；

图7是第二组实验仿真图，图7（a）和图7（b）的拍摄时间分别为1997.05和1997.08，大小均为290×350，图7（c）为参考图；

图8是对应图5的将本发明方法与Mumford和Shah提出的水平集模型（CVK参见D.MumfordandJ.Shah,“Optimal approximation by piecewisesmoothfunctionsandassociatedvariationalproblems,”Commun.Pure Appl.Math.,vol.42,no.5,pp.577–685,1989），以及MLSK方法进行比较的结果图。

具体实施方式

本发明提出了一种基于小波分解的多尺度水平集SAR图像变化检测方法，其特征是：（1）首先对两幅配准后的同一地区不同时相的SAR图像采用对数比操 作算子做出差异图；（2）通过平稳小波分解（SWT）对差异图进行多层小波分解，得到不同分辨率的图像；（3）对低分辨率的图像通过水平集算法进行初步分割，所得到的分割图像的轮廓作为更高一层分辨率图像水平集算法的初始化曲线；（4）按照步骤（3）逐层进行，直到对原高分辨率图像进行水平集分割得到最终的分割图像。

如图1所示。

主流程图步骤特征是：

步骤101：开始基于小波分解的多尺度水平集SAR图像变化检测方法；

步骤102：使用对数比方法构造两幅不同时间同一地区的图像I₁,I₂的差异影像图DI并对其进行平稳小波分解获得不同尺度的图像；

步骤103：在尺度为L的图像（即最低分辨率的图像）上运行水平集算法；

步骤104：由低分辨率图像到高分辨率图像上依次进行水平集演化分割，每一次的低分辨率图像的分割结果作为下一层高分辨率图像水平集演化的初始化曲线，直到对原始图像进行水平集演化；

步骤105：获得最终的分割结果；

步骤106：结束基于小波分解的多尺度水平集SAR图像变化检测方法。

如图2所示。

所述的步骤102，包括如下步骤：

步骤201：开始平稳小波分解SAR图像的方法；

步骤202：分别取两时相SAR图像I₁,I₂在同一位置x上的像素I₁(x)和I₂(x)，通过对数比方法构造差异图DI(x)：

DI(x)＝log(I₁(x)/I₂(x))

步骤203：设置尺度水平L，这里取L=3；

步骤204：利用平稳小波分解对所得到的对数比差异图DI(x)进行L-1层分解，从而获得不同分辨率的图像；

步骤205：结束平稳小波分解SAR图像的方法。

如图3所示。

所述的步骤103，包括如下步骤：

步骤301：开始运行水平集算法；

步骤302：使用Kittler阈值分割算法对差异图DI(x)进行初始分割，得到初始分割曲线φ⁰，令k=0，φ=φ^k；

步骤303：以曲线φ^k为初始化曲线进行水平集演化，最小化能量函数F^MS(I,C)，得到一个轮廓曲线C来分割图像，其中F^MS(I,C)定义如下：

$F^{\mathrm{MS}}(I,C)=\underset{\mathrm{\&Omega;}}{\&Integral;}{|X-I|}^2\mathrm{dxdy}+\mathrm{\&lambda;}\underset{\mathrm{\&Omega;}\backslash C}{\&Integral;}{\left|\mathrm{\&Delta;I}\right|}^2\mathrm{dxdy}+\mathrm{\&mu;}\left|C\right|$

其中λ和μ为正的常数，|C|是曲线C的长度，Ω表示图像区域，I为去除边界的图像区域，分割曲线若闭合曲线C没有位于两个同质区域的边界时，能量函数并不能达到能量最小，只有当闭合曲线C位于两个同质区域的边界时，能量函数才能达到能量最小。但是这个优化问题属于非凸问题，所以考虑一种特殊情况，I为分段常数函数，上述式子化为：

$\begin{array}{c}{F}^{\mathrm{MS}}(C,c_1,c_2)={\&Integral;}_{\mathrm{inside}\left(C\right)}{|X-c_1|}^2\mathrm{dxdy}+\\ {\&Integral;}_{\mathrm{outside}\left(C\right)}{|X-c_2|}^2\mathrm{dxdy}+\mathrm{\&mu;}\left|C\right|\end{array}$

式中式子中c₁和c₂为曲线C内部和外部的灰度均值，前两项称为全局能量项，最后一项作为调整项。为了便于上式的求解，将轮廓曲线C用函数φ刻画如下：

$\left\{\begin{array}{c}C=\&PartialD;\mathrm{\&omega;}=\{(x,y)\&Element;\mathrm{\&Omega;}:\mathrm{\&phi;}(x,y)=0\},\\ \mathrm{inside}\left(C\right)=\mathrm{\&omega;}=\{(x,y)\&Element;\mathrm{\&Omega;}:\mathrm{\&phi;}(x,y)>0\},\\ \mathrm{outside}\left(C\right)=\mathrm{\&Omega;}\backslash \stackrel{\&OverBar;}{\mathrm{\&omega;}}=\{(x,y)\&Element;\mathrm{\&Omega;}:\mathrm{\&phi;}(x,y)<0\}\end{array}\right.$

能量函数式子可以表达为：

$\begin{array}{c}{F}^{\mathrm{MS}}(\mathrm{\&phi;},c_1,c_2)={\&Integral;}_{\mathrm{\&Omega;}}{|X-c_1|}^{2}H\left(\mathrm{\&phi;}\right)\mathrm{dxdy}+\\ {\&Integral;}_{\mathrm{\&Omega;}}{|X-c_2|}^2(1-H\left(\mathrm{\&phi;}\right))\mathrm{dxdy}+\mathrm{\&mu;}{\&Integral;}_{\mathrm{\&Omega;}}|\&dtri;H\left(\mathrm{\&phi;}\right)|\mathrm{dxdy}\end{array}$

其中H为Heaviside函数，可近似表达为H_ε：

$H_{\mathrm{\&epsiv;}}=\frac{1}{2}(1+\frac{2}{\mathrm{\&pi;}}\arctan \left(\frac{z}{\mathrm{\&epsiv;}}\right))$

因此上述最小化问题的求解可通过下面的Euler-Lagrange方程:

$\frac{\mathrm{d\&phi;}}{\mathrm{dt}}=\&dtri;H\left(\mathrm{\&phi;}\right)[\mathrm{\&mu;}\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)-{(X-c_1)}^2+{(X-c_1)}^2]$

求解时，先将φ固定不变，求与c₁和c₂相关的最小化，得到c₁ ^k和c₂ ^k，再通过下式计算φ^k+1：

${\mathrm{\&phi;}}^{k+1}={\mathrm{\&phi;}}^k-\mathrm{\&Delta;t}\frac{{\&PartialD;F}^{\mathrm{MS}}}{\&PartialD;\mathrm{\&phi;}}({\mathrm{\&phi;}}^{k-1},c_1^k,c_2^k)$

其中Δt为时间步长；

步骤304：重新初始化φ使其为曲线的符号距离函数；

步骤305：达到收敛即停止迭代过程，否则转向步骤303；

步骤306：得到最终分割图；

步骤307：结束运行水平集算法。

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

1、仿真参数

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

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

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

③正确分类的概率PCC：PCC=(TP+TN)/(TP+FP+TN+FN)；

④衡量检测结果图与参考图一致性的Kappa系数：Kappa=(PCC-PRE)/(1-PRE),其中

$\mathrm{PRE}=\frac{(\mathrm{TP}+\mathrm{FP})\&CenterDot;\mathrm{Nc}+(\mathrm{FN}+\mathrm{TN})\&CenterDot;\mathrm{Nu}}{N^2}.$

2、仿真内容

①将本发明方法与MLSK方法进行比较，给出低分辨率到高分辨率图像上水平集的分割结果，通过第一组实验仿真图完成；

②将本发明方法与CVK方法和MLSK方法进行比较，通过第二组实验仿真图完成；

3、仿真实验结果及分析

①反映Bern城市的SAR图像如图4所示，图4（a）和图4（b）的拍摄时间分别为1999.04和1999.05，大小均为301×301，图4（c）为参考图。

采用本发明提出的MSLSK方法，从低分辨率到高分辨率图像上水平集演化分割的结果图如图5所示，其中图5（a），图5（b）和图5（c）分别表示低分辨率，中分辨率和高分辨率图像上的分割结果，可以看出该方法在低分辨率图像上已经可以得到较好的分割结果，在此基础上逐层进行的更高层分辨率图像的水平集演化进一步达到去噪的效果；由YakoubBazi提出的MSLK方法从低分辨率到高分辨率图像上水平集演化分割的结果图如图6所示，其中图6（a），图 6（b）和图6（c）分别表示低分辨率，中分辨率和高分辨率图像上的分割结果，与图5组图比较，该方法在低分辨率上细节保持的不太好，且噪声较多，从而造成在此基础上的更高分辨率图像的分割结果对噪声的鲁棒性比MSLSK方法较差。这是由于水平集演化方法的初始化曲线的重要性决定的。Bern组图的MLSK和MSLSK方法的变化检测结果性能指标比较如表1所示。

表1Bern组图的MLSK和MSLSK方法的比较

从表1中可以看出，本发明方法的漏检数较低，验证了该方法的理论正确性，虽然两种方法的正确检测率是相当的，但是本发明方法有更高的Kappa值，说明其与参考图的一致性更高。

②反映Ottawa地区水灾的SAR图像如图7所示，图7（a）和图7（b）的拍摄时间分别为1997.05和1997.08，大小均为290×350，图7（c）为参考图。

采用本发明方法与CVK和MLSK方法所得的变化检测结果图的比较如图8所示，其中图8（a），图8（b）和图8（c）分别表示采用本发明方法，CVK以及MSLK方法所得的变化检测结果图。通过与参考图像图7（c）对比，不太容易从主观上评价这三种方法的优差。但可以从三种方法的性能指标上来进行定量评价。这三种方法得到的变化检测结果分析如表2所示。

表2CVK，MSLK和MSLSK方法得到的变化检测结果分析

从表2可以看出本发明方法可以大大降低漏检率，这在变化检测的应用当中是相当重要的，同时该方法也保证了对于噪声的鲁棒性。本发明方法的正确检测率高达98.31%，由于另外两种方法，且其Kappa值为93.60%，高于CVK和MLSK方法的Kappa。因此，本发明方法漏检数FN和错检数FP相对较均衡，得到的变化检测结果图比较稳定。

Claims (2)

1.一种基于小波分解的多尺度水平集SAR图像变化检测方法，其特征是：包括如下步骤：

步骤101：开始基于小波分解的多尺度水平集SAR图像变化检测方法；

步骤102：使用对数比方法构造两幅不同时间同一地区的图像, 的差异影像图并对其进行平稳小波分解获得不同尺度的图像；

步骤103：在尺度为L的图像，即最低分辨率的图像，上运行水平集算法；

步骤104：由低分辨率图像到高分辨率图像上依次进行水平集演化分割，每一次的低分辨率图像的分割结果作为下一层高分辨率图像水平集演化的初始化曲线，直到对原始图像进行水平集演化；

步骤105：获得最终的分割结果；

步骤106：结束基于小波分解的多尺度水平集SAR图像变化检测方法；

所述的步骤103，包括如下步骤：

步骤301：开始运行水平集算法；

步骤302：使用Kittler阈值分割算法对差异图进行初始分割，得到初始分割曲线，令k=0，=；

步骤303：以曲线为初始化曲线进行水平集演化，最小化能量函数，得到一个轮廓曲线C来分割图像，其中定义如下：

其中λ 和μ为正的常数， 是曲线C的长度，Ω 表示图像区域，I为去除边界的图像区域，分割曲线C Ω；若闭合曲线C 没有位于两个同质区域的边界时，能量函数并不能达到能量最小，只有当闭合曲线C 位于两个同质区域的边界时，能量函数才能达到能量最小；但是这个优化问题属于非凸问题，所以考虑一种特殊情况，I为分段常数函数，上述式子化为：

式中式子中和为曲线C内部和外部的灰度均值，前两项称为全局能量项，最后一项作为调整项；为了便于上式的求解，将轮廓曲线C用函数刻画如下：

能量函数式子表达为：

其中H为Heaviside函数，表达为：

因此上述最小化问题的求解可通过下面的Euler-Lagrange方程:

求解时，先将固定不变，求与和相关的最小化，得到和，再通过下式计算：

其中Δt为时间步长；

步骤304：重新初始化使其为曲线的符号距离函数；

步骤305：达到收敛即停止迭代过程，否则转向步骤303；

步骤306：得到最终分割图；

步骤307：结束运行水平集算法。

2.根据权利要求1所述的基于小波分解的多尺度水平集SAR图像变化检测方法，其特征是：所述的步骤102，包括如下步骤：

步骤201：开始平稳小波分解SAR图像的方法；

步骤202：分别取两时相SAR图像, 在同一位置上的像素和，通过对数比方法构造差异图：

步骤203：设置尺度水平L，这里取L=3；

步骤204：利用平稳小波分解对所得到的对数比差异图进行L-1层分解，从而获得不同分辨率的图像；

步骤205：结束平稳小波分解SAR图像的方法。