基于自动差分聚类算法的遥感图像分割方法

摘要

本发明公开了一种基于自动差分聚类算法的遥感图像分割方法，主要解决现有图像分割技术中计算复杂度较高、分割效果不好的问题。其实现步骤为：(1)输入待分割图像，提取待分割图像特征；(2)产生聚类数据；(3)随机抽取聚类数据初始化种群；(4)根据个体的标签位激活聚类中心；(5)根据激活的聚类中心，计算个体适应度值；(6)用改进的差分进化方法进化种群；(7)对进化后种群进行类别数振荡操作；(8)利用FCM更新质心；(9)利用更新后质心判断终止条件，并记录最优个体；(10)对最优个体解码，分配类别标号并输出分割图像。本发明具有分割精度高、边缘定位准确的优点，可用于目标识别。

说明书

技术领域

本发明属于图像处理技术领域，更进一步涉及遥感图像分割方法，可用于目标识 别。

背景技术

图像分割就是把图像分成若干个特定的、具有独特性质的区域并提出感兴趣目标 的技术和过程。目前，人们更多采用基于聚类分析的方法来进行图像分割。用基于聚 类分析的方法分割图像是将图像空间中的像素用对应的特征空间点表示，根据它们在 特征空间的聚集对特征空间进行分割，然后将它们映射回原图像空间，从而达到图像 分割的目的。

为了更加准确、全面的获取图像分割信息，近几年出现了一些应用自动聚类方法 实现图像分割的技术，可以在图像分割类别数未知的情况下对图像进行分割，而且可 以达到更好的分割结果和分割精度。

U.Maulik在学术期刊“Automatic Fuzzy Clustering Using Modified Differential  Evolution for Image Classification”（IEEE Trans.Geosci.Remote Sens.VOL.48,NO.9,Sep. 2010）中公开了一种基于改进差分进化的自动模糊聚类的图像分割方法。该方法首先 对传统的差分进化方法进行了改进，提出了一种加入模糊策略的自动聚类方法，然后 将其用于图像分割，最后通过XB指标获得最终的图像分割结果。虽然该方法在区域一 致性和边缘保持方面有一定的改善效果，但是该方法仍然存在的不足是，由于采用了 过多的进化技术，增加了整个分割过程的计算复杂度，使得分割速度较慢，同时，该 方法选取的聚类指标不恰当，图像分割结果较差。

发明内容

本发明针对上述现有技术存在的不足，提出了一种新的基于自动差分聚类算法的 遥感图像分割方法。以解决现有图像分割技术中收敛速度慢、稳定性较差，计算复杂 度高、细节保持性能不好的缺点，提高图像分割的精度。

本发明实现上述目的的思路是：在对输入的图像进行特征提取和分水岭分割获得 聚类数据后，先随机抽取聚类数据初始化种群，再使用PBM指标评价聚类性能，然后 用改进的差分进化方法进化种群，最后对种群进行类别数振荡操作并用FCM方法更新 质心，通过不断的迭代更新获得最优个体，最优个体所对应的类别标号作为像素的灰 度值，得到图像分割结果。其实现步骤包括如下：

(1)输入待分割图像I，分别提取待分割图像I的小波特征向量和纹理特征向量， 并用小波特征向量和纹理特征向量表示待分割图像I的每一个像素点v；

(2)产生聚类数据：

用待分割图像I的膨胀变换减去其腐蚀变换，得到形态梯度图像I₁，计算形态 梯度图像I₁的灰度值矩阵的平方积，得到反映图像边缘的浮点活动图像I₂，用分水岭 方法初分割浮点活动图像I₂，得到不同的图像块，对每一个图像块的所有像素点的特 征取平均值，获得代表初始聚类数据的每一块的特征值；

(3)随机抽取聚类数据初始化种群：

种群中每个个体编码长度为L+L×d，每个个体分为两个部分，前L位为标签位， 后L×d位为L个待激活的聚类中心，d为聚类数据的维数，随机初始化每个个体的 标签位，并随机抽取L个聚类数据作为待激活的聚类中心，令当前迭代次数t=1；

