CN103927737A - 基于非局部均值的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于非均值滤波的SAR图像变化检测方法。包括对同一地域不同时间获取的两幅SAR图像进行预处理；利用预处理后两幅SAR图像构造比值差异影像图；遍历比值差异影像图每个像素，计算每个像素点的平滑指数矩阵；对预处理后两幅SAR图像分别进行非局部均值滤波后作比值运算得非局部均值滤波比值图；用平滑指数作为权重，将比值差异影像图和非局部均值滤波比值图像求和得到最终的差异影像图；运用模糊局部C均值聚类法分割该最终的差异影像图得到变化检测结果图。本发明利用图像平滑指数特性保持差异图边缘信息，在平滑指数小的匀质区域，用非局部均值修正像素，有效抑制了噪声，更好地表现真实变化信息，提高了变化检测结果精确度。

Description

基于非局部均值的SAR图像变化检测方法

技术领域

本发明属于遥感图像变化检测领域，具体地说，涉及一种基于非局部均值思想构造差异影像的SAR（Synthetic Aperture Radar，合成孔径雷达）图像变化检测方法。

背景技术

合成孔径雷达（SAR）图像变化检测是从不同时间获取同一地理区域的多时相遥感影像，定性地分析和确定地表变化过程和特征的技术，与光学成像相比，合成孔径雷达（SAR）具有全天候，全天时获取数据的能力，同时具有穿越一定的植被以及遮盖物的能力，与光学成像相比，其更容易辨别地面的目标，一般将其作为光学传感器的有力补充。SAR图像变化检测技术正广泛的应用于各个领域，例如环境监控，农业研究，城市区域研究，森林监控等方面。近年来逐渐出现了许多种多时相SAR图像变化检测方法，最常用的SAR图像变化检测技术主要由两步组成：首先对多时相SAR图像构造差异影像图，然后在差异影像图基础上提取变化区域。

差异影像构造过程是SAR图像变化检测的重要部分，通过对两幅不同时间获取的同一地域的SAR图像进行比较处理得到差异影像信息图，差异影像信息图的精度高低直接影响变化检测结果的性能。SAR图像变化检测中对经过配准和辐射校正的SAR图像进行比较构造差异影像图是至关重要且必不可少的一步，差异影像图的质量直接决定后续分析处理的精确度，进而影响到整个SAR图像变化检测***的性能。

现有的变化检测方法往往不能很好的去除SAR图像中固有的噪声问题，基于比值法的变化检测方法对变化区域不敏感，漏检率较高，基于均值比的变化检测方法易将未变化区域归于变化区域，误减率较高。用现有的差异影像构造方法构造出的差异影像图往往噪声较多，在使用一些去噪方法去噪后经常会造成图像信息丢失较多，构造出的差异影像往往质量不高，很难得到既包含大量图像细节变化信息又克服噪声影响的差异影像图，造成变化检测结果误差很大，影响对检测结果正确率的分析。

总而言之，目前需要本领域技术人员迫切解决的一个技术问题是：如何最大程度的突出变化区域信息而有效的抑制未变化区域的信息，构造出既包含有大量有效信息又基本不受噪声影响的差异影像图，从中提取变化信息，提高变化检测结果的精确度。

发明内容

有鉴于此，本发明所要解决的技术问题在于，针对现有的差异图生成方法产生的差异图质量不高，信息丢失较多的问题，提出了基于非局部均值构造差异影像图，用于SAR图像变化检测的方法，针对SAR图像的特点，用像素平滑指数作为权重，将非局部均值的思想引入差异影像构造过程，构造包含更多有效信息并能在一定程度上抑制噪声的差异影像图，提高SAR图像变化检测精度。

本发明提供了一种基于非局部均值的SAR图像变化检测方法，具体步骤包括如下：

(1)对同一地域不同时间获取的两幅SAR图像进行滤波去噪，辐射校正和几何配准的预处理，得到预处理后的两幅SAR图像I₁，I₂；

(2)利用预处理后的两幅SAR图像I₁和I₂，构造比值差异影像图

(3)遍历比值差异影像图D_R的每个像素，计算每个像素点的平滑指数矩阵其中，μ(x)为以像素点i为中心的邻域内的像素值均值，σ(x)为以像素点i为中心的邻域内像素值的方差；

(4)对预处理后的SAR图像I₁和I₂分别进行非局部均值滤波，得到图像NL(I₁)和NL(I₂)；

