CN103593855B - 基于粒子群优化和空间距离测度聚类的图像分割方法 - Google Patents

Abstract

本发明公开了一种基于粒子群优化和空间距离测度聚类的图像分割处理方法，主要解决现有聚类图像分割技术存在局部错分现象和区域杂点多的问题。其实现步骤为：(1)输入原始图像，提取像素特征，并进行分水岭分割；(2)根据分割后的区域计算邻接矩阵和产生聚类数据；(3)利用聚类数据随机初始化种群；(4)计算种群的隶属度矩阵和适应度值，升级个体最优和全局最优，并进化种群，(5)更新迭代次数，若达到预先设定的最大迭代次数，则输出最佳隶属度矩阵，否则继续执行步骤(4)；(6)根据最佳隶属度矩阵，按照最大概率原则进行标记，得到分割结果。本发明与现有技术相比，区域一致性好，分割正确率高，可用于SAR图像的目标识别。

Description

基于粒子群优化和空间距离测度聚类的图像分割方法

技术领域

本发明属于图像处理技术领域，特别是一种涉及纹理图像和SAR图像的分割方法，可应用于目标识别。

背景技术

图像分割是图像处理的关键技术之一。图像分割的结果就是将图像分成若干个部分，每部分代表图像中不一样的特征，并把同一部分像素标记为同一个值。现有的图像分割方法主要有基于区域的方法，基于边缘检测的方法，基于聚类的方法等。目前，人们更多采用基于聚类分析的方法来进行图像分割。用基于聚类分析的方法分割图像是利用已知的训练样本集在图像的特征空间找到决策分类点、线或面，然后将它们映射回原图像空间，以实现对图像的划分，从而达到图像分割的目的。

聚类，是一种重要的数据分析方法，是无监督分割。聚类分析把在某方面具有相同属性的一组模型聚成一类，在现有的聚类方法中，模糊C均值聚类是最重要的聚类算法之一，但是，众所周知，它存在对初始点的选取十分敏感，而且只能得到局部最优解的缺点。

为了克服上述缺点，很多学者将其他技术与该聚类算法相结合，得到了有效的聚类算法，其中进化算法就是有效的结合技术，主要包括遗传算法，粒子群优化算法等。其中，粒子群优化算法是一种崭新的随机优化方法，它具有程序实现简单、控制参数少的特点，所以得到了广泛应用。在进化算法与聚类算法结合技术中，为了降低计算复杂度，往往先对图像进行预处理，即采用分水岭算法对其预分割，将图像分割成互不重叠的不规则的小块。由此会产生过分割，而后为了降低过分割，可以对其进行优化，并将该优化过程转化为聚类问题。但是在该聚类过程中，由于只考虑区域的属性，如灰度，纹理等，而忽略了区域与区域之间的关系，缺乏空间信息的完整性，造成区域内杂点较多，区域一致性较差，分割不理想的不足。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种基于粒子群优化和空间距离测度聚类的图像分割方法，以保证图像分割信息的完整性，保持分割后的区域一致性，提高图像分割的质量。

实现本发明目的的技术方案是：将待分割图像分成互不重叠的区域，利用图像区域的特征信息进行聚类，从而达到分割的目的。在聚类过程中，不仅考虑图像区域纹理特征的影响，而且考虑区域空间邻域信息的影响。其实现步骤如下：

(1)输入待分割图像，提取该图像的特征，并计算该图像的梯度，得到梯度图像，对梯度图像进行分水岭分割，得到N个互不重叠的区域，N≥1000；

(2)在N个互不重叠的区域中，分别对每两个区域进行膨胀扩充，获得每个区域与其他各区域的邻接关系：若两个区域相交则它们的邻接关系记为1，否则，邻接关系记为0，根据这些邻接关系形成区域邻接矩阵；

(3)对每个区域中所有像素点特征取均值，获得每一个区域的特征向量，作为初始聚类数据点集合Z＝{z₁,z₂…z_N}；

(4)利用初始聚类数据点集合，随机初始化大小为M的种群:

