结合多阈值分割与模糊聚类的SAR图像变化检测方法

1.一种结合多阈值分割与模糊聚类的SAR图像变化检测方法，包括如下步骤：
(1)中值滤波
选取3×3中值滤波器对两幅待检测的SAR图像进行预处理，得到中值滤波后的两幅图像；
(2)求出归一化后的对数比差异图像
2a)采用对数比差异公式，由中值滤波后的两幅图像，求出对数比差异图像；
2b)采用归一化公式对对数比差异图像进行归一化处理后得到归一化后的对数比差异图像；
(3)多阈值分割
采用基于标准粒子群优化的大津法，对归一化后的对数比差异图像进行多阈值分割，分割后图像中的像素被分成三类：没有发生变化的像素、发生了变化的像素和不能确定是否发生了变化的像素；
所述的基于标准粒子群优化的大津法对归一化后的对数比差异图像进行多阈值分割的具体步骤如下：
第一步，把用来分割图像的n个阈值编码为标准粒子群优化中的一个n维粒子；
第二步，使用标准粒子群优化搜索使得大津法目标函数最大的n个阈值，大津法目标函数如下所示：
2 2 2 2
Maximize σ＝ω0(μ0-μT)+ω1(μ1-μT)+...+ωn(μn-μT)
2
归一化后的对数比差异图像可以被n个阈值分成n+1个部分，在上式中，σ表示这n+1个部分之间的方差，ω0，ω1，...，ωn分别表示各个部分出现的概率，μ0，μ1，...，μn分别表示各个部分的灰度平均值，μT表示整个图像的灰度平均值；
在标准粒子群优化中粒子第d维飞行速度的更新公式为：
其中， 表示种群中第i个粒子的第d维飞行速度，ω表示惯性权重，c1和c2表示学习因子， 表示种群中第i个粒子的第d维位置， 和 是两个相互独立的均匀分布于
0和1之间的随机数，pBesti是第i个粒子的历史最优位置，gBest是整个种群的历史最优位置；
在标准粒子群优化中粒子第d维位置的更新公式为：
其中， 表示种群中第i个粒子的第d维位置， 表示种群中第i个粒子的第d维飞行速度；
第三步，n个阈值把归一化后的对数比差异图像分成n+1个部分，按照灰度值从小到大的顺序记为：C0，C1，...，Cn；把C0中的像素归为没有发生变化的像素，把Cn中的像素归为发生了变化的像素，而把C1至Cn-1中的像素归为不能确定是否发生了变化的像素，从而把归一化后的对数比差异图像中的像素分成了三类；
(4)初始化隶属度矩阵
4a)初始化隶属度矩阵U0，U0存储像素隶属于非变化类的隶属度，隶属度用一个在区间[0，1]上取值的实数来表示，多阈值分割后属于没有发生变化的像素在U0中的隶属度取1，多阈值分割后属于发生了变化的像素在U0中的隶属度取0，多阈值分割后属于不能确定是否发生了变化的像素在U0中的隶属度则随机生成；
4b)初始化隶属度矩阵U1，U1存储像素隶属于变化类的隶属度，隶属度用一个在区间[0，1]上取值的实数来表示，多阈值分割后属于没有发生变化的像素在U1中的隶属度取0，多阈值分割后属于发生了变化的像素在U1中的隶属度取1，多阈值分割后属于不能确定是否发生了变化的像素在U1中的隶属度则等于1减去该像素在U0中的隶属度；
(5)模糊聚类
采用FLICM算法对多阈值分割后属于不能确定是否发生了变化的像素进行模糊聚类，迭代更新直至达到预定的终止条件，即迭代前后隶属度的最大改变量小于0.00001；
(6)去模糊
根据归一化后的对数比差异图像中每个像素隶属于非变化类的隶属度和隶属于变化类的隶属度，将像素判归到隶属度较大的那一类，从而把归一化后的对数比差异图像中所有的像素分成了两类：非变化类和变化类；
(7)输出变化检测结果。
2.根据权利要求1所述的结合多阈值分割与模糊聚类的SAR图像变化检测方法，其特征在于，步骤2a)中所述的对数比差异公式如下：
I3＝|log(I1+1)-log(I2+1)|
其中，I3表示对数比差异图像的像素灰度值，I1和I2分别表示中值滤波后的两幅图像的像素灰度值。
3.根据权利要求1所述的结合多阈值分割与模糊聚类的SAR图像变化检测方法，其特征在于，步骤2b)中所述的归一化公式如下：
其中，ID表示归一化后的对数比差异图像的像素灰度值，I3表示对数比差异图像的像素灰度值，Imax表示I3中最大的灰度值，Imin表示I3中最小的灰度值。
4.根据权利要求1所述的结合多阈值分割与模糊聚类的SAR图像变化检测方法，其特征在于，步骤(5)中所述的FLICM模糊聚类算法的具体实现如下：
FLICM算法使用隶属度表示像素可能属于某个类别的程度，它在聚类过程中同时考虑了图像灰度信息和空间邻域信息，FLICM的目标函数如下所示：
其中，J表示目标函数值，N表示多阈值分割后属于不能确定是否发生了变化的像素组成的集合，C是取值为2的类别数目，uki表示第i个像素在第k类上的隶属度，xi表示第i
2
个像素的灰度值，vk表示第k类的聚类中心，d (xi，vk)表示第i个像素到第k类的聚类中心的距离度量，m表示取值为2的模糊加权指数，Gki表示第i个像素和第k类的聚类中心之间的模糊系数，它的定义如下所示：
其中，Gki表示第i个像素和第k类的聚类中心之间的模糊系数，Ni表示以第i个像素为中心的大小为3×3的局部窗口，xj表示第j个像素的灰度值且该像素是落入局部窗口Ni中的非中心像素，vk表示第k类的聚类中心，ukj表示第j个像素在第k类上的隶属度，eij
2
表示第i个像素和第j个像素之间的空间欧式距离，d(xj，vk)表示第j个像素到第k类的聚类中心的距离度量，m表示取值为2的模糊加权指数；
FLICM算法采用迭代计算的方式求出使得其目标函数最小的一组隶属度值，FLICM算法中隶属度的迭代计算公式如下所示：
其中，uki表示第i个像素在第k类上的隶属度，C是取值为2的类别数目，xi表示第i个像素的灰度值，vk表示第k类的聚类中心，vt表示第t类的聚类中心，Gki表示第i个像素和第k类的聚类中心之间的模糊系数，Gti表示第i个像素和第t类的聚类中心之间的模糊系数，d2(xi，vk)表示第i个像素到第k类的聚类中心的距离度量，d2(xi，vt)表示第i个像素到第t类的聚类中心的距离度量，m表示取值为2的模糊加权指数；
FLICM算法中聚类中心的迭代计算公式如下所示：
其中，vk表示第k类的聚类中心，N表示多阈值分割后属于不能确定是否发生了变化的像素组成的集合，uki表示第i个像素在第k类上的隶属度，xi表示第i个像素的灰度值，m表示取值为2的模糊加权指数；
使用FLICM算法对多阈值分割后属于不能确定是否发生了变化的像素进行模糊聚类，直至达到预定的终止条件，即迭代前后隶属度的最大改变量小于0.00001。

