小波域局部信息约束的鲁棒模糊聚类算法

李昌兴，王佳烨，吴成茂，乔彩彩，郝思佳

(1.西安邮电大学 理学院，陕西 西安 710121；2.西安邮电大学 通信与信息工程学院，陕西 西安 710121；3.西安邮电大学 电子工程学院，陕西 西安 710121；4.河南工学院 经济学院，河南 新乡 453003)

图像分割是根据像素的特征将图像分成不重叠的子区域，提取特定的特征区域，以备图像理解与图像识别[1]。目前，图像分割已广泛应用于生物医学工程[2]、智能交通[3]、军事安防[4-5]、智慧农业工程[6]等诸多领域。图像的特征提取和相关聚类算法的选择是图像分割的关键[7]。

模糊C-均值聚类算法[8](Fuzzy C-means Clustering Algorithm，FCM)因其简单易实现被广泛应用，但FCM算法对噪声极敏感。结合空间约束的模糊C-均值聚类[9](Fuzzy C-means Dlustering with Spatial Information Constraints,FCM_S)算法在FCM算法的基础上引入邻域信息，提高了算法抗噪性，然而，FCM_S算法在每次迭代过程中都需要计算像素点的邻域信息。为此，提出了FCM_S1与FCM_S2算法[10]，分别利用均值滤波和中值滤波在迭代前计算像素的邻域信息，算法的运行效率和抗噪性均有提升，但丢失了部分图像细节信息。广义模糊C-均值聚类算法[11](Generalized Fuzzy C-Means Clustering Algorithm，GFCM)在FCM算法的基础上引入了一个新的模糊参数，多数情况下GFCM算法性能优于FCM算法，但抗噪性依然较差。在此基础上，有学者提出了含有空间约束的广义模糊C-均值聚类[12](Generalized Fuzzy C-Means Clustering with Spatial Information Constraints,GFCMS)算法，通过中值滤波提升了算法鲁棒性。然而，上述算法中都需要人为控制的正则化参数。为此，后续的研究提出了模糊局部信息C-均值聚类算法[13](Fuzzy Local Information C-Means Clustering Algorithm，FLICM)，该算法无需人为设置正则化参数，自适应能力有所提升，但当噪声强度较大时，算法的抗噪性依然不理想。后来有学者在FLICM的基础上用邻域点的像素信息修正模糊加权因子，提出了局部空间信息的模糊C-均值算法[14](Weighted Fuzzy Local Information C-Means Clustering Algorithm,WFLICM)，以及在WFLICM算法中融合了核函数的核加权模糊局部信息C-均值聚类算法[15](Kernel Weighted Fuzzy Local Information C-Means Clustering Algorithm，KWFLICM)，进一步提高了算法对噪声的抑制能力，但也导致图像细节信息丢失，分割精度下降，且加权因子与核函数的引入使得算法的复杂度增加。

针对上述模糊聚类算法对图像纹理细节信息保留不佳、抗噪性能较差的问题，在FCM_S算法中融入了纹理信息，提出了一种基于小波变换的含有局部信息约束的鲁棒模糊聚类彩色图像分割算法。利用小波变换提取图像纹理特征，对所得小波系数建立广义高斯分布(Generalized Gaussian Density,GGD)模型[16]，并以Kullback-Leibler(KL)散度作为两个GGD模型之间的差异度量。同时，在Lab颜色空间[17-18]计算样本与聚类中心之间的欧式平方距离。分别在Lab颜色空间与小波域引入邻域信息。利用拉格朗日乘子法构建优化模型，得到迭代的鲁棒模糊聚类算法，以期算法能达到满意的分割精度与抗噪鲁棒性。

1 预备知识

1.1 FCM算法

FCM算法[8]通过隶属度强度将像素划分成不同类别，以实现图像分割。FCM算法的目标函数定义为

(1)

相关约束条件为

