CN103745482A - 一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法 - Google Patents

Abstract

本发明提出了一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法，蝙蝠算法作为一种性能优良的演化算法，它既降低了传统单阈值图像分割方法对初始参数的敏感度；同时，蝙蝠算法作为一种先进的演化算法，可以直接采用十进制进行编码，相对于大部分演化算法采用二进制编码而言，计算过程更为方便，大大提高整个图像分割过程的效率，快速得到最优的图像分割效果。

Description

一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法

技术领域

本发明涉及图像处理技术领域，尤其涉及一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法。

背景技术

公知，图像分割就是把图像分成若干个特定的、具有独特性质的区域并提出感兴趣目标的技术和过程。它是由图像处理到图像分析的关键步骤。图像分割的任务是将输入图像分割为一些独立的区域，使同一区域具有相同的属性，而使不同的区域具有不同的属性。对于图像分割问题，现在已有大量的方法，但鉴于图像种类繁多、所包含的数据量大、细节变化多端的特点，迄今为止还没有一种图像分割方法适合于所有的情况。

阈值分割方法是一种基于区域的图像分割技术，它假定能够根据灰度值将不同的目标或者目标与背景区分开，通过设定不同的阈值，把图像像素点分为若干类。其中，合理的灰度级门限的选取是阈值法所要解决的关键问题。虽然阈值分割方法具有简单快速的优点，但是有两个主要的缺点限制其性能：

其一是传统的单阈值分割方法过于依赖初始参数。现已引进了多阈值的方法来解决图像分割问题[1]，多阈值的方法是通过适当增加参数的个数，以得到更加合理的灰度级门限。参数个数越多，门限值越合理，但整个算法的时间复杂度也会随之上升。

其二是传统的多阈值分割方法时间复杂度很高。由于传统的多阈值分割方法是一个不断修正的过程，当超过两个阈值时，最优阈值可能的组合次数呈灰度级的指数上升。已有研究者引进了差分演化算法到多阈值图像分割中来解决此问题[2]。演化算法作为一种并行搜索过程，只要给定算法的初值，它就会依据适应度函数不断进行搜索并找到最优解，而不需要很多的迭代过程，也不用每次对隶属度函数进行计算。同时，蝙蝠算法作为一种先进的演化算法，可以直接采用十进制进行编码，相对于大部分演化算法采用二进制编码而言，计算过程更为方便，能在很大程度上提高图像分割的效率。

参考文献：

[1].Wang Lei,Duan Huichuan,Application of Otsu'method in multi-threshold imagesegmentation,Computer Engineering and Design,2008,29（11）,P:2844-2846.

[2].Erik Cuevas,Daniel Zaldivar,Marco Pérez-Cisneros,A novel multi-threshold segmentationapproach based on differential evolution optimization,Expert Systems with Applications,vol：37,Issue：7，Issue date：2010,P:5265–5271.

[3].Tao Wenbing,Jin Hai,Liu Liman,Object segmentation using ant colony optimizationalgorithm and fuzzy entropy,Pattern Recognition Letters,2007,28,P:788-796.

[4].Ming-Huwi Horng,Multilevel thresholding selection based on the artificial bee colonyalgorithm for image segmentation,Expert Systems with Applications,2011,38,P:13785-13791.

[5].Xin-She Yang,Bat algorithm for multi-objective optimisation,International Journal ofBio-Inspired Computation,2011,3，P:267-274.

[6].Valentin Osuna-Enciso,Erik Cuevas,Hunberto Sossa,A comparison of nature inspiredalgorithms for multi-threshold image segmentation,Expert Systems with Applications,2011,40,P:1213-1219.

[7].张小红,宁红梅,基于混沌粒子群和模糊聚类的图像分割算法,计算机应用研究,2011,28,P:4786-4789.

发明内容

本发明为了解决上述的技术问题，提出了一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法。

本发明的技术方案是：一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法，包括下述步骤：

步骤1，输入待分割图像，种群规模S以及停机条件；

步骤2，依据设定的种群规模S，运行蝙蝠算法，初始种群的各变量，产生初始种群G_k，对于第i个个体而言，其位置向量x_i（这里i值范围是1≤i≤S，S表示种群规模），其取值范围为（0，255），速度v_i、声波频率f_i、脉冲响度A_i、发射速率r_i五个变量随机产生，代数n表示蝙蝠算法的运行代数，在初始条件下，令n=0；

步骤3，对不满足大小关系的个体进行扰动；

每个个体需满足0≤a₁＜b₁＜c₁＜a₂＜b₂＜c₂≤255的大小关系，各参数的取值范围均为（0，255），在初始化的过程中需要对不满足大小关系的个体进行如下的顺序扰动：

3.1.对a₁产生扰动，得到a₁ ¹，如果a₁ ¹＜0,则令a₁ ¹＝0，如果a₁ ¹＞250，则令a₁ ¹＝250；

3.2.对b₁产生扰动，得到b₁ ¹，如果b₁ ¹＜a₁ ¹+1,则令b₁ ¹＝a₁ ¹+1，如果b₁ ¹＞251，则令b₁ ¹＝251；

