CN102938071A - 基于非局部均值的sar图像变化检测模糊聚类分析方法 - Google Patents

Abstract

本发明公开一种基于非局部均值的SAR图像变化检测模糊聚类分析方法。实现过程主要包括：输入一幅由两幅不同时间相同地域的SAR图像的差异图；按全局快速模糊C均值聚类(FGFCM)算法中的相似性度量指标修正差异图像素值，得到局部空间信息像素值矩阵；对差异图作非局部均值处理生成非局部滤波的像素值矩阵；对上述两个矩阵加权求和生成完整的像素值矩阵；运用FGFCM算法对其进行聚类，生成变化检测二值结果图，整体完成对两幅SAR图像的变化检测。本发明兼顾了图像局部空间信息和非局部均值信息，并将其有机结合，使图像分析的聚类过程中既有效克服噪声影响并保留图像细节，得到更为精确的差异图分析结果。

Description

基于非局部均值的SAR图像变化检测模糊聚类分析方法

技术领域

本发明属于SAR图像变化检测技术领域，涉及SAR图像变化检测中的差异图分析技术。具体地说是提出了一种基于非局部均值的SAR图像变化检测模糊聚类分析方法，用来进行SAR图像变化检测中对差异图的分类分析，克服原有方法变化区域检测错误率较高的问题，提高SAR图像变化检测中的检测精度及速度。

背景技术

随着合成孔径雷达(synthetic aperture radar，SAR)技术的快速发展，SAR系统可以全天候、全天时获取图像数据，是较好的变化检测图像源。SAR图像变化检测是通过对同一地区的不同时期的两幅SAR图像进行比较分析，根据图像之间的差异得到所需的地物变化信息。SAR图像变化检测技术在森林覆盖变化、土地利用与覆盖、城市环境变化等领域具有广泛的应用前景。

SAR图像变化检测可以分为三个阶段：首先，对获得的SAR图像进行预处理，包括图像配准和辐射校正；然后，对校正过的图像进行比较，生成包含变化信息的差异图像；最后，通过分析差异图，得到最终的变化检测结果图像(即用二值图像代表变化和未变化两类)。差异图分析技术是SAR图像变化检测中的关键技术之一。常见的差异图分析理论可以分为两类，阈值理论和聚类理论。

在SAR图像变化检测中，许多不同的阈值分割方法被应用于差异图分析。KI阈值分割方法是Kittler和Illingworth提出的经典的最小误差阈值方法，该方法假设变化区域与非变化区域的直方图统计分布符合高斯分布，通过最小化惩罚函数来自动确定阈值。但该方法中基于高斯模型的假设并不准确，实际SAR图像数据分布并不符合高斯分布，所以其检测错误率仍较高。意大利G.Moser，S.B.Serpico等人在最小误差阈值基础上提出了一种广义最小误差阈值方法，该方法以三种概率分布模型：Nakagami分布、Log-Normal分布和Weibull分布为基础扩展了最小误差阈值方法。该阈值方法由于使用了更加适合SAR图像数据分布的模型，所以呈现了较好的检测结果。但是，这三种模型仍不能完全拟合差异图直方图，所以其检测错误率仍然较高。阈值理论方法简单，时间迅速但需要估计差异图变化类和未变化类的统计模型，这样大大限制了它的应用范围。聚类理论可以克服阈值理论的限制性，无需对变化类和未变化类的分布模型进行估计，可以提高分类精度，进而大大改善变化检测的正确率。

聚类理论是指将物理或抽象对象的集合分组成由类似的对象组成的多个类的过程。由聚类所生成的簇是一组数据对象的集合，这些对象与同一簇中的对象彼此相似，与其它簇中的对象相异。其中K均值聚类算法和模糊C均值聚类算法是理论发展最成熟和最为常用的聚类算法。K均值算法是典型的硬聚类算法，其主要缺陷是依赖初始聚类中心的选择，一旦初始值选的不好，可能无法得到有效地聚类结果。另外，该算法需要不断地进行样本分类调整，不断地计算调整后新的聚类中心，因此，K均值聚类算法的时间开销是非常大的。由于硬聚类把每个待识别的对象严格的划分到某类中，具有非此即彼的性质，而模糊聚类能够描述样本类属的中介性，能够客观地反映现实世界，模糊聚类已逐渐成为聚类分析的主流。在众多模糊聚类算法中，模糊C均值聚类(FCM)算法应用最为广泛，发展最为成熟，并有一系列基于模糊C均值聚类算法的改进算法，克服了原始FCM处理图像不考虑图像空间信息及容易陷入局部最优的缺陷，对原始FCM作出了改进。

