CN105205807A - 基于稀疏自动编码机的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于稀疏自动编码机的遥感图像变化检测方法，主要针对现有技术的变化检测结果漏检与误检的像素较多，且直接对差异图进行处理无法利用差异图隐含信息的缺点。其实现步骤为：(1)读入图像；(2)构造差异图；(3)训练稀疏自动编码机；(4)提取特征；(5)对特征进行模糊聚类；(6)去模糊处理；(7)输出变化检测结果。本发明先采用稀疏自动编码机提取差异图，再利用提取特征对遥感图像进行变化检测，从而减少了变化检测结果中像素的漏检与误检的个数，提高了变化检测的精度。

Description

基于稀疏自动编码机的遥感图像变化检测方法

技术领域

本发明属于计算机技术领域，更进一步涉及图像处理技术领域中的基于稀疏自动编码机的遥感图像变化检测方法。本发明由两幅不同时相的遥感图像得到差异图像，利用稀疏自动编码机对差异图像进行特征提取，再对特征进行分类完成遥感图像的变化检测。本发明可以应用于自然灾害检测与救援时灾区遥感图像的变化检测，城市发展规划，地质研究等领域，完成对特定地区的遥感图像的变化进行检测。

背景技术

变化检测是通过分析同一地区不同时相的遥感图像，检测出该地区中地物随时间发生变化的信息。随着近年来遥感和信息技术的发展，变化检测已经成为当前遥感图像分析研究的一个重要方向，广泛地应用于许多领域，如资源保护中森林覆盖的动态监测、城市规划中土地利用与覆盖的变化信息获取、军事侦察中战略目标的动态监视以及自然灾害评估等。

随着现代社会日益增长的需求和技术的进步，极大地促进了变化检测理论和方法的研究，多样的空间数据产品极大丰富了变化检测理论和方法的研究内容。

夏晨阳，石爱业和吴国宝共同发表的论文“基于SVM混合核的遥感图像变化检测”(《信息技术》，2014年第8期，P38-P41)中提出了一种基于SVM混合核的遥感图像变化检测的方法。该方法首先利用CVA算法构造差异影像了，其次利用灰度共生矩阵提取差异影像的纹理特征与差异影像的灰度特征组成特征向量，接着利用差异影像的直方图选择置信度高的训练样本，并利用构造的SVM混合核进行训练得到分类超平面，最后利用SVM混合核函数对差异影像进行二分类得到最后的变化检测结果。该方法存在的不足之处是，对差异图进行了多种处理，使改方法流程较为复杂。另外该方法没有考虑任何像元的空间信息，以至最终变化检测结果对图像噪声和成像伪影十分敏感，最终的变化检测结果误检和漏检的像素较多。

西安电子科技大学申请的专利“基于区域和Kmeans聚类的遥感图像变化检测方法”(专利申请号201310114150，公开号CN103198480A)公开了一种基于区域的遥感图像变化检测方法。该方法通过对差异图像进行最大熵阈值分析提取感兴趣区域和肯定非变化区域，并利用这两个区域的特征采用Kmeans算法对两个区域进行分类，完成对遥感图像的变化检测。但是该方法仍然存在的不足之处是，直接对差异图像素的辉度值进行处理，对于像素的辉度值的无监督分类并不能利像素的特征，细节丢失严重，最终的变化检测结果效果不佳。

发明内容

本发明的目的在于克服上述现有技术存在的不足，提供一种基于稀疏自动编码机的遥感图像变化检测方法，该方法在传统变化检测的基础上加入稀疏自动编码机，通过训练好的稀疏自动编码机，可以提取差异图的特征，在对特征进行模糊聚类从而进一步获得变化检测结果。该方法思路简单明确，在传统方法上对差异图进行了特征提取，有效的提高了遥感图像变化检测的精度。

本发明包括如下步骤：

(1)输入图像：

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

(2)构造差异图：

按照下式，构造已配准的同一地区不同时刻的两幅遥感图像的差异图：

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

其中，X表示已配准的同一地区不同时刻的两幅遥感图像的差异图，|·|表示绝对值操作，log表示以10为底的对数操作，X₁和X₂分别表示已配准的同一地区不同时刻的两幅遥感图像；

(3)训练稀疏自动编码机：