3.3.对c₁产生扰动，得到c₁ ¹，如果c₁ ¹＜b₁ ¹+1,则令c₁ ¹＝b₁ ¹+1，如果c₁ ¹＞252，则令c₁ ¹＝252；

3.4.对a₂产生扰动，得到a₂ ¹，如果a₂ ¹＜c₁ ¹+1,则令a₂ ¹＝c₁ ¹+1，如果a₂ ¹＞253，则令a₂ ¹＝253；

3.5.对b₂产生扰动，得到b₂ ¹，如果b₂ ¹＜a₂ ¹+1,则令b₂ ¹＝a₂ ¹+1，如果b₂ ¹＞254，则令b₂ ¹＝254；

3.6.对c₂产生扰动，得到c₂ ¹，如果c₂ ¹＜b₂ ¹+1,则令c₂ ¹＝b₂ ¹+1，如果c₂ ¹＞255，则令c₂ ¹＝255；

步骤4，计算灰度值为k的像素属于暗、灰、亮的隶属度值和条件概率；

${\mathrm{\&mu;}}_d\left(k\right)=\left\{\begin{array}{cc}1& k\&le;a_1\\ 1-\frac{{(k-a_1)}^2}{(c_1-a_1)*(b_1-a_1)}& a_1<k\&le;b_1\\ \frac{{(k-c_1)}^2}{(c_1-a_1)*(c_1-b_1)}& b_1<k\&le;c_1\\ 0& k>c_1\end{array}\right.$

${\mathrm{\&mu;}}_m\left(k\right)=\left\{\begin{array}{cc}0& k\&le;a_1\\ \frac{{(k-a_1)}^2}{(c_1-a_1)*(c_1-b_1)}& a_1<k\&le;b_1\\ 1-\frac{{(k-c_1)}^2}{(c_1-a_1)*(c_1-b_1)}& b_1<k\&le;c_1\\ 1& \\ 1-\frac{{(k-a_2)}^2}{(c_2-a_2)*(b_2-a_2)}& a_2<k\&le;b_2\\ \frac{{(k-c_2)}^2}{(c_2-a_2)*(c_2-b_2)}& b_2<k\&le;c_2\\ 0& k>c_2\end{array}\right.$

${\mathrm{\&mu;}}_b\left(k\right)=\left\{\begin{array}{cc}0& k\&le;a_2\\ \frac{{(k-a_2)}^2}{(c_2-a_2)*(c_2-b_2)}& a_2<k\&le;b_2\\ 1-\frac{{(k-c_2)}^2}{(c_2-a_2)*(c_2-b_2)}& b_2<k\&le;c_2\\ 1& k>c_2\end{array}\right.$

$p_d=\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}{p_k*\mathrm{\&mu;}}_d\left(k\right)$

$p_m=\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}{p}_k*{\mathrm{\&mu;}}_m\left(k\right)$

$p_b=\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}{p}_k*{\mathrm{\&mu;}}_b\left(k\right)$

其中，μ_d(k)、μ_m(k)、μ_b(k)分别表示灰度值为k(0≤k≤255)的像素属于暗、灰、亮的隶属度；p_k为图像灰度直方图的概率分布，p_d、p_m、p_b分别表示表示灰度值为k(0≤k≤255)的像素属于暗、灰、亮的条件概率；

步骤5,计算种群中个体的适应度并找出其全局最优解；

根据每个个体对应的像素属于暗、灰、亮隶属度和条件概率，通过模糊熵公式计算每个个体的适应度值，将每个个体适应度值H_s(i)（这里i值范围是1≤i≤S，S表示种群规模）表中最大的值fitness_generation作为种群的全局最优解，并存入最优解集；

步骤6，采用蝙蝠算法对种群中各变量进行更新；

调整声波频率f_i产生新的解并更新速度v_i和位置x_i，更新公式如下：

f_i＝f_min+(f_max-f_min)^*β

v_i ⁿ＝v_i ^n-1+(x_i ⁿ-x_globalbest)^*f_i

x_i ⁿ＝x_i ^n-1+v_i ⁿ

其中，n表示蝙蝠算法的运行代数，f_min和f_max分别表示声波频率的范围，β∈[0,1]是一个随机向量，保证f_i在[f_min,f_max]范围内，x_globalbest为全局最优解的位置向量，速度v_i可为正值或负值，位置向量x_i可任意方向移动；

步骤7，从最优解集中选择一个解，并在该最优解附近形成一个局部解，进而在该局部解附近形成一个新解；

若n=1,局部解取步骤5中全局最优解所对应的位置变量，若n≠1，局部解直接取前一代个体适应度的全局最优解所对应的位置变量x_globalbest，新解的计算公式如下：

x_new＝x_globalbest+ε^*avgA^k

其中，x_new表示得到的新解，ε∈[-1,1]是一个任意的实数，avgA^k是个体在当前代的平均脉冲响度，由脉冲响度A_i计算所得；

步骤8，计算新解中个体的适应度；

步骤9，判断是否对个体参数进行更新；