全局快速模糊C均值聚类(FGFCM)算法在原始FCM中引入空间相似性度量指标，既考虑了邻域的灰度值信息，也考虑了邻域的空间信息，能综合图像中的空间信息得到更为精确的聚类结果，且聚类过程在灰度统计直方图上进行，图像的灰度值往往远小于图像大小，这样，在时间复杂度上FGFCM比原始FCM更有优势。

综上，现有技术中包括现有的差异图分析方法对差异图的质量依赖性高，且对差异图中噪声敏感，现有的常用差异图分析方法检测错误率较高，以致SAR图像变化检测系统精度不高。

发明内容

本发明的目的在于：运用改进模糊C均值聚类将SAR图像变化检测中的差异信息图分为变化/未变化两类，克服现有差异图分析方法对差异图中的噪声敏感，检测错误率较高的缺陷。

SAR图像变化检测可以分为三个阶段，预处理阶段，比较产生差异图阶段和差异图分析阶段，对实现变化检测而言，这三个阶段是不可或缺但又相对独立的部分，差异图分析作为SAR图像变化检测的终极阶段，其任务是分析由上一阶段生成的差异图像，即把上一阶段生成的差异信息图作为图像源，进行分类分析处理，得到能显示变化/未变化区域的二值图像，完成SAR图像变化检测，本发明正是处于SAR图像变化检测的终极阶段。

本发明的技术方案是：首先输入由两幅不同时间相同地域的SAR图像构造的差异图，然后按全局快速模糊C均值聚类(FGFCM)算法中的相似性度量指标修正像素值，得到考虑到局部空间信息的像素值矩阵，接着对差异图进行非局部均值滤波处理，得到非局部滤波后的像素值矩阵，再将局部空间信息矩阵与非局部信息矩阵加权求和，生成完整的像素值矩阵，最后运用FGFCM算法对其进行聚类，再通过FGFCM聚类结果生成变化检测结果图，完成对两幅SAR图像中变化区域的最终检测，本发明其具体实现步骤如下：

步骤1通过星载合成孔径雷达获取两幅不同时间相同地域的SAR图像，将两幅不同时间相同地域的图像，输入到安装有矩阵实验室和Visual C++6.0软件的计算机中，使SAR图像中每个像素点由其像素灰度值，即像素值表示；利用相关软件处理经过滤波去噪，辐射校正与几何配准的预处理SAR图像X₁，图像X₂构造的差异图像X，

使得在差异图像上低灰度级呈现为无变化区域，高灰度级呈现为变化区域；

步骤2通过Matlab或C++编程，或Matlab和C++混合编程实现对差异图X中每个像素点x_i按相似性度量方法进行处理，得到局部空间信息像素矩阵，差异图X中每个像素点按公式

计算，其中S_ij表示以像素点i为中心的邻域内第j个像素值和中心像素的相似度，x_j为邻域内第j个像素点的像素值，x′_i为修正后新的像素矩阵中第i个像素点的像素值；

步骤3对差异图X中每个像素点x_i作非局部均值处理，得到非局部滤波后像素矩阵，非局部滤波后像素矩阵中每个像素点通过公式求得，其中

是指以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，

是像素点i和在搜索窗口内像素点p的相似度权重，且满足0≤w_ip≤1和

为非局部滤波后的像素矩阵中第i个像素点的像素值；

步骤4对差异图X中每个像素点计算给予权重λ_i，

其中

A_i，k，A_p，k分别代表以像素点i和像素点p为中心的第k个像素点的像素值，U_t代表在搜索窗口内相似度指数U(x_i，x_p)由大到小排列的第t个，m表示要取相似度指数U(x_i，x_p)的个数；

步骤5对步骤2生成的局部空间信息像素矩阵和步骤3非局部均值滤波后的像素矩阵进行加权求和，得到新的完整像素值矩阵η，

其中η_i为像素值矩阵η中第i个像素点的像素值；

步骤6设置FGFCM算法中聚类个数c，模糊指数m，终止条件ε，最大迭代次数T，搜索窗口半径r，邻域窗口半径s；