(5)将经过非局部均值滤波后的图像NL(I₁)和NL(I₂)作比值运算，得到非局部均值滤波比值图D_NR；

(6)将平滑指数作为权重，将比值差异影像图D_R和非局部均值滤波比值图像D_NR进行求和运算，得到最终的差异影像图

（7）运用模糊局部C均值聚类方法对最终的差异影像图DI进行图像分割，生成变化检测结果图，完成对两幅SAR图像变化信息的最终检测。

所述步骤（4）中对预处理后的SAR图像I₁和I₂分别进行非局部均值滤波，步骤如下：

对SAR图像I₁进行非局部均值运算，遍历图像I₁每个像素点，计算其中是指在SAR图像I₁中以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，是像素点i和在搜索窗口内像素点p的相似度权重，且满足0≤w_ip≤1和 ${\mathrm{\Σ}}_{p\∈W_i^r}{w}_{\mathrm{ip}}=1,w_{\mathrm{ip}}(p\∈W_i^r)$ 由公式 $w_{\mathrm{ip}}=\frac{1}{Z_i}\exp (-\underset{k=1}{\overset{{(2s+1)}^2}{\mathrm{\Σ}}}\frac{1}{h}\log (\frac{A_{i,k}}{A_{p,k}}+\frac{A_{p,k}}{A_{i,k}}))$ 求得，其中s为邻域窗口半径且s＝3，参数h用于控制指数函数的衰减，理论上讲，非局部均值要取遍图像内每个点的7×7邻域块，由于在图像较大的情况下，这样的时间复杂度太高，故通常只选取像素点附近较大的一块区域（即搜索窗口）进行非局部均值运算，在本发明中令r＝10，即在一个21×21的区域内进行非局部运算，A_i,k，A_p,k分别代表以像素点i和像素点p为中心的第k个像素点的像素值。为对SAR图像I₁非局部滤波后的像素矩阵中第i个像素点的像素值，得到SAR图像I₁的非局部均值滤波图像NL(I₁)。

对SAR图像I₂进行非局部均值运算，遍历图像I₂每个像素点，计算其中是指在SAR图像I₂中以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，是像素点i和在搜索窗口内像素点p的相似度权重，且满足0≤w_ip≤1和 ${\mathrm{\Σ}}_{p\∈W_i^r}{w}_{\mathrm{ip}}=1,w_{\mathrm{ip}}(p\∈W_i^r)$ 由公式 $w_{\mathrm{ip}}=\frac{1}{Z_i}\exp (-\underset{k=1}{\overset{{(2s+1)}^2}{\mathrm{\Σ}}}\frac{1}{h}\log (\frac{A_{i,k}}{A_{p,k}}+\frac{A_{p,k}}{A_{i,k}}))$ 求得，其中s为邻域窗口半径且s＝3，参数h用于控制指数函数的衰减，令r＝10，即在一个21×21的区域内进行非局部运算，A_i,k，A_p,k分别代表以像素点i和像素点p为中心的第k个像素点的像素值。为对SAR图像I₂非局部滤波后的像素矩阵中第i个像素点的像素值，得到SAR图像I₁的非局部均值滤波图像NL(I₂)。

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

1、本发明采用平滑指数作为权重，在构造差异影响图的过程中，引入非局部均值修正过的像素值，结合本来的比值信息和非局部思想对图像进行修正，贡献一定的权重，构造出了包含更多细节信息的较高质量的差异影像图。

2、本发明在构造差异影像图的过程中结合了非局部均值运算的方法，较好地去除了差异影像图中的噪点，提高了差异影像图的性能，增加了变化检测结果的精确度。

3、本发明采用模糊局部C均值的方法分割差异影像图，可以更大程度的抑制变化区域的背景信息，增强变化区域的变化信息，从而获得更高的检测精度。

4、本发明引入非局部均值滤波的方法构造既保持变化细节信息又充分抑制噪声的差异信息图，便于后期差异图分析，提高SAR图像变化检测中的检测精度，减少错误率。

附图说明

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

图2是本发明仿真使用的Bern地区两幅SAR图像和一幅标准参考图；

图3是本发明基于非局部均值方法构造的Bern地区两幅SAR图像差异影像图；

图3是用本发明和现有的基于比值法的变化检测方法、基于均值比的变化检测方法分别对Bern地区两幅SAR图像进行变化检测的仿真结果图。