1)0≤uki≤1,i=1,2,…,n,k=1,2,…,c

其中：m为模糊因子；uki表示第i(i=1,2,…,n)个像素xi相对于第k(k=1,2,…,c)个聚类中心vk的隶属度；‖xi-vk‖2表示第i个像素xi到第k个聚类中心vk的欧式平方距离。

依据拉格朗日乘子法得到隶属度uki与聚类中心vk的迭代表达式为

(2)

(3)

FCM算法采用如下步骤。

步骤1设置聚类数c、模糊因子m、迭代终止条件ε以及最大迭代次数tmax，并令初始迭代次数t=0。

步骤6根据隶属度最大原则对像素分类。

1.2 FCM_S算法

FCM_S算法[9]将空间邻域信息约束作为惩罚项加入FCM算法的目标函数，其目标函数可以表示为

(4)

相关约束条件为

1)0≤uki≤1,i=1,2,…,n,k=1,2,…,c

其中：参数ξ控制邻域像素对聚类结果的贡献；NR表示以xi为中心的邻域窗口内含有的像素数；Ni为邻域窗口内像素的集合。

依据拉格朗日乘子法得到隶属度uki与聚类中心vk的迭代表达式为

(5)

(6)

2 小波域局部约束的模糊聚类算法

FCM_S算法虽考虑了邻域信息，使得算法抗噪性能有所提升。但其分割精度较低，且当图像受噪声污染严重时，算法的分割效果仍不理想[19-20]。为了进一步提升FCM_S算法的抗噪鲁棒性，在其目标函数中引入纹理信息，并提出一种小波域局部信息约束的模糊聚类算法。改进算法的原理示意图如图1所示。

图1 改进算法的原理示意图

2.1 图像的特征表示

Lab颜色空间包含3个通道，即一个亮度通道l，以及两个颜色通道a和b[17]。l分量取值范围为[0,100]，其值越大亮度越高。a分量取值范围为[-128,127]，其值由负数到正数，对应颜色从绿色到红色；b分量取值范围为[-128,127]，其值从负数到正数，对应颜色从蓝色到黄色。Lab颜色空间具有感知均匀性，其色域宽阔，更符合人类视觉感知[18]。

令xi=(li,ai,bi)与xj=(lj,aj,bj)分别表示像素i处与像素j处的颜色信息，则像素i与像素j的平方欧氏距离[19]可以表示为

‖xi-xj‖2=(li-lj)2+(ai-aj)2+(bi-bj)2

(7)

设Ii表示以第i个像素为中心的21×21的图像窗口，并对图像在窗口上进行N级小波变换，分解得到N个低频子带和3N个高频子带，建立广义高斯分布模型来表征图像纹理特征变化的分布[20]，图像纹理特征变化的分布为

(8)

其中：Γ(z)为伽马函数；w为n个一系列独立的小波系数；α为广义高斯分布的尺度参数，与标准差有关；β为形状参数，描述分布的衰减速度。当β=1和β=2时，广义高斯分布进化为拉普拉斯分布和高斯分布[19-20]。

为了快速准确地估计形状参数，利用凸变换[21]求解参数α与β。求解形状参数β的表达式为[19,21]

(9)

利用牛顿-拉夫森迭代法[22]求解形状参数

(10)

以初始值βint=1.0可从半开区间中的任意起点收敛到唯一的全局最优解，获得精确的形状估计值[19-21]，再由式(8)计算尺度参数

(11)

KL散度是衡量两个概率分布差异的非对称性的度量[23]，其计算式可表示为

(12)

将式(8)代入式(12)可得

(13)

则像素i与像素j处的局部窗口上小波分解所得高频子带GGD模型的KL散度之和[24]可表示为

(14)

其中，αir与βir分别表示在以像素i为中心的窗口上小波变换所得第r(r=1,2,3,…,3N)个高频子带GGD模型的尺度参数与形状参数。

2.2 算法构建与最优解推导