步骤7对差异图完整像素值矩阵η按设定聚类个数随机初始化隶属度矩阵μ_ki，μ_ki表示第i个像素点属于第k类的隶属度，满足

和0≤μ_ki≤1，令迭代计数b＝1；

步骤8更新差异图完整像素值矩阵η聚类中心

其中，v_k代表第k类的聚类中心，M表示差异图中的像素值统计个数，γ_i表示像素值为i的像素点个数，且满足 ${\mathrm{\Σ}}_{i=1}^M{\mathrm{\γ}}_i=N;$

步骤9更新差异图完整像素值矩阵η模糊隶属度矩阵其中v_k代表第k类的聚类中心，m为模糊指数，c为聚类个数；

步骤10根据FGFCM聚类算法目标函数定义，计算目标函数

当

重复步骤8和步骤9，直到

成立，聚类结束；

步骤11根据聚类结果，命完整像素值矩阵η中像素值较大的类输出像素值为255，完整像素值矩阵η中像素值较小的类输出像素值为0，得到变化检测二值结果图，差异信息图分为变化/未变化两类，二值结果图表明实现了SAR图像变化检测的分类。

本发明将非局部均值思想引入到FGFCM算法中，通过设置权重，让像素点的局部空间信息与非局部均值有机结合，得到新的包含更多有效正确信息的像素值矩阵，聚类进行在新产生的像素值矩阵统计直方图上，在降低时间复杂度的同时得到更为精确的分析结果，提高SAR图像变化检测差异图分析正确率，得到更佳的变化检测结果。

本发明的实现还在于：步骤4中计算差异图X中每个像素点的给予权重λ_i，具体实现步骤包括有：

4.1计算索窗口内相似度指数U(x_i，x_p)，按公式

进行计算，其中A_i，k，A_p，k分别代表以像素点i和像素点p为中心的第k个像素点的像素值。

4.2计算给与权重，为对相似度指数U(x_i，x_p)的部分项求和，U_t代表在搜索窗口内相似度指数U(x_i，x_p)由大到小排列的第t个，m表示要取相似度指数U(x_i，x_p)的个数，根据大量实验，m＝2r时可得最佳效果，r为以像素点i为中心的搜索窗口半径。

由于相似度指数是评价搜索窗口内其余像素对中心像素的影响的自然指标，相似度指数大表示邻域块之间的非局部信息对中心像素影响较大，可以较多地使用非局部分量，反之，像素点邻域内的像素对本身影响更重。这样的权重表达式可以有效地提高去噪效率，降低误检率。本发明通过适当设置权重，让像素点的局部空间信息与非局部均值有机结合，可以得到新的包含更多有效正确信息的像素值矩阵，让聚类进行在新产生的像素值矩阵统计直方图上，在降低时间复杂度的同时得到更为精确的分割结果，可以明显提高SAR图像变化检测差异图分析正确率，得到更佳的变化检测结果。

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

1、本发明在全局快速模糊C均值聚类算法的基础上加入非局部均值信息，使得模糊聚类应用于斑点噪声很多的SAR图像变化检测差异图上，在降低时间复杂度的同时提高最终分割的精度，非局部均值思想的加入结合原有的局部空间修正可以最大程度的抑制乘性噪声，从而使得SAR图像变化检测得到更好的结果。通过与几种SAR图像变化检测聚类方法对比，本发明方法得到的SAR图像变化检测结果的性能最佳。

2、本发明在向全局快速模糊C均值聚类算法中引入非局部均值时，巧妙地利用了相似性较大的邻域块之间的非局部信息对中心像素影响较大，反之，像素点邻域内的像素对本身影响更重，得到的权重表达式可以有效地提高去噪效率，降低误检率。

3、本发明与其它经典的SAR图像变化检测分析方法相比，获得了较低的检测错误率，所获得的最终二值检测结果图最为贴合实际变化结果，验证了本发明的有效性。

附图说明

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

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

图3是本发明方法和三种对比方法对Bern地区SAR影像的实验结果图；

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

图5是本发明方法和三种对比方法对黄河入海口地区SAR影像的实验结果图。

具体实施方式

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