用新解x_new中个体的适应度用H_{s_new}(i)与H_s(i)进行比较，（若n=1,H_s(i)取步骤5中全局最优解的个体适应度H_s(i)，若n≠1，H_s(i)取算法上一代计算产生的最优解集的个体适应度H_s(i)），若H_s(i)＜H_{s_new}(i)(1≤i≤S)且脉冲响度A_i大于人为设定阈值R_A，则进入步骤10；否，则进入步骤11；

步骤10，对原解x中第i个个体的位置变量用新解x_new中第i个个体的位置变量进行替换并对脉冲响度A_i和发射速率r_i进行更新后，进入步骤11，脉冲响度A_i和发射速率r_i的更新公式如下；

$A_i^{n+1}=\mathrm{\&alpha;}{A}_i^m,r_i^{n+1}=r_i^0[1-\exp (-\mathrm{\&gamma;n})]$

其中n表示蝙蝠算法的运行代数，exp表示以自然对数e为底的指数函数，α和β是两个常量；

步骤11，计算种群全局最优解，通过与先前全局最优解的比较，对全局最优解进行更新；

计算种群中个体的适应度，得到当前的全局最优解fitness_{new_gen}与先前的fitness_generation进行比较，若fitness_{new_gen}＞fitnessgeneration，则将其更新为种群新的全局最优解；否则，将保留原先的全局最优解fitness_generation；

步骤12，判断是否满足停机条件，若是，则进入步骤13，否则返回步骤6；

步骤13，根据最终得到的全局最优解fitnessgeneration对应的6个参数a₁、b₁、c₁、a₂、b₂、c₂，得到最终分割阈值；

阈值t₁和阈值t₂的公式如下：

若(a₁+c₁)/2≤b₁≤c₁,则 $t_1=a_1+\sqrt{(c_1-a_1)*(b_1-a_1)/2}$

若a₁≤b₁＜(a₁+c₁)/2,则 $t_1=c_1-\sqrt{(c_1-a_1)*(c_1-b_1)/2}$

若(a₂+c₂)/2≤b₂≤c₂,则 $t_2=a_2+\sqrt{(c_2-a_2)*(b_2-a_2)/2}$

若a₂≤b₂＜(a₂+c₂)/2,则 $t_2=c_2-\sqrt{(c_2-a_2)*(c_2-b_2)/2}$

使用t₁和t₂对图像进行分割，即将图像中灰度值在[0,t₁)范围内的像素作为第一类，图像中灰度值在[t₁,t₂)范围内的像素作为第二类，图像中灰度值在[t₂,255]范围内的像素作为第三类，并将最终分割图像输出。

而且，步骤1中，群体规模S=50，停机条件为：连续二十代运算得到的全局最优解保持不变。

而且，步骤6中，蝙蝠算法频率范围的取值为f_min＝0、f_max＝10。

而且，步骤9中，阈值R_A=1.5。

而且，步骤10中，α＝β＝0.9。

而且，步骤5、步骤8和步骤11中所述的适应度值计算公式如下：

$H_d=-\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}\frac{p_k*{\mathrm{\&mu;}}_d\left(k\right)}{p_d}*\mathrm{In}\left(\frac{p_k*{\mathrm{\&mu;}}_d\left(k\right)}{p_d}\right)$

$H_m=-\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}\frac{p_k*{\mathrm{\&mu;}}_m\left(k\right)}{p_m}*\mathrm{In}\left(\frac{p_k*{\mathrm{\&mu;}}_m\left(k\right)}{p_m}\right)$

$H_b=-\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}\frac{p_k*{\mathrm{\&mu;}}_b\left(k\right)}{p_b}*\mathrm{In}\left(\frac{p_k*{\mathrm{\&mu;}}_b\left(k\right)}{p_b}\right)$

H＝H_d+H_m+H_b

其中，H为图像的总模糊熵，H_d、H_m、H_b分别表示暗、灰、亮这三个类别中各类的类内模糊熵，μ_d(k)、μ_m(k)、μ_b(k)和p_d、p_m、p_b分别表示灰度值为k(0≤k≤255)的像素属于暗、灰、亮的隶属度值和条件概率，p_k为图像灰度直方图的概率分布,H_s(i)(1≤i≤S)表示每个个体适应度值，也是每个个体的模糊熵值。

附图说明

图1是本发明的流程图

图2是本发明和对比方法在仿真实验中应用的Lena图

图3-1是采用本发明的方法（以下简称BAT-Multithreshold）在Lena图像进行仿真实验时得到的结果图

图3-2是采用基于粒子优化算法的双阈值模糊熵方法（以下简称PSO-Multithreshold）在Lena图像进行仿真实验时得到的结果图

图3-3是采用基于遗传算法的双阈值模糊熵方法（以下简称GA-Multithreshold）方法在Lena图像进行仿真实验时得到的结果图

图4是BAT-Multithreshold、PSO-Multithreshold、GA-Multithreshold三种方法在Lena图像中适应度优化曲线的比较图

图5是本发明和对比方法在仿真实验中应用的Peppers图

图6-1是采用本发明的方法在Peppers图像进行仿真实验时得到的结果图

图6-2是采用PSO-Multithreshold方法在Peppers图像进行仿真实验时得到的结果图

图6-3是采用GA-Multithreshold方法在Peppers图像进行仿真实验时得到的结果图