(3a)将稀疏自动编码机的网络层数设定为输入层、隐含层和输出层的3层；

(3b)将稀疏自动编码机的隐含层设定为20个节点，输出层设定为9个节点；

(3c)初始化稀疏自动编码机网络中所设定节点的权值和偏置；

(3d)将已配准的同一地区不同时刻的两幅遥感图像的差异图矩阵作为训练数据，利用反向传导方法，计算稀疏自动编码机的隐含层与输出层的权值与偏置，完成稀疏自动编码机的训练；

(4)提取特征：

将已配准的同一地区不同时刻的两幅遥感图像的差异图，输入到训练好的稀疏自动编码机中，将稀疏自动编码机隐含层输出的已配准的同一地区不同时刻的两幅遥感图像的差异图的每个像素的20组数据，作为已配准的同一地区不同时刻的两幅遥感图像的差异图的每个像素的特征；

(5)对特征进行模糊聚类：

利用已配准的同一地区不同时刻的两幅遥感图像的差异图的每个像素的特征，采用模糊C均值方法,对已配准的同一地区不同时刻的两幅遥感图像的差异图进行模糊聚类，迭代更新直至迭代前后隶属度的最大改变量小于0.00001；

(6)去模糊处理：

采用模糊C均值方法，根据对已配准的同一地区不同时刻的两幅遥感图像的差异图进行模糊聚类得到的隶属度值，将已配准的同一地区不同时刻的两幅遥感图像的差异图中的每个像素判归到隶属度值大的一类，从而将已配准的同一地区不同时刻的两幅遥感图像的差异图中的所有像素分为变化类和非变化类两类；

(7)输出变化检测结果。

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

第一，由于本发明采用稀疏自动编码机，对差异图进行特征提取，提取的特征中包括了像元空间信息，克服了现有技术采用传统方法变化检测结果对图像噪声和成像伪影十分敏感，最终变化检测结果漏检与误检的像素较多的缺点，使得本发明具有受噪声影响小，变化检测结果正确率高的优点。

第二，由于本发明采用模糊c均值方法对差异图的特征进行模糊聚类，克服了现有技术采用传统方法直接对差异图像的素辉度值进行处理，无法利用差异图像素中的隐含信息这一缺点，使得本发明具有保留变化检测结果细节较为丰富，精度较高的优点。

附图说明

图1是本发明的流程图；

图2是仿真实验中Ottawa地区水灾的遥感图像；

图3是仿真实验中对Ottawa地区水灾遥感图像进行变化检测图及其对比图。

具体实施方式

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

参照附图1，对本发明的实现步骤做进一步的描述。

步骤1,输入图像。

输入已配准的同一地区不同时刻的两幅遥感图像。

步骤2,构造差异图。

按照下式，构造已配准的同一地区不同时刻的两幅遥感图像的差异图：

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

其中，X表示已配准的同一地区不同时刻的两幅遥感图像的差异图，|·|表示绝对值操作，log表示以10为底的对数操作，X₁和X₂分别表示已配准的同一地区不同时刻的两幅遥感图像。

步骤3，训练稀疏自动编码机。

将稀疏自动编码机的网络层数设定为3层，包含输入层、隐含层和输出层。将稀疏自动编码机的隐含层设定为20个节点，输出层设定为9个节点。

在[-0.015,0.015]范围内随机选取360个随机数，用所选取的360个随机数初始化稀疏自动编码机网络中节点的权值，在[-0.015,0.015]范围内随机选取29个随机数初始化稀疏自动编码机网络中节点的偏置。

将已配准的同一地区不同时刻的两幅遥感图像的差异图矩阵作为训练数据，利用反向传导方法，计算稀疏自动编码机的隐含层与输出层的权值与偏置，完成稀疏自动编码机的训练。

反向传传导方法具体步骤如下：

第一步，按照下式，计算稀疏自动编码机的损失值：

$J=\[\frac{1}{m}\mathrm{\Σ}\left(\frac{1}{2}\right|\left|h\right(x)-x|{|}^2)\]+0.5\×{10}^{-5}\×\mathrm{\Σ}{\left(W\right)}^2+\mathrm{\Σ}\left(0.05\right|\left|\widehat{\mathrm{\ρ}}\right)$