本发明是一种基于非局部均值的SAR图像变化检测模糊聚类分析方法，本发明在原有的考虑到图像空间信息的FGFCM算法上进一步引入了非局部均值思想，使非局部均值与局部信息有机结合，在抑制差异图噪声的同时极大程度地保持了图像的边缘和细节，使聚类进行在包含更多有效信息的新的像素值矩阵上，得到更为精确的分析结果，降低了整个SAR图像变化检测错误率，提高了SAR图像变化检测系统性能。本发明用于分析SAR图像变化检测差异信息图，得到最终二值结果图。本发明现阶段可以运行在32位XP系统(及以上)计算机，matlab7.0及以上平台，Visual C++6.0平台。

实施例1

本发明是一种基于非局部均值的SAR图像变化检测模糊聚类分析方法，参见图1，首先对两幅不同时间相同地域的SAR图像构造差异图，然后按全局快速模糊C均值聚类(FGFCM)算法中的相似性度量指标修正像素值，得到考虑到局部空间信息的像素值矩阵，接着对差异图进行非局部均值滤波处理，得到非局部滤波后的像素值矩阵，再将局部空间信息矩阵与非局部信息矩阵加权结合，生成完整的像素值矩阵，最后运用FGFCM算法对其进行聚类，再通过该聚类结果生成变化检测结果图，完成对两幅SAR图像中变化区域的最终检测，基于非局部均值的SAR图像变化检测模糊聚类分析具体实现步骤包括有：

步骤1通过星载合成孔径雷达获取两幅不同时间相同地域的SAR图像，将两幅不同时间相同地域的图像，输入到安装有矩阵实验室和Visual C++6.0软件的计算机中，利用相关软件进行处理，使SAR图像中每个像素点由其像素灰度值即像素值表示。相关软件主要指矩阵实验室和Visual C++6.0软件和常用计算机工具软件，本例中相关软件是利用Matlab7.0(及更高版本)和Visual C++6.0软件进行处理，处理过程包括滤波去噪，辐射校正与几何配准的预处理SAR图像X₁，图像X₂构造的差异图像X，每个像素点通过公式

处理，其中，x_1i和x_2i分别为图像X₁，图像X₂的第i个像素点的像素值，x_i为差异图像X的第i个像素点的像素值。

步骤2通过Matlab编程实现对差异图X中每个像素点x_i按相似性度量方法进行处理。