(4)根据每个个体的标签位激活相应个体的聚类中心；

判断个体标签T每一位值的大小是否大于0.5，如果大于0.5，则激活相应的聚 类数据点，否则不予激活；所有被激活的聚类数据点组合成为该个体的聚类中心；

(5)根据激活的聚类中心，采用PBM有效性指标公式计算每个个体适应度值；

(6)用改进的差分进化方法对种群进行变异和交叉，并利用个体适应度值进行种 群更新；

(7)对更新后的种群进行类别数振荡操作：

7a)对种群中的每个个体依照振荡规则进行类别数振荡操作，每个个体获得新的 类别数；

7b)计算每个个体中L个数据点的密度，并依照密度大小进行排序；

7c)将新的类别数与旧的类别数进行比较，如果小于，则依照密度排名从该个体 中选择密度较大的聚类数据点加入聚类中心，否则，从已有聚类中心中淘汰密度排名 较小的聚类数据点。

(8)利用FCM方法对振荡操作后的聚类中心进行更新；

(9)判断种群迭代的次数是否达到使用者设定的最大进化代数T=100，或者相邻 两代最优个体的适应度差值是否小于规定数值E=10^-4，如果是，则保存最优个体并执 行步骤(10)，否则，返回步骤(4)；

(10)对获得的最优个体进行解码，求得每个聚类数据的类别标号，并用这些类别 标号作为像素的灰度值，输出图像的分割结果。

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

第一，本发明由于在图像分割预处理中，采用了小波分解方法和灰度共生矩阵方 法提取了小波特征向量和纹理特征向量，获得了更多的图像细节信息，解决了现有技 术仅提取小波特征丢失图像细节信息的缺点，使得本发明提高了整体分割精度。

第二，本发明在图像分割的前处理中采用了分水岭初分割图像获得聚类图像块， 用块的平均特征向量代表块特征，克服了现有技术直接对像素点聚类计算量大的缺 点，使得本发明具有更快的分割速度。

第三，本发明在图像分割的聚类过程中采用了自动聚类思想，克服了现有技术对 聚类类别数的依赖性，使得本发明可以在获得图像最优分割结果的前提下同时自动获 得图像的正确聚类类别数，这样更加符合实际问题的研究。

第四，本发明在图像分割的聚类过程中采用了改进的差分进化方法以及聚类类别 数振荡策略，能够在复杂的解空间中进行有效的搜索并保持种群的多样性，克服了已 有技术容易陷入局部最优影响分割结果的缺点，使得本发明可以获得更加准确的区域 一致性和更好的边缘保持性能。

附图说明

图1为本发明的流程图；

图2为本发明与现有技术在一幅光学遥感图像上的分割结果对比图；

图3为本发明与现有技术在一幅Ku波段SAR图像上的分割结果对比图。

具体实施方式

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

步骤1，输入待分割图像I，分别提取待分割图像I的小波特征向量和纹理特征 向量，并用小波特征向量和纹理特征向量表示待分割图像I的每一个像素点v。

1a)利用小波分解方法获取小波特征向量：

小波分解方法采用了对图像进行窗口大小为16×16的三层小波变换，得到由子带 系数所构成的10维小波特征向量。

1b)利用灰度共生矩阵方法提取纹理特征向量，即先将待处理图像量化为16个灰 度级，再依次令两个像素点连线与横轴的方向夹角为0°、45°、90°和135°，按照下式 分别计算四个方向的灰度共生矩阵：

P(i,j)＝#{(x₁,y₁),(x₂,y₂)∈M×N|f(x₁,y₁)＝r,f(x₂,y₂)＝s}

其中，P(i，j)为灰度共生矩阵在坐标(i，j)位置上的元素，#为集合{}的元素个数，(x₁，y₁) 和(x₂，y₂)为距离等于1的两个像素点坐标，∈为集合中的属于符号，M×N为待处理 图像的大小，|为概率论中的条件符号，r为(x₁，y₁)处像素点矢量化后的灰度值，s为 (x₂，y₂)处像素点矢量化后的灰度值；然后根据灰度共生矩阵分别选取该矩阵四个方向 上的对比度、同质性和能量值，获得像素的12维纹理特征向量；