其中，J表示稀疏自动编码机的损失值,m表示已配准的同一地区不同时刻的两幅遥感图像的差异图矩阵总列数,Σ表示求和操作，||·||表示范数操作，x表示已配准的同一地区不同时刻的两幅遥感图像的差异图矩阵，h(x)表示稀疏自动编码机输出的矩阵，W表示稀疏自动编码机网络中所有节点的权值，表示相对熵，表示稀疏自动编码机隐含层输出数据的平均数。

第二步，按照下式，计算稀疏自动编码机输出层的残差值：

${\mathrm{\δ}}_2=\frac{\∂}{\∂z_2}\frac{1}{2}\left|\right|x-h\left(x\right)|{|}^2$

其中，δ₂表示稀疏自动编码机输出层的残差值，表示偏导操作，z₂表示稀疏自动编码机输出层的输入数据，||·||表示范数操作，x表示已配准的同一地区不同时刻的两幅遥感图像的差异图矩阵，h(x)表示稀疏自动编码机输出的矩阵。

第三步，按照下式，计算稀疏自动编码机隐含层的残差值：

${\mathrm{\δ}}_1=(\mathrm{\Σ}W\×{\mathrm{\δ}}_2)\×\frac{e^{(-z_2)}}{{(1+e^{(-x)})}^2}$

其中，δ₁表示稀疏自动编码机隐含层的残差值，Σ表示求和操作，W表示稀疏自动编码机网络中所有节点的权值，δ₁表示稀疏自动编码机输出层的残差值，e表示自然常数，z₁表示稀疏自动编码机隐含层的输入数据。

第四步，按照下式，更新稀疏自动编码机输出层的权值和偏置：

$\stackrel{\‾}{W_2}=W_2-a_2{\mathrm{\δ}}_2$

$\stackrel{\‾}{b_2}=b_2-{\mathrm{\δ}}_2$

其中，表示更新后的稀疏自动编码机输出层节点的权值，W₂表示稀疏自动编码机输出层节点的权值，表示更新后的稀疏自动编码机输出层节点的偏置，b₂表示稀疏自动编码机输出层节点的偏置，a₂表示稀疏自动编码机输出层的输出数据，δ₂表示稀疏自动编码机输出层的残差值。

第五步，按照下式，更新稀疏自动编码机隐含层的权值和偏置：

$\stackrel{\‾}{W_1}=W_1-a_1{\mathrm{\δ}}_1$

$\stackrel{\‾}{b_1}=b_1-{\mathrm{\δ}}_1$

其中，表示更新后的稀疏自动编码机隐含层节点的权值，W₁表示稀疏自动编码机隐含层节点的权值，表示跟新后稀疏自动编码机隐含层节点的偏置，b₁表示稀疏自动编码机隐含层节点的偏置，a₁表示稀疏自动编码机隐含层的输出数据，δ₁表示稀疏自动编码机隐含层的残差值。

第六步，重复执行第一步、第二步、第三步、第四步、第五步，得到稀疏自动编码机的损失值，用所获得的损失值更新稀疏自动编码机的权值与偏置。

第七步，判断相邻两次稀疏自动编码机的损失值的变化量是否小0.001，若是，则停止更新稀疏自动编码机的隐含层与输出层的权值与偏置，否则，执行第一步、第二步、第三步、第四步、第五步、第七步。

步骤4，提取特征。

将已配准的同一地区不同时刻的两幅遥感图像的差异图，输入到训练好的稀疏自动编码机中，将稀疏自动编码机隐含层将输出的已配准的同一地区不同时刻的两幅遥感图像的差异图的每个像素的20组数据，作为已配准的同一地区不同时刻的两幅遥感图像的差异图的每个像素的特征，完成特征的提取。

步骤5，对特征进行模糊聚类。

利用已配准的同一地区不同时刻的两幅遥感图像的差异图的每个像素的特征，采用模糊C均值方法,对已配准的同一地区不同时刻的两幅遥感图像的差异图进行模糊聚类，迭代更新直至迭代前后隶属度的最大改变量小于0.00001。

模糊C均值方法具体步骤如下：

第一步，按照下式，计算模糊C均值方法的目标函数值：

$P=\underset{i=1}{\overset{N}{\Σ}}\underset{k=1}{\overset{2}{\Σ}}\[u_{ki}\×d(x_i,v_k)\]$