2.1遍历差异图X每个像素点，计算以像素点i为中心的邻域内第j个像素值和中心像素i的相似度S_ij， $S_{\mathrm{ij}}=\left\{\begin{array}{cc}\exp (-\max \left(\right|p_i-p_j|,|q_i-q_j\left|\right)/{\mathrm{\λ}}_s-{\left|\right|x_i-x_j\left|\right|}^2/{\mathrm{\λ}}_g{\mathrm{\σ}}_i^2),& i\≠j\\ 0& i=j\end{array}\right.,$ 其中(p_i，q_i)，(p_j，q_j)分别为像素点i和邻域内第j个像素点的空间坐标，x_i，x_j分别为像素点i和邻域内第j个像素点的像素值，λ_s和λ_g分别为S_ij的两个延展因子，其中N_i代表以像素点i为中心的邻域，N_R代表此邻域内除像素点i以外的像素点个数，在本发明运用于SAR图像变化检测中，取7×7邻域，即邻域半径为3。

2.2遍历差异图X每个像素点，即i遍历差异图中每个像素点进行相似性度量方法计算，

其中S_ij代表以像素点i为中心的邻域内第j个像素值和中心像素i的相似度，x_j为邻域内第j个像素点的像素值，x′_i为修正后，即考虑到局部空间信息的新的像素矩阵中第i个像素点的像素值，得到新的考虑到空间邻域信息的像素值矩阵。

步骤3对差异图X中每个像素点x_i作非局部均值处理。

3.1计算在搜索窗口内，以像素点p为中心的邻域块与以像素点i为中心的邻域块的相似度权重， $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用于控制指数函数的衰减，

是指以像素点i为中心的半径为r的搜索窗口，理论上讲，非局部均值要取遍图像内每个点的7×7邻域块，由于在图像较大的情况下，这样的时间复杂度太高，故通常只选取像素点附近较大的一块区域(即搜索窗口)进行非局部均值运算，在本发明中令r＝10，即在一个21×21的区域内进行非局部运算，A_i，k，A_p，k分别代表以像素点i和像素点p为中心的第k个像素点的像素值。

3.2遍历差异图X每个像素点，即i遍历差异图中每个像素点进行非局部均值计算，即其中

是指以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，

是像素点i和在搜索窗口内像素点p的相似度权重，且满足0≤w_ip≤1和

为非局部滤波后的像素矩阵中第i个像素点的像素值。

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

步骤4对差异图X中每个像素点计算给予权重λ_i。

4.1计算索窗口内相似度指数U(x_i，x_p)，按公式

进行计算，其中A_i，k，A_p，k分别代表以像素点i和像素点p为中心的第k个像素点的像素值。

4.2计算给与权重，为对相似度指数U(x_i，x_p)的部分项求和，U_t代表在搜索窗口内相似度指数U(x_i，x_p)由大到小排列的第t个，m表示要取相似度指数U(x_i，x_p)的个数，根据大量实验，m＝2r时可得最佳效果。

步骤5i遍历差异图中每个像素点，对考虑到局部空间信息的新的像素矩阵和非局部均值滤波后的像素矩阵进行加权求和，

得到新的完整像素值矩阵η，其中η_i为像素值矩阵η中第i个像素点的像素值。为了降低SAR图像变化检测的错误率克服现有方法较高的检测错误率，提高检测精度，针对SAR图像变化检测的特点，本发明将非局部均值思想引入到SAR图像变化检测聚类分析中，与其它现有的方法相比能够获得较低的检测错误率。

步骤6设置FGFCM算法中聚类个数c，模糊指数m，终止条件ε，最大迭代次数T，搜索窗口半径r，邻域窗口半径s，本发明用于分析SAR图像变化检测差异图，只需把图像分为变化和未变化两部分，因此聚类个数c＝2。令m＝2，通常模糊指数均取值为2，ε＝0.00001，T＝200，r＝10，s＝3。

步骤7对完整像素值矩阵η按设定聚类个数随机初始化隶属度矩阵μ_ki，表示第i个像素点属于第k类的隶属度，满足

和0≤μ_ki≤1，令迭代计数b＝1。

步骤8更新完整像素值矩阵η聚类中心i遍历差异图中每个像素点，其中，v_k代表第k类的聚类中心，M表示差异图中的像素值统计个数，γ_i表示像素值为i的像素点个数，且满足

步骤9更新完整像素值矩阵η模糊隶属度矩阵

i遍历差异图中每个像素点，其中v_k代表第k类的聚类中心，m＝2，c＝2。

步骤10计算目标函数

其中M表示差异图中的像素值统计个数，γ_i表示像素值为i的像素点个数，v_k代表第k类的聚类中心，μ_ki表示第i个像素点属于第k类的隶属度m＝2，c＝2，当

重复第(8)，(9)步，直到 $|J_m^{b+1}-J_m^b|\≤\mathrm{\ϵ},$ 聚类结束。

步骤11根据聚类结果，命像素值较大的类输出像素值为255，像素值较小的类输出像素值为0，得到变化检测二值结果图。完整像素值矩阵η中像素值较小的类输出像素值为0，差异信息图分为变化/未变化两类，二值结果图表明实现了SAR图像变化检测的分类。

本发明在全局快速模糊C均值聚类算法的基础上加入非局部均值信息，使得模糊聚类应用于斑点噪声很多的SAR图像变化检测差异图上，在降低时间复杂度的同时提高最终分割的精度，非局部均值思想的加入结合原有的局部空间修正可以最大程度的抑制乘性噪声，从而使得SAR图像变化检测得到更好的结果。通过与几种SAR图像变化检测聚类方法对比，本发明方法得到的SAR图像变化检测结果的性能最佳。

实施例2

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

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

步骤2对差异图X中每个像素点x_i按相似性度量方法进行运算处理，详细实施步骤如实施例2中的步骤2。

步骤3对差异图X中每个像素点x_i作非局部均值处理，详细实施步骤如实施例1中的步骤3。

步骤4对差异图X中每个像素点计算给予权重λ_i，详细实施步骤如实施例1中的步骤4。

步骤5对考虑到局部空间信息的新的像素矩阵和非局部均值滤波后的像素矩阵进行加权求和，详细实施步骤如实施例1中的步骤5。

步骤6设置FGFCM算法中聚类个数c，模糊指数m，终止条件ε，最大迭代次数T，搜索窗口半径r，邻域窗口半径s，详细实施步骤如实施例1中的步骤6。

步骤7对像素值矩阵η按设定聚类个数随机初始化隶属度矩阵μ_ki，详细实施步骤如实施例1中的步骤7。

步骤8更新聚类中心v_k，详细实施步骤如实施例1中的步骤8。

步骤9更新模糊隶属度矩阵μ_ki，详细实施步骤如实施例1中的步骤9。

步骤10计算目标函数，并判定是否满足终止条件，详细实施步骤如实施例1中的步骤10。

步骤11根据聚类结果，得到变化检测二值结果图，如图3(d)，详细实施步骤如实施例1中的步骤11。

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

实施例3

基于非局部均值的SAR图像变化检测模糊聚类分析方法同实施例1-2，

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

1、仿真参数

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

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

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

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

④正确率：

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

为了验证基于非局部均值的SAR图像变化检测模糊聚类分析方法的优越性，将本发明方法与在差异信息图分析阶段比较常见的模糊C均值聚类(FCM)方法和基于邻域信息的模糊C均值聚类(FLICM)方法及全局快速模糊C均值聚类(FGFCM)方法性能做出比较。对带参考图的真实SAR图像数据Bern地区进行了实验。这里把模糊C均值聚类方法简称为FCM，基于邻域信息的模糊C均值聚类方法简称为FLICM，全局快速模糊C均值聚类方法简称为FGFCM，本发明方法简称为NLFCM.

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

对在真实Bern区域SAR图像数据上产生的差异信息图X进行差异图分析的实验结果图如图3所示。

在图3中图3(a)代表模糊C均值聚类(FCM)方法分析结果图像，3(b)代表基于邻域信息的模糊C均值聚类(FLICM)方法分析结果图像，3(c)代表全局快速模糊C均值聚类(FGFCM)方法分析结果图像，3(d)代表本发明(NLFCM)方法分析差异图的二值结果图。

从图3中可见，图3(a)中含有较多噪声点，由于原始的FCM方法没有考虑图像的空间信息，仅对单个像素点作聚类处理，不能抑制差异图中的噪声点，使得结果总错误数较高，图3(b)和图3(c)的分析方法涉及到空间信息，在杂点抑制方面有较好表现，但是由于只考虑局部信息，使得图像边缘信息弱化，失去了原有图像的边缘细节信息，使得误检率较高，图3(d)中在保持边缘细节的同时很好的抑制了噪声点，分析效果明显优于图3(a)，图3(b)和图3(c)。

差异信息图用四种分析方法作分析性能指标如下表所示：

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

从表1中可以看出，本发明在变化检测总错误数上，相比其他对比算法是最少的，而且在漏检数上，本发明与原始FCM相比有明显提升，这是由于本发明在FGFCM的基础上利用空间信息的指引作用，克服噪声的同时也能保持图像细节，在误检数方面，本发明(NLFCM)相对RFLICM和FGFCM有显著提高，这是因为本发明充分结合非局部均值，保留了差异信息图中的细节，边缘表达更佳；这样在误检数与漏检数上本发明相对都比较均衡，数量上也处于中间水平，总错误数最少，说明了本发明方法的优越性。总体来看，通过表1性能结果和图3的视觉效果均说明了本发明方法对差异信息图分析更为适用，在聚类的过程中抑制了差异信息图中的噪声和最大程度上保留了图像细节，最终提高了SAR图像变化检测的精确度，使得总错误数更少，误检数与漏检数上相对都比较均衡。

实施例4

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

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

步骤2对差异图X中每个像素点x_i按相似性度量方法进行运算处理，详细实施步骤如实施例1中的步骤2。

步骤3对差异图X中每个像素点x_i作非局部均值处理，详细实施步骤如实施例1中的步骤3。

步骤4对差异图X中每个像素点计算给予权重λ_i，详细实施步骤如实施例1中的步骤4。

步骤5对考虑到局部空间信息的新的像素矩阵和非局部均值滤波后的像素矩阵进行加权求和，详细实施步骤如实施例1中的步骤5。

步骤6设置FGFCM算法中聚类个数c，模糊指数m，终止条件ε，最大迭代次数T，搜索窗口半径r，邻域窗口半径s，详细实施步骤如实施例1中的步骤6。

步骤7对像素值矩阵η按设定聚类个数随机初始化隶属度矩阵μ_ki，详细实施步骤如实施例1中的步骤7。

步骤8更新聚类中心v_k，详细实施步骤如实施例1中的步骤8。

步骤9更新模糊隶属度矩阵μ_ki，详细实施步骤如实施例1中的步骤9。

步骤10计算目标函数，并判定是否满足终止条件，详细实施步骤如实施例1中的步骤10。

步骤11根据聚类结果，得到变化检测二值结果图，如图5(d)，详细实施步骤如实施例1中的步骤11。

实施例5

基于非局部均值的SAR图像变化检测模糊聚类分析方法同实施例1-4，

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

1、仿真参数同实施例3

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

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

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

各种差异信息图分析方法在真实黄河入海口区域SAR图像数据的实验结果图如图5所示，其中图5(a)代表模糊C均值聚类(FCM)方法分析结果图像，图5(b)代表基于邻域信息的模糊C均值聚类(FLICM)方法分析结果图像，图5(c)代表全局快速模糊C均值聚类(FGFCM)方法分析结果图像，图5(d)代表本发明(NLFCM)方法分析差异图的二值结果图。

图5(a)-5(d)反映出受噪声污染比较严重的SAR图像生成的含噪差异信息图的分析情况，从视觉上分辨可见，图5(a)中噪声杂点极其严重，变化区域边缘有所保持，但出现大量成片的不规则噪声点，图5(b)中变化区域的边缘较为规整，对未变区域的噪声抑制不明显，图5(c)中变化区域边缘模糊，未变区域噪声较多，图5(d)为本发明分析差异结果图得到的二值结果图，由于NLFCM在FGFCM的基础上加入了非局部信息，可以利用大量图像信息修正像素点，抑制噪声影响，并尽可能地保持图像细节，然后使聚类进行在修正过的新的像素矩阵上，图5(d)效果明显优于图5(a)，图5(b)和图5(c)。

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

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

从表2中可以明显看出，本发明在变化检测总错误数上，相比其他对比算法是最少的，说明了本发明方法的优越性。

从附图5中看到，原始FCM聚类方法产生结果图像噪声很多，这是由于原始FCM聚类方法没有考虑到任何图像空间信息，对图像噪声完全没有抑制作用，从表2中的误检数也可以给出同样的结论，FLICM聚类方法对变化区域的边缘保留比较突出，是由于聚类过程中考虑到图像的邻域空间信息但并未对图像进行修正，而对于这种受噪声影响很大的图像，空间邻域信息的指导作用是很有限的，而FGFCM聚类方法产生的结果边缘被弱化，噪声也未得到有效控制，本发明的NLFCM聚类方法是总错误率最小的。总体来看，通过表2性能结果和参考图5的视觉效果说明了本发明方法在聚类的过程中抑制了差异信息图中的噪声和最大程度上保留了图像细节，最终提高了SAR图像变化检测的精确度，使得总错误数更少，误检数与漏检数上相对都比较均衡。

综上所述本发明方法能够得到更好的SAR图像变化检测结果.本发明与原始FCM和FGFCM及FLICM相比，兼顾了图像局部空间信息和非局部均值，有机结合以上两者，使聚类过程中既可以有效克服噪声影响并保留图像细节，得到更为精确的差异图分析结果。

通过两个实验的分析可以看出，本发明方法的结果与其他方法相比更能抑制噪声，保持图像细节，尤其在处理本身含噪情况严重的图像，优势更为明显，从而得到更佳的SAR图像变化检测的结果。

综上，本发明的基于非局部均值的SAR图像变化检测模糊聚类分析方法，其实现过程主要包括：首先对两幅不同时间相同地域的SAR图像构造差异图，然后按全局快速模糊C均值聚类(FGFCM)算法中的相似性度量指标修正像素值，得到考虑到局部空间信息的像素值矩阵，接着对差异图进行非局部均值滤波处理，得到非局部滤波后的像素值矩阵，再将局部空间信息矩阵与非局部信息矩阵加权求和，生成完整的像素值矩阵，最后运用FGFCM算法对其进行聚类，再通过该聚类结果生成变化检测结果图，完成对两幅SAR图像中变化区域的最终检测。本发明在原有的考虑到图像空间信息的FGFCM算法上进一步引入了非局部均值思想，使非局部均值与局部信息有机结合，在抑制差异图噪声的同时极大程度地保持了图像的边缘和细节，使聚类进行在包含更多有效信息的新的像素值矩阵上，得到更为精确的分析结果，降低了整个SAR图像变化检测错误率，提高了SAR图像变化检测系统性能。

Claims (2)

1.一种基于非局部均值的SAR图像变化检测模糊聚类分析方法，其特征在于：基于非局部均值的SAR图像变化检测模糊聚类分析具体实现步骤包括有：

步骤1通过星载合成孔径雷达获取两幅不同时间相同地域的SAR图像，将两幅不同时间相同地域的图像，输入到安装有矩阵实验室和Visual C++6.0软件的计算机中，使SAR图像中每个像素点由其像素灰度值即像素值表示；利用相关软件处理经过滤波去噪，辐射校正与几何配准的预处理SAR图像X₁，图像X₂构造的差异图像X，使得在差异图像上低灰度级呈现为无变化区域，高灰度级呈现为变化区域；

步骤2通过Matlab或C++编程，或Matlab和C ++混合编程实现对差异图X中每个像素点x_i按相似性度量方法进行处理，得到局部空间信息像素矩阵，差异图X中每个像素点按公式

计算，其中S_ij表示以像素点i为中心的邻域内第j个像素值和中心像素的相似度，x_j为邻域内第j个像素点的像素值，x′_i为修正后新的像素矩阵中第i个像素点的像素值；

步骤3对差异图X中每个像素点x_i作非局部均值处理，得到非局部滤波后像素矩阵，非局部滤波后像素矩阵中每个像素点通过公式

求得，其中

是指以像素点i为中心的半径为r的搜索窗口，x_p是像素点p的像素值，

是像素点i和在搜索窗口内像素点p的相似度权重，且满足0≤w_ip≤1和

为非局部滤波后的像素矩阵中第i个像素点的像素值；

步骤4对差异图X中每个像素点计算给予权重λi，

其中

A_i，k，A_p，k分别代表以像素点i和像素点p为中心的第k个像素点的像素值，U_t代表在搜索窗口内相似度指数U(x_i，x_p)由大到小排列的第t个，m表示要取相似度指数U(x_i，x_p)的个数；

步骤5对步骤2生成的局部空间信息像素矩阵和步骤3非局部均值滤波后的像素矩阵进行加权求和，得到新的完整像素值矩阵η，

其中η_i为像素值矩阵η中第i个像素点的像素值；

步骤6设置FGFCM算法中聚类个数c，模糊指数m，终止条件ε，最大迭代次数T，搜索窗口半径r，邻域窗口半径s；

步骤7对差异图完整像素值矩阵η按设定聚类个数随机初始化隶属度矩阵μ_ki，μ_ki表示第i个像素点属于第k类的隶属度，满足

和0≤μ_ki≤1，令迭代计数b＝1；

步骤8更新差异图完整像素值矩阵η聚类中心

其中，v_k代表第k类的聚类中心，M表示差异图中的像素值统计个数，γ_i表示像素值为i的像素点个数，且满足 ${\mathrm{\Σ}}_{i=1}^M{\mathrm{\γ}}_i=N;$

步骤9更新差异图完整像素值矩阵η模糊隶属度矩阵

其中v_k代表第k类的聚类中心，m为模糊指数，c为聚类个数；

步骤10根据FGFCM聚类算法目标函数定义，计算目标函数

当

重复步骤8和步骤9，直到

成立，聚类结束；

步骤11根据聚类结果，命聚类后完整像素值矩阵η中像素值较大的类输出像素值为255，完整像素值矩阵η中像素值较小的类输出像素值为0，得到变化检测二值结果图，差异信息图分为变化/未变化两类，二值结果图表明实现了SAR图像变化检测的分类。

2.根据权利要求1所述的基于非局部思想的SAR图像变化检测模糊聚类分析方法，其特征在于：步骤4中计算差异图X中每个像素点的给予权重λ_i，具体实现步骤包括有：

4.1计算搜索窗口内相似度指数U(x_i，x_p)，按公式

进行计算，其中A_i，k，A_p，k分别代表以像素点i和像素点p为中心的第k个像素点的像素值；

4.2计算给与权重，对相似度指数U(x_i，x_p)的部分项求和，U_t代表在搜索窗口内相似度指数U(x_i，x_p)由大到小排列的第t个，m表示要取相似度指数U(x_i，x_p)的个数，m取值范围为1.5r≤m≤2r，r为以像素点i为中心的搜索窗口半径。

Images