CN104200472A - 基于非局部小波信息的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开一种基于非局部小波信息的遥感图像变化检测方法，主要解决遥感图像的变化检测方法存在的检测精度不高的问题。其实现步骤为：(1)读入数据；(2)构造差异图；(3)小波分解；(4)对高频部分应用基于非局部信息的方法进行去噪；(5)小波逆变换；(6)使用基于局部模糊C均值聚类的方法进行分割。本发明考虑到噪声主要分布在高频细节信息部分，将图像的去噪过程在高频细节部分进行，这既保护了低频信息的完整性，又能有效地去除噪声。本发明将单个像素变为结合邻域信息的向量进行处理，保护了图像本身的结构，提高了图像对噪声的鲁棒性。

Description

基于非局部小波信息的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，更进一步涉及一种基于非局部小波信息的遥感图像变化检测方法。该方法可应用于环境变化中的湖泊水位动态监测、农作物生长状态的动态监测、军事侦察等领域，能够快速的检测出两时相遥感图像的变化信息。

背景技术

遥感图像的变化检测是指通过分析在不同时间来自同一地区的两幅或多幅遥感图像，检测出该地区的地物随时间发生的变化信息。遥感图像的变化检测已经广泛地应用于如森林资源的动态监测、土地覆盖和利用的变化监测、农业资源调查、城市规划布局、环境监测分析、自然灾害评估、地理数据更新以及军事侦察中战略目标(如道路、桥梁、机场)等的动态监视等许多领域。

SAR图像变化检测中，主要分为三个步骤：1.差异图生成，2.对差异图进行分割，3.分析变化检测结果图。其中，差异图的生成在变化检测中起着至关重要的作用，它的质量好坏直接影响变化检测结果的优劣。目前，应用较为广泛的基本的差异图生成算法主要有差值法，比值法，对数比法等。早期的研究已经证实，差值操作在高灰度级区域会比低灰度级区域产生的误差更多，这一点不利于构造差异影像直方图的统计模型。另外，比值操作对于乘性噪声是非常敏感的，对数比值操作能够降低乘性斑点噪声的影响，因此得到了学者的广泛应用。

西安电子科技大学在其专利申请“基于图像融合的光学遥感图像变化检测”(专利申请号：201210234076.1，公开号：CN 102750705A)中提出了一种基于图像融合的光学遥感图像变化检测方法。其实现过程为：对同一地域不同时间获取的两幅光学遥感图像预处理后，分别构造差值法差异图和比值法差异图；对这两幅差异图分别进行N层小波分解后获取各分解层高频带和低频带的小波系数；对高频带和低频带的小波系数分别用不同的融合算子融合处理，得到高频带和低频带的融合小波系数；对高频带和低频带的融合小波系数进行逆变换得到融合后差异图；运用模糊局部C均值聚类法分割融合后的差异图，得到变化检测结果。该方法虽然能都解决现有变化检测技术检测结果精度不高的问题，但仍存在的不足是，该方法在图像受到噪声干扰时，图像的变化检测结果迅速变差，抗噪性能差。

西北工业大学在其专利申请“一种基于边缘和灰度的遥感图像变化检测方法”(专利申请号：201010292888.2，公开号：CN 101968885A)中提出了一种基于边缘和灰度的遥感图像变化检测方法。该技术方案是利用基于双边滤波的Canny算法进行多时相图像边缘特征提取，然后对灰度差值图像进行OSTU阈值分割和边缘提取，获得灰度特征。再将所提取的边缘和灰度特征进行综合，检测遥感图像的变化区域。该方法由于充分利用了图像中的线性特性的同时，利用灰度差值图像弥补了因配准误差造成的断线，提高了检测准确性，但仍存在的不足是，不能有效地衡量去除噪声和保留图像细节信息之间的关系，不能有效地保留图像的细节信息。

发明内容

本发明针对上述现有技术存在的不足，提出了一种基于非局部小波变换遥感图像变化检测方法。本发明既可以较好的保留图像的细节信息，又对噪声具有很好的鲁棒性，能够在去除噪声的同时保留细节信息。

本发明实现上述目的的思路是：在对读入的数据进行对数比操作后，先对构造的差异图像进行小波分解，再对分解后图像高频细节部分进行基于非局部信息的去噪处理，接着进行小波逆变换，最后使用基于局部模糊C均值聚类算法进行分割，得到最终的变化检测结果图。