其中，P表示目标函数值，Σ表示求和操作，N表示差异图的像素的个数，i表示差异图的像素的序号，k表示对差异图模糊聚类后的类别，u_ki表示差异图的第i个像素在第k类上的隶属度，x_i表示差异图的第i个像素的特征，v_k表示差异图的第k类的聚类中心，d(x_i,v_k)表示差异图的第i个像素到第k类的聚类中心的欧氏距离。

第二步，按照下式，计算差异图的隶属度：

$u_{ki}=\frac{m(x,v)}{d(x_i,v_k)}$

其中，u_ki表示差异图的第i个像素在第k类上的隶属度，m(x,v)表示差异图的每一个像素到每一类的聚类中心的所有欧氏距离的和，x表示差异图像素的特征，v表示差异图的聚类中心，d(x_i,v_k)表示差异图的第i个像素到第k类的聚类中心的欧氏距离，x_i表示差异图的第i个像素的特征，v_k表示差异图的第k类的聚类中心，v_t表示差异图的第t类的聚类中心。

第三步，按照下式，更新差异图的聚类中心：

$v_k=\frac{\underset{i=1}{\overset{N}{\Σ}}{u}_{ki}{x}_i}{\underset{i=1}{\overset{N}{\Σ}}{u}_{ki}}$

其中，v_k表示差异图的第k类的聚类中心，k表示对差异图模糊聚类后的类别，Σ表示求和操作，N表示差异图的像素的个数，i表示差异图的像素的序号，u_ki表示差异图的第i个像素在第k类上的隶属度，k表示对差异图模糊聚类后的类别，x_i表示差异图的第i个像素的特征。

第四步，重复执行第一步，第二步，第三步，判断迭代前后隶属度的最大改变量是否小于0.00001，若是，停止迭代，否则，执行第一步。

步骤6，去模糊处理。

根据对已配准的同一地区不同时刻的两幅遥感图像的差异图进行模糊聚类得到的隶属度值，将已配准的同一地区不同时刻的两幅遥感图像的差异图中的每个像素判归到隶属度值大的一类，从而将已配准的同一地区不同时刻的两幅遥感图像的差异图中的所有像素分为变化类和非变化类两类。

步骤7，输出变化检测结果。

下面结合附图2和附图3对本发明的效果做进一步说明。

1.仿真实验条件：

本实验的硬件测试平台是：处理器为IntelCore2CPU，主频为2.33GHz，内存2GB，软件平台为：WindowsXP操作系统和MatlabR2012a。

本发明处理的图像为两幅Ottawa地区水灾的遥感图像，如图2所示，图2(a)和图2(b)的拍摄时间分别为2004年4月和2004年5月，大小均为290×350。

对于具有参考图的实验仿真图，可进行定量的变化检测分析：

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

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

③正确分类的概率PCC:PCC＝(TP+TN)/(TP+FP+TN+FN)；

④衡量检测结果图与参考图一致性的Kappa系数：

Kappa＝(PCC-PRE)/(1-PRE),其中：

PRE＝[(TP+FP)×Nc+(FN+TN)×Nu]/N2

这里，N表示总像素个数，Nc和Nu分别表示实际的变化像素和未变化像素数。

2.仿真内容：

在实验中分别对Ottawa地区水灾的遥感图像做基于SVM的变化检测、基于Kmeans聚的变化检测和基于稀疏自动编码机的变化检测，并对比了3中种方法变化检测图的效果。

3.仿真结果分析：

反映Ottawa地区水灾的遥感图像如图2所示，图2(a)和图2(b)的拍摄时间分别为2004年4月和2004年5月，大小均为290×350。其中，图3(a)是Ottawa地区水灾的变化检测参考图，图3(b)是基于SVM的变化检测图，图3(c)是基于Kmeans聚类的变化检测图和基于稀疏自动编码机的变化检测图，图3(d)是基于稀疏自动编码机的变化检测图。通过与变化检测参考图3(a)对比，可以看出基于SVM的变化检测方法，变化检测图的漏检与误检像素点较多，变化检测结果较差。基于Kmeans聚类的变法检测方法噪声点非常多，细节丢失较为严重，有部分多变化的区域没有检测出来。而本方法无论是最终变化检测结果的正确率还是对噪声的鲁棒性都比前两种方法要好。两种方法的变化检测结果对比如表1所示。

