CN102930519B - 基于非局部均值的sar图像变化检测差异图生成方法 - Google Patents

Abstract

本发明公开一种基于非局部均值的SAR图像变化检测差异图生成方法。实现过程主要包括：先对经过预处理的两幅不同时间相同地域的SAR图像构造比值图，然后求出比值图的平滑指数矩阵，接着对两幅SAR图像运用非局部均值滤波，再将非局部均值处理后的两幅图像作比值运算，再将平滑指数作为权重对比值和非局部均值滤波比值图求和，最后生成差异图像。本发明在生成差异图阶段利用图像平滑指数特性，在平滑指数大的图像边缘，充分发挥像素值本身的决定性作用，保持差异图边缘信息，在平滑指数小的匀质区域，用非局部均值对其像素进行修正，可以有效抑制噪声，更好地表示变化真实情况，得到质量更佳的差异信息图，保证了后续分析能有较好的结果。

Description

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

技术领域

本发明属于SAR图像变化检测技术领域，涉及SAR图像变化检测中的差异图生成阶段。具体地说是提出了一种基于非局部思想的SAR图像变化检测差异图生成方法，用于生成包含更多有效信息，并能一定程度上抑制噪声的差异信息图，便于后续分析处理，提高SAR图像变化检测精度。

背景技术

SAR图像变化检测是从不同时间获取同一地理区域的多时相遥感影像，定性地分析和确定地表变化过程和特征的技术。由于与光学遥感系统相比，SAR系统具有全天时、全天候获取数据的能力，所以SAR图像变化检测技术正广泛的应用于各个领域，例如环境监控，农业研究，城市区域研究，森林监控等方面。

差异图生成是SAR图像变化检测的重要部分，通过对两幅不同时间在同一地域上的SAR图像进行比较得到差异信息图，以供后续分析产生变化/未变化二值图像，差异信息图精度高低也直接影响变化检测的性能。SAR图像变化检测中，差异图的生成是对经过预处理(包括图像配准和辐射校正)的SAR图像进行算术运算，主要经过差值运算，比值运算，对数比值运算得到初步的差异信息图，通过对信息图进行图像分割就可以得到最终的二值图像。

非局部均值思想常用于图像去噪，是对双边滤波的一个推广，图像中往往包含有许多冗余信息，充分利用这些冗余信息为去除图像噪声服务，这是非局部平均滤波模型的主要思想。冗余信息即是指图像中部分区域灰度的相似程度，根据相似度来进行平滑去噪是非局部均值图像去噪的一个优点。非局部均值(NL-means)模型的主要特点是：该方法不是用图像中单个像素的灰度值进行比较，而是对该像素周围的整个灰度的分布状况进行比较，根据灰度分布的相似性来贡献权值。

SAR图像变化检测中对经过配准和辐射校正的SAR图像进行比较生成差异图是至关重要且必不可少的一步，差异图的质量直接决定后续分析处理的精确度，进而影响到整个SAR图像变化检测系统的性能。

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

发明内容

本发明的目的在于：针对现有的差异图生成方法产生的差异图质量不高，信息丢失较多的问题，提出了基于非局部均值生成差异图方法，针对SAR图像变化检测的特点，以像素平滑指数作权重将非局部均值引入差异图生成过程，对初步差异图进行修正，提高检测精度，生成了包含有效信息更多的差异图，与其它现有的方法相比能够产生更便于后续处理的差异图像。一句话，本发明要解决的技术问题就是提高差异图精度和质量，使生成的SAR差异信息图更能体现变化区域的信息。

本发明的技术方案是：基于非局部均值SAR图像变化检测差异图生成具体实现步骤包括有：

步骤1通过星载合成孔径雷达获取两幅不同时间相同地域的SAR图像，将获取的两幅SAR图像，输入到安装有矩阵实验室和Visual C++6.0软件的计算机中，利用相关软件处理，相关软件包括矩阵实验室和Visual C++6.0软件和计算机配置的软件工具，处理包括：滤波去噪，辐射校正与几何配准的预处理，得到SAR图像I₁和图像I₂。