本发明的步骤包括如下：

(1)读入数据：

读入同一地区不同时刻获取的已配准的两幅遥感图像；

(2)构造差异图：

对读入的两幅遥感图像，使用对数比公式，获得两幅遥感图像的差异图；

(3)小波分解：

对两幅遥感图像的差异图，使用离散小波变换DWT进行分解，得到一幅低频轮廓图像和三幅高频细节图像，三幅高频细节图像分别为水平方向，竖直方向和对角方向的高频细节图像；

(4)去噪：

(4a)从水平方向的高频细节图像中任选一点k，以k为中心，选取一个大小为l×l的邻域窗口，其中，l取值为(3，5，7)奇数；

(4b)将邻域窗口的全部像素点的灰度值以列优先的方式排成一列，得到该邻域窗口的中心像素点k的特征向量；

(4c)重复步骤(4a)、步骤(4b)，直至得到水平方向高频细节图像中的全部像素点的特征向量；

(4d)按照下式，计算水平方向高频细节图像去噪时的权重值：

$w_{\mathrm{ij}}=e^{-\frac{{||R_i-R_j||}^2}{{\mathrm{\σ}}^2}}$

其中，w_ij表示水平方向高频细节图像中的点j对点i去噪时的影响权重值，e表示自然常数，R_i表示水平方向高频细节图像中的点i的特征向量，R_j表示水平方向高频细节图像中的点j的特征向量，σ表示水平方向高频细节图像的方差；

(4e)按照下式，对水平方向高频细节图像去噪，得到水平方向高频细节图像中的像素点的真实灰度值：

$t_i=\frac{1}{W_i}\underset{j\∈s_i}{\mathrm{\Σ}}{w}_{\mathrm{ij}}{x}_j$

其中，t_i表示水平方向高频细节图像中的点i的真实灰度值，W_i表示水平方向高频细节图像中的点i的去噪权重值的归一化常量，w_ij表示水平方向高频细节图像中的点j对点i去噪时的影响权重值，x_j表示水平方向高频细节图像中的点j的像素灰度值，S_i表示水平方向高频细节图像中的点i的搜索窗口；

(4f)重复步骤(4d)、步骤(4e)，直至得到水平方向高频细节图像中的全部像素点的真实灰度值；

(4g)采用步骤(4a)、步骤(4b)、步骤(4c)、步骤(4d)、步骤(4e)、步骤(4f)的相同方法，对步骤(3)中得到的竖直方向高频细节图像和对角方向高频细节图像分别进行去噪处理，得到去噪后竖直方向和对角方向的高频细节图像；

(5)小波逆变换：

对水平、竖直和对角方向形成的三幅去噪后高频细节图像和步骤(3)中得到的低频轮廓图像共四幅图像，使用小波逆变换IDWT，得到去噪后的两幅遥感图像的差异图；

(6)分割差异图：

使用基于局部模糊C均值聚类算法FLICM，对去噪后的两幅遥感图像的差异图进行分割，得到最终的变化检测结果图。

本发明与现有技术相比具有以下优点：

第一，本发明利用基于非局部信息的去噪处理方法，克服了现有技术在图像受到噪声干扰时，图像的变化检测结果迅速变差，抗噪性能差的问题，使得本发明对噪声的鲁棒性得到了提高。

第二，本发明通过对对数比差异图进行小波分解，仅在高频细节部分进行基于非局部信息的去噪处理，保留低频轮廓部分信息，克服了现有技术不能有效地衡量去除噪声和保留图像细节信息之间的关系和不能有效地保留图像的细节信息的问题，使得本发明的变化检测精度得到了提高。

附图说明

图1为本发明的流程图；

图2为本发明的第一组两时相遥感数据图像及其变化检测参考图；

图3为本发明与现有对比实验方法的差异图；

图4为本发明与现有对比实验方法的变化检测结果图；

图5为本发明的第二组两时相遥感数据图像及其变化检测参考图；

图6为本发明与现有对比实验方法的差异图；

图7为本发明与现有对比实验方法的变化检测结果图；

具体实施方式

下面结合附图对本发明做进一步的详细描述。

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