(4a)随机初始化各粒子的位置X＝{x₁,x₂…x_M}、速度V＝{v₁,v₂…v_M}和隶属度矩阵U＝[u_ij]，i＝1,2,…c,j＝1,2,…N，u_ij表示第j个聚类点隶属于第i类的隶属度，c为聚类类别数目；

(4b)将各粒子的当前位置作为各粒子的最优位置，从所有粒子的最优位置中随机选择一个位置作为种群的全局最优位置，并初始化各粒子最优位置的适应度值和全局最优位置的适应度值；

(5)计算各粒子中聚类点集合分别到各聚类中心的距离矩阵：

D＝[d_ij]，i＝1,2,…c，j＝1,2…N，

其中，d_ij表示为第j聚类点到聚类中心i的距离，式中为第j个聚类点到聚类中心i的纹理特征距离，为第j个聚类点到聚类中心i的空间特征距离，α为权重参数，取值为0.3；

(6)利用第j个聚类点隶属于第i类的隶属度计算各粒子的隶属度矩阵U＝[u_ij]，式中m为模糊系数，取值为2。

(7)根据各粒子的隶属度矩阵，计算种群中的各粒子的适应度值式中，u_ij表示第j个聚类点隶属于第i类的隶属度，d_ij表示为第j个聚类点到聚类中心i的距离。

(8)升级每个粒子的最优位置和种群全局最优位置：将粒子的适应度值与该粒子最优位置的适应度值比较，如果小于最优位置的适应度值，则用粒子的位置代替该粒子最优位置，否则粒子最优位置不变；将粒子最优位置的适应度值与全局最优位置的适应度值比较，如果小于全局最优位置的适应度值，则用粒子最优位置代替全局最优位置，否则全局最优位置不变。

(9)更新各粒子的速度和位置，产生下一代种群；

(10)更新迭代次数t，直到迭代次数达到预先设定的最大值max gen，则输出最佳隶属度矩阵，否则返回到(5)进行下一代迭代；

(11)根据最佳隶属度矩阵计算各个区域邻域内的区域属于不同类的概率，按照最大概率原则依次对每个区域进行标记类别号；

(12)对不同类别号的像素点赋不同的灰度值，得到分割后的图像，并输出。

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

1、本发明由于考虑了区域的邻域空间信息，提高了区域一致性，保证了图像信息的完整性；

2、本发明将最大隶属度原则改为有邻域信息的最大概率原则进行区域块的标记进一步提高了分割精确度。

附图说明

图1是本发明的实现流程图；

图2是用本发明与现有方法在一幅两类纹理图像上的分割结果对比图；

图3是用本发明与现有方法在一幅四类纹理图像上的分割结果对比图；

图4是用本发明与现有方法在一幅分为田地和森林的机载SAR图像上的分割结果对比图；

图5是用本发明与现有方法在一幅分为山体和平原的机载SAR图像上的分割结果对比图。

具体实施方式

参照图1，本发明的具体实现步骤如下：

步骤一、输入待分割图像，提取该图像的特征。

(1a)对于任意像素点i′，利用小波分解，提取图像的10维小波特征向量；

(1b)对于任意像素点i′，计算0°,45°,90°,135°四个方向上的灰度共生矩阵，选取该四个矩阵上的三个统计量，即对比度、同质性和角二阶，获得像素点i的12维纹理特征向量；

(1c)将上述10维小波特征向量和12维纹理特征向量合并成22维特征向量，作为第i′个像素点的纹理特征；

(1d)对图像中的所有像素点重复步骤(1a)-(1c),得到原始图像所有像素点的特征。

步骤二、计算待分割图像的梯度。

对待分割图像分别进行膨胀变换和腐蚀变换，再用膨胀变换后的结果减去其腐蚀变换后的结果，由此得到待分割图像形态梯度图像。

步骤三、对形态梯度图像进行分水岭分割，得到N个互不重叠的区域，N≥1000。

步骤四、在N个互不重叠的区域中，分别对每两个区域进行膨胀扩充，膨胀扩充后，若两个区域相交，则将它们的邻接关系记为1，否则，将邻接关系记为0，由此获得每个区域与其他各区域的邻接关系，根据这些邻接关系得到区域邻接矩阵。

