CN104680536B - 利用改进的非局部均值算法对sar图像变化的检测方法 - Google Patents

Abstract

本发明公开了一种利用改进的非局部均值算法对SAR图像变化的检测方法，主要解决现有技术去噪漏检像素个数多和变化检测精度不高的问题。其实现过程是：1)输入变化前后的SAR图像；2)对于两幅图像中的每个像素点，计算其与搜索区域内所有像素点的基于对数比的傅里叶核加权距离并引入二维高斯函数计算该点与搜寻区域内所有像素点的权重；3)根据权重求该点新的灰度值并替换其原灰度值，得到去噪后两幅图像；4)利用比值算子将去噪后的两幅图像生成差异图；4)利用FLICM算法分割差异图，输出变化检测结果图。本发明克服了去噪过程中的过平滑现象，更好的保留了细节信息，降低了错误率，提高了图像变化检测的精确度，可用于环境监测。

Description

利用改进的非局部均值算法对SAR图像变化的检测方法

技术领域

本发明属于图像处理技术领域，特别涉及对SAR图像变化的检测方法，可用于环境监测、土地利用、森林采伐的监测。

背景技术

合成孔径雷达SAR图像具有高分辨率、全天候、全天时的特点，可方便地获得同一地区不同时间的图像。SAR图像的变化检测是对同一地区在不同时间所得的图像进行比较得到差异图，然后对差异图进行分析得到变化区域和不变区域的过程。SAR图像的变化检测在环境监测、农作物生长状况监测、土地利用、森林采伐监测等方面有着非常广泛的应用。

在最近几十年，关于SAR图像变化检测的研究得到越来越多的关注，许多变化检测方法被提出。在目前常用的方法中，SAR图像的变化检测分为以下三个过程：

1)对得到的图像进行预处理，主要包括辐射校正、几何校正和图像配准；

2)对预处理后的图像进行比较，获得差异图；

3)对差异图进行分析，获得变化检测结果图像。

因为SAR图像的乘性斑点噪声处理起来比较困难，所以它的去噪处理和抑制斑点噪声显得尤为重要。由此可知，第一步和第二步对SAR图像的变化检测至关重要。对于相关的处理，目前比较优秀的算法有Yaoguo Zheng等提出的使用结合差异图和k-均值聚类的SAR图像变化检测的算法CDI、公茂果提出的使用小波变换融合基于均值比和对数比的差异图的方法FDWT和Buades等提出的基于非局部均值NLM的图像去噪方法。在这几种方法中，基于小波变换融合的方法虽然降低了图像斑点噪声的影响，却对一些边缘和细小部分的图像信息检测不明显，而且去噪能力也不是足够好。基于结合差异图和k-均值聚类的SAR图像变化检测虽然具有很强的去噪能力，但是在去噪过程中出现过平滑现象，在一定程度上抑制了图像的变化信息，因此造成图像的漏检率过大。而非局部均值去噪的方法仅仅考虑了图像中像素点的灰度属性，而忽略了它的空间结构属性，因此依然会丢失部分细节信息，出现模糊过平滑现象，尤其是当含噪图像块的相似块较少时，去噪性能会大大降低。

发明内容

本发明针对上述已有技术的不足，提出一种利用改进的非局部均值算法对SAR图像变化的检测方法，以减小图像的漏检率，提高图像变化检测的精确度。

为了实现上述目标，本发明的技术方案包括如下步骤：

(1)分别输入变化前的SAR原图像Q₁和变化后的SAR图像Q₂；

(2)对变化前的原图像Q₁和变化后的图像Q₂分别进行去噪处理，得到不含噪的原图像Q₃和不含噪的变化后图像Q₄：

2a)分别选取图像中以第i像素点和第j个像素点为中心，大小为5×5的局部分析窗N_i和N_j，这两个分析窗N_i和N_j内的像素点灰度值向量分别为x(N_i)和x(N_j)；对比和的比值，选取其中较大者进行取对数操作，得到对数比距离d₁(i,j)为：

2b)将对数比距离d₁(i,j)引入到傅里叶核加权距离的求解公式中，得到基于对数比值的傅里叶核的加权距离：