步骤1，读入数据。

读入中同一地区不同时刻获取的已经过配准和几何校正的两幅遥感图像。

本发明实施例中输入第一组的原始遥感图像如图2所示，图2(a)、图2(b)源自欧洲遥感图像库中的Bern城市的不同时刻的两幅原始的遥感图像，分别拍摄于1999.04和1999.05，其大小均为256×256个像素。图2(c)为Bern城市变化检测参考图。

本发明实施例中输入第二组的原始遥感图像如图5所示，图5(a)、图5(b)源自加拿大遥感图像库中的Ottawa地区的不同时刻的两幅原始的遥感图像，分别拍摄于1997.05和1997.08，其大小均为290×350个像素。图5(c)为Ottawa地区变化检测参考图。

步骤2，构造差异图。

对读入的两幅遥感图像，使用对数比公式，获得两幅遥感图像的差异图，对数比公式为：

$X=\left|\log \frac{X_2}{X_1}\right|$

其中，X表示两幅遥感图像的差异图，X₁和X₂分别表示同一地区不同时刻已配准的两幅遥感图像，log表示自然对数操作。

参照图3，图3(a)是对图2(a)和图2(b)使用对数比公式得到的两幅遥感图像的差异图。

参照图6，图6(a)是对图5(a)和图5(b)使用对数比公式得到的两幅遥感图像的差异图。

步骤3，小波分解。

对两幅遥感图像的差异图，使用离散小波变换DWT进行分解，得到一幅低频轮廓图像和三幅高频细节图像，三幅高频细节图像分别为水平方向，竖直方向和对角方向的高频细节图像。

步骤4，去噪。

(4a)从水平方向的高频细节图像中任选一点k，以k为中心，选取一个大小为l×l的邻域窗口，其中，l取值为(3，5，7)奇数。

(4b)将邻域窗口的全部像素点的灰度值以列优先的方式排成一列，得到该邻域窗口的中心像素点k的特征向量。

(4c)重复步骤(4a)、步骤(4b)，直至得到水平方向高频细节图像中的全部像素点的特征向量。

(4d)按照下式，计算水平方向高频细节图像去噪时的权重值：

$w_{\mathrm{ij}}=e^{-\frac{{||R_i-R_j||}^2}{{\mathrm{\σ}}^2}}$

其中，w_ij表示水平方向高频细节图像中的点j对点i去噪时的影响权重值，e表示自然常数，R_i表示水平方向高频细节图像中的点i的特征向量，R_j表示水平方向高频细节图像中的点j的特征向量，σ表示水平方向高频细节图像的方差。

(4e)按照下式，对水平方向高频细节图像去噪，得到水平方向高频细节图像中的像素点的真实灰度值。

$t_i=\frac{1}{W_i}\underset{j\∈s_i}{\mathrm{\Σ}}{w}_{\mathrm{ij}}{x}_j$

其中，t_i表示水平方向高频细节图像中的点i的真实灰度值，W_i表示水平方向高频细节图像中的点i的去噪权重值的归一化常量，w_ij表示水平方向高频细节图像中的点j对点i去噪时的影响权重值，x_j表示水平方向高频细节图像中的点j的像素灰度值，S_i表示水平方向高频细节图像中的点i的搜索窗口。

(4f)重复步骤(4d)、步骤(4e)，直至得到水平方向高频细节图像中的全部像素点的真实灰度值。

(4g)采用步骤(4a)、步骤(4b)、步骤(4c)、步骤(4d)、步骤(4e)、步骤(4f)的相同方法，对步骤(3)中得到的竖直方向高频细节图像和对角方向高频细节图像分别进行去噪处理，得到去噪后竖直方向和对角方向的高频细节图像。

步骤5，小波逆变换。

对水平、竖直和对角方向形成的三幅去噪后高频细节图像和步骤(3)中得到的低频轮廓图像共四幅图像，使用小波逆变换IDWT，得到去噪后的两幅遥感图像的差异图。

参照图3，图3(b)是对图3(a)去噪后的两幅遥感图像的差异图。

参照图6，图6(b)是对图6(a)去噪后的两幅遥感图像的差异图。

步骤6，分割差异图。