图7是BAT-Multithreshold、PSO-Multithreshold、GA-Multithreshold三种方法在Peppers图像中适应度优化曲线的比较图

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

参照图1，本发明的流程图，一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法，具体实现过程包括下述步骤：

步骤1，输入待分割图像，种群规模S以及停机条件；

其中，群体规模S是种群中所包含个体的具体数目，在本发明中，取S=50。本发明设定的停机条件为：连续二十代运算得到的全局最优解保持不变。

步骤2，依据设定的种群规模S，运行蝙蝠算法，初始种群的各变量，

依据设定的种群规模S，产生初始种群G_k。在本方法中，由于种群表示的是数字图像的灰度值，而灰度值的取值范围是（0，255），故对于种群中任意一个个体中的一位，其取值范围只能取（0，255），同时对蝙蝠算法中的其它变量进行初始化。

对于蝙蝠算法而言，需要用到位置、速度、声波频率、脉冲响度、发射速率这五个参数。其中，对于第i个个体而言，其位置向量x_i（1≤i≤S），其中S表示种群规模，在本发明中表示数字图像的灰度值，其取值范围是(0,255)，速度v_i、声波频率f_i、脉冲响度A_i、发射速率r_i这4个参数依据种群规模S随机产生。最后，用代数n表示蝙蝠算法的运行代数。在初始条件下，令n=0。

步骤3，对不满足大小关系的个体进行扰动；

根据模糊熵图像分割方法，每3个参数确定一个阈值，为了确定两个阈值，种群中每个个体由6个参数a₁、b₁、c₁、a₂、b₂、c₂组成，每个参数的取值范围均为（0，255），且每个个体必须满足0≤a₁＜b₁＜c₁＜a₂＜b₂＜c₂≤255的大小关系，在初始化的过程中需要对不满足大小关系的个体进行如下的顺序扰动：

对a₁产生扰动，得到a₁ ¹，如果a₁ ¹＜0,则令a₁ ¹＝0，如果a₁ ¹＞250，则令a₁ ¹＝250；

对b₁产生扰动，得到b₁ ¹，如果b₁ ¹＜a₁ ¹+1,则令b₁ ¹＝a₁ ¹+1，如果b₁ ¹＞251，则令b₁ ¹＝251；

对c₁产生扰动，得到c₁ ¹，如果c₁ ¹＜b₁ ¹+1,则令c₁ ¹＝b₁ ¹+1，如果c₁ ¹＞252，则令c₁ ¹＝252；

对a₂产生扰动，得到a₂ ¹，如果a₂ ¹＜c₁ ¹+1,则令a₂ ¹＝c₁ ¹+1，如果a₂ ¹＞253，则令a₂ ¹＝253；

对b₂产生扰动，得到b₂ ¹，如果b₂ ¹＜a₂ ¹+1,则令b₂ ¹＝a₂ ¹+1，如果b₂ ¹＞254，则令b₂ ¹＝254；

对c₂产生扰动，得到c₂ ¹，如果c₂ ¹＜b₂ ¹+1,则令c₂ ¹＝b₂ ¹+1，如果c₂ ¹＞255，则令c₂ ¹＝255；

步骤4，计算灰度值为k的像素属于暗、灰、亮的隶属度值和条件概率；

${\mathrm{\&mu;}}_d\left(k\right)=\left\{\begin{array}{cc}1& k\&le;a_1\\ 1-\frac{{(k-a_1)}^2}{(c_1-a_1)*(b_1-a_1)}& a_1<k\&le;b_1\\ \frac{{(k-c_1)}^2}{(c_1-a_1)*(c_1-b_1)}& b_1<k\&le;c_1\\ 0& k>c_1\end{array}\right.$

${\mathrm{\&mu;}}_m\left(k\right)=\left\{\begin{array}{cc}0& k\&le;a_1\\ \frac{{(k-a_1)}^2}{(c_1-a_1)*(c_1-b_1)}& a_1<k\&le;b_1\\ 1-\frac{{(k-c_1)}^2}{(c_1-a_1)*(c_1-b_1)}& b_1<k\&le;c_1\\ 1& \\ 1-\frac{{(k-a_2)}^2}{(c_2-a_2)*(b_2-a_2)}& a_2<k\&le;b_2\\ \frac{{(k-c_2)}^2}{(c_2-a_2)*(c_2-b_2)}& b_2<k\&le;c_2\\ 0& k>c_2\end{array}\right.$

${\mathrm{\&mu;}}_b\left(k\right)=\left\{\begin{array}{cc}0& k\&le;a_2\\ \frac{{(k-a_2)}^2}{(c_2-a_2)*(c_2-b_2)}& a_2<k\&le;b_2\\ 1-\frac{{(k-c_2)}^2}{(c_2-a_2)*(c_2-b_2)}& b_2<k\&le;c_2\\ 1& k>c_2\end{array}\right.$

$p_d=\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}{p_k*\mathrm{\&mu;}}_d\left(k\right)$

$p_m=\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}{p}_k*{\mathrm{\&mu;}}_m\left(k\right)$

$p_b=\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}{p}_k*{\mathrm{\&mu;}}_b\left(k\right)$