步骤2采用Matlab或C++中其中之一编程，或使用Matlab和C++混合编程实现对两幅SAR图像I₁和图像I₂构造比值差异图采用Matlab或C++中其中之一单独编程，或者采用Matlab和C++混合编程均可实现构造比值差异图。

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

步骤4对两幅不同时间相同地域的SAR图像I₁、图像I₂分别进行非局部均值修正，得到经过非局部均值滤波过的SAR图像NL(I₁)和NL(I₂)，在计算中，NL(I₁)和NL(I₂)分别为SAR图像I₁、图像I₂每个像素点经过非局部均值处理后生成的新的像素矩阵，每个像素点的非局部均值修正像素值按以下公式计算，其中是指以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，是像素点i和在搜索窗口内像素点p的相似度权重，且满足0≤w_ip≤1和 为非局部滤波后的像素矩阵中第i个像素点的像素值。

步骤5将经过非局部均值滤波过的SAR图像NL(I₁)和NL(I₂)作比值运算得到非局部均值滤波比值图D_NR，其中NL(I₁)，NL(I₂)分别为SAR图像I₁和图像I₂经过非局部均值滤波后的图像。经过非局部均值滤波后的图像包含更多有效变化细节，本发明充分利用非局部均值对图像噪声的抑制作用，将非局部均值引入差异图生成阶段，通过非局部均值修正初步差异图，获得包含更多有效变化细节并克服噪声的高质量差异信息图率。

步骤6将平滑指数作为权重对比值差异图D_R和非局部修正后的比值图像D_NR求和，得到最终的差异图像即图像DI为SAR图像I₁和图像I₂的差异信息图，保存数据，作为下一步差异图分析的图像源。

本发明首先对两幅不同时间相同地域的SAR图像构造比值差异图，然后求出比值差异图的各个像素点的平滑指数，接着对两幅不同时间相同地域的SAR图像利用非局部均值方法修正各个像素值，再将非局部均值滤波后的两幅图像作比值，再以平滑指数作为权重对比值图和非局部滤波比值图求和，最后得到差异信息图像。

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

本发明的实现还在于：步骤4中对于SAR图像进行非局部均值修正具体实现步骤包括有：

4.1遍历SAR图像I每个像素点i，计算像素点i和在搜索窗口内像素点p的相似度权重且相似度权重满足0≤w_ip≤1和其中，s为邻域窗口半径且s＝3，参数h用于控制指数函数的衰减，理论上讲，非局部均值要取遍图像内每个点的7×7邻域块，由于在图像较大的情况下，这样的时间复杂度太高，故通常只选取像素点附近较大的一块区域(即搜索窗口)进行非局部均值运算，在本发明中令r＝10，即在一个21×21的区域内进行非局部运算，A_i，k，A_p，k分别代表以像素点i和像素点p为中心的第k个像素点的像素值。

4.2遍历SAR图像I每个像素点i，进行非局部均值运算，其中，为对SAR图像I非局部滤波后的像素矩阵中第i个像素点的像素值，得到SAR图像I的非局部均值滤波图像NL(I)。

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

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

1、本发明在生成差异图阶段利用图像平滑指数特点，在平滑指数大的地方一般是图像边缘，差异图像素值本身对像素值有决定性影响；平滑指数小的地方作为匀质区域含有较多冗余信息，用非局部均值对其像素进行修正后能更好的表示真实情况，与现有方法相比，本发明方法得到的SAR图像变化检测结果的性能最佳。

2、本发明在生成差异图过程中结合非局部均值运算较好地去除差异图噪点，增加了差异图的精确度，从而保证了后续分析能有较好的结果。

3、本发明与其它经典的SAR图像变化检测阈值方法相比，获得了较低的检测错误率，同时ROC曲线下的面积最大，验证了本发明的有效性。

附图说明

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