其中，i1是傅里叶核在以像素点i为中心，大小为5×5的局部分析窗N_i内构造的权重，是由傅里叶核构造的以i1为权重的加权距离；

2c)利用基于对数比值的傅里叶核的加权距离，引入二维高斯核函数，计算第i像素点与其搜索区域S_i内第j个像素点的权重w(i,j)：

其中，i2是傅里叶核在以像素点i为中心，大小为7×7的局部分析窗N_i内构造的权重，是由傅里叶核构造的以i2为权重的加权距离，Z(i)的表达式如下所示：

其中，0≤w(i,j)≤1且Z(i)表示归一化常数，S_i是以像素点i为中心，大小为5×5的搜索区域，h是平滑参数，取值为1/3，N_i和N_j分别表示以第i个像素点和第j个像素点为中心的5×5的局部分析窗，x(N_i)和x(N_j)分别是局部分析窗N_i和N_j内像素点的灰度值向量；

2d)对以第i个像素点为中心的搜索区域S_i内所有像素点的灰度值进行加权平均，求得图像中第i个像素点的新的灰度值INL(i)：

INL(i)＝∑w(i,j)X(j),j∈S_i，

其中，X(j)是去噪声前原图像中第j个像素点的灰度值，w(i,j)是与X(j)对应的权值，表示第i个像素点和第j个像素点之间的相似性；

2e)对变化前的SAR原图像Q₁和变化后的SAR图像Q₂分别重复步骤2a)到步骤2d)，直到两个图像中每个像素点都获得新的灰度值；

2f)对变化前的SAR原图像Q₁和变化后的SAR图像Q₂分别用获得的每个像素点的新的灰度值替代该像素点的原灰度值，即可得到不含噪的变化前的SAR原图像Q₃和不含噪的变化后的SAR图像Q₄；

(3)利用比值算子得到不含噪的原图像Q₃和不含噪的变化后图像Q₄的差异图Q₅；

(4)利用FLICM聚类算法对差异图Q₅的灰度值进行类别数为2的聚类分割，其中一类为变化类，另一类为非变化类；

(5)输出由步骤(4)得到的变化类图像，即为变化前SAR的图像Q₁和变化后SAR图像Q₂的变化检测结果图。

本发明与上述处理SAR图像变化检测的三个优秀算法相比有以下优点：

1.本发明中使用的是傅里叶核加权的距离，能够更好地分辨、提取并放大SAR图像的有用特征信息；

2.本发明使用一种基于对数比值距离的方法，可以将乘性斑点噪声转换为加性斑点噪声，对图像的乘性斑点噪声有更好的鲁棒性；

3.本发明两次利用傅里叶核挖掘图像的空间结构信息，然后利用二维高斯函数求权重，使得去噪过程在后续边缘检测中不会偏向任一方向，提高了在检测边缘细小部分的性能，有效地减少了漏检像素点的个数，较大程度地提高了图像变化检测的精确度。

附图说明

图1是本发明的流程图；

图2是本发明实验使用的Bern地区的SAR图像的数据集；

图3是现有的Bern地区的SAR图像变化检测的参考图；

图4是分别用CDI算法、FDWT算法、NLM算法和本发明对图2进行变化检测的结果图；

图5是本发明实验使用的Ottawa地区的SAR图像的数据集；

图6是现有的Ottawa地区的SAR图像变化检测的参考图；

图7是分别用CDI算法、FDWT算法、NLM算法和本发明对图5进行变化检测的结果图。

具体实施方式

参照图1，本发明的实现步骤如下。

步骤1：分别输入变化前的SAR原图像Q₁和变化后的SAR图像Q₂。

针对同一地区不同时间拍摄的两幅SAR图像，将时间较早拍摄的SAR图像定义为变化前的SAR原图像，记为Q₁，将时间较晚拍摄的SAR图像定义为变化后的SAR图像，记为Q₂。

步骤2：对变化前的原图像Q₁和变化后的图像Q₂分别进行去噪处理，得到不含噪的原图像Q₃和不含噪的变化后图像Q₄。