1c)用小波特征向量和纹理特征向量表示待分割图像的每一个像素点：

每个像素点用22维特征向量表示，其中小波特征向量10维，纹理特征向量12 维。

步骤2,产生聚类数据

2a)用待分割图像I的膨胀变换减去其腐蚀变换，得到形态梯度图像I₁，形态 梯度图像I₁用每个坐标点的灰度值表示，灰度值计算公式为：

I₁(x,y)=δ(I)(x,y)-ε(I)(x,y)

其中，I₁(x,y)为形态梯度图像I₁每个点坐标的灰度值，δ(I)和ε(I)分别为图像的 膨胀算子和腐蚀算子，(x,y)为距离等于1的像素点坐标；

2b)计算形态梯度图像I₁的灰度值矩阵的平方积，得到反映图像边缘的浮点活 动图像I₂；

2c)用分水岭方法初分割浮点活动图像I₂，得到不同的图像块；

2d)对每一个图像块的所有像素点的特征取平均值，获得代表初始聚类数据的 每一块的特征值。

步骤3，随机抽取聚类数据初始化种群。

设置种群规模N＝20、最大迭代次数T＝100和当前迭代次数t=0,种群中每个个体 编码长度为L+L×d，每个个体分为两个部分，前L位为标签位，后L×d位为L个待 激活的聚类中心，d为聚类数据的维数，初始化时为每一位赋值[0,1]之间的数。

步骤4，根据每个个体的标签位的值激活相应个体的聚类中心，如果标签位每一 位的值大于0.5，则激活相应的聚类数据点，否则不予激活；将所有被激活的聚类数 据点组合成为该个体的聚类中心。

步骤5，根据激活的聚类中心，采用PBM有效性指标公式计算每个个体适应度 值：

$\mathrm{PBM}={(\frac{1}{K}\times \frac{E_1}{E_K}\times D_K)}^2$

其中，PBM为待测染色体的适应度值；K为该染色体当前的类别数；E1为样本 点到聚类中心的距离总和，其值为常量；E_K为聚类数据集的类内距离和，D_K为聚类 中心的最大类间距离；

$E_K=\underset{k=1}{\overset{K}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}{u}_{k,j}\left|\right|x_j-z_k\left|\right|$

$D_K=\underset{i,j=1}{\overset{K}{\mathrm{mak}}}\left|\right|z_i-z_j\left|\right|$

其中，n为数据点的个数；u_k,j为第j个数据点在第k类的隶属度；x_j为第j 个数据点；z_k为第k个聚类中心；D_K为聚类中心的最大类间距离，z_i为第i个聚类 中心；z_j为第j个聚类中心，聚类中心为一个d维矢量；||||表示求矢量z_i-z_j的范数。

步骤6，用改进的差分进化方法对种群进行变异和交叉，并利用个体适应度值进 行种群更新；

6a）用改进的差分进化方法对种群进行变异操作，变异操作的公式为：

$V_i(t+1)=\left(\begin{array}{c}{Z}_i\left(t\right)+F(Z_n\left(t\right)-Z_m\left(t\right)),\mathrm{rand}\left(\mathrm{0,1}\right)>\alpha \\ Z_{\mathrm{GBest}}\left(t\right)+\alpha (Z_{\mathrm{LBest}}\left(t\right)-Z_j\left(t\right)),\mathrm{rand}\left(\mathrm{0,1}\right)<\alpha \end{array}\right)$

其中，V_i(t+1)为第i条染色体经过变异操作后的得到的个体；t为当前的迭代次 数；Z_i(t)为当代种群中第i条染色体；i，m，n为随机产生的自然数，其值不能大于 种群规模；rand(0,1)为[0,1]之间的随机数；Z_GBest(t)为全局最优个体；Z_LBest(t)为 当代种群中的局部最优个体；Z_j(t)为当代随机个体；F为传统差分进化的变异因子，