图2是Bern地区两幅SAR影像图、标准图和变化参考图；

图3是本发明变化检测差异图生成方法和现有技术两种方法对Bern地区SAR影像的实验结果图；

图4采用本发明获得的ROC曲线与其他方法的ROC曲线的比较图；

图5黄河入海口地区两幅SAR影像图、标准图和变化参考图；

图6本发明变化检测差异图生成方法和现有技术两种方法对黄河入海口地区SAR影像的实验结果图；

图7采用本发明获得的ROC曲线与其他方法的ROC曲线的比较图。

具体实施方式

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

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

本发明现阶段可以运行在32位XP系统(及以上)计算机，matlab7.0及以上平台，Visual C++6.0平台。

实施例1

本发明是一种基于非局部均值的SAR图像变化检测差异图生成方法，参见图1，首先输入两幅经由预处理的不同时间相同地域的SAR图像，对两幅SAR图像作比值运算，然后计算比值差异图的各个像素点的平滑指数，接着对两幅SAR图像进行非局部均值滤波，修正各个点像素值，再将非局部修正后的两幅图像作比值运算，接着将平滑指数作为权重对比值和非局部修正后的图构成的新比值图求和，最后得到差异图像，基于非局部均值的SAR图像变化检测差异图生成具体实现步骤包括有：

步骤1通过星载合成孔径雷达获取两幅不同时间相同地域的SAR图像，输入获取的两幅不同时间相同地域的SAR图像，将获取的两幅SAR图像输入32位及以上XP系统(及以上)计算机，利用Matlab7.0(及更高版本)和Visual C++6.0软件处理经过滤波去噪，辐射校正与几何配准的预处理SAR图像I₁和图像I₂。

经过滤波去噪，辐射校正与几何配准的预处理SAR图像I₁和图像I₂。

步骤2采用Matlab或C++中其中之一编程，或使用Matlab和C++混合编程实现对两幅SAR图像I₁和图像I₂构造比值差异图，本例中采用Matlab单独编程实现对两幅SAR图像I₁和图像I₂构造比值差异图。

比值差异图构造方法通过公式生成比值差异图像D_R，使得在比值差异图像上低灰度级呈现为无变化区域，高灰度级呈现为变化区域。

步骤3遍历比值差异图像D_R的每个像素，计算图像上每个像素点的平滑指数矩阵其中 其中，x_i代表像素点为中心的半径为n的邻域内第i个像素点的像素灰度值，据平滑指数特性，在平滑指数大的地方一般是图像边缘，平滑指数小的地方为匀质区域。平滑指数矩阵是平滑指数的集合。

步骤4对两幅不同时间相同地域的SAR图像I₁、图像I₂分别进行非局部均值修正，得到经过非局部均值滤波过的SAR图像NL(I₁)和NL(I₂)，在计算中，NL(I₁)和NL(I₂)分别为SAR图像I₁、图像I₂每个像素点经过非局部均值处理后生成的新的像素矩阵，每个像素点的非局部均值修正像素值按以下公式计算，其中是指以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，是像素点i和在搜索窗口内像素点p的相似度权重，且满足0≤w_ip≤1和 为非局部滤波后的像素矩阵中第i个像素点的像素值。

步骤5将非局部均值处理过的SAR图像NL(I₁)，NL(I₂)中每个像素点作比值运算得到非局部滤波比值图，使得在非局部滤波比值图像上低灰度级呈现为无变化区域，高灰度级呈现为变化区域。本发明利用非局部均值对图像噪声的抑制作用，通过非局部均值修正初步差异图，得到包含更多有效变化细节并克服噪声的高质量差异信息图，保证后续分析的正确率。

步骤6将平滑指数作为权重对比值差异图和非局部滤波比值图像求和，对两幅图像上每个对于点作加权求和运算得到最终的差异图像DI，即SAR图像I₁，图像I₂的差异信息图，保存数据，作为下一步差异图分析的图像源。