步骤五、对每个区域中所有像素点特征取均值，获得每一个区域的特征向量，作为初始聚类数据点集合：Z＝{z₁,z₂…z_N}。

步骤六、利用初始聚类数据点集合Z，随机初始化大小为M的种群。

(6a)从数据点集合Z中随机选取c个数据点作为种群中当前粒子的位置，这c个数据点中的每个数据点代表一个聚类中心，由此得到种群中各粒子的当前位置X＝{x₁,x₂…x_M}，其中c为聚类类别数目；

(6b)随机初始化种群中速度V＝{v₁,v₂…v_M}和隶属度矩阵U＝[u_ij]，i＝1,2,…c,j＝1,2,…N，u_ij表示第j个聚类点隶属于第i类的隶属度，c为聚类类别数目。

步骤七、将各粒子的当前位置作为各粒子的最优位置P＝[p_r]，r＝1,2…M，其中，p_r表示为第r个粒子的最优位置，再从所有粒子的最优位置P中随机选择一个位置作为种群的全局最优位置p_g,并初始化各粒子最优位置的适应度值Pbest＝[pbest_r]，r＝1,2…M为无穷大，全局最优位置的适应度值gbest为无穷大，其中pbest_r表示为第r个粒子最优位置p_r的适应度值。

步骤八、计算第j个聚类点到聚类中心i的距离d_ij：

(8a)计算第j个聚类点到聚类中心i的纹理特征距离

$d_{ij}^F=\left|\right|j-i|{|}_2,$

其中，||·||₂为2范数；

(8b)根据步骤四得到的邻接矩阵，选取第j行所有为1的点作为第j个聚类点的邻域；

(8c)计算第j个聚类点到聚类中心i的空间特征距离

$d_{ij}^S=1-\frac{\underset{k\∈{\mathrm{\η}}_j}{\Σ}{u}_{ik}{\mathrm{\β}}_{ik}}{\underset{i=1}{\overset{c}{\Σ}}\underset{k\∈{\mathrm{\η}}_j}{\Σ}{u}_{it}{\mathrm{\β}}_{ik}},$

其中，u_ik表示第k个聚类点属于第i类的隶属度，η_j为第j个聚类点的邻域,β_ik是邻域因子， 表示第k个聚类点到聚类中心i的纹理特征距离；

(8d)根据聚类点j到聚类中心i的纹理特征距离和聚类点j到聚类中心i的空间特征距离计算第j个聚类点到聚类中心i的距离d_ij：

$d_{ij}=\mathrm{\α}*d_{ij}^F+(1-\mathrm{\α})*d_{ij}^S,$

式中，α为权重参数，取值为0.3。

步骤九、根据第j个聚类点到聚类中心i的距离d_ij，计算种群中的各粒子的适应度值。

(9a)根据第j个聚类点到聚类中心i的距离d_ij,得到每个粒子中聚类点集到各聚类中心的距离矩阵D＝[d_ij]，i＝1,2,…c，j＝1,2…N。

(9b)根据第j个聚类点到聚类中心i的距离d_ij，计算第j个聚类点隶属于第i类的隶属度式中m为模糊系数，取值为2，由此得到各粒子的隶属度矩阵U＝[u_ij]，i＝1,2,…c，j＝1,2…N。

(9c)根据各粒子的隶属度矩阵U和距离矩阵D，计算种群中的各粒子的适应度值Fitness＝[fitness_r]，r＝1,2…M，其中

步骤十、根据各粒子的适应度值，升级每个粒子的最优位置和种群全局最优位置。

(10a)将各粒子的适应度值Fitness与各自粒子最优位置的适应度值Pbest比较：如果粒子r的适应度值fitness_r小于最优位置的适应度值pbest_r，则用该粒子的位置x_r作为该粒子的最优位置p_r，否则该粒子的最优位置p_r不变；