其中，μ_d(k)、μ_m(k)、μ_b(k)分别表示灰度值为k(0≤k≤255)的像素属于暗、灰、亮的隶属度；p_k为图像灰度直方图的概率分布，p_d、p_m、p_b分别表示表示灰度值为k(0≤k≤255)的像素属于暗、灰、亮的条件概率；

步骤5,计算种群中个体的适应度并找出其全局最优解；

当得到每个个体对应的像素属于暗、灰、亮隶属度和条件概率后，即可通过模糊熵公式计算每个个体的适应度值。模糊熵计算公式如下式：

$H_d=-\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}\frac{p_k*{\mathrm{\&mu;}}_d\left(k\right)}{p_d}*\mathrm{In}\left(\frac{p_k*{\mathrm{\&mu;}}_d\left(k\right)}{p_d}\right)$

$H_m=-\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}\frac{p_k*{\mathrm{\&mu;}}_m\left(k\right)}{p_m}*\mathrm{In}\left(\frac{p_k*{\mathrm{\&mu;}}_m\left(k\right)}{p_m}\right)$

$H_b=-\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}\frac{p_k*{\mathrm{\&mu;}}_b\left(k\right)}{p_b}*\mathrm{In}\left(\frac{p_k*{\mathrm{\&mu;}}_b\left(k\right)}{p_b}\right)$

H＝H_d+H_m+H_b

其中，H为图像的总模糊熵，H_d、H_m、H_b分别表示暗、灰、亮这三个类别中各类的类内模糊熵，μ_d(k)、μ_m(k)、μ_b(k)和p_d、p_m、p_b分别表示灰度值为k(0≤k≤255)的像素属于暗、灰、亮的隶属度值和条件概率，p_k为图像灰度直方图的概率分布,H_s(i)(1≤i≤S)（这里i值范围是1≤i≤S，其中S表示种群规模）表示每个个体适应度值，也是每个个体的模糊熵值，将H_s(i)中最大的值fitness_generation作为种群的全局最优解，并存入最优解集；

步骤6，采用蝙蝠算法对种群中各变量进行更新；

调整声波频率f_i产生新的解并更新速度v_i和位置x_i，每产生一组新解，就需要对表示该组解的各个变量进行更新，声波频率f_i、速度v_i和位置x_i的更新公式如下：

f_i＝f_min+(f_max-f_min)^*β

v_i ⁿ＝v_i ^n-1+(x_i ⁿ-x_globalbest)^*f_i

x_i ⁿ＝x_i ^n-1+v_i ⁿ

其中，n表示蝙蝠算法的运行代数，v_i ⁿ、x_i ⁿ分别表示蝙蝠算法的速度变量和位置变量，f_min和f_max分别表示声波蝙蝠算法频率的范围，频率范围的取值为f_min＝0、f_max＝10，β∈[0,1]是一个随机向量，保证f_i在[f_min,f_max]范围内，x_globalbest为全局最优解的位置向量，速度v_i可为正值或负值，位置向量x_i可任意方向移动；

步骤7，从最优解集中选择一个解，并在该最优解附近形成一个局部解，进而在该局部解附近形成一个新解；

若n=1,局部解取步骤5中全局最优解所对应的位置变量。若n≠1，局部解直接取前一代个体适应度的全局最优解所对应的位置变量x_globalbest，新解的计算公式如下：

x_new＝x_globalbest+ε^*avgA^k

其中，x_new表示得到的新解，ε∈[-1,1]是一个任意的实数，avgA^k是个体在当前代的平均脉冲响度，由脉冲响度A_i计算所得；

步骤8，计算新解中个体的适应度，新解x_new中个体的适应度为J_{q_new}(i_b)；

步骤9，判断是否对个体参数进行更新；

用新解x_new中个体的适应度用H_{s_new}(i)与H_s(i)进行比较，（若n=1,H_s(i)取步骤5中全局最优解的个体适应度H_s(i)。若n≠1，H_s(i)取算法上一代计算产生的最优解集的个体适应度H_s(i)），若H_s(i)＜H_{s_new}(i)(1≤i≤S)且脉冲响度A_i大于人为设定阈值R_A，则进入步骤10；否，则进入步骤11；

步骤10，对原解x中第i个个体的位置变量用新解x_new中第i个个体的位置变量进行替换并对脉冲响度A_i和发射速率r_i进行更新后，进入步骤11，脉冲响度A_i和发射速率r_i的更新公式如下；

$A_i^{n+1}=\mathrm{\&alpha;}{A}_i^m,r_i^{n+1}=r_i^0[1-\exp (-\mathrm{\&gamma;n})]$

其中n表示蝙蝠算法的运行代数，exp表示以自然对数e为底的指数函数，α和β是两个常量；在本发明中，取α＝β＝0.9。

步骤11，计算种群全局最优解，通过与先前全局最优解的比较，对全局最优解进行更新；

计算种群中个体的适应度，得到当前全局最优解fitness_{new_gen}与先前的fitness_generation进行比较，若fitness_{new_gen}＞fitness_generation，则将其更新为种群新的全局最优解；否则，将保留原先的全局最优解fitness_generation；