具体实施方式

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

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

步骤1，对同一地域不同时间获取的两幅SAR图像进行滤波去噪，辐射校正和几何配准的预处理，得到预处理后的两幅SAR图像I₁，I₂。

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

步骤2，利用预处理后的两幅SAR图像I₁和I₂，构造比值差异影像图

将图像I₁中位于i行j列的像素点(i,j)的灰度值I₁(i,j)和相对应的图像I₂中位于i行j列的像素点的灰度值I₂(i,j)，通过比值运算D_R(i,j)＝I₁(i,j)I₂(i,j)，得到比值法差异图D_R中位于i行j列的像素点(i,j)的灰度值D_R(i,j)；对图像I₁和图像I₂中每个位于i行j列的像素点的灰度值从左到右，从上到下都进行比值运算，构造出比值差异影像图D_R。

步骤3，遍历比值差异影像图D_R的每个像素，计算每个像素点的平滑指数矩阵其中，μ(x)为以像素点i为中心的邻域内的像素值的均值，σ(x)为以像素点i为中心的邻域内像素值的方差。

这里，其中，x_i代表像素点为中心的半径为n的邻域内第i个像素点的像素灰度值，平滑指数矩阵为与比值差异影像图规模一样的矩阵。

据平滑指数特性，在平滑指数大的地方一般是图像边缘，平滑指数小的地方为匀质区域。本发明在差异影像图构造阶段利用图像平滑指数特性，在平滑指数大的地方一般是图像边缘，差异影像图的像素值权重较大；平滑指数小的地方作为匀质区域常含有较多冗余信息，用非局部均值对其像素进行修正后能更好的表示真实情况，这样得到的SAR图像变化检测结果的性能最佳。

步骤4，对预处理后的SAR图像I₁和I₂分别进行非局部均值滤波，得到图像NL(I₁)和NL(I₂)。

4.1对SAR图像I₁进行非局部均值运算，遍历图像I₁每个像素点，计算其中是指在SAR图像I₁中以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，是像素点i和在搜索窗口内像素点p的相似度权重，且满足0≤w_ip≤1和 ${\mathrm{\Σ}}_{p\∈W_i^r}{w}_{\mathrm{ip}}=1,w_{\mathrm{ip}}(p\∈W_i^r)$ 由公式 $w_{\mathrm{ip}}=\frac{1}{Z_i}\exp (-\underset{k=1}{\overset{{(2s+1)}^2}{\mathrm{\Σ}}}\frac{1}{h}\log (\frac{A_{i,k}}{A_{p,k}}+\frac{A_{p,k}}{A_{i,k}}))$ 求得，其中s为邻域窗口半径且s＝3，参数h用于控制指数函数的衰减，理论上讲，非局部均值要取遍图像内每个点的7×7邻域块，由于在图像较大的情况下，这样的时间复杂度太高，故通常只选取像素点附近较大的一块区域（即搜索窗口）进行非局部均值运算，在本发明中令r＝10，即在一个21×21的区域内进行非局部运算，A_i,k，A_p,k分别代表以像素点i和像素点p为中心的第k个像素点的像素值。为对SAR图像I₁非局部滤波后的像素矩阵中第i个像素点的像素值，得到SAR图像I₁的非局部均值滤波图像NL(I₁)。

4.2对SAR图像I₂进行非局部均值运算，遍历图像I₂每个像素点，计算其中是指在SAR图像I₂中以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，是像素点i和在搜索窗口内像素点p的相似度权重，且满足0≤w_ip≤1和 ${\mathrm{\Σ}}_{p\∈W_i^r}{w}_{\mathrm{ip}}=1,w_{\mathrm{ip}}(p\∈W_i^r)$ 由公式 $w_{\mathrm{ip}}=\frac{1}{Z_i}\exp (-\underset{k=1}{\overset{{(2s+1)}^2}{\mathrm{\Σ}}}\frac{1}{h}\log (\frac{A_{i,k}}{A_{p,k}}+\frac{A_{p,k}}{A_{i,k}}))$ 求得，其中s为邻域窗口半径且s＝3，参数h用于控制指数函数的衰减，令r＝10，即在一个21×21的区域内进行非局部运算，A_i,k，A_p,k分别代表以像素点i和像素点p为中心的第k个像素点的像素值。为对SAR图像I₂非局部滤波后的像素矩阵中第i个像素点的像素值，得到SAR图像I₁的非局部均值滤波图像NL(I₂)。