(10b)将各粒子最优位置的适应度值Pbest与全局最优位置的适应度值gbest比较：如果粒子r最优位置的适应度值pbest_r小于全局最优位置的适应度值gbest，则用该粒子最优位置p_r代替全局最优位置p_g，否则全局最优位置p_g不变。

步骤十一、根据各粒子最优位置P和全局最优位置p_g，利用如下公式更新各粒子的速度和位置：

$v_i^{t+1}=\mathrm{\δ}*\[p_i^t+{\mathrm{\λ}}_1*{\mathrm{rand}}_1\·(p_i^t-x_i^t)+{\mathrm{\λ}}_2*{\mathrm{rand}}_2\·(p_g^t-x_i^t)\],$

$x_i^{t+1}=x_i^t+v_i^t,$

其中，表示t+1代第i个粒子的速度，表示t代第i个粒子的速度，表示t代第i个粒子的位置；表示t+1代第i个粒子的位置，表示t代第i个粒子的最好位置，表示t代全局最优位置，λ₁，λ₂表示两个不同的正的加速常数，rand₁，rand₂均表示[O，1]之间的均匀分布随机数组，δ取值为0.729，·表示点乘运算。

步骤十二、更新迭代次数t，直到迭代次数达到预先设定的最大值maxgen＝100，输出最佳隶属度矩阵，执行步骤十三，否则返回到步骤八进行下一代迭代。

步骤十三、根据最佳隶属度矩阵计算各个区域的邻域属于不同类的概率，按照最大概率原则依次对每个区域进行标记，得到分割结果。

(13a)计算各个区域的邻域属于不同类的概率，其计算公式如下：

$p_{ij}=\frac{N_j^i}{\underset{i=1}{\overset{c}{\Σ}}{N}_j^i},$

其中，p_ij表示为区域j的邻域属于第i类的概率，为区域j的邻域按照最大隶属度原则属于第i类的个数；

(13b)在每个区域j中，根据区域j的邻域属于第i类的概率p_ij值的大小，选择最大的概率；利用最大概率对每个区域j进行类别标号，然后由每个区域的类别标号得到每个像素点的类别标号；对不同的类别标号赋于不同的灰度值，得到最终的分割结果。

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

1.仿真条件：本发明的仿真在windows XP，SPI，CPU Inter Core 2 Duo，基本频率2.33Ghz，软件平台为Matlab2007b运行。

2.仿真内容：应用本发明方法和模糊C均值方法，粒子群聚类方法分别对两幅合成纹理图像和两幅SAR图像进行分割实验，并从细节信息的完整性、同质区域的内部一致性和边缘边界的清晰性评价这些方法的各自性能。合成纹理图像还可以用平均分割正确率结果来评价各种方法的性能。

仿真1，用本发明方法以及模糊C均值方法，粒子群聚类方法分别对两类合成纹理图像image1进行分割，其结果如图2所示。其中图2(a)为两类合成纹理图像image1；图2(b)为图像image1的理想分割结果；图2(c)为用现有模糊C均值方法对图2(a)进行分割得到的结果；图2(d)为用现有粒子群聚类方法对图2(a)进行分割得到的结果；图2(e)为用本发明方法对图2(a)进行分割得到的结果。

仿真2，用本发明方法以及模糊C均值方法，粒子群聚类方法分别对,四类合成纹理图像image2进行分割，其结果如图3所示。其中图3(a)为三类合成纹理图像image2；图3(b)为对图像image2的理想分割结果；图3(c)为用现有模糊C均值方法对图3(a)进行分割得到的结果；图3(d)为用现有粒子群聚类方法对图3(a)进行分割得到的结果；图3(e)为用本发明方法对图3(a)进行分割得到的结果。

仿真3，用本发明方法以及模糊C均值方法，粒子群聚类方法分别对机载SAR图像，其结果如图4所示。其中图4(a)为两类机载SAR原图像，分为田地和森林；图4(b)为用现有模糊C均值方法对图4(a)进行分割得到的结果；图4(c)为用现有粒子群聚类方法对图4(a)进行分割得到的结果；图4(d)为用本发明方法对图4(a)进行分割得到的结果。