步骤12，判断是否满足停机条件，若是，则进入步骤13，否则返回步骤6；

若不满足停机条件，则返回步骤6，直到算法运行终止，得到全局最优解。停机条件按如下方式设置：设置一个最优解集和一个计数变量，将每次得到的全局最优解fitness_generation存入到该最优解集中。在程序运行中，如果当前代全局最优解与上一代全局最优解保持不变的话，计数变量加1。根据多次实验运行结果，计数变量为20时算法具有最好的搜索性能，因此在本发明中计数变量设定为20，即如果连续二十代全局最优解保持不变的话，此时满足停机条件，程序终止，进入步骤13。如果程序运行中，计数变量未达到20，而当前代全局最优解与上一代全局最优解不同，则计数变量清零，从零开始计数，蝙蝠算法运行代数n=n+1，并返回步骤6。

步骤13，根据最终得到的全局最优解fitness_generation对应的6个参数a₁、b₁、c₁、a₂、b₂、c₂，得到最终分割阈值，并对图像进行分割输出；阈值t₁和阈值t₂的公式如下：

若(a₁+c₁)/2≤b₁≤c₁,则 $t_1=a_1+\sqrt{(c_1-a_1)*(b_1-a_1)/2}$

若a₁≤b₁＜(a₁+c₁)/2,则 $t_1=c_1-\sqrt{(c_1-a_1)*(c_1-b_1)/2}$

若(a₂+c₂)/2≤b₂≤c₂,则 $t_2=a_2+\sqrt{(c_2-a_2)*(b_2-a_2)/2}$

若a₂≤b₂＜(a₂+c₂)/2,则 $t_2=c_2-\sqrt{(c_2-a_2)*(c_2-b_2)/2}$

使用t₁和t₂对图像进行分割，即将图像中灰度值在[0,t₁)范围内的像素作为第一类，图像中灰度值在[t₁,t₂)范围内的像素作为第二类，图像中灰度值在[t₂,255]范围内的像素作为第三类，并将最终分割图像输出。

为了验证本发明的有效性，将本发明与基于粒子算法优化的双阈值模糊熵方法、基于遗传算法优化的双阈值模糊熵方法进行了仿真比较。在仿真实验过程中，采用演化算法的分割方法停机条件都是最优解保持连续二十代不变。本发明的参数设置为：S=50，R_A=1.5。基于上述设定的前提下，在不同的图像上分别进行仿真。

在仿真实验中，我们采用了三个已知指标来评价分割结果的好坏。他们分别是文献[7]中提出的搜索失败次数Fail_Search，文献[3-4]中定义的模糊熵和算法迭代次数Iter_Num。

搜索失败次数[7]Fail_Search是每10次独立运算算法最终无法得到最优解的次数，它主要用来判定提出的方法能否很好地应用的图像分割过程中；根据文献[3-4]中对模糊熵的定义，模糊熵越大时对应的图像分割结果越好；而算法迭代次数Iter_Num是整个算法最终需要完成迭代的次数，它主要是判定算法运行效率的标准，当算法迭代次数Iter_Num越小时表示算法运行效率越高。

本发明中，仿真实验是在Lena图和Peppers图上进行的。参照图3-1、3-2、3-3、4和图6-1、6-2、6-3和图7，他们是本发明和对比方法在Lena图和Peppers图上的仿真实验示例图。对比图3-1、3-2、3-3表示的各个分割结果可以得出，本发明的Lena图分割结果中可以基本将人物从背景图像中分割出来，且具有更加平滑的区域特征。对比图6-1、6-2、6-3表示的各个分割结果可以得出，本发明的Peppers图分割结果中各个辣椒具有更加连续的边缘特征，在一定程度上抑制了噪声对最后结果的影响。对比图4和图7，可以看出本发明与PSO-Multithreshold和GA-Multithreshold方法相比，收敛速度明显更快。

表1给出了本发明和对比方法对Lena图和Peppers图的100次独立分割后各项指标的平均值。可以看出本发明中的搜索失败次数Fail_Search相对较低。同时，本发明中算法迭代次数Iter_Num与PSO-Multithreshold和GA-Multithreshold方法相比明显偏小，而且本发明得到的模糊熵值也高于PSO-Multithreshold和GA-Multithreshold方法，说明本发明在性能和效率方面都优于PSO-Multithreshold和GA-Multithreshold方法。

表1Lena与Peppers图像上有效性验证指标

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其它的任何未违背本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (6)

1.一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法，其特征在于，包括下述步骤：

步骤1，输入待分割图像，种群规模S以及停机条件；

步骤2，依据设定的种群规模S，运行蝙蝠算法，初始种群的各变量，产生初始种群G_k，对于第i个个体而言，其位置向量x_i（这里i值范围是1≤i≤S，S表示种群规模），其取值范围为（0，255），速度v_i、声波频率f_i、脉冲响度A_i、发射速率r_i五个变量随机产生，代数n表示蝙蝠算法的运行代数，在初始条件下，令n=0；

步骤3，对不满足大小关系的个体进行扰动；