使用基于局部模糊C均值聚类算法FLICM，对去噪后的两幅遥感图像的差异图进行分割，得到最终的变化检测结果图，基于局部模糊C均值聚类算法FLICM的步骤如下：

(6a)初始化：

设定聚类数目c＝2，模糊指数m＝2和迭代截止条件阈值ε＝10^-6，随机初始化模糊隶属度矩阵U⁽⁰⁾，设置循环计数器b＝0，按照下式，进行初迭代截止条件阈值ε＝10^-6始化：

U^(b)＝[u_ij](i＝1,2；j＝1,2,…,N)

其中，U^(b)表示模糊隶属度矩阵，b表示循环计数器，u_ij表示两幅遥感图像的差异图中像素点j属于第i类的模糊隶属度，N表示两幅遥感图像的差异图中像素点的总个数。

(6b)利用下式计算聚类中心：

$v_k=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ki}}^m{x}_i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ki}}^m}$

其中，ν_k表示第k类的聚类中心，μ_ki表示两幅遥感图像的差异图中像素点i属于第k类的模糊隶属度，x_i表示两幅遥感图像的差异图中像素点i的灰度值，m表示模糊指数，N表示两幅遥感图像的差异图中像素点的总个数；

(6c)按照下式，计算空间距离项：

$G_{\mathrm{ki}}=\underset{j\∈N_i}{\mathrm{\Σ}}\frac{1}{d_{\mathrm{ij}}+1}{(1-{\mathrm{\μ}}_{\mathrm{kj}})}^m{\left|\right|x_j-v_k\left|\right|}^2$

其中，G_ki表示两幅遥感图像的差异图中像素点i到第k类聚类中心的空间距离，N_i表示点两幅遥感图像的差异图中像素点i的空间邻域，d_ij表示点i的空间邻域中的点j到点i的欧式距离，d_ij＝||x_j-x_i||²，x_j表示两幅遥感图像的差异图中像素点i的空间邻域中的点j的灰度值，ν_k表示第k类的聚类中心，m表示模糊指数；

(6d)按照下式，计算模糊隶属度矩阵的元素，得到更新后的模糊隶属度矩阵U^(b+1)：

${\mathrm{\μ}}_{\mathrm{ki}}=\frac{1}{\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{\left(\frac{{\left|\right|x_i-v_k\left|\right|}^2+G_{\mathrm{ki}}}{{\left|\right|x_i-v_j\left|\right|}^2+G_{\mathrm{ji}}}\right)}^{1/(m-1)}}$

其中，μ_ki表示两幅遥感图像的差异图中像素点i属于第k类的模糊隶属度，G_ki表示像素点i到第k个聚类中心的空间距离，x_i表示两幅遥感图像的差异图中像素点i的灰度值，ν_k表示第k类的聚类中心，m表示模糊指数，c表示聚类数目；

(6e)判断max{U^(b)-U^(b+1)}＜ε是否成立，如果成立，结束该聚类算法，得到最终的模糊隶属度矩阵U＝U^(b+1)，由模糊隶属度矩阵U得到各像素点的分类结果；否则，令b＝b+1,执行(6b)；其中，ε表示迭代截止条件阈值，ε＝10^-6，U^(b)表示更新前的模糊隶属度矩阵，U^(b+1)表示更新后的模糊隶属度矩阵。

参照图4，图4(a)、图4(b)是对图3(a)、图3(b)分别进行分割后得到的变化检测结果图。

参照图7，图7(a)、图7(b)是对图6(a)、图6(b)分别进行分割后得到的变化检测结果图。

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

1.仿真条件：

本发明的仿真是在Intel(R)Core(TM)i3 CPU 3.19GHz 3.05G内存WindowsXP系统下，Matlab r2012a运行平台上，完成的仿真实验。

2.仿真实验结果及分析：

从图4和图7这两组变化检测结果对比可以发现，本发明在视觉上均能得到较好的结果，为了进一步比较客观的说明问题，分别比较了变化检测图像的正确率及kappa系数。对于第一组实验对数比方法得到的变化监测结果图4(a)的正确率为99.49％，本发明方法得到的变化检测结果图4(b)的正确率为99.57％，图4(a)的kapaa系数为0.8331，本发明方法图4(b)的kapaa系数为0.8653，本发明方法的漏检数由306减少到225；第二组实验变化检测结果如图7所示，图7(a)的正确率和kappa系数分别为0.9733和0.8935，本发明图7(b)所示的结果正确率为0.9796，kappa系数为0.9196。通过这些指标可以定量的看出本发明应用于SAR图像变化检测产生了较好的效果。