为了提高图像的分割性能，在提取颜色信息的基础上，增加图像的纹理信息，并引入小波域局部邻域信息。改进算法的目标函数可表示为

(15)

需要满足约束条件为

1)0≤uki≤1,i=1,2,…,n,k=1,2,…,c

其中：λ和τ为正则化参数；Dki表示以像素i为中心的局部窗口和以聚类中心像素k为中心的局部窗口上小波变换所得高频子带所对应的GGD模型KL散度之和；uki表示像素i对聚类中心像素k的隶属度；vk表示包含所有颜色特征的聚类中心。

为使目标函数服从于隶属度和唯一性，利用拉格朗日乘子将其转化为无约束极值问题，即

(16)

其中，μ=(μi)是一系列拉格朗日乘子。

首先对式(16)关于隶属度uki求偏导

(17)

(18)

将式(17)带入式(18)可得

(19)

则隶属度uki的迭代更新表达式为

(20)

对式(16)关于vk求偏导，可得

(21)

求解式(21)可得颜色聚类中心vk的迭代表达式

(22)

求解GGD模型参数，对式(16)关于αkr求偏导可得

(23)

求解式(23)，并令

Hkir=Ψ((βkr+1)/βir)/βir+lgαir

其中，Ψ(z)为双伽马函数。则尺度参数αkr的迭代式和形状参数βkr的超越方程可以分别表示为

(24)

(25)

考虑到超越方程的求解较为复杂，将求解方程f(βkr)=0时βkr的取值问题，转化为求解f(βkr)2=0的极小值问题。利用二分法[24]对非线性方程进行一维搜索，根据使用牛顿-拉弗森迭代法求解结果，β的取值绝大多数在[0.5,3]区间内，为了快速更新形状参数聚类中心βkr，设置二分法的搜索区间为[0.5,3]。

各聚类中心的更新与邻域信息密切相关，且聚类中心在各空间独立更新。隶属度的更新与像素颜色信息、GGD模型参数以及邻域信息都相关。改进算法具体的实现步骤如下。

步骤1设置聚类数c，正则化参数λ和τ，阈值ε，最大迭代次数tmax以及Ni与Nr，令初始迭代次数t=0。

步骤2对局部窗口Ii进行级小波分解，建立GGD模型，并根据式(10)与式(11)求解GGD模型参数αir和βir。

步骤3随机初始化隶属度U0。

步骤5根据式(20)计算Ut+1，若满足|Jt+1-Jt|<ε或t=tmax，则输出(Ut,Vt,At,Bt)并执行步骤6；否则，返回执行步骤4，并令t=t+1。

3 实验结果及分析

使用数据集BSD500，纹理图像集Describable Textures Dataset和遥感图像集NWPU-RESISC45中的图像进行实验。将FCM、FCM_S、FCM_S1、FCM_S2、GFCM、GFCMS、FLICM、WFLICM以及KWFLICM算法等9种算法作为对比算法。实验测试环境为Matlab R2018a。

设置参数λ=0.5，τ=0.01，ε=10-1，FCM_S的ξ1=3.8，FCM_S1的ξ2=3.8，FCM_S2的ξ3=3.8，GFCM的ξ4=0.9，GFCMS的ξ5=0.9、ξ5=6，邻域窗口为5×5，则邻域内除中心像素外像素数Nr=24。

采用划分系数Vpc、划分熵Vpe、F值、误分率以及峰值信噪比(Peak Signal-to-Noise Ratio,PSNR)作为图像分割效果的衡量指标[26]。其中：划分系数Vpc与划分熵Vpe都反映分割矩阵的模糊性，Vpc值越接近1，Vpe越接近0，聚类越明确，算法性能越好[27]；F值为精确率与召回率的调和平均值，综合分割算法的查准率与查全率指标，当F值较高时则能说明分割方法具有较高的有效性[21]；误分率描述错误分类的像素占比，该值越小，分割效果越好[28]；PSNR值越大，算法抑制噪声的能力越强，图像分割的效果越好[28]。