2a)分别选取图像中以第i像素点和第j个像素点为中心，大小为5×5的局部分析窗N_i和N_j，这两个分析窗N_i和N_j内的像素点灰度值向量分别为x(N_i)和x(N_j)；对比和的比值，选取其中较大者进行取对数操作，得到对数比距离d₁(i,j)为：

2b)将对数比距离d₁(i,j)引入到傅里叶核加权距离的求解公式中，得到基于对数比值的傅里叶核的加权距离：

其中，i1是傅里叶核在以像素点i为中心，大小为5×5的局部分析窗N_i内构造的权重，是由傅里叶核构造的以i1为权重的加权距离。

2c)利用基于对数比值的傅里叶核的加权距离，引入二维高斯核函数，计算第i像素点与其搜索区域S_i内第j个像素点的权重w(i,j)：

其中，i2是傅里叶核在以像素点i为中心，大小为7×7的局部分析窗N_i内构造的权重，是由傅里叶核构造的以i2为权重的加权距离，Z(i)的表达式如下所示：

其中，0≤w(i,j)≤1且Z(i)表示归一化常数，S_i是以像素点i为中心，大小为5×5的搜索区域，h是平滑参数，取值为1/3，N_i和N_j分别表示以第i个像素点和第j个像素点为中心的5×5的局部分析窗，x(N_i)和x(N_j)分别是局部分析窗N_i和N_j内像素点的灰度值向量；

2d)对以第i个像素点为中心的搜索区域S_i内所有像素点的灰度值进行加权平均，求得图像中第i个像素点的新的灰度值INL(i)：

INL(i)＝∑w(i,j)X(j),j∈S_i，

其中，X(j)是去噪声前原图像中第j个像素点的灰度值，w(i,j)是与X(j)对应的权值，表示第i个像素点和第j个像素点之间的相似性；

2e)对变化前的SAR原图像Q₁和变化后的SAR图像Q₂分别重复步骤2a)到步骤2d)，直到两个图像中每个像素点都获得新的灰度值；

2f)对变化前的SAR原图像Q₁和变化后的SAR图像Q₂分别用获得的每个像素点的新的灰度值替代该像素点的原灰度值，即可得到不含噪的变化前的SAR原图像Q₃和不含噪的变化后的SAR图像Q₄。

步骤3：利用比值算子构造差异图Q₅。

用不含噪的原图像Q₃和不含噪的变化后图像Q₄构造比值算子，并选取和中较小者为1的被减对象，将取得的结果作为构造的差异图Q₅：

步骤4：利用FLICM聚类算法对差异图Q₅的灰度值进行类别数为2的聚类分割。

4a)设置模糊化参数m＝2，聚类中心个数r＝2，循环计数器b＝0，并将步骤3获得的差异图的灰度矩阵转化为一个列向量X＝[x_i]，x_i为第i个像素点的灰度值，i＝1,2,…N，N为总像素点数，将第i个像素点属于第k类的隶属度记为u_ki，u_ki∈[0,1]且k表示聚类的类标，随机初始化一个隶属度矩阵U^b＝[u_ki]；

4b)按照如下公式计算第k类的聚类中心c_k：

其中，是模糊化参数为m的第i个像素点属于第k类的隶属度，k＝1,2…,r，r取值为2；

4c)按照如下公式计算隶属度矩阵U^(b+1)：

U^(b+1)＝{u_ki}，

其中， ||x_i-c_k||²表示目标点i和聚类中心c_k之间的欧几里德距离，d_ij表示i像素点和j像素点的欧几里德距离；

4d)设定聚类停止条件为ε＝0.001，若max{U^(b)-U^(b+1)}<ε，则停止迭代，否则，令b＝b+1，跳转至步骤4b)继续迭代，直到满足聚类的停止条件；

4e)获得聚类后的图像，其中一类为变化类，另一类为非变化类。

步骤5：输出由步骤(4)得到的变化类图像，即为变化前的SAR原图像Q₁和变化后的SAR图像Q₂的变化检测结果图。

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

1.实验条件

实验环境为：windows XP，SPI，CPU Pentium(R)4，基本频率2.4GHZ，软件平台为MatlabR2010a。