SAR图像的差异信息图包含两个时刻之间的变化信息，是作下一步分析的基础，经过差异信息图分析等后续处理，能到得到质量有保证的SAR图像变化检测结果。差异信息图的质量直接影响整个SAR图像变化检测系统的精度。

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

实施例2

基于非局部均值的SAR图像变化检测差异图生成方法同实施例1，参照图1，实现本发明SAR图像变化检测差异图生成，首先对两幅不同时间相同地域的SAR图像构造比值差异图，然后求出比值差异图的各个像素点的平滑指数，接着对两幅不同时间相同地域的SAR图像利用非局部均值方法修正各个像素值，再将非局部修正后的两幅图像作比值，再将平滑指数作为权重对比值和非局部修正后的图构成的新比值图求和，最后得到差异图像。下边通过本例对该发明的实现过程进行详细说明：

步骤1通过星载合成孔径雷达获取两幅不同时间相同地域的SAR图像，将获取的两幅SAR图像输入到32位及以上XP系统(及以上)计算机中，经由Matlab7.0(及更高版本)和Visual C++6.0软件处理，将这两幅SAR图像进行滤波去噪，辐射校正与几何配准的预处理，得到处理后的两幅图像I₁，I₂。

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

步骤2使用Matlab和C++混合编程实现对SAR图像I₁，I₂构造比值差异图像。

构造方法通过公式生成比值差异图像D_R，使得在比值差异图像上低灰度级呈现为无变化区域，高灰度级呈现为变化区域，其中x₁(l，t)，x₁(l，t)分别表示在SAR图像I₁，I₂中坐标为(l，t)的像素点的像素灰度值，D_R(l，t)为比值差异图像D_R中坐标为(l，t)的像素点的像素灰度值。

步骤3遍历比值差异图像D_R的每个像素点，即(k，t)遍历比值差异图像D_R的每个像素点计算图像上每个像素点的平滑指数其中 x_j代表以像素点(l，t)为中心的半径为n的邻域内第j个像素点的像素灰度值，平滑指数指标代表此像素点是否为匀值区域，n＝3。为比值差异图像D_R的平滑指数集合，为与比值差异图像规模一样的矩阵。

步骤4分别对两幅不同时间相同地域的SAR图像I₁，图像I₂进行非局部均值滤波。

4.1对SAR图像I₁进行非局部均值运算，遍历图像I₁每个像素点，计算其中是指在SAR图像I₁中以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，是像素点i和在搜索窗口内像素点p的相似度权重，且满足0≤w_ip≤1和 由公式 $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}}))$ 求得，其中 $Z_i=\underset{p\∈W_i^r}{\mathrm{\Σ}}\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和 由公式 $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}}))$ 求得，其中 $Z_i=\underset{p\∈W_i^r}{\mathrm{\Σ}}\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将非局部均值处理过的SAR图像NL(I₁)，NL(I₂)中每个坐标对应点作比值运算 $D_{\mathrm{NR}}(l,t)=\frac{\min (\mathrm{NL}\left(I_1\right)(l,t),\mathrm{NL}\left(I_2\right)(l,t))}{\max (\mathrm{NL}\left(I_1\right)(l,t),\mathrm{NL}\left(I_2\right)(l,t))}\×255,$ 其中NL(I₁)(l，t)和NL(I₂)(l，t)分别为SAR图像NL(I₁)，NL(I₂)中坐标为(l，t)的像素点的像素灰度值，D_NR(l，t)代表非局部滤波比值图D_NR中坐标为(l，t)的像素点的像素灰度值，这样使得在非局部滤波比值图像D_NR上低灰度级呈现为无变化区域，高灰度级呈现为变化区域。

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

本发明运用引入非局部均值，修正初步差异图，生成保持变化细节并充分抑制噪声的新的差异信息图，便于后期差异图分析，提高SAR图像变化检测中的检测精度，减少错误率。

实施例3

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

