CN102800074B - 基于轮廓波变换的sar图像变化检测差异图生成方法 - Google Patents

Abstract

本发明公开一种基于轮廓波变换的SAR图像变化检测差异图生成方法。实现过程主要包括：先对两幅经预处理的不同时间相同地域的SAR图像构造对数比值图像和均值比值图像；用Contourlet变换处理使其生成对应的Contourlet系数；依不同融合规则分别对两幅图像系数分高频和低频进行计算；对融合后的Contourlet系数进行逆Contourlet变换生成变化差异图。本发明利用均值比值图像和对数比值图像分别提取其高频和低频的不同特点，基于Contourlet变换的图像融合充分利用了各源图像互补信息，使SAR图像变化检测得到更好结果，获得了较低的检测错误率，抑制图像噪声，提高分析精度。与其它差异图生成方法对比，本发明噪声抑制力强，边缘保持佳，变化信息得以最大程度保留。

Description

基于轮廓波变换的SAR图像变化检测差异图生成方法

技术领域

本发明属于SAR图像变化检测技术领域，主要涉及SAR图像变化检测中的差异图生成阶段。具体是一种基于轮廓波(Contourlet)变换的SAR图像变化检测差异图生成方法。用于SAR图像变化检测。

背景技术

遥感图像变化检测是指通过对不同时期同一区域的遥感图像进行比较分析，根据图像之间的差异得到变化信息。与可见光和红外遥感相比，SAR系统具有全天时、全天候获取数据的能力，并对地物有一定穿透能力，所以SAR图像变化检测技术正广泛的应用于各个领域，例如环境监测、土地利用/覆盖、农业调查、城市变化分析、军事侦察和打击效果评估等方面。

SAR图像变化检测中，差异图的生成是重要的步骤，对经过预处理(包括图像配准和辐射校正)的SAR图像进行算术运算，主要有差值运算，比值运算，对数比值运算得到初步的差异信息图，通过对信息图进行图像分割就可以得到最终的二值图像。

在现有的差异图生成图像方法里，算术方法最为常用，差值法，比值法及其衍生的对数比值法，均值比值法，其中，均值比值法和最为广泛使用的对数比值法对原始的比值法有了很大改进。E.E.Kuruoglu和J.Zerubia在2004年提出的对数比值法将SAR图像系统里的乘性斑点噪声转换为加性噪声，并将图像像素值的数量级转换到对数尺度上，便于后续处理里的噪声抑制，是最为常用的差异图生成方法；2007年J.Inglad和G.Mercier提出的均值比值法是将图像每个像素点取邻域平均值，有利于去除图像噪声，提高分类精度。然而，均值比值图在利用均值抑制噪声的同时模糊了图像的边缘，减少了图像细节信息，降低了变化检测正确率；对数比值图在抑制背景区域的同时也减小了变化区域的像素值，降低了变化检测的精度。

图像融合在自动目标识别、遥感、机器人视觉、智能系统、医学图像处理以及军事领域有着广泛的应用。它具有可以充分利用各源图像互补信息；更适合人的视觉感受；适合进一步分析等优点而广受关注。

小波变换是近年来最为著名的变换，因其具有多尺度和多分辨率的特点而被广泛应用于图像，语音等信号处理领域。2012年，公茂果等将小波融合方法用在SAR图像变化检测生成差异图信息图上并取得很好的效果。由于二维小波变换方向选择性差，在图像表达中不可避免地在图像边缘细节纹理处引入一定程度的模糊，在变化检测领域，也有一定的局限性。

Contourlet变换作为一种新型的多尺度，多分辨率分析工具，不但具备小波的多分辨率，局部性和临界采样等特性，还具备小波欠缺的多方向性和各向异性，所以Contourlet变换比小波变换更能稀疏地表达二维图像。遥感图像具有信息量大，图像细节丰富，Contourlet变换比小波变换有更丰富的系数来表达图像并且能用较少的系数表示图像的细节，所以本发明将Contourlet变换应用于SAR图像差异图的生成过程中。

在常用的SAR图像变化检测差异图生成方法中，对数比值法对变化区域不敏感，漏检率较高，均值比值法易将未变区域归于变化区域，误检率较高，都直接导致检测总错误数较大。而差异图像的生成至关重要，如果差异图中本来含有较多的噪声成分，对后续的分析正确率有直接影响，进而决定整个变化检测系统的正确率。

发明内容

本发明的目的在于：为了降低SAR图像变化检测的错误率，克服现有方法较高的检测错误率，提供一种将图像有效信息最大化，减少错误率，便于后期差异图分析，改善最终分析处理精度进而提高检测精度的基于轮廓波(Contourlet)变换的SAR图像变化检测差异图生成方法。