每个个体需满足0≤a₁＜b₁＜c₁＜a₂＜b₂＜c₂≤255的大小关系，各参数的取值范围均为（0，255），在初始化的过程中需要对不满足大小关系的个体进行如下的顺序扰动：

3.1.对a₁产生扰动，得到a₁ ¹，如果a₁ ¹＜0,则令a₁ ¹＝0，如果a₁ ¹＞250，则令a₁ ¹＝250；

3.2.对b₁产生扰动，得到b₁ ¹，如果b₁ ¹＜a₁ ¹+1,则令b₁ ¹＝a₁ ¹+1，如果b₁ ¹＞251，则令b₁ ¹＝251；

3.3.对c₁产生扰动，得到c₁ ¹，如果c₁ ¹＜b₁ ¹+1,则令c₁ ¹＝b₁ ¹+1，如果c₁ ¹＞252，则令c₁ ¹＝252；

3.4.对a₂产生扰动，得到a₂ ¹，如果a₂ ¹＜c₁ ¹+1,则令a₂ ¹＝c₁ ¹+1，如果a₂ ¹＞253，则令a₂ ¹＝253；

3.5.对b₂产生扰动，得到b₂ ¹，如果b₂ ¹＜a₂ ¹+1,则令b₂ ¹＝a₂ ¹+1，如果b₂ ¹＞254，则令b₂ ¹＝254；

3.6.对c₂产生扰动，得到c₂ ¹，如果c₂ ¹＜b₂ ¹+1,则令c₂ ¹＝b₂ ¹+1，如果c₂ ¹＞255，则令c₂ ¹＝255；

步骤4，计算灰度值为k的像素属于暗、灰、亮的隶属度值和条件概率；

${\mathrm{\&mu;}}_d\left(k\right)=\left\{\begin{array}{cc}1& k\&le;a_1\\ 1-\frac{{(k-a_1)}^2}{(c_1-a_1)*(b_1-a_1)}& a_1<k\&le;b_1\\ \frac{{(k-c_1)}^2}{(c_1-a_1)*(c_1-b_1)}& b_1<k\&le;c_1\\ 0& k>c_1\end{array}\right.$

${\mathrm{\&mu;}}_m\left(k\right)=\left\{\begin{array}{cc}0& k\&le;a_1\\ \frac{{(k-a_1)}^2}{(c_1-a_1)*(c_1-b_1)}& a_1<k\&le;b_1\\ 1-\frac{{(k-c_1)}^2}{(c_1-a_1)*(c_1-b_1)}& b_1<k\&le;c_1\\ 1& \\ 1-\frac{{(k-a_2)}^2}{(c_2-a_2)*(b_2-a_2)}& a_2<k\&le;b_2\\ \frac{{(k-c_2)}^2}{(c_2-a_2)*(c_2-b_2)}& b_2<k\&le;c_2\\ 0& k>c_2\end{array}\right.$

${\mathrm{\&mu;}}_b\left(k\right)=\left\{\begin{array}{cc}0& k\&le;a_2\\ \frac{{(k-a_2)}^2}{(c_2-a_2)*(c_2-b_2)}& a_2<k\&le;b_2\\ 1-\frac{{(k-c_2)}^2}{(c_2-a_2)*(c_2-b_2)}& b_2<k\&le;c_2\\ 1& k>c_2\end{array}\right.$

$p_d=\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}{p_k*\mathrm{\&mu;}}_d\left(k\right)$

$p_m=\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}{p}_k*{\mathrm{\&mu;}}_m\left(k\right)$

$p_b=\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}{p}_k*{\mathrm{\&mu;}}_b\left(k\right)$

其中，μ_d(k)、μ_m(k)、μ_b(k)分别表示灰度值为k(0≤k≤255)的像素属于暗、灰、亮的隶属度；p_k为图像灰度直方图的概率分布，p_d、p_m、p_b分别表示表示灰度值为k(0≤k≤255)的像素属于暗、灰、亮的条件概率；

步骤5,计算种群中个体的适应度并找出其全局最优解；

根据每个个体对应的像素属于暗、灰、亮隶属度和条件概率，通过模糊熵公式计算每个个体的适应度值，将每个个体适应度值H_s(i)（这里i值范围是1≤i≤S，S表示种群规模）表中最大的值fitness_generation作为种群的全局最优解，并存入最优解集；

步骤6，采用蝙蝠算法对种群中各变量进行更新；

调整声波频率fi产生新的解并更新速度v_i和位置x_i，更新公式如下：

f_i＝f_min+(f_max-f_min)^*β

v_i ⁿ＝v_i ^n-1+(x_i ⁿ-x_globalbest)^*f_i

x_i ⁿ＝x_i ^n-1+v_i ⁿ

其中，n表示蝙蝠算法的运行代数，f_min和f_max分别表示声波频率的范围，β∈[0,1]是一个随机向量，保证f_i在[f_min,f_max]范围内，x_globalbest为全局最优解的位置向量，速度v_i可为正值或负值，位置向量x_i可任意方向移动；

步骤7，从最优解集中选择一个解，并在该最优解附近形成一个局部解，进而在该局部解附近形成一个新解；