F＝0.5*(1+rand(0,1))，式中，α为改进后差分进化的变异因子，

α＝1/[1+exp(-(1/t))]，式中，exp(·)为以自然对数e为底的指数函数。

6b）用改进的差分进化方法对种群进行交叉操作，交叉操作的公式为：

其中，Q_k,j(t+1)为第k条染色体第j位经过交叉操作后得到的值；k为第k条染 色体；j为染色体的第j位；v_k,j(t+1)为第k条染色体第j位经过变异操作后得到的值； Z_k,j(t)为当代种群中第k条染色体第j位的值；rand_j(0,1)为[0,1]之间的随机数； rand(k)为[1，k]之间的随机数，CR为变异概率，公式为：

CR＝(C_rmax-C_rmin)*(T-t)/T，

式中，Cr_max为最大变异概率，取值为1.0；Cr_min为最小变异概率，取值为0.5；

6c）利用个体适应度值进行种群更新，种群更新公式如下：

其中，Z_i(t+1)为第i个个体经过差分进化后得到的新个体；Q_i(t+1)为第i个个 体经过变异交叉后得到的个体；f(·)为适应度值函数。

步骤7，对更新后的种群进行类别数振荡操作。

7a)对种群中的每个个体依照振荡规则进行类别数振荡操作：

7a1)修正类别数：在种群进化的前期，将当前迭代次数t与最大迭代次数T 进行比较，若其比值小于0.33，则将当前染色体类别数k_old加上[-3,3]之间的一个随机 数，作为染色体新的类别数k_new；若比值处于[0.33,0.66]之间，则将当前染色体的类别 数k_old加上[-2,2]之间的一个随机数，作为新的类别数k_new；若比值大于0.66，则将当 前染色体的类别数k_old加上[-1,1]之间的一个随机数，作为新的类别数k_new；

7a2)按照修正后的类别数对聚类数据进行重新聚类。

7b)计算每个个体中L个数据点的密度，并依照密度大小进行排序；

7c)将新的类别数k_new与当前染色体类别数k_old进行比较，如果k_new＜k_old，则依照 密度排名从该个体中选择密度较大的聚类数据点加入聚类中心，否则，从已有聚类中 心中淘汰密度排名较小的聚类数据点。

步骤8，利用FCM方法对振荡操作后的聚类中心进行更新；

FCM方法的聚类中心更新公式为：

$z_k=\frac{{\Sigma }_{i=1}^n{u}_{k,i}^m{x}_i}{{\Sigma }_{i=1}^n{u}_{k,i}^m},1\le k\le K$

其中，z_k为更新操作后的聚类中心，n为数据点的个数，m为模糊指数，取 值m=2；x_i为第i个数据点，K为当前个体的类别数，u_k,j为第j个数据点在第k类 的隶属度，

$u_{k,j}=\frac{{\left(\frac{1}{D(z_k,x_i)}\right)}^{\frac{1}{m-1}}}{{\Sigma }_{j=1}^K{\left(\frac{1}{D(z_j,x_i)}\right)}^{\frac{1}{m-1}}},1\le k\le K,1\le i\le n$

其中，D(z_k,x_i)为第i个数据点与第k个聚类中心的欧氏距离，D(z_j,x_i)为第i 个数据点与第j个聚类中心的欧氏距离，z_k为第k个聚类中心，z_j为第j个聚类中心。

步骤9，判断种群迭代的次数是否达到使用者设定的最大进化代数T=100，或者 相邻两代最优个体的适应度差值是否小于规定数值E=10^-4，如果是，则保存最优个体 并执行步骤(10)，否则，返回步骤(4)。

步骤10，对获得的最优个体进行解码，求得每个聚类数据的类别标号，并用这 些类别标号作为像素的灰度值，输出图像的分割结果。

解码的具体步骤是：

首先，根据最优个体标签确定最优个体的聚类中心；