实验的第一组图像是瑞士Bern地区的SAR图像数据集。如图2所示，其中，图2(a)是Bern地区变化前的原始图像，图2(b)是Bern地区变化后的SAR图像。两幅图分别拍摄于1999年4月和1999年5月，反映了Bern郊区附近的水灾情况，两幅图像的尺寸均为301×301。Bern地区的SAR图像变化检测的标准结果如图3所示。

实验的第二组是加拿大Ottawa地区的SAR图像数据集。如图5所示，5(a)是变化前的原始图像，5(b)是变化后的SAR图像。两幅图分别拍摄于1997年5月和1997年8月，两图之间发生的变化部分反映是加拿大Ottawa地区夏季雨季来临，洪水淹没渥太华部分陆地区域的情况，两幅大小均为290×350。Ottawa地区的SAR图像变化检测的标准结果如图6所示。

2.实验内容和实验结果

实验一：用本发明和上述现有CDI算法、FDWT算法、NML算法分别对图2的变化进行检测，检测结果如图4所示，其中，4(a)为CDI算法对图2变化的检测结果图，4(b)是FDWT算法对图2变化的检测结果图，4(c)是NML算法对图2变化的检测结果图，而4(d)是本发明对图2变化的检测结果图。

实验二：用本发明和现有CDI算法、FDWT算法、NML算法分别对图5的变化进行检测，检测结果如图7所示，其中，7(a)为CDI算法对图5变化的检测结果图，7(b)是FDWT算法对图5变化的检测结果图，7(c)是NML算法对图2变化的检测结果图，7(d)是本发明对图5变化的检测结果图。

从图4(a)和图7(a)可以看出，CDI算法对图像处理出现了过平滑现象，在抑制斑点噪声的同时造成了变化区域信息的流失。

从图4(b)和图7(b)可以看出，FDWT算法在图像的边缘存在斑点噪声，丢失了一些边缘细节信息。

从图4(c)和图7(c)可以看出，NML算法在背景区域存在一些噪声点，而且没有很好的考虑到边缘信息，所以对乘性斑点噪声不具有很好的鲁棒性。

从图4(d)和图7(d)可以看出，本发明将斑点抑制到了很小的程度，并且在边缘细节部分的变化信息也很完整，可以把图像的变化区域和非变化区域区分的更加明显。

本发明和现有CDI算法、FDWT算法、NML算法对图2和图5的变化检测数据列如表1。

表1 检测结果数据表

表1中展示了四种评价指标：漏检个数FN、误检个数FP、总错误数OE和Kappa系数。其中，漏检个数FN是参考图中发生变化但是实验结果图中未发生变化的像素个数，误检个数FP是参考图中未发生变化但是实验结果图中发生变化的像素个数，总错误数OE是漏检个数和误检个数的总和，Kappa系数是一种衡量图像分割精度的指标，其结果越接近于1，表示所得的结果精确度越高。

由表1可得，本发明较现有CDI算法、FDWT算法、NLM算法，可以获得最低的总错误数和最好的卡帕系数，且本发明检测的结果图最接近参考图。

Claims (2)

1.一种利用改进的非局部均值算法对SAR图像变化的检测方法，包括如下步骤：

(1)分别输入变化前的SAR原图像Q₁和变化后的SAR图像Q₂；

(2)对变化前的原图像Q₁和变化后的图像Q₂分别进行去噪处理，得到不含噪的原图像Q₃和不含噪的变化后图像Q₄：

2a)分别选取图像中以第i像素点和第j个像素点为中心，大小为5×5的局部分析窗N_i和N_j，这两个分析窗N_i和N_j内的像素点灰度值向量分别为x(N_i)和x(N_j)；对比和的比值，选取其中较大者进行取对数操作，得到对数比距离d₁(i,j)为：

$d_1(i,j)=log\left(max\right(\frac{x\left(N_i\right)}{x\left(N_j\right)},\frac{x\left(N_j\right)}{x\left(N_i\right)}\left)\right);$

2b)将对数比距离d₁(i,j)引入到傅里叶核加权距离的求解公式中，得到基于对数比值的傅里叶核的加权距离：