本发明的技术方案是，其具体实现步骤如下：

(1)输入两幅经过滤波去噪，辐射校正与几何配准的预处理的不同时间相同地域的SAR图像X_A、X_B。

(2)对该两幅SAR图像构造对数比值图其中X_l为生成的对数比值差异图，其中X_l(i，j)，X_A(i，j)，X_B(i，j)分别为对数比值图X_l，图像X_A，图像X_B在图像中坐标为(i，j)的像素点的像素值。

(3)对该两幅图像构造均值比值图其中μ_A(i，j)，μ_B(i，j)分别为图像X_A，X_B中以坐标为(i，j)的像素点为中心邻域像素平均均值，X_m(i，j)为均值比值图X_m图像中坐标为(i，j)的像素点的像素值。

(4)分别对对数比值图像X_l、均值比值图像X_m进行Contourlet变换，得到对数比值图的Contourlet系数Y_l和均值比值图的Contourlet系数Y_m。

(5)对得到的对数比值图的Contourlet系数Y_l和均值比值图的Contourlet系数Y_m分为低频部分和高频部分进行融合。

5.1获取融合后图像低频系数，按取平均规则融合得到低频系数，即：其中和分别为Y_m和Y_l的低频部分，由于低频系数代表图像的轮廓信息，包含图像的变化部分，本发明对低频系数利用平均规则旨在增强融合图像低频部分，即增强变化部分的表现。