其次，计算聚类数据中每一个数据点到各个聚类中心的欧式距离，其中最短欧式 距离所在的聚类中心标号就是该数据点的类别标号；

最后，将类别标号赋给相应图像块的每一个像素。

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

1.仿真条件

本发明的仿真是在主频2.5GHZ的Pentium Dual_Core CPU E5200、内存2GB的 硬件环境和MATLAB R2009a的软件环境下进行的。

2.仿真内容

仿真一：用本发明与现有的方法对来自圣地亚哥地区牛尾洲的港口图的部分光学 遥感图像数据进行测试分割，仿真结果如图2，其中：

图2(a)为原始测试遥感图像，该图像有两种类标，一种为陆地，另一种为港口， 图像大小为256×256；

图2(b)为用本发明的分割结果图；

图2(c)为现有技术中的An Evolutionary Approach to Multiobjective Clustering的 分割结果图；

图2(d)为现有技术中的遗传算法聚类技术的分割结果；

图2(e)为现有技术中的K均值方法的分割结果图。

仿真二：用本发明与现有的方法对一副具有1米分辨率的Ku波段的SAR图像数 据进行测试分割，仿真结果如图3，其中：

图3(a)为原始测试SAR图像，此SAR图像有三种地物类别，即：河流，植被和 农作物,图像大小为256×256；

图3(b)为本发明的分割结果图；

图3(c)为现有技术中的An Evolutionary Approach to Multiobjective  Clustering(MOCK)的分割结果图；

图3(d)为现有技术中的遗传算法聚类技术的分割结果图；

图3(e)为现有技术中的K均值方法的分割结果图。

3.仿真效果分析

仿真一：用本发明与现有的方法对来自圣地亚哥地区牛尾洲的港口图的部分光学 遥感图像数据进行测试分割，仿真结果如图2，对其仿真效果进行分析，其中：

图2(a)为原始测试遥感图像，该图像有两种类标，一种为陆地，另一种为港口， 图像大小为256×256；

图2(b)为用本发明的分割结果图，从图2(b)中可以看出,本发明将水域部分从图 像中清晰的分割出来，边缘完整准确，区域一致性好；

图2(c)为现有技术中的An Evolutionary Approach to Multiobjective Clustering的 分割结果图，从图2(c)中可以看出，该技术将中间部分的一块陆地错分成了水域；

图2(d)为现有技术中的遗传算法聚类技术的分割结果，从图2(d)中可以看出， 该技术将左边的一些陆地错分成了水域，同时在中间偏右的陆地区域有严重的错分；

图2(e)为现有技术中的K均值方法的分割结果图，从图2(e)可以看出，图2(e) 相比图2(d)效果要好一些，但也把右上方、左下方、中间部分的陆地错分成了水域。

仿真二：用本发明与现有的方法对一副具有1米分辨率的Ku波段的SAR图像数 据进行测试分割，仿真结果如图3，对其仿真效果进行分析，其中：

图3(a)为原始测试SAR图像，此SAR图像有三种地物类别，即：河流，植被和 农作物,图像大小为256×256；

图3(b)为本发明的分割结果图，从图3(b)中可以看出，利用本发明，水域的堤 岸分的更清楚，比图3(c)产生了更多的均匀的植被区域，同时本发明相比于其它三 个方法，对三类地物的分割结果更接近真实分布；

图3(c)为现有技术中的An Evolutionary Approach to Multiobjective  Clustering(MOCK)的分割结果，从图3(c)中可以看出，该技术在水域的分割一致性上 有所提高，同时农作物的分割结果也有很大的改进，然而一些小的植被区域没有检测 出来；

图3(d)和图3(e)为现有技术中的遗传算法聚类技术的分割结果图和K均值方法 的分割结果图，从图3(d)和图3(e)中可以看出，两种方法对水域和植被的边缘定位 的均不准确，同时还将水域附近的植被错分为了农作物，二者均没有找到水上的桥梁 同时把右上方的大片农作物错分成了水域。