$dis\left(x\right(N_i),x(N_j\left)\right)=\left|\right|\log \left(\max \right(\frac{x\left(N_i\right)}{x\left(N_j\right)},\frac{x\left(N_j\right)}{x\left(N_i\right)}\left)\right)|{|}_{i1}^2$

其中，i1是傅里叶核在以像素点i为中心，大小为5×5的局部分析窗N_i内构造的权重，是由傅里叶核构造的以i1为权重的加权距离；

2c)利用基于对数比值的傅里叶核的加权距离，引入二维高斯核函数，计算第i像素点与其搜索区域S_i内第j个像素点的权重w(i,j)：

$w(i,j)=\frac{1}{Z\left(i\right)}\frac{\exp (-\frac{\left|\right|log\left(max\right(\frac{x\left(N_i\right)}{x\left(N_j\right)},\frac{x\left(N_j\right)}{x\left(N_i\right)}\left)\right)|{|}_{i1}^2+|\left|log\left(max\right(\frac{x\left(N_i\right)}{x\left(N_j\right)},\frac{x\left(N_j\right)}{x\left(N_i\right)}\left)\right)\right|{|}_{i2}^2}{4*h^2})}{4*\mathrm{\&pi;}*h^2}$

其中，i2是傅里叶核在以像素点i为中心，大小为7×7的局部分析窗N_i内构造的权重，是由傅里叶核构造的以i2为权重的加权距离，Z(i)的表达式如下所示：

$Z\left(i\right)=\underset{j\&Element;S_i}{\&Sigma;}\frac{\exp (-\frac{\left|\right|log\left(max\right(\frac{x\left(N_i\right)}{x\left(N_j\right)},\frac{x\left(N_j\right)}{x\left(N_i\right)}\left)\right)|{|}_{i1}^2+|\left|log\left(max\right(\frac{x\left(N_i\right)}{x\left(N_j\right)},\frac{x\left(N_j\right)}{x\left(N_i\right)}\left)\right)\right|{|}_{i2}^2}{4*h^2})}{4*\mathrm{\&pi;}*h^2},$

其中，0≤w(i,j)≤1且Z(i)表示归一化常数，S_i是以像素点i为中心，大小为5×5的搜索区域，h是平滑参数，取值为1/3，N_i和N_j分别表示以第i个像素点和第j个像素点为中心的5×5的局部分析窗，x(N_i)和x(N_j)分别是局部分析窗N_i和N_j内像素点的灰度值向量；

2d)对以第i个像素点为中心的搜索区域S_i内所有像素点的灰度值进行加权平均，求得图像中第i个像素点的新的灰度值INL(i)：

INL(i)＝∑w(i,j)X(j),j∈S_i,

其中，X(j)是去噪声前原图像中第j个像素点的灰度值，w(i,j)是与X(j)对应的权值，表示第i个像素点和第j个像素点之间的相似性；

2e)对变化前的SAR原图像Q₁和变化后的SAR图像Q₂分别重复步骤2a)到步骤2d)，直到两个图像中每个像素点都获得新的灰度值；

2f)对变化前的SAR原图像Q₁和变化后的SAR图像Q₂分别用获得的每个像素点的新的灰度值替代该像素点的原灰度值，即可得到不含噪的变化前的SAR原图像Q₃和不含噪的变化后的SAR图像Q₄；

(3)利用比值算子得到不含噪的原图像Q₃和不含噪的变化后图像Q₄的差异图Q₅；

(4)利用FLICM聚类算法对差异图Q₅的灰度值进行类别数为2的聚类分割，其中一类为变化类，另一类为非变化类；

(5)输出由步骤(4)得到的变化类图像，即为变化前SAR的图像Q₁和变化后SAR图像Q₂的变化检测结果图。

2.根据权利要求书中1中所提及的方法，其中步骤(3)所述的利用比值算子得到不含噪的原图像Q₃和不含噪的变化后图像Q₄的差异图Q₅，按照如下公式计算：

$Q_5=1-min(\frac{Q_3}{Q_4},\frac{Q_4}{Q_3}).$