由图4和图7的仿真结果可以说明，针对SAR图像变化检测，本发明存在一定的优势，将图像使用小波变换只对高频部分进行处理，具有针对性，能够有效的去除噪声而不影响图像的低频信息；在去除噪声的过程中，使用图像的邻域信息构成的特征向量替代单个像素，用全局像素点的特征向量与所求像素点特征向量灰度值的相似性确定出全局像素对于所求像素点灰度值的贡献，提高了对噪声的鲁棒性，保护了图像本身的结构信息。

综上所述，本发明能够更加有效的应用于SAR图像变化检测中。

Claims (5)

1.一种基于非局部小波信息的遥感图像变化检测方法，包括如下步骤：

(1)读入数据：

读入同一地区不同时刻获取的已配准的两幅遥感图像；

(2)构造差异图：

对读入的两幅遥感图像，使用对数比公式，获得两幅遥感图像的差异图；

(3)小波分解：

对两幅遥感图像的差异图，使用离散小波变换DWT进行分解，得到一幅低频轮廓图像和三幅高频细节图像，三幅高频细节图像分别为水平方向，竖直方向和对角方向的高频细节图像；

(4)去噪：

(4a)从水平方向的高频细节图像中任选一点k，以k为中心，选取一个大小为l×l的邻域窗口，其中，l取值为(3，5，7)奇数；

(4b)将邻域窗口的全部像素点的灰度值以列优先的方式排成一列，得到该邻域窗口的中心像素点k的特征向量；

(4c)重复步骤(4a)、步骤(4b)，直至得到水平方向高频细节图像中的全部像素点的特征向量；

(4d)按照下式，计算水平方向高频细节图像去噪时的权重值：

$w_{\mathrm{ij}}=e^{-\frac{{||R_i-R_j||}^2}{{\mathrm{\σ}}^2}}$

其中，w_ij表示水平方向高频细节图像中的点j对点i去噪时的影响权重值，e表示自然常数，R_i表示水平方向高频细节图像中的点i的特征向量，R_j表示水平方向高频细节图像中的点j的特征向量，σ表示水平方向高频细节图像的方差；

(4e)按照下式，对水平方向高频细节图像去噪，得到水平方向高频细节图像中的像素点的真实灰度值：

$t_i=\frac{1}{W_i}\underset{j\∈s_i}{\mathrm{\Σ}}{w}_{\mathrm{ij}}{x}_j$

其中，t_i表示水平方向高频细节图像中的点i的真实灰度值，W_i表示水平方向高频细节图像中的点i的去噪权重值的归一化常量，w_ij表示水平方向高频细节图像中的点j对点i去噪时的影响权重值，x_j表示水平方向高频细节图像中的点j的像素灰度值，S_i表示水平方向高频细节图像中的点i的搜索窗口；

(4f)重复步骤(4d)、步骤(4e)，直至得到水平方向高频细节图像中的全部像素点的真实灰度值；

(4g)采用步骤(4a)、步骤(4b)、步骤(4c)、步骤(4d)、步骤(4e)、步骤(4f)的相同方法，对步骤(3)中得到的竖直方向高频细节图像和对角方向高频细节图像分别进行去噪处理，得到去噪后竖直方向和对角方向的高频细节图像；

(5)小波逆变换：

对水平、竖直和对角方向形成的三幅去噪后高频细节图像和步骤(3)中得到的低频轮廓图像共四幅图像，使用小波逆变换IDWT，得到去噪后的两幅遥感图像的差异图；

(6)分割差异图：

使用基于局部模糊C均值聚类算法FLICM，对去噪后的两幅遥感图像的差异图进行分割，得到最终的变化检测结果图。

2.根据权利要求1所述的基于非局部小波信息的遥感图像变化检测方法，其特征在于，步骤(2)所述的对数比公式如下：

$X=\left|\log \frac{X_2}{X_1}\right|$

其中，X表示两幅遥感图像的差异图，X₁和X₂分别表示同一地区不同时刻已配准的两幅遥感图像，log表示自然对数操作。