步骤1在Bern地域的不同时间获取两幅合成孔径雷达(SAR)图像，将这两幅SAR图像进行滤波去噪，辐射校正与几何配准的预处理，预处理后的两幅SAR图像为SAR图像I₁，SAR图像I₂，其中预处理后得到的图像I₁如图2(a)所示，图2(a)是1999年4月Bern地区的地貌信息，预处理后得到的图像I₂如图2(b)所示，图2(b)是1999年5月Bern地区的地貌信息。

步骤2利用预处理后的如图2(a)所示的图像I₁和图2(b)所示的图像I₂，用比值运算构造比值图D_R，详细实施步骤如实施例2中的步骤2。

步骤3计算比值图D_R上每个像素点的平滑指数矩阵详细实施步骤如实施例2中的步骤3。

步骤4分别对两幅不同时间相同地域的SAR图像I₁，I₂进行非局部均值滤波，得到非局部均值滤波图像NL(I₁)，NL(I₂)，详细实施步骤如实施例2中的步骤4。

步骤5将非局部均值处理过的SAR图像NL(I₁)，NL(I₂)中每个坐标对应点作比值运算，得到非局部滤波比值图像D_NR，详细实施步骤如实施例2中的步骤5。

步骤6将平滑指数作为权重对比值差异图和非局滤波比值图像求和，得到最终的差异图像DI，保存数据，作为下一步差异图分析的图像源。详细实施步骤如实施例2中的步骤6。

本实施例步骤由Matlab和C++混合编程实现，其中，步骤1及步骤4-6是通过Matlab实现，步骤2-3由C++语言实现，该实施例给出了一个具体实施方案，本领域专业人员可以通过本发明中的步骤指引实现本发明中的基于非局部均值的SAR图像变化检测差异图生成，可以直接沿用本发明各步骤中的编程语言，也可根据自身习惯选用其他高级语言实现本发明。

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

1、仿真参数

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

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

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

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

④正确率：

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

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

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

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

在图3中图3(a)代表均值比值(MR)图像，图3(b)代表对数比值(LR)图像，图3(c)代表本发明(NLR)方法产生的差异图像。

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

各种方法在真实Bern区域SAR图像数据的实验结果图如参考图4所示。

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

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

从表1中可以看出，本发明在变化检测总错误数上，相比其他对比算法是最少的，尤其漏检测数与其他方法相比有显著优势，说明了本发明方法的优越性。从图3中可以看到，MR方法对图像模糊较多，导致较高的漏检率，这种方法的漏检测数是三种方法中最高的，从表1中的漏检数也可以给出同样的结论，LR方法对变化区域有一定的缩小。而NLR方法是总错误率最小的。总体来看，通过表1性能结果和图3的视觉效果及参考图4ROC曲线形状说明了本发明方法生成的差异图性能良好，便于后续分析处理，提高最终变化检测的精确度，使得总错误数更少，误检数与漏检数上相对都比较均衡。综上所述本发明方法能够得到更好的SAR图像变化检测结果。

实施例4

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

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

步骤2利用预处理后的如图5(a)所示的图像I₁和图5(b)所示的图像I₂，用比值运算构造比值图D_R。

步骤3计算比值图D_R上每个像素点的平滑指数矩阵

步骤4分别对两幅不同时间相同地域的SAR图像I₁，I₂进行非局部均值滤波，得到非局部均值滤波图像NL(I₁)，AL(I₂)。

步骤5将非局部均值处理过的SAR图像NL(I₁)，NL(I₂)中每个坐标对应点作比值运算，得到非局部滤波比值图像D_NR。

步骤6将平滑指数作为权重对比值差异图和非局滤波比值图像求和，得到最终的差异图像DI，保存数据，作为下一步差异图分析的图像源。

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

1、仿真参数同实施例3

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

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

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

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

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

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

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