非局部均值是近年来常用的图像去噪方法，由于传统局部均值只考虑到像素周边部分的像素值对本身的影响，会导致弱化边缘，而双边滤波中只涉及图像中像素点之间的距离和相似度对像素值的作用，没有顾及像素点周边的整体信息。非局部较好的综合以上两种滤波方式的优点，既考虑到图像邻域块对本身的指导作用，又考虑到全图中其他非邻域的像素值对本身的影响。通过像素点邻域块与其他像素点邻域块之间的相似性修正像素值，可以保留图像细节，充分抑制图像噪声。非局部均值结合图像的平滑指标，可以在平滑指数较大的边缘区域保留较多本身像素值，在平滑指数较小的匀质区域取非局部均值成分更多，这样，既可以利用非局部均值对噪声的有效滤除，又能充分保留差异图像边缘信息。

步骤5，将经过非局部均值滤波后的图像NL(I₁)和NL(I₂)作比值运算，得到非局部均值滤波比值图像D_NR。

将图像NL(I₁)中位于i行j列的像素点(i,j)的灰度值NL（I₁(i,j)）和相对应的图像NL(I₂)中位于i行j列的像素点的灰度值NL（I₂(i,j)），通过比值运算D_NR(i,j)＝NL（I₁(i,j)）/NL（I₁(i,j)），得到非局部均值滤波比值图像D_NR中位于i行j列的像素点(i,j)的灰度值D_NR(i,j)；对图像NL(I₁)和图像NL(I₂)I₁中每个位于i行j列的像素点的灰度值从左到右，从上到下都进行比值运算，构造出非局部均值滤波比值图像D_NR。

步骤6，将平滑指数作为权重对比值差异图和非局滤波比值图像求和，对两幅图像上每个坐标对应点作加权求和运算 $\mathrm{DI}(i,j)=\∂(i,j)*D_R(i,j)+(1-\∂(i,j))*D_{\mathrm{NR}}(i,j),$ 其中DI(i,j)为求和差异图DI中坐标为(i,j)的像素点的像素灰度值，得到最终的差异影像图DI，即SAR图像I₁和图像I₂的差异信息图。

平滑指数是评价图像的重要指标，是对每个像素点邻域内的方差和均值的比值，像素点的平滑指数越大，代表此像素点为图像边缘部分；像素点的平滑指数越小，代表此像素点为处于图像中非边缘的匀质区域。相对而言，匀质区域内冗余信息较多，可以通过以平滑指数为权重，将非局部均值修正过的像素值引入到差异图构造过程中，即在生成差异图的过程中，可以结合本来的比值信息和非局部思想对图像进行修正，贡献一定的权重，产生更为合理的差异图。

,步骤7，运用模糊局部C均值聚类方法对最终的差异影像图DI进行图像分割，生成变化检测结果图，完成对两幅SAR图像变化信息的最终检测。

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

1、仿真条件：

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

2、仿真参数

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

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

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

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

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

为了验证基于非局部均值滤波的SAR图像变化检测方法的优越性，将本发明方法与传统的基于比值法的变化检测方法和基于均值比值法算法的变化检测方法的性能做出对比和定量分析。这里把基于比值法的变化检测方法简称为LR，基于均值比值法算法的变化检测方法简称为MR，本发明方法简称为NLR。

将本发明中的变化检测方法应用在如图2所示真实Bern地区的SAR图像数据集上，并与现有的基于比值法的变化检测方法和基于均值比值法算法的变化检测方法进行对比分析。

如如图2所示的Bern地区三幅图像，其中图2(a)表示1999年4月Bern地区的地貌信息，即为第一时间获取图像。图2(b)表示1999年5月Bern地区的地貌信息即为第二时间获取图像。图2(c)表示变化检测的标准参考图。

用图2（a）和图2（b）所示的两幅SAR图像分别使用基于比值法的变化检测方法、基于均值比值法算法的变化检测方法和本发明方法得到变化检测结果图，如图3所示，其中图3(a)为基于比值法的变化检测方法(LR)结果图，图3(b)为基于均值比的变化检测方法(MR)结果图，图3(c)代表本发明(NLR)方法的结果图。

将图3所示的三种方法的性能指标进行对比，如表1所示：

表1Bern地区三种变化检测方法结果性能比较