3.1 无噪图像的测试与分析

为了验证改进的算法对图像纹理的分割精度，对自然图像“海星”“鳞片”以及遥感图像“立交”进行实验。3种原始图像如图2所示。

图2 原始图像

其中，自然图像“海星”与“鳞片”包含大量纹理细节信息，遥感图像“立交”目标类别与细节繁多，并且颜色相近，能够进一步地测试算法的分割精度。10种算法对不同图像分割的具体结果分别如图3至图8所示，分割性能指标具体数据如表1所示。

图3 10种算法对图像海星的分割结果

表1 10种算法对图2(a)-图(c)的分割性能指标

从视觉效果来看，图3和图5显示，FCM、FCM_S以及GFCM算法与改进算法在一定程度上均保留了图像轮廓与细节部分，但通过布局放大图像图4与图6来看，改进算法保留的细节信息更加清晰完整，并且去除了部分杂点。其余算法均丢失了图像的细节信息，大多区域边缘光滑呈块状分布。这是因为，中值滤波、均值滤波或者权重因子的引入使得算法在聚类过程中考虑了邻域像素信息，这虽有利于对去除杂点和噪声点，但对图像的细节信息处理效果不理想。

图4 图3分割结果的局部放大

图5 10种算法对图像鳞片的分割结果

图6 图5中分割结果的局部放大

另外，由图7与图8可见，改进算法的分割结果中，公路边界完整清晰，且有效地分割出了公路上目标较小的汽车。相对于其余算法，改进算法能较好地去除图像中的杂点，并保留了部分细节信息。由表1可以看出，改进算法的误分率相对于其余算法均明显减小，F值明显提升。结合图3至图8与表1可以发现，改进算法对于图像细节信息保留较完整，分割结果较为理想。

图7 10种算法对图像立交的分割结果

3.2 噪声图像的测试与分析

为了进一步验证改进算法的抗噪鲁棒性，选取“蝴蝶”“建筑”和“机场”等3幅图像进行测试。用于测试的原始图像与加噪图像分别如图9所示。其中：对图像蝴蝶添加强度为10%的椒盐噪声；对图像建筑添加均值为0，归一化方差为0.1的高斯噪声；由于遥感图像所含噪声水平大小未知，对图像机场添加归一化方差为0.05的斑点噪声。各算法分割结果分别如图10至图13所示，分割性能各指标具体数据如表2所示。

图8 图7中分割结果的局部放大

图9 原始图像与加噪图像

图10 10种算法对图9(d)的分割结果

由图10可知，在椒盐噪声干扰下，FCM_S2、GFCMS、KWFLICM以及改进算法均能获得较理想的分割结果。由图11中对图10分割结果的局部放大可见，FCM_S2算法与GFCMS算法仍有部分噪声点，KWFLICM算法虽有效地抑制了椒盐噪声的影响，但丢失了“蝴蝶”翅膀上花纹等细节信息。相较于其他算法，改进算法残留的噪声点较少，具有出明显的抗噪性能优势，且所得分割目标细节信息较清晰完整。其他算法的分割结果中仍存在大量噪声点，且出现了错分现象。分割无噪纹理图像非常有效的FCM和GFCM算法严重失真。这是因为，这些算法没有利用局部约束信息，导致其在没有杂声干扰情况下，对图像的细节信息保留较好，但不具有抗噪鲁棒性。由图12与图13可见，在高斯噪声与斑点噪声的干扰下，对比算法的分割结果都存在大量噪声点，改进算法对噪声的抑制能力较好，获得的同质区域合理清晰。从表2中可以看出，相对于其他算法，改进算法能获得较大的峰值信噪比与较小的误分率，进一步表明改进算法分割精度更高，对噪声的抑制能力更强。