仿真4，用本发明方法以及模糊C均值方法，粒子群聚类方法分别对机载SAR图像，其结果如图5所示。其中图5(a)为两类机载SAR图像，分为山体和平原；图5(b)为用现有模糊C均值方法对图5(a)进行分割得到的结果；图5(c)为用现有粒子群聚类方法对图5(a)进行分割得到的结果；图5(d)为用本发明方法对图5(a)进行分割得到的结果。

3.仿真结果分析：

从图2(c)和图3(c)的分割结果中可以看出，模糊C均值聚类方法存在一定的杂点，且分割的比较凌乱。

从图2(d)和3(d)的分割结果中可以看出，粒子群聚类算法比模糊C均值边缘分得较好一些，但是由于没有考虑区域之间的空间信息还是存在一定的杂点。

从图2(e)和3(e)的分割结果中可以看出,本发明方法提高了区域一致性，有比较好的分割结果。

从图4(b)，5(b)的分割结果中可以看出，模糊C均值聚类方法存较多的杂点。

从图4(c)，5(c)的分割结果中可以看出,粒子群聚类算法由于没有考虑区域之间的空间信息还是存在一定的杂点。

从图4(d)，5(d)的分割结果中可以看出,本发明方法去除了杂点,保持了较高的区域一致性，有比较好的分割结果。

表1给出了不同的方法对Image1和Image2的独立分割运行20次后的平均分割正确率。

表1 不同的方法对Image1和Image2的分割正确率

图像	模糊C均值	粒子群聚类	本发明
				Image1	0.9356	0.9413	0.9641
Image2	0.9591	0.9607	0.9710

从表1可见，本发明与其他两种算法的对比，得到了最高的平均正确率。

Claims (6)

1.一种基于粒子群优化和空间距离测度聚类的图像分割处理方法，包括如下步骤：

(1)输入待分割图像，提取该图像的特征，并计算该图像的梯度，得到梯度图像，对梯度图像进行分水岭分割，得到N个互不重叠的区域，N≥1000；

(2)在N个互不重叠的区域中，分别对每两个区域进行膨胀扩充，获得每个区域与其他各区域的邻接关系：若两个区域相交则它们的邻接关系记为1，否则，邻接关系记为0，根据这些邻接关系形成区域邻接矩阵；

(3)对每个区域中所有像素点特征取均值，获得每一个区域的特征向量，作为初始聚类数据点集合Z＝{z₁,z₂…z_N}；

(4)利用初始聚类数据点集合，随机初始化大小为M的种群:

(4a)随机初始化各粒子的位置X＝{x₁,x₂…x_M}、速度V＝{v₁,v₂…v_M}和隶属度矩阵U＝[u_ij]，i＝1,2,…c,j＝1,2,…N，u_ij表示第j个聚类点隶属于第i类的隶属度，c为聚类类别数目；

(4b)将各粒子的当前位置作为各粒子的最优位置，从所有粒子的最优位置中随机选择一个位置作为种群的全局最优位置，并初始化各粒子最优位置的适应度值和全局最优位置的适应度值；

(5)计算各粒子中聚类点集合分别到各聚类中心的距离矩阵：

D＝[d_ij]，i＝1,2,…c，j＝1,2…N，

其中，d_ij表示为第j聚类点到聚类中心i的距离，式中为第j个聚类点到聚类中心i的纹理特征距离，为第j个聚类点到聚类中心i的空间特征距离，α为权重参数，取值为0.3；

(6)利用第j个聚类点隶属于第i类的隶属度计算各粒子的隶属度矩阵U＝[u_ij]，式中m为模糊系数，取值为2；

(7)根据各粒子的隶属度矩阵，计算种群中的各粒子的适应度值式中，u_ij表示第j个聚类点隶属于第i类的隶属度，d_ij表示为第j个聚类点到聚类中心i的距离；