从表1中可以看出，本发明在变化检测总错误数上，相比其他对比算法是最少的，黄河图像本身受噪声影响程度较大，通过非局部方法参与修正，误检数与漏检数上都有很大程度抑制，说明了本发明方法的优越性。从图6中可以看到，MR方法对图像模糊较多，导致较高的错检率，从表2中的漏检数也可以给出同样的结论，LR方法对变化区域有一定的缩小，这两种方法的总错误数是四种方法中最高的。而本发明(NLR)方法对以上两种缺陷有一定的综合改善，NLR方法是总错误率最小的。总体来看，通过表2性能结果和参考图6的视觉效果及参考图7中ROC曲线形状说明了本发明方法生成的差异图性能良好，便于后续分析处理，提高最终变化检测的精确度，使得总错误数更少，误检数与漏检数上相对都比较均衡。综上所述本发明方法能够得到更好的SAR图像变化检测结果.

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

综上，本发明的基于非局部均值的SAR图像变化检测差异图生成方法，其实现过程主要包括：首先对两幅不同时间相同地域的SAR图像构造比值差异图，然后求出比值差异图的各个像素点的平滑指数，接着对两幅不同时间相同地域的SAR图像利用非局部均值方法修正各个像素值，再将非局部修正后的两幅图像作比值，再将平滑指数作为权重对比值和非局部修正后的图构成的新比值图求和，最后得到差异图像。本发明在生成差异图阶段利用图像平滑指数特点，在平滑指数大的地方一般是图像边缘，差异图像素值本身对像素值有决定性影响；平滑指数小的地方作为匀质区域含有较多冗余信息，用非局部均值对其像素进行修正后能更好的表示真实情况，本发明方法得到的SAR图像变化检测结果的性能最佳。本发明在生成差异图过程中引入非局部均值运算较好地去除差异图噪点，增加了差异图的精确度，从而保证了后续分析能有较好的结果。

Claims (1)

1.一种基于非局部均值SAR图像变化检测差异图生成方法，其特征在于：基于非局部均值SAR图像变化检测差异图生成具体实现步骤包括有：

步骤1通过星载合成孔径雷达获取两幅不同时间相同地域的SAR图像，将获取的两幅SAR图像，输入到安装有矩阵实验室软件的计算机中，利用相关软件处理，经过滤波去噪，辐射校正与几何配准的预处理，得到SAR图像I₁和图像I₂；

步骤2采用Matlab或C++中其中之一编程，或使用Matlab和C++混合编程实现对两幅SAR图像I₁和图像I₂构造比值差异图

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

步骤4对两幅不同时间相同地域的SAR图像I₁、图像I₂分别进行非局部均值修正，得到经过非局部均值滤波过的SAR图像NL(I₁)和NL(I₂)，在计算中，NL(I₁)和NL(I₂)分别为SAR图像I₁、图像I₂每个像素点经过非局部均值处理后生成的新的像素矩阵，每个像素点的非局部均值修正像素值按以下公式计算，其中W_i ^r是指以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，w_ip是像素点i和在搜索窗口内像素点p的相似度权重，其中p∈W_i ^r，且满足0≤w_ip≤1和 为非局部滤波后的像素矩阵中第i个像素点的像素值；对于SAR图像进行非局部均值修正具体实现步骤包括有：

4.1遍历SAR图像I每个像素点i，计算像素点i和在搜索窗口内像素点p的相似度权重且相似度权重满足0≤w_ip≤1和 ${\mathrm{\Σ}}_{p\∈{W_i}^r}{w}_{\mathrm{ip}}=1,$ 其中， $Z_i=\underset{p\∈W_i^r}{\mathrm{\Σ}}\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个像素点的像素值；

4.2遍历SAR图像I每个像素点i，进行非局部均值运算，其中，为对SAR图像I非局部滤波后的像素矩阵中第i个像素点的像素值，得到SAR图像I的非局部均值滤波图像NL(I)；

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

步骤6将平滑指数作为权重对比值差异图D_R和非局部修正后的比值图像D_NR求和，得到最终的差异图像即图像DI为SAR图像I₁和图像I₂的差异信息图，保存数据，作为下一步差异图分析的图像源。