图11 图10中分割结果的局部放大

图12 10种算法对图9(e)的分割结果

图13 10种算法对图9(f)的分割结果

表2 10种算法对图9(d)-图9(f)的分割性能指标

3.3 小波变换窗口测试

在实验中，对图像在窗口Ii上进行小波变换，为了研究窗口的尺寸大小对于改进算法分割性能的影响，将Ii的大小分别设置为13×13、15×15、17×17、19×19、21×21、23×23、25×25、27×27和29×29，对图2(a)-图2(c)进行实验，所得的误分率如表3所示。

表3 改进算法在不同尺寸小波变换窗口时的误分率

若窗口Ii过小，则窗口内图像信息过少，也可能导致窗口内几乎全部为低频信息，不利于牛顿-拉夫森算法快速求解GGD模型的参数；若窗口过大，则窗口有可能跨越多个特征区域，使得GGD模型参数不能合理表征该像素点周围的像素纹理特征，导致分割性能下降和运行时间的增加。由表8可知，选取窗口Ii的大小为21×21有利于获得最佳分割结果。

3.4 复杂度分析与时间开销测试

为了比较上述对比算法与改进算法的执行效率，进行计算复杂度分析。对于一幅像素为n且可分为c类的图像，各算法计算复杂度如表4所示。其中：l表示迭代次数；q表示当前像素所在邻域窗口大小。从表4可以看出，FCM算法和GFCM算法的计算复杂度最低，因为这两种算法仅包含最简单的迭代。FCM_S算法和FLICM算法需要计算每次迭代中邻域信息，FCM_S1、GFCMS与FCM_S2算法需要计算每个像素均值或中值信息，计算复杂度较高。WFLICM和KWFLICM算法的计算复杂度包含相似性度量因子的计算和算法的迭代，这两种算法的计算复杂度最高。改进算法的复杂度为算法迭代复杂度，技术复杂度居中。

表4 各算法计算复杂度分析

为了直观地展现各算法的时间开销，将测得各算法对不同图像的运行时间绘制成柱状图，不同算法对各测试图像的运行时间对比如图14所示。从图14中可以看出，WFLICM算法与KWFLICM算法较为耗时，这是由于加权因子的引入和邻域信息的多次计算。改进算法的运行时间较长，仅小于FLICM、WFLICM和KWFLICM算法，这是因为，改进算法的计算复杂度虽然不高，但其特征提取过程较复杂，尤其是小波变换、建立GGD模型以及求解模型参数的过程耗时较长。图14中改进算法的时间开销不包括特征提取过程，仅为聚类算法时间开销。改进算法对不同大小的图像特征提取的时间开销如表5所示。可以看出，改进算法对图像特征提取过程的时间消耗与图像大小有关。

图14 10种算法对各测试图像的运行时间对比

表5 改进算法处理不同图像时特征提取的时间消耗

改进算法在图像特征提取过程中耗时较长，但其聚类过程可在较短时间内收敛，且能达到更理想的分割效果。

4 结语

提出了一种基于小波变换的局部信息约束的鲁棒模糊聚类算法。将纹理信息嵌入目标函数，并在Lab颜色空间与小波域分别引入了邻域信息。实验结果表明，相较于其他对比算法，对于纹理图像与噪声图像，改进算法能够获得满意的分割结果，具有更高的分割精度与更强的抗噪鲁棒性。特别是针对受不同强度及类型的噪声污染的图像，改进算法能获得更合理的同质区域，有效地抑制噪声，又一定程度地保留了图像的细节信息。另外，改进算法能获得更低的误分率与更高的峰值信噪比。但是，仍存在一些待优化的问题，例如，改进算法中的特征提取与参数求解过程导致算法时间开销上升，若能以一个高效的模型来表征图像纹理特征，则算法的效率将大大提升。另外，若能使算法自适应地确定聚类数目，则将能够有效地提升算法的实用性。