若n=1,局部解取步骤5中全局最优解所对应的位置变量，若n≠1，局部解直接取前一代个体适应度的全局最优解所对应的位置变量x_globalbest，新解的计算公式如下：

x_new＝x_globalbest+ε^*avgA^k

其中，x_new表示得到的新解，ε∈[-1,1]是一个任意的实数，avgA^k是个体在当前代的平均脉冲响度，由脉冲响度A_i计算所得；

步骤8，计算新解中个体的适应度；

步骤9，判断是否对个体参数进行更新；

用新解x_new中个体的适应度用H_{s_new}(i)与H_s(i)进行比较，（若n=1,H_s(i)取步骤5中全局最优解的个体适应度H_s(i)，若n≠1，H_s(i)取算法上一代计算产生的最优解集的个体适应度H_s(i)），若H_s(i)＜H_{s_new}(i)(1≤i≤S)且脉冲响度A_i大于人为设定阈值R_A，则进入步骤10；否，则进入步骤11；

步骤10，对原解x中第i个个体的位置变量用新解x_new中第i个个体的位置变量进行替换并对脉冲响度A_i和发射速率r_i进行更新后，进入步骤11，脉冲响度A_i和发射速率r_i的更新公式如下；

$A_i^{n+1}=\mathrm{\&alpha;}{A}_i^m,r_i^{n+1}=r_i^0[1-\exp (-\mathrm{\&gamma;n})]$

其中n表示蝙蝠算法的运行代数，exp表示以自然对数e为底的指数函数，α和β是两个常量；

步骤11，计算种群全局最优解，通过与先前全局最优解的比较，对全局最优解进行更新；

计算种群中个体的适应度，得到当前的全局最优解fitness_{new_gen}与先前的fitness_generation进行比较，若fitness_{new_gen}＞fitness_generation，则将其更新为种群新的全局最优解；否则，将保留原先的全局最优解fitness_generation；

步骤12，判断是否满足停机条件，若是，则进入步骤13，否则返回步骤6；

步骤13，根据最终得到的全局最优解fitness_generation对应的6个参数a₁、b₁、c₁、a₂、b₂、c₂，得到最终分割阈值；

阈值t₁和阈值t₂的公式如下：

若(a₁+c₁)/2≤b₁≤c₁,则 $t_1=a_1+\sqrt{(c_1-a_1)*(b_1-a_1)/2}$

若a₁≤b₁＜(a₁+c₁)/2,则 $t_1=c_1-\sqrt{(c_1-a_1)*(c_1-b_1)/2}$

若(a₂+c₂)/2≤b₂≤c₂,则 $t_2=a_2+\sqrt{(c_2-a_2)*(b_2-a_2)/2}$

若a₂≤b₂＜(a₂+c₂)/2,则 $t_2=c_2-\sqrt{(c_2-a_2)*(c_2-b_2)/2}$

使用t₁和t₂对图像进行分割，即将图像中灰度值在[0,t₁)范围内的像素作为第一类，图像中灰度值在[t₁,t₂)范围内的像素作为第二类，图像中灰度值在[t₂,255]范围内的像素作为第三类，并将最终分割图像输出。

2.根据权利要求1所述的一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法，其特征在于：步骤1中，群体规模S=50，停机条件为：连续二十代运算得到的全局最优解保持不变。

3.根据权利要求1所述的一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法，其特征在于：步骤6中，蝙蝠算法频率范围的取值为f_min＝0、f_max＝10。

4.根据权利要求1所述的一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法，其特征在于：步骤9中，阈值R_A=1.5。

5.根据权利要求1所述的一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法，其特征在于：步骤10中，α＝β＝0.9。

6.根据权利要求1所述的一种基于蝙蝠算法优化模糊熵的双阈值图像分割方法，其特征在于：步骤5、步骤8和步骤11中所述的适应度值计算公式如下：

$H_d=-\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}\frac{p_k*{\mathrm{\&mu;}}_d\left(k\right)}{p_d}*\mathrm{In}\left(\frac{p_k*{\mathrm{\&mu;}}_d\left(k\right)}{p_d}\right)$

$H_m=-\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}\frac{p_k*{\mathrm{\&mu;}}_m\left(k\right)}{p_m}*\mathrm{In}\left(\frac{p_k*{\mathrm{\&mu;}}_m\left(k\right)}{p_m}\right)$

$H_b=-\underset{k=0}{\overset{255}{\mathrm{\&Sigma;}}}\frac{p_k*{\mathrm{\&mu;}}_b\left(k\right)}{p_b}*\mathrm{In}\left(\frac{p_k*{\mathrm{\&mu;}}_b\left(k\right)}{p_b}\right)$

H＝H_d+H_m+H_b

其中，H为图像的总模糊熵，H_d、H_m、H_b分别表示暗、灰、亮这三个类别中各类的类内模糊熵，μ_d(k)、μ_m(k)、μ_b(k)和p_d、p_m、p_b分别表示灰度值为k(0≤k≤255)的像素属于暗、灰、亮的隶属度值和条件概率，p_k为图像灰度直方图的概率分布,H_s(i)(1≤i≤S)表示每个个体适应度值，也是每个个体的模糊熵值。

Images