从表1中可以看出，本发明方法的变化检测结果对图像中像素点的漏检数最少，总错误数最低，提高了变化检测结果的精度，显示了本发明方法的优越性。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (3)

1.一种基于非局部均值滤波的SAR图像变化检测方法，其特征在于，包括如下步骤：

(1)对同一地域不同时间获取的两幅SAR图像进行滤波去噪，辐射校正和几何配准的预处理，得到预处理后的两幅SAR图像I₁，I₂；

(2)利用预处理后的两幅SAR图像I₁和I₂，构造比值差异影像图

(3)遍历比值差异影像图D_R的每个像素，计算每个像素点的平滑指数矩阵其中，μ(x)为以像素点i为中心的邻域内的像素值均值，σ(x)为以像素点i为中心的邻域内像素值的方差；

(4)对预处理后的SAR图像I₁和I₂分别进行非局部均值滤波，得到图像NL(I₁)和NL(I₂)；

(5)将经过非局部均值滤波后的图像NL(I₁)和NL(I₂)作比值运算，得到非局部均值滤波比值图D_NR；

(6)将平滑指数作为权重，将比值差异影像图D_R和非局部均值滤波比值图像D_NR进行求和运算，得到最终的差异影像图

（7）运用模糊局部C均值聚类方法对最终的差异影像图DI进行图像分割，生成变化检测结果图，完成对两幅SAR图像变化信息的最终检测。

2.如权利要求1所述的步骤（3）中遍历比值差异影像图D_R的每个像素，计算每个像素点的平滑指数矩阵，其特征在于，所述平滑指数矩阵为与比值差异影像图规模一样的矩阵

3.如权利要求1所述步骤（4）中对预处理后的SAR图像I₁和I₂分别进行非局部均值滤波，得到图像NL(I₁)和NL(I₂)其特征是，所述对预处理后的SAR图像I₁和I₂分别进行非局部均值滤波的步骤为：

对SAR图像I₁进行非局部均值运算，遍历图像I₁每个像素点，计算其中是指在SAR图像I₁中以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，是像素点i和在搜索窗口内像素点p的相似度权重，且满足0≤w_ip≤1和 ${\mathrm{\Σ}}_{p\∈W_i^r}{w}_{\mathrm{ip}}=1,w_{\mathrm{ip}}(p\∈W_i^r)$ 由公式 $w_{\mathrm{ip}}=\frac{1}{Z_i}\exp (-\underset{k=1}{\overset{{(2s+1)}^2}{\mathrm{\Σ}}}\frac{1}{h}\log (\frac{A_{i,k}}{A_{p,k}}+\frac{A_{p,k}}{A_{i,k}}))$ 求得，其中s为邻域窗口半径且s＝3，参数h用于控制指数函数的衰减，理论上讲，非局部均值要取遍图像内每个点的7×7邻域块，由于在图像较大的情况下，这样的时间复杂度太高，故通常只选取像素点附近较大的一块区域（即搜索窗口）进行非局部均值运算，在本发明中令r＝10，即在一个21×21的区域内进行非局部运算，A_i,k，A_p,k分别代表以像素点i和像素点p为中心的第k个像素点的像素值。为对SAR图像I₁非局部滤波后的像素矩阵中第i个像素点的像素值，得到SAR图像I₁的非局部均值滤波图像NL(I₁)；

对SAR图像I₂进行非局部均值运算，遍历图像I₂每个像素点，计算其中是指在SAR图像I₂中以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，是像素点i和在搜索窗口内像素点p的相似度权重，且满足0≤w_ip≤1和 ${\mathrm{\Σ}}_{p\∈W_i^r}{w}_{\mathrm{ip}}=1,w_{\mathrm{ip}}(p\∈W_i^r)$ 由公式 $w_{\mathrm{ip}}=\frac{1}{Z_i}\exp (-\underset{k=1}{\overset{{(2s+1)}^2}{\mathrm{\Σ}}}\frac{1}{h}\log (\frac{A_{i,k}}{A_{p,k}}+\frac{A_{p,k}}{A_{i,k}}))$ 求得，其中s为邻域窗口半径且s＝3，参数h用于控制指数函数的衰减，令r＝10，即在一个21×21的区域内进行非局部运算，A_i,k，A_p,k分别代表以像素点i和像素点p为中心的第k个像素点的像素值。为对SAR图像I₂非局部滤波后的像素矩阵中第i个像素点的像素值，得到SAR图像I₁的非局部均值滤波图像NL(I₂)。