3.根据权利要求1所述的基于非局部小波信息的遥感图像变化检测方法，其特征在于，步骤(4e)所述的水平方向高频细节图像中的点i的去噪权重值的归一化常量W_i计算公式如下：

$W_i=\underset{j\∈S_i}{\mathrm{\Σ}}{w}_{\mathrm{ij}}$

其中，W_i表示水平方向高频细节图像中的点i的去噪权重值的归一化常量，S_i表示水平方向高频细节图像中的点i的搜索窗口，w_ij表示水平方向高频细节图像中的点j对点i去噪时的影响权重值，i表示水平方向高频细节图像中的点i，j表示点i的搜索窗口中的点j。

4.根据权利要求1所述的基于非局部小波信息的遥感图像变化检测方法，其特征在于，步骤(4e)所述的水平方向高频细节图像中的点i的搜索窗口为：以水平方向高频细节图像中的点i为中心的21×21个像素点的窗口。

5.根据权利要求1所述的基于非局部小波信息的遥感图像变化检测方法，其特征在于，步骤(6)所述的基于局部模糊C均值聚类算法FLICM的步骤如下：

第一步，初始化：

设定聚类数目c＝2，模糊指数m＝2和迭代截止条件阈值ε＝10^-6，随机初始化模糊隶属度矩阵U⁽⁰⁾，设置循环计数器b＝0，按照下式，进行初迭代截止条件阈值ε＝10^-6始化：

U^(b)＝[u_ij](i＝1,2；j＝1,2,…,N)

其中，U^(b)表示模糊隶属度矩阵，b表示循环计数器，u_ij表示两幅遥感图像的差异图中像素点j属于第i类的模糊隶属度，N表示两幅遥感图像的差异图中像素点的总个数；

第二步，利用下式计算聚类中心：

$v_k=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ki}}^m{x}_i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ki}}^m}$

其中，ν_k表示第k类的聚类中心，μ_ki表示两幅遥感图像的差异图中像素点i属于第k类的模糊隶属度，x_i表示两幅遥感图像的差异图中像素点i的灰度值，m表示模糊指数，N表示两幅遥感图像的差异图中像素点的总个数；

第三步，按照下式，计算空间距离项：

$G_{\mathrm{ki}}=\underset{j\∈N_i}{\mathrm{\Σ}}\frac{1}{d_{\mathrm{ij}}+1}{(1-{\mathrm{\μ}}_{\mathrm{kj}})}^m{\left|\right|x_j-v_k\left|\right|}^2$

其中，G_ki表示两幅遥感图像的差异图中像素点i到第k类聚类中心的空间距离，N_i表示点两幅遥感图像的差异图中像素点i的空间邻域，d_ij表示点i的空间邻域中的点j到点i的欧式距离，d_ij＝||x_j-x_i||²，x_j表示两幅遥感图像的差异图中像素点i的空间邻域中的点j的灰度值，ν_k表示第k类的聚类中心，m表示模糊指数；

第四步，按照下式，计算模糊隶属度矩阵的元素，得到更新后的模糊隶属度矩阵U^(b+1)：

${\mathrm{\μ}}_{\mathrm{ki}}=\frac{1}{\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{\left(\frac{{\left|\right|x_i-v_k\left|\right|}^2+G_{\mathrm{ki}}}{{\left|\right|x_i-v_j\left|\right|}^2+G_{\mathrm{ji}}}\right)}^{1/(m-1)}}$

其中，μ_ki表示两幅遥感图像的差异图中像素点i属于第k类的模糊隶属度，G_ki表示像素点i到第k个聚类中心的空间距离，x_i表示两幅遥感图像的差异图中像素点i的灰度值，ν_k表示第k类的聚类中心，m表示模糊指数，c表示聚类数目；

第五步，判断max{U^(b)-U^(b+1)}＜ε是否成立，如果成立，结束该聚类算法，得到最终的模糊隶属度矩阵U＝U^(b+1)，由模糊隶属度矩阵U得到各像素点的分类结果；否则，令b＝b+1,执行第二步；其中，ε表示迭代截止条件阈值，ε＝10^-6，U^(b)表示更新前的模糊隶属度矩阵，U^(b+1)表示更新后的模糊隶属度矩阵。