5.2获取融合后图像高频系数，按取邻域能量最小原则规则融合得到高频系数，即： $Y_{N,t}^f\left\{k\right\}(i,j)=\left\{\begin{array}{cc}{Y}_{N,t}^l\left\{k\right\}(i,j),& D_{N,t}^l\left\{k\right\}(i,j)\≤D_{N,t}^m\left\{k\right\}(i,j)\\ Y_{N,t}^m\left\{k\right\}(i,j),& D_{N,t}^l\left\{k\right\}(i,j)\≥D_{N,t}^m\left\{k\right\}(i,j)\end{array}\right.,$ 其中代表坐标位于(i，j)像素的contourlet系数M×N邻域内的能量，和分别为Y_l和Y_m的高频部分，代表Contourlet分解中第N层分解的第t个方向的系数，其中k＝2，3，4，5，6，N＝1，2，3，4，5，t＝1，2，...，2^N，因为高频系数代表图像的细节信息，包括图像中的边缘和线条，本发明选取源图像中邻域能量较小的系数作为融合后图像的系数可以抑制背景区域(未变化区域)的信息，这样变化区域和背景区域就可以呈现出较大的差别，便于后续的分类处理。

(6)由Contourlet生成融合后图像的contourlet系数Y_f，Y_f由低频系数Y_f{1}和高频系数构成。

(7)将融合后得到的Contourlet系数Y_f进行逆Contourlet变换，也称为图像重构，得到目标图像X_F，即SAR图像X_A、X_B的差异信息图，该差异信息图在抑制噪声的同时极大地增强了变化的有效信息，可以用于SAR图像变化检测的后续分析处理。

在常见的SAR图像变化检测体系中，可以分为三个阶段：预处理阶段，包括图像配准和辐射校正；比较阶段，即通过比较方法生成包含变化信息的差异图像；分析阶段，指对比较产生的差异图像进行分析得到最终的二值结果图像。其中，对经过配准和辐射校正的SAR图像进行比较生成差异图是至关重要且必不可少的一步，差异图的质量直接决定后续分析处理的精确度，进而影响到整个SAR图像变化检测系统的性能。近年来，由于SAR图像变化检测在环境监测、土地利用/覆盖、农业调查、城市变化分析、军事侦察和打击效果评估等领域的应用，在实际领域的运用过程中对其精确度的要求日益提高，实际运用的需要促进了SAR图像变化检测的理论发展和研究，本发明正是在这样的背景下进行的创新和研究。

在现有的SAR图像变化检测方法里，像素级分析最为普遍，然而像素级分析不能充分体现图像的相关信息及表达细节。Contourlet变换作为一种新型的多尺度，多分辨率分析工具，具备多分辨率，局部性临界采样特性，多方向性和各向异性，可以稀疏地表达二维图像。遥感图像具有信息量大，图像细节丰富的特点，这样Contourlet变换具有丰富的系数表达图像并且能用较少的系数表示图像的细节。结合SAR图像变化检测的目的，即扩大变化区域与未变化区域的差别，进而将变化区域从未变化区域中分割出来，本发明将Contourlet变换引入并应用于SAR图像差异图的生成过程中，本发明运用基于Contourlet变换的融合可以利用到多幅源图像的信息，并且在适合的融合规则指引下，能在最大程度上从源图像中获取有效变化信息。

本发明运用图像融合技术(采用Contourlet变换)将多种差异图融合，可以充分利用各源图像的优点，将图像的有效信息最大化，对最终的分析处理精度有所改善，提高检测精度。

本发明首先对两幅不同时间相同地域的SAR图像构造对数比值图像和均值比值图像，然后用Contourlet变换处理使其生成对应的Contourlet系数，接着根据设计的融合入规则对两幅图像的系数分高频和低频进行计算，最后对融合后的Contourlet系数进行逆Contourlet变换产生差异图。

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

1、本发明在利用现有的差异图像生成方法中最为常见的对数比值法，均值比值法，克服均值比值图的误检率较高，对数比值图漏检率较高的缺陷，本发明经过大量的分析和实验，对高频和低频部分分别设计融合规则，通过在高频部分利用对数比值图充分抑制图像的背景信息，在低频部分利用均值比值图增强图像的变化区域，较好的综合它们各自的优点，即均值比值图漏检率较低，对数比值图误检率较低，综合可以得到更能代表变化区域的差异图像。

2、本发明与常见的小波融合方法相比，克服了离散二维小波变换方向选择性差，在图像表达中不可避免地在图像边缘细节纹理处引入一定程度的模糊的缺陷，更好的保留了图像的有效信息，可以得到包含更多有用信息的差异信息图。

3、本发明与其它经典的SAR图像变化检测差异图生成方法相比，获得了较低的检测错误率，抑制图像噪声和提高分析精度得到了改善。

附图说明

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

图2是本发明涉及的Bern地区两幅SAR影像图、标准图和变化参考图，图2(a)表示1999年4月Bern地区的地貌信息；图2(b)表示1999年5月Bern地区的地貌信息；图2(c)表示变化检测的标准结果图；

图3是本发明方法和三种对比方法对Bern地区SAR影像的实验结果图，其中，图3(a)代表对数比值图像，图3(b)代表均值比值图像，图3(c)代表小波融合图像，图3(d)代表本发明方法产生的差异图像；

图4是本发明的ROC曲线与其他方法的比较，ROC曲线代表在阈值遍历的情况下对应的分析正确率，曲线下的面积越大，表示差异图质量越高，适应于后续的分析；

图5是本发明涉及的黄河入海口地区两幅SAR影像图、标准图和变化参考图，图5(a)表示1999年4月黄河入海口地区的地貌信息，图5(b)表示1999年5月黄河入海口地区的地貌信息，图5(c)表示变化检测的标准结果图；

图6是本发明方法和三种对比方法对黄河入海口地区SAR影像的实验结果图，其中，图6(a)代表对数比值图像，图6(b)代表均值比值图像，图6(c)代表小波融合图像，图6(d)代表本发明方法产生的差异图像；

图7是本发明的ROC曲线与其他方法的比较，ROC曲线代表在阈值遍历的情况下对应的分析正确率，曲线下的面积越大，表示差异图质量越高，适应于后续的分析。

具体实施方式

下面结合附图对本发明详细说明

本发明是一种基于Contourlet变换的SAR图像变化检测差异图生成方法，用于生成SAR图像变化检测的差异信息图。本发明现阶段可以运行在32位XP系统(及以上)计算机，matlab7.0及以上平台，Visual C++6.0平台。

实施例1

本发明是一种基于Contourlet变换的SAR图像变化检测差异图生成方法，参见图1，首先对两幅经过预处理的不同时间相同地域的SAR图像构造对数比值图像和均值比值图像，然后用Contourlet变换处理使其生成对应的Contourlet系数，接着根据设计的融合规则对两幅图像的系数分高频和低频进行计算，最后对融合后的Contourlet系数进行逆Contourlet变换产生差异图，基于融合的SAR图像变化检测差异图生成具体实现步骤包括有：

步骤1输入两幅不同时间相同地域的经过滤波去噪，辐射校正与几何配准的预处理SAR图像X_A，X_B。

步骤2对该两幅图像X_A，X_B构造对数比值图其中X_l为生成的对数比值差异图，其中X_l(i，j)，X_A(i，j)，X_B(i，j)分别为对数比值图X_l，图像X_A，图像X_B在图像中坐标为(i，j)的像素点的像素值。

步骤3对该两幅图像X_A、X_B，构造均值比值图其中μ_A(i，j)，μ_B(i，j)分别为图像X_A，X_B中以坐标为(i，j)的像素点为中心邻域像素平均均值，X_m(i，j)为均值比值图X_m图像中坐标为(i，j)的像素点的像素值。

步骤4分别对对数比值图像X_l、均值比值图像X_m进行Contourlet变换，得到对数比值图的Contourlet系数Y_l和均值比值图的Contourlet系数Y_m。

步骤5对得到的对数比值图的Contourlet系数Y_l和均值比值图的Contourlet系数Y_m分为低频部分和高频部分进行融合。

5.1获取融合后图像低频系数，按取平均规则融合得到低频系数，即：其中和分别为Y_m和Y_l的低频部分。

5.2 获取融合后图像高频系数，按取邻域能量最小原则规则融合得到高频系数，即： $Y_{N,t}^f\left\{k\right\}(i,j)=\left\{\begin{array}{cc}{Y}_{N,t}^l\left\{k\right\}(i,j),& D_{N,t}^l\left\{k\right\}(i,j)\≤D_{N,t}^m\left\{k\right\}(i,j)\\ Y_{N,t}^m\left\{k\right\}(i,j),& D_{N,t}^l\left\{k\right\}(i,j)\≥D_{N,t}^m\left\{k\right\}(i,j)\end{array}\right.,$ 其中代表坐标位于(i，j)像素的contourlet系数M×N邻域内的能量，和分别为Y_m和Y_l的第N层第t个方向的的高频部分，k＝N+1。

步骤6由Contourlet生成融合后图像的contourlet系数Y_f，Y_f由低频系数Y_f{1}和高频系数构成。

步骤7将融合后得到的Contourlet系数Y_f进行逆Contourlet变换，也称为图像重构，得到目标图像X_F，即SAR图像X_A、X_B的差异信息图。

本发明利用现有的差异图生成图像方法中最为常见的对数比值法，均值比值法，并克服均值比值图的误检率较高，对数比值图漏检率较高的缺陷，较好的综合利用它们各自的优点，即均值比值图漏检率较低，对数比值图误检率较低，融合后得到更能代表变化区域的差异图像。

实施例2

基于Contourlet变换的SAR图像变化检测差异图生成方法同实施例1，参照图1，实现本发明Contourlet变换用于SAR图像变化检测差异图，首先输入经过预处理的两幅不同时间相同地域的SAR图像，对该两幅图像分别构造对数比值图像和均值比值图像，然后用Contourlet变换处理使其生成对应的Contourlet系数，接着根据设计的融合规则对两幅图像的Contourlet系数分高频和低频进行计算，最后对融合后的Contourlet系数进行逆Contourlet变换产生差异图。下边通过本例对该发明的实现过程进行详细说明：

步骤1在相同地域的不同时间获取两幅合成孔径雷达(SAR)图像，将这两幅SAR图像进行滤波去噪，辐射校正与几何配准的预处理，得到处理后的两幅图像X_A，X_B，其中预处理后得到的图像X_A，预处理后得到的图像X_B。

通过预处理可以消除图像的几何误差，以达到对同一区域不同图像的地理坐标的匹配，消除传感器自身引起的噪声和大气辐射引起的辐射噪声。

步骤2用预处理后的两幅图像X_A，X_B构造对数比值图X_l，构造方式为：将图像X_A中位于i行j列的像素点(i，j)的灰度值X_A(i，j)和相对应的图像X_B中位于i行j列的像素点的灰度值X_B(i，j)，通过对数比值运算得到对数比值差异图X_l中位于i行j列的像素点(i，j)的灰度值X_l(i，j)，若X_l(i，j)为0，则表示图像X_A中该像素点没有随时间的变化而发生变化，否则，认为图像X_A中的该像素点发生了变化；对图像X_A和图像X_B中每个位于i行j列的像素点的灰度值从左到右，从上到下都进行差值运算，得到对数比值差异图X_l。

步骤3对该两幅图像X_A、X_B，构造均值比值图，构造方式为：将图像X_A中位于i行j列的像素点(i，j)的灰度值X_A(i，j)和相对应的图像X_B中位于i行j列的像素点的灰度值X_B(i，j)，通过均值比值运算其中μ_A(i，j)，μ_B(i，j)分别为图像A，B的邻域像素平均值，得到均值比值法差异图X_m中位于i行j列的像素点(i，j)的灰度值X_m(i，j)，如果X_m(i，j)接近于0，表示图像X_A中位于i行j列的像素点没有随时间的变化而发生变化，反之，若X_m(i，j)远大于0，则认为该像素点发生了变化；对图像X_A和图像X_B中每个位于i行j列的像素点的灰度值从左到右，从上到下都进行此运算，构造比值法差异图X_m。

步骤4对对数比值图像X_l、均值比值图像X_m进行N层离散Contourlet变换，得到每幅差异图在N个分解层上高频带和低频带上的多方向多尺度分解的Contourlet系数，N＝5，本发明经过大量实验，选定N＝5时效果较佳。

4.1在第一层，即N＝1，将对数比值图X_l分解成一个低频子带图像和一个带通子带图像将带通子带图像进行2^N级方向分解，在此为2个方向分解，得到2个高频带Contourlet系数其中t＝1，2；将均值比值图X_m分解成一个低频子带图像和一个带通子带图像将带通子带图像进行2^N级方向分解，在此为2个方向分解，得到2个高频带Contourlet系数其中t＝1，2。

4.2 在第二层，即N＝2，将对数比值图X_l经步骤4.1生成的低频子带图像分解成一个低频子带图像和一个带通子带图像得到4个高频带Contourlet系数其中t＝1，2，3，4；将对数比值图X_m经步骤4.1生成的低频子带图像分解成一个低频子带图像和一个带通子带图像得到4个高频带Contourlet系数其中t＝1，2，3，4。

4.3在第三层，即N＝3，将对数比值图X_l经步骤4.2生成的低频子带图像分解成一个低频子带图像和一个带通子带图像得到8个高频带Contourlet系数其中t＝1，2，…，8；将对数比值图X_m经步骤4.2生成的低频子带图像分解成一个低频子带图像和一个带通子带图像得到8个高频带Contourlet系数其中t＝1，2，…，8。

4.4在第四层，即N＝4，将对数比值图X_l经步骤4.3生成的低频子带图像分解成一个低频子带图像和一个带通子带图像得到16个高频带Contourlet系数其中t＝1，2，…，16；将对数比值图X_m经步骤4.3生成的低频子带图像分解成一个低频子带图像和一个带通子带图像得到16个高频带Contourlet系数其中t＝1，2，…，16。

4.5在第五层，即N＝5，将对数比值图X_l经步骤4.4生成的低频子带图像分解成一个低频子带图像和一个带通子带图像得到一个低频带Contourlet系数和32个高频带Contourlet系数其中t＝1，2，…，32；将对数比值图X_m经步骤4.4生成的低频子带图像分解成一个低频子带图像和一个带通子带图像得到一个低频带Contourlet系数和16个高频带Contourlet系数其中t＝1，2，…，32。本发明通过5层Contourlet分解为进一步的系数处理作好了数据准备。

步骤5将均值比值图和对数比值图的Contourlet系数分为高频和低频处理，在Contourlet变换中，低频系数代表图像轮廓信息，高频代表图像细节信息。

5.1低频系数处理如下：

$Y_f\left\{1\right\}=\left\{Y_5^m\right\{1\}+Y_5^l\{1\left\}\right\}/2---\left(1\right)$

其中，m为均值比值图的子图像，l为对数比值图的子图像，f为融合后图像的子图像，(i，j)表示图像中像素点的第i行j列。

由于Contourlet变换分解后的低频信息反映的是图像的概貌信息，能表达变化部分的信息，为使融合图像中保留的更多的两幅图像中的有用变化信息，故对变换后两幅图的低频信息取加权平均分量作为融合图像的低频分量。

5.2高频系数处理如下：

$Y_{N,t}^f\left\{k\right\}(i,j)=\left\{\begin{array}{cc}{Y}_{N,t}^l\left\{k\right\}(i,j),& D_{N,t}^l\left\{k\right\}(i,j)\≤D_{N,t}^m\left\{k\right\}(i,j)\\ Y_{N,t}^m\left\{k\right\}(i,j),& D_{N,t}^l\left\{k\right\}(i,j)\≥D_{N,t}^m\left\{k\right\}(i,j)\end{array}\right.---\left(2\right)$

其中，代表对数比值图中以(i，j)为中心第N层中第t个方向的Contourlet系数的邻域内能量，代表均值比值图中以(i，j)为中心内第N层中第t个方向的Contourlet系数的邻域能量，其中k＝2，3，4，5，6，将N取遍1到5，将t取遍每层的1到2^N个方向，即可得到所有需要的Contourlet系数。

由于高频信息代表图像的细节信息，比如具有跳跃性的边缘和线条特征，为了最大程度的抑制背景信息的，本发明选择局部区域能量最小的Contourlet系数作为融合后图像的Contourlet系数。

本发明在差异图像生成方法中利用了现有的最常见的对数比值法，均值比值法，并克服均值比值图的误检率较高，对数比值图漏检率较高的缺陷，本发明经过大量的分析和实验，对高频和低频部分分别设计融合规则，通过在高频部分利用对数比值图充分抑制图像的背景信息，在低频部分利用均值比值图增强图像的变化区域，较好的综合它们各自的优点，即均值比值图漏检率较低，对数比值图误检率较低，取长补短，在频域里优化了数据。

步骤6由Contourlet生成融合后图像的Contourlet系数Y_f，Y_f由低频系数Y_f{1}和高频系数构成，其中k＝2，3，4，5，6，N＝1，2，3，4，5，t＝1，2，…，2^N，让k，N，t遍历这些取值，即可得到整体Contourlet系数Y_f。

步骤7将融合后得到的Contourlet系数Y_f进行逆Contourlet变换，也称为图像重构，得到目标图像X_F，即SAR图像X_A、X_B的差异信息图。本发明经融合后得到的差异图像更能代表变化区域的变化信息。用作SAR图像变化检测下一步的继续分析处理。

本发明运用Contourlet变换处理差异图以生成新的差异图，便于后期差异图分析，提高SAR图像变化检测中的检测精度，减少错误率。

实施例3

基于Contourlet变换的SAR图像变化检测差异图生成方法同实施例1-2，参照图1，本实例中采用本发明对Bern地域的不同时间获取两幅合成孔径雷达(SAR)图像及参考图进行仿真，实现步骤如下：

步骤1在Bern地域的不同时间获取两幅合成孔径雷达(SAR)图像，将这两幅SAR图像进行滤波去噪，辐射校正与几何配准的预处理，处理后的两幅图像X_A，X_B，其中预处理后得到的图像X_A如图2(a)所示，预处理后得到的图像X_B如图2(b)所示。

步骤2利用预处理后的如图2(a)所示的图像X_A和图2(b)所示的图像X_B，用对数比值运算构造如图3(a)所示的对数比值图X_l，详细实施步骤如实施例2中的步骤2。

步骤3利用预处理后的如图2(a)所示的图像X_A和图2(b)所示的图像X_B，用均值比值运算构造如图3(b)所示的均值比值法图X_m，详细实施步骤如实施例2中的步骤3。

步骤4对如图3(a)所示的对数比值图X_l和如图3(b)所示的均值比值图X_m分别进行N层Contourlet分解，得到每幅差异图在N个分解层上多个方向高频带和低频带上的Contourlet系数，N＝5，详细实施步骤如实施例2中的步骤4。

步骤5对各分解层各个方向的高频带和低频带的Contourlet系数用不同的融合规则进行融合处理，即对低频带Contourlet系数采用取平均值的方法进行融合，得到低频带融合Contourlet系数，对高频带Contourlet系数采用选能量最小的方法进行融合处理，得到高频带融合Contourlet系数，详细实施步骤如实施例2中的步骤5。

步骤6由Contourlet生成融合后图像的Contourlet系数Y_f，详细实施步骤如实施例1中的步骤6。

步骤7将融合后得到的Contourlet系数Y_f进行逆Contourlet变换，所得到的目标图像即如图3(d)所示的融合后差异图X_F，详细实施步骤如实施例2中的步骤7。

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

1、仿真参数

对于具有参考图的实验仿真图组，可进行定量的变化检测结果分析，主要评价指标有：

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

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

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

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

为了验证基于Contourlet融合的SAR图像变化检测差异图生成方法的优越性，将本发明方法与在差异图生成阶段比较常见的对数比值法和均值比值法算法及小波融合方法性能做出比较。对带参考图的真实SAR图像数据Bern地区进行了实验。这里把对数比值法简称为LR，均值比值法简称为MR，小波融合法简称为WT，Contourlet变换融合法简称为CT.

将本发明和现有技术中的三种方法分别应用在真实Bern SAR图像数据上，进行差异图生成。实验相关图像如图2所示。其中图2(a)表示1999年4月Bern地区的地貌信息，即为第一时间获取图像X_A，图2(b)表示1999年5月Bern地区的地貌信息，即为第二时间获取图像X_B，图2(c)表示变化检测的参考标准结果图。

各种差异图在真实Bern区域SAR图像数据的实验结果图如图3所示。

在图3中图3(a)代表对数(LR)比值图像，3(b)代表均值比值(MR)图像，3(c)代表小波融合(WT)图像，3(d)代表本发明(CT)方法产生的差异图像，CT为本发明所用Contourlet变化的缩写。

在图4中，ROC曲线代表在阈值遍历的情况下对应的分析正确率，曲线下的面积越大，表示差异图质量越高，适应于后续的分析。

各种方法的差异图用FLICM聚类方法作分析性能指标如下表所示：

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

从表1中可以看出，本发明在变化检测总错误数上，相比其他对比算法是最少的，而且在误检数与漏检数上相对都比较均衡，数量上也处于中间水平，说明了本发明方法的优越性。从附图3中可以看到，MR方法对图像模糊较多，导致较高的误检率，从表1中的误检数也可以给出同样的结论，LR方法对变化区域有一定的缩小，这两种方法的总错误数是四种方法中最高的。而WT方法对以上两种缺陷有一定的综合改善，经由融合后的图像可以结合源图像的优点，综合出较好的结果，而CT方法是总错误率最小的。总体来看，通过表1性能结果和图3的视觉效果及图4ROC曲线形状均说明了本发明方法生成的差异图性能良好，质量高的差异图像便于后续分析处理，最终提高了SAR图像变化检测的精确度，使得总错误数更少，误检数与漏检数上相对都比较均衡。

实施例4

基于Contourlet变换的SAR图像变化检测差异图生成方法同实施例1-3，利用本发明对黄河入海口地域的不同时间获取两幅合成孔径雷达(SAR)图像仿真进行差异图生成。

步骤1在黄河入海口地域的不同时间获取两幅合成孔径雷达(SAR)图像，将这两幅SAR图像进行滤波去噪，辐射校正与几何配准的预处理，得到处理后的两幅图像X_A，X_B，其中预处理后得到的图像X_A如图5(a)所示，预处理后得到的图像X_B如图5(b)所示。

步骤2利用预处理后的如图5(a)所示的图像X_A和图5(b)所示的图像X_B，用对数比值运算构造如图6(a)所示的对数比值图X_l。

步骤3利用预处理后的如图5(a)所示的图像X_A和图5(b)所示的图像X_B，用均值比值运算构造如图6(b)所示的均值比值法图X_m。

步骤4对如图6(a)所示的对数比值图X_l和如图6(b)所示的均值比值图X_m分别进行N层Contourlet分解，得到每幅差异图在N个分解层上多个方向高频带和低频带上的Contourlet系数，N＝5。

步骤5对各分解层个方向的高频带和低频带的Contourlet系数用不同的融合规则进行融合处理，即对低频带Contourlet系数采用取平均值的方法进行融合，得到低频带融合Contourlet系数，对高频带Contourlet系数采用选能量最小的方法进行融合处理，得到高频带融合Contourlet系数。

步骤6由Contourlet生成融合后图像的Contourlet系数Y_f。

步骤7将融合后得到的Contourlet系数Y_f进行逆Contourlet变换，所得到的目标图像即如图6(d)所示的融合后差异图X_F。

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

1、仿真参数同实施例3

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

将本发明和现有三种差异图生成方法分别应用于真实黄河入海口区域SAR图像数据上，进行差异图生成。实验相关图像如参考图5所示。

图5(a)表示2008年6月入海口区域地貌信息，即为第一时间获取图像X_A，图5(b)表示2009年6月入海口区域地貌信息，即为第二时间获取图像X_B，图5(b)中可见矩形区域为新生成的农田，图5(c)表示变化检测的参考标准结果图。

各种差异图生成方法在真实黄河入海口区域SAR图像数据的实验结果图如图6所示，其中，图6(a)代表对数(LR)比值图像，6(b)代表均值比值(MR)图像，6(c)代表小波融合(WT)图像，6(d)代表本发明(CT)方法产生的差异图像。

图6(a)-6(d)均反映出了变化，从视觉上分辨可见，图6(a)在抑制背景未变化区域的同时也抑制了变化区域的表现，图6(b)中变化区域的边缘被模糊严重，对未变区域的噪声抑制也不明显，图6(c)为小波融合后的到的差异图，图6(d)为本发明Contourlet融合后的到的差异图，由于融合策略可以综合源图像的有用信息，图6(c)和图6(d)效果明显优于图6(a)和图6(b)。因为本发明所采用的Contourlet变换比小波变换有更多的方向选择，从图6(c)和图6(d)比较可见，Contourlet融合后的图像的边缘及细节保持能力强于小波融合后的图像。

对于低分辨率的SAR图像差异图，人眼感知是有局限的，要得到精确的数据必须借助于ROC曲线判定和常用评价指标进行比较分析。

在图7中，ROC曲线代表在阈值遍历的情况下对应的分析正确率，曲线下的面积越大，表示差异图质量越高，适应于后续的分析。

各种方法生成的差异图用FLICM聚类方法作分析性能指标如下表所示：

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

从表2中可以看出，本发明在变化检测总错误数上，相比其他对比算法是最少的，而且在误检数与漏检数上相对都比较均衡，数量上也处于中间水平，说明了本发明方法的优越性。

从附图6中可以看到，MR方法产生的图像比较模糊，对边缘信息保留不够明显，导致误检率较高，从表2中的漏检数也可以给出同样的结论，LR方法对变化区域的保留不够突出，从后续分析结果来看，导致很多变化区域不能被判断为已经变化，使得漏检像素较多，这两种方法的总错误数是四种方法中最高的。而小波融合方法对以上两种缺陷略有改善，本发明的CT方法是总错误率最小的。总体来看，通过表2性能结果和参考图6的视觉效果及参考图7中ROC曲线形状说明了本发明方法生成的差异图性能良好，便于后续分析处理，提高最终变化检测的精确度，使得总错误数更少，误检数与漏检数上相对都比较均衡。

综上所述本发明方法能够得到更好的SAR图像变化检测结果.本发明与常见的小波融合方法相比，克服了离散二维小波变换方向选择性差，在图像表达中不可避免地在图像边缘细节纹理处引入一定程度的模糊的缺陷，更好的保留了图像的有效信息，可以得到包含更多有用信息的差异信息图。

通过两个实验的分析可以看出，本发明方法的结果与其他方法相比更为精确，能够生成效果更佳的差异图便于后续分析处理，从而得到更佳的SAR图像变化检测的结果。

综上，本发明的基于轮廓波变换的SAR图像变化检测差异图生成方法，其实现过程主要包括：先对两幅经预处理的不同时间相同地域的SAR图像构造对数比值图像和均值比值图像；用Contourlet变换处理使其生成对应的Contourlet系数；依不同融合规则分别对两幅图像系数分高频和低频进行计算；对融合后的Contourlet系数进行逆Contourlet变换生成变化差异图。本发明利用均值比值图像和对数比值图像分别提取其高频和低频的不同特点，基于Contourlet变换的图像融合充分利用了各源图像互补信息，使SAR图像变化检测得到更好结果，获得了较低的检测错误率，抑制图像噪声，提高分析精度。与其它差异图生成方法对比，本发明噪声抑制力强，边缘保持佳，变化信息得以最大程度保留。

Claims (1)

1.一种基于轮廓波变换的SAR图像变化检测差异图生成方法，其特征在于：基于轮廓波变换的SAR图像变化检测差异图生成具体实现步骤包括有：

(1)输入两幅不同时间相同地域的经过滤波去噪，辐射校正与几何配准的预处理SAR图像X_A，X_B；

(2)对该两幅图像X_A，X_B构造对数比值图其中X_l为生成的对数比值差异图，其中X_l(i，j)，X_A(i，j)，X_B(i，j)分别为对数比值图X_l，图像X_A，图像X_B在图像中坐标为(i，j)的像素点的像素值；

(3)对该两幅图像X_A、X_B，构造均值比值图其中μ_A(i，j)，μ_B(i，j)分别为图像X_A，X_B中以坐标为(i，j)的像素点为中心的邻域像素平均均值，X_m(i，j)为均值比值图X_m图像中坐标为(i，j)的像素点的像素值；

(4)分别对对数比值图像X_l、均值比值图像X_m进行Contourlet变换，得到对数比值图的Contourlet系数Y_l和均值比值图的Contourlet系数Y_m；

(5)对得到的对数比值图的Contourlet系数Y_l和均值比值图的Contourlet系数Y_m分为低频部分和高频部分进行融合：

5.1获取融合后图像低频系数，按取平均规则融合得到低频系数，即：其中和分别为Y_m和Y_l的低频部分；

5.2获取融合后图像高频系数，按取邻域能量最小原则融合得到高频系数，即： $Y_{N,t}^f\left\{k\right\}(i,j)=\left\{\begin{array}{cc}{Y}_{N,t}^l\left\{k\right\}(i,j),& D_{N,t}^l\left\{k\right\}(i,j)\≤D_{N,t}^m\left\{k\right\}(i,j)\\ Y_{N,t}^m\left\{k\right\}(i,j),& D_{N,t}^l\left\{k\right\}(i,j)>D_{N,t}^m\left\{k\right\}(i,j)\end{array}\right.,$ 其中代表坐标位于(i，j)像素的contourlet系数在M×N邻域内的能量，和分别为Y_l和Y_m的第N层第t个方向的高频部分，k＝N+1；代表均值比值图中以(i，j)为中心内第N层中第t个方向的Contourlet系数的邻域能量；

(6)由Contourlet生成融合后图像的contourlet系数Y_l，Y_f由低频系数Y_f{1}和高频系数构成；

(7)将融合后得到的Contourlet系数Y_f进行逆Contourlet变换，也称为图像重构，得到目标图像X_F，即SAR图像X_A、X_B的差异信息图。