表1三种变化检测方法最终结果的对比

从表1中可以看出本发明方法的变化检测结果无论是漏检数还是错检数都远远低于前两种对比方法，正确检测率和Kappa系数也较高，可见本发明方法提高了遥感图像变化检测的效果。

Claims (3)

1.基于稀疏自动编码机的遥感图像变化检测方法，包括如下步骤：

(1)输入图像：

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

(2)构造差异图：

按照下式，构造已配准的同一地区不同时刻的两幅遥感图像的差异图：

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

其中，X表示已配准的同一地区不同时刻的两幅遥感图像的差异图，|·|表示绝对值操作，log表示以10为底的对数操作，X₁和X₂分别表示已配准的同一地区不同时刻的两幅遥感图像；

(3)训练稀疏自动编码机：

(3a)将稀疏自动编码机的网络层数设定为输入层、隐含层和输出层的3层；

(3b)将稀疏自动编码机的隐含层设定为20个节点，输出层设定为9个节点；

(3c)初始化稀疏自动编码机网络中所设定节点的权值和偏置；

(3d)将已配准的同一地区不同时刻的两幅遥感图像的差异图矩阵作为训练数据，利用反向传导方法，计算稀疏自动编码机的隐含层与输出层的权值与偏置，完成稀疏自动编码机的训练；

(4)提取特征：

将已配准的同一地区不同时刻的两幅遥感图像的差异图，输入到训练好的稀疏自动编码机中，将稀疏自动编码机隐含层输出的已配准的同一地区不同时刻的两幅遥感图像的差异图的每个像素的20组数据，作为已配准的同一地区不同时刻的两幅遥感图像的差异图的每个像素的特征；

(5)对特征进行模糊聚类：

利用已配准的同一地区不同时刻的两幅遥感图像的差异图的每个像素的特征，采用模糊C均值方法,对已配准的同一地区不同时刻的两幅遥感图像的差异图进行模糊聚类，迭代更新直至迭代前后隶属度的最大改变量小于0.00001；

(6)去模糊处理：

根据对已配准的同一地区不同时刻的两幅遥感图像的差异图进行模糊聚类得到的隶属度值，将已配准的同一地区不同时刻的两幅遥感图像的差异图中的每个像素判归到隶属度值大的一类，从而将已配准的同一地区不同时刻的两幅遥感图像的差异图中的所有像素分为变化类和非变化类两类；

(7)输出变化检测结果。

2.根据权利要求1所述的基于稀疏自动编码机的遥感图像变化检测方法，其特征在于，步骤(3d)中所述反向传导方法的具体步骤如下：

第1步，按照下式，计算稀疏自动编码机的损失值：

$J=\[\frac{1}{m}\mathrm{\Σ}\left(\frac{1}{2}\right|\left|h\right(x)-x|{|}^2)\]+0.5\×{10}^{-5}\×\mathrm{\Σ}{\left(W\right)}^2+\mathrm{\Σ}\left(0.05\right|\left|\widehat{\mathrm{\ρ}}\right)$

其中，J表示稀疏自动编码机的损失值,m表示已配准的同一地区不同时刻的两幅遥感图像的差异图矩阵总列数,Σ表示求和操作，||·||表示范数操作，x表示已配准的同一地区不同时刻的两幅遥感图像的差异图矩阵，h(x)表示稀疏自动编码机输出的矩阵，W表示稀疏自动编码机网络中所有节点的权值，表示相对熵，表示稀疏自动编码机隐含层输出数据的平均数；

第2步，按照下式，计算稀疏自动编码机输出层的残差值：

${\mathrm{\δ}}_2=\frac{\∂}{\∂z_2}\frac{1}{2}\left|\right|x-h\left(x\right)|{|}^2$

其中，δ₂表示稀疏自动编码机输出层的残差值，表示偏导操作，z₂表示稀疏自动编码机输出层的输入数据，||·||表示范数操作，x表示已配准的同一地区不同时刻的两幅遥感图像的差异图矩阵，h(x)表示稀疏自动编码机输出的矩阵；

第3步，按照下式，计算稀疏自动编码机隐含层的残差值：

${\mathrm{\δ}}_1=(\mathrm{\Σ}W\×{\mathrm{\δ}}_2)\×\frac{e^{(-z_2)}}{{(1+e^{(-x)})}^2}$