结合多阈值分割与模糊聚类的SAR图像变化检测方法\n技术领域\n[0001] 本发明属于计算机技术领域，更进一步涉及图像处理技术领域中的结合多阈值分割与模糊聚类的SAR图像变化检测方法。本发明由两幅不同时相的SAR图像得到差异图像，再对差异图像进行多阈值分割和模糊聚类，实现SAR图像变化检测，可用于地物覆盖与利用、自然灾害监测与评估，城市规划，地图更新等领域。\n背景技术\n[0002] 合成孔径雷达(Synthetic Aperture Radar，SAR)技术近年来应用得越来越广泛。\nSAR图像与普通光学遥感图像相比，具有可以全天候获取的特点，并且随着SAR图像分辨率不断提高，基于它的图像处理技术也越来越多。图像的变化检测技术是指通过对同一地区在不同时间的两幅遥感图像的对比和分析，得到该地区在这两个时期的地貌变化信息。目前基于SAR图像的变化检测技术已经得到了比较多的应用，比如在森林覆盖变化、自然灾害监测与评估等方面。\n[0003] 对于如何从两幅SAR图像中找出变化的区域，大致有两种方法：一种是分类后比较法，就是先对不同的SAR图像各自进行分类，然后再比较分类后的结果，找出变化的区域；另外一种方法则是从不同的SAR图像中先得到它们的差异图像，再对差异图像进行分析，从而找到变化区域，后一种方法目前应用得多一些。对差异图像的处理实际就是要对它进行一种两类别的划分，因此聚类算法可以用来解决对差异图像的分类问题。\n[0004] 模糊C均值(Fuzzy C-Means，FCM)算法是最流行的聚类算法之一，它采用模糊理论建立样本类别的不确定性描述，通过优化目标函数得到样本对所有聚类中心的隶属度，进而判定出每个样本的类别。但是由于FCM没有考虑样本的空间信息，所以对噪声比较敏感，影响了聚类的效果。模糊局部信息C均值(Fuzzy Local Information C-Means，FLICM)算法是一种比较新的聚类算法，它针对FCM对噪声比较敏感的缺点作了改进，在聚类的目标函数中引入了样本的局部空间信息，取得了比FCM更好的聚类效果，但由于对所有样本都增加了邻域信息的计算，与FCM相比，其运算量较大。\n[0005] 西安电子科技大学申请的专利“基于量子免疫克隆的SAR图像变化检测方法”(专利申请号201010230980.6，公开号CN101908213A)，该方法通过量子比特定义聚类中心，搜索最优的聚类中心并得到全局阈值，但该方法存在的不足是，只使用了FCM的目标函数来构造抗体亲合度函数，并没有考虑像素的邻域信息。\n发明内容\n[0006] 本发明的目的在于克服上述现有技术存在的不足，把多阈值分割和模糊聚类结合起来，用于SAR图像变化检测，以降低变化检测时间，且提高变化检测的精度。\n[0007] 为实现上述目的，本发明在进行变化检测时，首先利用多阈值分割判别出“简单”像素的类别，再将其余“困难”像素留给FLICM模糊聚类算法来判别。\n[0008] 本发明包括如下步骤：\n[0009] (1)中值滤波\n[0010] 选取常用的3×3中值滤波器对两幅待检测的SAR图像进行预处理，得到中值滤波后的两幅图像；\n[0011] (2)求出归一化后的对数比差异图像\n[0012] 2a)采用对数比差异公式，由中值滤波后的两幅图像，求出对数比差异图像；\n[0013] 2b)采用归一化公式对对数比差异图像进行归一化处理后得到归一化后的对数比差异图像；\n[0014] (3)多阈值分割\n[0015] 采用基于标准粒子群优化的大津法，对归一化后的对数比差异图像进行多阈值分割，分割后图像中的像素被分成三类：没有发生变化的像素、发生了变化的像素和不能确定是否发生了变化的像素；\n[0016] (4)初始化隶属度矩阵\n[0017] 4a)初始化隶属度矩阵U0，U0存储像素的隶属于非变化类的隶属度，隶属度用一个在区间[0，1]上取值的实数来表示，多阈值分割后属于没有发生变化的像素在U0中的隶属度取1，多阈值分割后属于发生了变化的像素在U0中的隶属度取0，多阈值分割后属于不能确定是否发生了变化的像素在U0中的隶属度则随机生成；\n[0018] 4b)初始化隶属度矩阵U1，U1存储像素的隶属于变化类的隶属度，隶属度用一个在区间[0，1]上取值的实数来表示，多阈值分割后属于没有发生变化的像素在U1中的隶属度取0，多阈值分割后属于发生了变化的像素在U1中的隶属度取1，多阈值分割后属于不能确定是否发生了变化的像素在U1中的隶属度则等于1减去该像素在U0中的隶属度；\n[0019] (5)模糊聚类\n[0020] 采用FLICM算法对多阈值分割后属于不能确定是否发生了变化的像素进行模糊聚类，迭代更新直至达到预定的终止条件，即迭代前后隶属度的最大改变量小于0.00001；\n[0021] (6)去模糊\n[0022] 根据归一化后的对数比差异图像中每个像素的隶属于非变化类和变化类的隶属度值，将像素判归到隶属度值较大的那一类，从而把归一化后的对数比差异图像中所有的像素分成了两类：非变化类和变化类；\n[0023] (7)输出变化检测结果。\n[0024] 本发明与现有的技术相比具有以下优点：\n[0025] 1.本发明首先利用多阈值分割判别出“简单”像素的类别，再将其余“困难”像素留给FLICM算法来判别，与现有的FLICM算法直接对像素进行聚类和判别相比，本发明降低了变化检测的运算量和运行时间。\n[0026] 2.本发明在使用FLICM算法进行模糊聚类时，局部窗口Ni中的某些非中心像素的类别已经在多阈值分割阶段确定，这些确定的邻域信息改善了现有技术进行模糊聚类的效果，使得本发明进一步提高了变化检测的精度。\n附图说明\n[0027] 图1是本发明的流程图；\n[0028] 图2是仿真实验中对Ottawa地区水灾SAR图像进行变化检测的结果；\n[0029] 图3是仿真实验中对墨西哥地区火灾SAR图像进行变化检测的结果；\n[0030] 图4是仿真实验中对Bern城市水灾SAR图像进行变化检测的结果。\n具体实施方式\n[0031] 参照图1，本发明中结合多阈值分割与模糊聚类的SAR图像变化检测方法包括如下步骤：\n[0032] (1)中值滤波\n[0033] 选取常用的3×3中值滤波器对两幅待检测的SAR图像进行预处理，得到中值滤波后的两幅图像；\n[0034] (2)求出归一化后的对数比差异图像\n[0035] 2a)采用如下的对数比差异公式，由中值滤波后的两幅图像，求出对数比差异图像：\n[0036] I3＝|log(I1+1)-log(I2+1)|\n[0037] 其中，I3表示对数比差异图像的像素灰度值，I1和I2分别表示中值滤波后的两幅图像的像素灰度值；\n[0038] 2b)采用如下的归一化公式对对数比差异图像进行归一化处理后得到归一化后的对数比差异图像：\n[0039] \n[0040] 其中，ID表示归一化后的对数比差异图像的像素灰度值，I3表示对数比差异图像的像素灰度值，Imax表示I3中最大的灰度值，Imin表示I3中最小的灰度值；\n[0041] (3)多阈值分割\n[0042] 采用基于标准粒子群优化的大津法对归一化后的对数比差异图像进行多阈值分割，分割后图像中的像素被分成三类：没有发生变化的像素、发生了变化的像素和不能确定是否发生了变化的像素；\n[0043] 多阈值分割的具体步骤如下：\n[0044] 3a)把用来分割图像的n个阈值编码为标准粒子群优化中的一个n维粒子；\n[0045] 3b)使用标准粒子群优化搜索使得大津法目标函数最大的n个阈值，大津法目标函数如下所示：\n[0046] Maximize σ2＝ω0(μ0-μT)2+ω1(μ1-μT)2+...+ωn(μn-μT)2[0047] 归一化后的对数比差异图像可以被n个阈值分成n+1个部分，在上式中，σ2表示这n+1个部分之间的方差，ω0，ω1，...，ωn分别表示各个部分出现的概率，μ0，μ1，…，μn分别表示各个部分的灰度平均值，μT表示整个图像的灰度平均值；\n[0048] 在标准粒子群优化中粒子第d维飞行速度的更新公式为：\n[0049] \n[0050] 其中，表示种群中第i个粒子的第d维飞行速度，ω表示惯性权重，c1和c2表示学习因子， 表示种群中第i个粒子的第d维位置， 和 是两个相互独立的均匀分布于0和1之间的随机数，pBesti是第i个粒子的历史最优位置，gBest是整个种群的历史最优位置；\n[0051] 在标准粒子群优化中粒子第d维位置的更新公式为：\n[0052] \n[0053] 其中， 表示种群中第i个粒子的第d维位置， 表示种群中第i个粒子的第d维飞行速度；\n[0054] 3c)通过标准粒子群优化的搜索，找到使得大津法目标函数最大的n个阈值，这n个阈值把归一化后的对数比差异图像分成n+1个部分，按照灰度值从小到大的顺序记为：\nC0，C1，...，Cn。把C0中的像素归为没有发生变化的像素，把Cn中的像素归为发生了变化的像素，而把C1至Cn-1中的像素归为不能确定是否发生了变化的像素，从而把归一化后的对数比差异图像中的像素分成了三类；\n[0055] (4)初始化隶属度矩阵\n[0056] 4a)初始化隶属度矩阵U0，U0存储像素的隶属于非变化类的隶属度，隶属度用一个在区间[0，1]上取值的实数来表示，多阈值分割后属于没有发生变化的像素在U0中的隶属度取1，多阈值分割后属于发生了变化的像素在U0中的隶属度取0，多阈值分割后属于不能确定是否发生了变化的像素在U0中的隶属度则随机生成；\n[0057] 4b)初始化隶属度矩阵U1，U1存储像素的隶属于变化类的隶属度，隶属度用一个在区间[0，1]上取值的实数来表示，多阈值分割后属于没有发生变化的像素在U1中的隶属度取0，多阈值分割后属于发生了变化的像素在U1中的隶属度取1，多阈值分割后属于不能确定是否发生了变化的像素在U1中的隶属度则等于1减去该像素在U0中的隶属度；\n[0058] (5)模糊聚类\n[0059] 采用FLICM算法对多阈值分割后属于不能确定是否发生了变化的像素进行模糊聚类，迭代更新直至达到预定的终止条件，即迭代前后隶属度的最大改变量小于0.00001；\n[0060] 5a)FLICM算法使用隶属度表示像素可能属于某个类别的程度，它在聚类过程中同时考虑了图像灰度信息和空间邻域信息，FLICM的目标函数如下所示：\n[0061] \n[0062] 其中，J表示目标函数值，N表示多阈值分割后属于不能确定是否发生了变化的像素组成的集合，C是取值为2的类别数目，uki表示第i个像素在第k类上的隶属度，xi表示\n2\n第i个像素的灰度值，vk表示第k类的聚类中心，d (xi，vk)表示第i个像素到第k类的聚类中心的距离度量，m表示取值为2的模糊加权指数，Gki表示第i个像素和第k类的聚类中心之间的模糊系数，它的定义如下所示：\n[0063] \n[0064] 其中，Gki表示第i个像素和第k类的聚类中心之间的模糊系数，Ni表示以第i个像素为中心的大小为3×3的局部窗口，xj表示第j个像素的灰度值且该像素是落入局部窗口Ni中的非中心像素，vk表示第k类的聚类中心，ukj表示第j个像素在第k类上的隶属\n2\n度，eij表示第i个像素和第j个像素之间的空间欧式距离，d (xj，vk)表示第j个像素到第k类的聚类中心的距离度量，m表示取值为2的模糊加权指数；\n[0065] 5b)FLICM算法采用迭代计算的方式求出使得其目标函数最小的一组隶属度值，FLICM算法中隶属度的迭代计算公式如下所示：\n[0066] \n[0067] 其中，uki表示第i个像素在第k类上的隶属度，C是取值为2的类别数目，xi表示第i个像素的灰度值，vk表示第k类的聚类中心，vt表示第t类的聚类中心，Gki表示第i个像素和第k类的聚类中心之间的模糊系数，Gti表示第i个像素和第t类的聚类中心之间的\n2 2\n模糊系数，d(xi，vk)表示第i个像素到第k类的聚类中心的距离度量，d(xi，vt)表示第i个像素到第t类的聚类中心的距离度量，m表示取值为2的模糊加权指数；\n[0068] 5c)FLICM算法中聚类中心的迭代计算公式如下所示：\n[0069] \n[0070] 其中，vk表示第k类的聚类中心，N表示多阈值分割后属于不能确定是否发生了变化的像素组成的集合，uki表示第i个像素在第k类上的隶属度，xi表示第i个像素的灰度值，m表示取值为2的模糊加权指数；\n[0071] (6)去模糊\n[0072] 根据归一化后的对数比差异图像中每个像素的隶属于非变化类和变化类的隶属度值，将像素判归到隶属度值较大的那一类，从而把归一化后的对数比差异图像中所有的像素分成了两类：非变化类和变化类；\n[0073] (7)输出变化检测结果。\n[0074] 本发明的效果可以通过以下仿真实验进一步说明：\n[0075] 1.仿真实验条件\n[0076] 仿真实验环境：操作系统为Windows XP，CPU为AMD Athlon1.60GHz，内存为\n1.75GB，编程平台为Visual C++6.0。\n[0077] 仿真实验一是Ottawa地区水灾的SAR图像变化检测，图像大小是290×350，图\n2(a)是1997年5月的SAR图像，图2(b)是1997年8月的SAR图像，图2(c)是相应的变化检测参考图，图2(d)是FCM算法的变化检测结果，图2(e)是FLICM算法的变化检测结果，图2(f)是本发明的变化检测结果。\n[0078] 仿真实验二是墨西哥地区火灾的SAR图像变化检测，图像大小是512×512，图\n3(a)是2000年4月的SAR图像，图3(b)是2002年5月的SAR图像，图3(c)是变化检测参考图，图3(d)是FCM算法的变化检测结果，图3(e)是FLICM算法的变化检测结果，图3(f)是本发明的变化检测结果。\n[0079] 仿真实验三是Bern城市水灾的SAR图像变化检测，图像大小为301×301，图4(a)是1999年4月的SAR图像，图4(b)是1999年5月的SAR图像，图4(c)是变化检测参考图，图4(d)是FCM算法的变化检测结果，图4(e)是FLICM算法的变化检测结果，图4(f)是本发明的变化检测结果。\n[0080] 2.仿真实验内容和结果\n[0081] 仿真实验中的参数设置如下：本发明中多阈值分割的阈值数目n设置为3(选择的依据在后面的阈值数目参数分析仿真实验中给出)；本发明中标准粒子群优化的两个学习因子均设为2.0，惯性权重设为随着搜索过程从0.9至0.4线性降低，且标准粒子群优化的终止条件为目标函数评价次数达到5000次；仿真实验中三种方法里的模糊聚类过程的终止条件都设为迭代前后隶属度的最大改变量小于0.00001。以下仿真实验中的数据均为20次独立运行后的平均结果。\n[0082] 仿真实验一：三种方法对这个问题的变化检测结果如图2(d)至(f)所示，对检测结果的具体对比分析见表1。从图2可以看出，本发明的检测结果与参考图最为接近。由表\n1可以看出，本发明漏检的像素数分别比FCM和FLICM少了1246个和1309个，而总的错误像素数也分别比两者少了1003个和372个。同时从运行时间上可以看出，虽然本发明增加了多阈值分割的过程，但也因此降低了后面进行模糊聚类的任务量，所以运行时间比直接使用FLICM算法进行聚类仍然快了很多。\n[0083] 表1Ottawa地区水灾变化检测结果\n[0084] \n[0085] 仿真实验二：三种方法对这个问题的变化检测结果如图3(d)至(f)所示，对检测结果的具体对比分析见表2。由表2可以看出，本发明漏检的像素数分别比FCM和FLICM少了623个和530个，而总的错误像素数也分别比两者少了149个和383个。本发明将只考虑灰度信息的多阈值分割与重点考虑邻域信息的FLICM算法有机地结合起来，提升了抗噪性能，大大降低了漏检数，使得总体检测性能优于两种所对比的方法。同时从表2可以看出本发明的运行时间约为FLICM算法的30％。\n[0086] 表2墨西哥地区火灾变化检测结果\n[0087] \n[0088] 仿真实验三：三种方法对这个问题的变化检测结果如图4(d)至(f)所示，对检测结果的具体对比分析见表3。由表3可以看出，对于该变化检测问题，本发明在漏检数方面仍然具有比较明显的优势，因而总的错误像素数也较少，只有112个，错误检测率只有\n0.12％，而且运行时间也很短，不到FLICM算法的25％。\n[0089] 表3Bern城市水灾变化检测结果\n[0090] \n[0091] 使用下式可以计算得出本发明在上述三组仿真实验中的运行时间比FLICM算法平均降低了大约70％：\n[0092] \n[0093] 其中，T表示本发明与FLICM算法相比运行时间的减少比例，T1表示本发明的运行时间，T2表示FLICM算法的运行时间。\n[0094] 3.阈值数目的参数分析仿真实验\n[0095] 本发明中进行多阈值分割时阈值数目n的选择对变化检测的效果有较大影响。为了合理地选择阈值数目n的大小，对阈值数目n进行了参数分析，在不同的阈值数目下独立运行10次本发明方法，平均后的检测结果如表4所示。从表4可以看出，阈值数目n过小会明显增加变化检测的错误检测率；阈值数目n过大则会引起运行时间的增加。综合考虑后，阈值数目n取3是比较合理的选择。\n[0096] 表4本发明中的阈值数目n的参数分析仿真实验结果\n[0097]