(8)升级每个粒子的最优位置和种群全局最优位置：将粒子的适应度值与该粒子最优位置的适应度值比较，如果小于最优位置的适应度值，则用粒子的位置代替该粒子最优位置，否则粒子最优位置不变；将粒子最优位置的适应度值与全局最优位置的适应度值比较，如果小于全局最优位置的适应度值，则用粒子最优位置代替全局最优位置，否则全局最优位置不变；

(9)更新各粒子的速度和位置，产生下一代种群；

(10)更新迭代次数t，直到迭代次数达到预先设定的最大值maxgen，则输出最佳隶属度矩阵，否则返回到(5)进行下一代迭代；

(11)根据最佳隶属度矩阵计算各个区域的邻域属于不同类的概率，按照最大概率原则依次对每个区域进行标记类别号；

(12)对不同类别号的像素点赋不同的灰度值，得到分割后的图像，并输出。

2.根据权利要求书1所述的基于粒子群优化和空间距离测度聚类的图像分割处理方法，其中步骤(1)中提取待分割图像的特征，按照以下步骤进行：

(2a)对于任意像素点i′，利用小波分解，提取图像的10维小波特征；

(2b)对于任意像素点i′，计算0°,45°,90°,135°四个方向上的灰度共生矩阵，选取该四个矩阵上的三个统计量，分别为对比度、同质性和角二阶，获得像素点i′的12维纹理特征；

(2c)将上述10维小波特征向量和12维纹理特征向量合并成22维纹理特征向量，作为第i′个像素点的特征；

(2d)对图像中的所有像素点重复步骤(2a)-(2c)，得到原始图像所有像素点的纹理特征。

3.根据权利要求书1所述的基于粒子群优化和空间距离测度聚类的图像分割处理方法，其中所述步骤(5)中第j个聚类点到聚类中心i的纹理特征距离通过如下公式计算：其中，||·||₂为2范数。

4.根据权利要求书1所述的基于粒子群优化和空间距离测度聚类的图像分割处理方法，其中所述步骤(5)中第j个聚类点到聚类中心i的空间特征距离通过如下公式计算：

$d_{ij}^S=1-\frac{\underset{k\∈{\mathrm{\η}}_j}{\Σ}{u}_{ik}{\mathrm{\β}}_{ik}}{\underset{i=1}{\overset{c}{\Σ}}\underset{k\∈{\mathrm{\η}}_j}{\Σ}{u}_{it}{\mathrm{\β}}_{ik}},$

其中，u_ik表示第k个聚类点属于第i类的隶属度，η_j为第j个聚类点的邻域,β_ik是邻域因子， 表示第j个聚类点到聚类中心i的纹理特征距离，表示第k个聚类点到聚类中心i的纹理特征距离，θ为参数，取值为0.7。

5.根据权利要求书1所述的基于粒子群优化和空间距离测度聚类的图像分割处理方法，其中步骤(9)所述的更新各粒子的速度和位置，按如下公式进行：

$v_i^{t+1}=\mathrm{\δ}*\[p_i^t+{\mathrm{\λ}}_1*{\mathrm{rand}}_1\·(p_i^t-x_i^t)+{\mathrm{\λ}}_2*{\mathrm{rand}}_2\·(p_g^t-x_i^t)\]$

$x_i^{t+1}=x_i^t+v_i^t$

其中，表示t+1代第i个粒子的速度，表示t代第i个粒子的速度，表示t代第i个粒子的位置，表示t+1代第i个粒子的位置，表示t代第i个粒子的最好位置，表示t代全局最优位置，λ₁，λ₂表示两个不同的正的加速常数，rand₁，rand₂均表示[O，1]之间的均匀分布随机数组，δ取值为0.729，·表示点乘运算。

6.根据权利要求书1所述的基于粒子群优化和空间距离测度聚类的图像分割处理方法，其中步骤(11)所述的根据最佳隶属度矩阵计算各个区域的邻域属于不同类的概率，按如下公式计算：

$p_{ij}=\frac{N_j^i}{\underset{i=1}{\overset{c}{\Σ}}{N}_j^i},$

其中，p_ij表示区域j的邻域属于第i类的概率，为区域j的邻域按照最大隶属度原则属于第i类的个数。