其中，δ₁表示稀疏自动编码机隐含层的残差值，Σ表示求和操作，W表示稀疏自动编码机网络中所有节点的权值，δ₁表示稀疏自动编码机输出层的残差值，e表示自然常数，z₁表示稀疏自动编码机隐含层的输入数据；

第4步，按照下式，更新稀疏自动编码机输出层的权值和偏置：

$\stackrel{\‾}{W_2}=W_2-a_2{\mathrm{\δ}}_2$

$\stackrel{\‾}{b_2}=b_2-{\mathrm{\δ}}_2$

其中，表示更新后的稀疏自动编码机输出层节点的权值，W₂表示稀疏自动编码机输出层节点的权值，表示更新后的稀疏自动编码机输出层节点的偏置，b₂表示稀疏自动编码机输出层节点的偏置，a₂表示稀疏自动编码机输出层的输出数据，δ₂表示稀疏自动编码机输出层的残差值；

第5步，按照下式，更新稀疏自动编码机隐含层的权值和偏置；

$\stackrel{\‾}{W_1}=W_1-a_1{\mathrm{\δ}}_1$

$\stackrel{\‾}{b_1}=b_1-{\mathrm{\δ}}_1$

其中，表示更新后的稀疏自动编码机隐含层节点的权值，W₁表示稀疏自动编码机隐含层节点的权值，表示跟新后稀疏自动编码机隐含层节点的偏置，b₁表示稀疏自动编码机隐含层节点的偏置，a₁表示稀疏自动编码机隐含层的输出数据，δ₁表示稀疏自动编码机隐含层的残差值；

第6步，重复执行第1步、第2步、第3步、第4步、第5步，得到稀疏自动编码机的损失值，用所获得的损失值更新稀疏自动编码机的权值与偏置；

第7步，判断相邻两次稀疏自动编码机的损失值的变化量是否小于0.001，若是，则停止更新稀疏自动编码机的隐含层与输出层的权值与偏置，否则，执行第1步、第2步、第3步、第4步、第5步、第7步。

3.根据权利要求1所述的基于稀疏自动编码机的遥感图像变化检测方法，其特征在于，步骤(5)中所述模糊C均值方法的具体步骤如下：

第1步，按照下式，计算模糊C均值方法的目标函数值：

$P=\underset{i=1}{\overset{N}{\Σ}}\underset{k=1}{\overset{2}{\Σ}}\[u_{ki}\×d(x_i,v_k)\]$

其中，P表示目标函数值，Σ表示求和操作，N表示差异图的像素的个数，i表示差异图的像素的序号，k表示对差异图模糊聚类后的类别，u_ki表示差异图的第i个像素在第k类上的隶属度，x_i表示差异图的第i个像素的特征，v_k表示差异图的第k类的聚类中心，d(x_i,v_k)表示差异图的第i个像素到第k类的聚类中心的欧氏距离；

第2步，按照下式，计算差异图的隶属度：

$u_{ki}=\frac{m(x,v)}{d(x_i,v_k)}$

其中，u_ki表示差异图的第i个像素在第k类上的隶属度，m(x,v)表示差异图的每一个像素到每一类的聚类中心的所有欧氏距离的和，x表示差异图像素的特征，v表示差异图的聚类中心，d(x_i,v_k)表示差异图的第i个像素到第k类的聚类中心的欧氏距离，x_i表示差异图的第i个像素的特征，v_k表示差异图的第k类的聚类中心，v_t表示差异图的第t类的聚类中心；

第3步，按照下式，更新差异图的聚类中心：

$v_k=\frac{\underset{i=1}{\overset{N}{\Σ}}{u}_{ki}{x}_i}{\underset{i=1}{\overset{N}{\Σ}}{u}_{ki}}$

其中，v_k表示差异图的第k类的聚类中心，k表示对差异图模糊聚类后的类别，Σ表示求和操作，N表示差异图的像素的个数，i表示差异图的像素的序号，u_ki表示差异图的第i个像素在第k类上的隶属度，k表示对差异图模糊聚类后的类别，x_i表示差异图的第i个像素的特征；

第4步，重复执行第1步、第2步、第3步，判断迭代前后隶属度的最大改变量是否小于0.00001，若是，则停止迭代，否则，执行第1步。