开放科学（资源服务）标志码（OSID）：

0 概述

图像分割是根据低层特征的相似性或相异性将图像分割成若干个连续的子区域，使同一区域内的特征相似，而不同区域内的特征相异。图像分割方法可大致分为基于边缘^[1]、基于阈值^[2]、基于聚类^[3-4]、基于区域^[5]、基于图论^[6]的5类方法。其中聚类是一种无监督分类方法^[7]，在像素没有标签的情况下，利用类间差别尽可能大而类内差别尽可能小的原则将原像素划分成不同的类，分为硬聚类^[8]和模糊聚类^[9]。硬聚类强调每个像素数据只能被归为一类，模糊聚类通过隶属度来强调每个像素数据属于每一类的程度，而不是硬性地将数据归到某一类。聚类分析的核心是利用数学方法将目标样本集按某种准则自动归为若干个互不重叠的类，使同一类中的样本相似，不同类中的样本互不相同，这一思想与图像分割相似。模糊聚类应用于分割问题时，引入隶属度表示样本属于某个特定集合的确定程度，体现了描述的模糊性。

在众多的模糊聚类算法中，模糊C-均值（Fuzzy C-Means，FCM）聚类算法^[10]是较为经典的无监督聚类算法，考虑了像素对各个聚类中心的隶属度，但并未考虑像素本身和邻域像素之间的关联性，对于含噪图像的聚类难以获得满意的分割结果。为了增强FCM算法的抗噪性，AHMED等^[11]将局部空间限制项引入FCM目标函数，提出一种基于空间信息的FCM算法（FCM_S），由于该算法在每次迭代计算过程中都需要计算像素的邻域项，因此时效性较差。CHEN等^[12]提出两种FCM_S的改进方法：FCM_S1和FCM_S2算法，这两种算法分别采用邻域像素的灰度均值和灰度中值来表示像素的局部空间信息，从而有效提高了算法的计算效率。KRINIDIS等^[13]利用像素点的局部信息构造空间模糊因子，并引入FCM提出模糊局部信息C-均值（Fuzzy Local Information C-Means，FLICM）聚类算法，该算法相比于FCM_S无需人为设置算法中的参数，能自适应控制噪声的干扰和保留图像细节信息间的平衡。ZHANG等^[14]引入包含非局部信息的像素相关度代替FLICM算法中的空间约束并改进模糊因子，提出非局部信息模糊C-均值聚类算法（Non-Local information Fuzzy C-Means，NLFCM），该算法有较好的抗噪性，但却抹去了图像的部分细节信息。

传统FCM算法利用欧式距离来计算像素到聚类中心的距离，原理是将像素空间中较为相近的像素聚在一起，仅适合聚类像素数目相差不悬殊的团状凸数据集，对于复杂的非凸数据集，聚类性能急剧下降。为此，很多学者将核函数^[15]的思想引入FCM算法，提出核模糊C-均值（Kernal Fuzzy C-Means，KFCM）聚类算法^[16]，该算法将像素数据通过非线性函数映射至高维特征空间，使得原空间线性不可分的像素变得线性可分，以此来改善聚类性能。由于高斯核函数的控制参数较少且性能较稳定，因此相比其他核函数应用更广泛。GONG等^[17]在FLICM算法的基础上进行改进，提出一种基于核度量的加权模糊局部信息C-均值（Kernel metric-based Weighted Fuzzy Local Information C-Means，KWFLICM）聚类算法。该算法利用空间差异性定义一个加权平衡模糊因子，包含了图像邻域像素的空间距离约束信息和空间灰度约束信息，并采用高斯核空间距离代替像素点间的欧式距离，提高了算法对噪声的鲁棒性。但由于引入过多的邻域信息，因此使得图像的部分细节信息被抹去。

BUADES等^[18]提出非局部均值（Non-Local Means，NLM）图像去噪方法，采用图像块间的欧式距离确定像素点间的相似度权重，通过像素点间的加权平均实现均值滤波效果。引导滤波器^[19]在一定程度上可以起到保边去噪的作用，且计算效率高，十分适合弱过滤工作^[20]，因此众多学者将引导滤波与模糊聚类算法相结合来提高模糊聚类算法的分割性能。引导滤波利用引导图像的信息对输入的图像进行滤波处理得到输出图像，相当于一个局部线性转换器使得输出和输入呈线性关系^[21]。考虑到引导滤波保边平滑的特点，GUO等^[22]提出图像引导模糊C-均值（Image Guided Fuzzy C-Means，IGFCM）聚类算法。在该算法中的引导滤波作为正则项加入模糊C-均值算法的目标函数，从而根据引导图像的信息对隶属度进行约束，但由于模型复杂难以计算，因此未得到很好的应用。为了简化IGFCM算法，GUO等^[23-24]提出基于引导滤波的改进模糊C-均值（modified FCM with Guided Filtering，FCM+GF）方法，在该方法中将需要分割的图像作为引导图像，对模糊C-均值聚类算法得到的隶属度进行处理，但该方法对强噪声图像分割结果并不理想，且仅限于灰度图像分割。对于彩色图像，执行多通道的引导滤波是一项困难的工作，鉴于此问题，XU等^[25]提出一种加权多通道引导滤波（Weighting Multi-channel Guided image Filtering，WMGF）算法来提高边缘保护，降低计算复杂度，基于该算法又提出一种新的彩色图像模糊C-均值聚类算法（FCM_WMGF），为了进一步提高不同噪声的抗噪性和边缘保护性，将多元形态重建方法^[26-27]引入改进的模糊C-均值聚类算法（MMR_FCM_WMGF）来获得更高的分割精度，该算法对于彩色图像具有较好的分割效果。LIU等^[28]提出一种鲁棒引导图像滤波（Robust Guided Image Filtering，RGIF）技术，将一个数据项和平滑项组成的优化框架融入引导滤波。SHEN等^[29]考虑输入图与引导图结构的不一致性，提出一种相互结构联合滤波算法并得到了较好的结果。

上述方法对于引导滤波的利用都是通过像素对模糊聚类算法得到的隶属度进行引导，然而该方法仅适用于处理低噪声图像。为了提高基于引导滤波的模糊聚类算法的分割性能和抑噪能力，本文将像素引导隶属度和隶属度引导像素的交替引导滤波引入KWFLICM算法，提出一种新的基于交替引导滤波的鲁棒模糊聚类算法_GFKWFLICM_GF。

1 相关工作 1.1 引导滤波

引导滤波是一种线性滤波器，基本思想是：在图像梯度比较大的地方实现保边，要求尽量不平滑，最好是输出与输入一样；在图像梯度较小的地方尽量平滑，输出与输入可以稍有不同。假设输入图为$ \boldsymbol{X} $，引导图为$ \boldsymbol{I} $，输出图为$ {\boldsymbol{X}}^{\text{'}} $，引导滤波的目的是使得输出图$ {\boldsymbol{X}}^{\text{'}} $和输入图$ \boldsymbol{X} $尽可能相似，同时纹理部分和引导图$ \boldsymbol{I} $尽可能相似。引导滤波器假设滤波后的结果和引导图像在滤波窗口内呈线性关系，即：

$ {X}_{j}^{\text{'}}={a}_{i}{I}_{j}+{b}_{i} \text{，} \forall j\in {w}_{i} $

(1)

其中：$ {w}_{i} $是以$ i $为中心像素的邻域窗口；$ {a}_{i} $和$ {b}_{i} $为该窗口所对应的线性系数，该线性系数保证了输出图像$ {\boldsymbol{X}}^{\text{'}} $和引导图$ \boldsymbol{I} $梯度的相对不变性，即$ \nabla {\boldsymbol{X}}^{\text{'}}=a\nabla \boldsymbol{I} $，在引导图有梯度的情况下，输出图才会有梯度。如果图像中任何一个地方的线性系数$ {a}_{i} $和$ {b}_{i} $都相同，则相当于对整个图像做了一个线性变换，并没有起到滤波的作用。若要达到保边滤波的效果，则需要在图像的任何地方$ {a}_{i} $和$ {b}_{i} $都是随着图像的内容变化的。为了求解引导滤波的线性系数$ {a}_{i} $和$ {b}_{i} $，假设所需要的系数能够使滤波输出图像的像素$ {X}_{j}^{\text{'}} $和输入图像的像素$ {X}_{j} $的差异最小，那么窗口$ {w}_{i} $的代价函数可以表示如下：

$ E({a}_{i}, {b}_{i})=\sum\limits _{j\in {w}_{i}}\left(\right({a}_{i}{I}_{j}+{b}_{i}-{X}_{j}{)}^{2}+\varepsilon {a}_{i}^{2}) $

(2)

其中：$ \varepsilon $是正则化系数，作用是防止$ {a}_{i} $过大。根据最小二乘法可以得出其线性系数的解：

$ {a}_{i}=\left(\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{I}_{j}{X}_{j}-{\stackrel{-}{I}}_{i}{\stackrel{-}{X}}_{i}\right)/(\varepsilon +{s}_{i}^{2}) $

(3)

$ {b}_{i}={\stackrel{-}{x}}_{i}-{a}_{i}{\stackrel{-}{I}}_{i} $

(4)

其中：$ {\stackrel{-}{I}}_{i} $和$ {s}_{i}^{2} $表示引导图$ \boldsymbol{I} $在窗口内的均值和方差；$ \left|{n}_{w}\right| $表示窗口内像素总数；$ {\stackrel{-}{X}}_{i} $表示窗口内输入图像$ \boldsymbol{X} $的均值。由于窗口是有尺寸的，输出图像中的一个像素点$ {X}_{j}^{\text{'}} $可以被不同窗口中的线性系数求得，且不同窗口得到的输出值不同，因此需要对这些值求均值，最终得到的$ {X}_{j}^{\text{'}} $表示如下：

$ {X}_{j}^{\text{'}}={\stackrel{-}{a}}_{i}{I}_{j}+{\stackrel{-}{b}}_{i} $

(5)

通过式（5）可以看出$ ({a}_{i}, {b}_{i}) $在空间上是有变化的，所以输出图像$ {\boldsymbol{X}}^{\text{'}} $不仅是引导图像$ \boldsymbol{I} $简单的尺度缩放，$ ({a}_{i}, {b}_{i}) $的值也取决于输入图像$ \boldsymbol{X} $。由于式（5）中对线性系数求平均值，$ ({a}_{i}, {b}_{i}) $的梯度变化是明显小于$ \boldsymbol{I} $的，因此$ {\boldsymbol{X}}^{\text{'}} $的像素值变化主要取决于引导图$ \boldsymbol{I} $，并且可以近似认为$ {\boldsymbol{X}}^{\text{'}}=a\nabla \boldsymbol{I} $。

当输入图像自身作为引导图时，$ \varepsilon =0 $、$ {a}_{j}=1 $、$ {b}_{j}=0 $，此时输入输出之间不存在变化。当$ \varepsilon > 0 $时，考虑两种极端情况。第一种情况是引导图像较平整，即在窗口内$ \boldsymbol{I} $为固定常数，那么$ {a}_{j}=0 $、$ {b}_{j}={\stackrel{-}{x}}_{i} $，此时模型退化为均值滤波。第二种情况是引导图像极其不平整，即在窗口内$ \boldsymbol{I} $的方差接近无穷大，那么$ {a}_{j}\approx 1 $、$ {b}_{j}\approx 0 $，此时输出与输入几乎无变化。通过方差评价图像是否平滑，方差大小由参数$ \varepsilon $表示，$ \varepsilon $越大，引导图中被判定为平整的区域越多，输出图像越平滑，反之，输出图像的平滑效果越差。

1.2 KWFLICM算法

FLICM算法通过模糊因子$ {G}_{ik} $平衡了局部邻域窗内各像素点对应的隶属度值，使得非噪声点像素及邻域窗口内的噪声像素的隶属度收敛到相似的值，因此算法鲁棒性得到了提升。KWFLICM算法对模糊因子$ {G}_{ik} $做出了改进，将局部方差系数构成的空间权重影响因子引入$ {G}_{ik} $，可以获得更多的空间相互作用信息，得到了鲁棒性更强的分割结果。KWFLICM算法的目标函数如下：

$ {J}_{m}(\boldsymbol{U}, \boldsymbol{V})=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}{\mu }_{ik}^{m}\left[\right(1-K({x}_{i}, {v}_{k}))+{G}_{ik}^{\text{'}}] $

(6)

其中：$ \boldsymbol{U}=({\mu }_{ik}{)}_{n\times c} $代表大小为$ n\times c $的隶属度矩阵，$ n $代表样本总数，$ c $代表样本被划分的类别数，$ {\mu }_{ik} $代表样本$ {x}_{i} $属于第$ k $类的程度；$ m $是模糊指数，本文默认m=2。

权重系数$ {G}_{ik}^{\text{'}} $定义如下：

$ {G}_{ik}^{\text{'}}=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}{\mu }_{ik}^{m}\sum\limits _{\underset{j\in {w}_{i}}{j\ne i}}{w}_{ij}(1-{\mu }_{jk}{)}^{m}(1-K({x}_{j}, {v}_{k})) $

(7)

权重因子$ {w}_{ij} $由空间位置信息$ {w}_{\mathrm{s}\mathrm{c}} $和邻域灰度信息$ {w}_{\mathrm{g}\mathrm{c}} $组成，详见文献[17]。利用最小二乘法优化其目标函数，得到隶属度和聚类中心的更新公式分别如下：

$ \begin{array}{l}{\mu }_{ik}=\\ \frac{1}{\sum\limits _{l=1}^{c}{\left[\frac{\left(\right(1-K({x}_{i}, {v}_{k}))\mathrm{ }+\sum\limits _{j\in {N}_{i}, j\ne i}{w}_{ij}(1-{\mu }_{jk}{)}^{m}(1-K({x}_{j}, {v}_{k}))}{\left(\right(1-K({x}_{i}, {v}_{l}))\mathrm{ }+\sum\limits _{j\in {N}_{i}, j\ne i}{w}_{ij}(1-{\mu }_{jl}{)}^{m}(1-K({x}_{j}, {v}_{l}))}\right]}^{1/(m-1)}}\end{array} $

(8)

$ {v}_{k}=\sum\limits _{i=1}^{n}\left({\mu }_{ik}^{m}K\right({x}_{i}, {v}_{k}\left){x}_{i}\right)/\sum\limits _{i=1}^{n}\left({\mu }_{ik}^{m}K\right({x}_{i}, {v}_{k}\left)\right) $

(9)

2 鲁棒模糊聚类算法 2.1 利用像素引导隶属度滤波的引导滤波

像素引导隶属度滤波的目的是得到新的隶属度，并使用新的隶属度矩阵作为划分矩阵来对图像进行聚类，目标函数如下：

$ E(\boldsymbol{U}, \boldsymbol{p}, \boldsymbol{q})=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}\sum\limits _{j\in {w}_{i}}\left[\right({p}_{ik}{I}_{j}+{q}_{ik}-{\mu }_{jk}{)}^{2}+\varepsilon {p}_{ik}^{2}] $

(10)

其中：$ {p}_{ik} $和$ {q}_{ik} $是当像素中心为$ i $时引导滤波的线性系数；$ \varepsilon $是防止$ {p}_{ik} $过大的正则化参数，本文默认$ \varepsilon =0.01 $。

该目标函数分别对$ {p}_{ik} $和$ {q}_{ik} $求偏导数并令其等于0，于是获得：

$ {q}_{ik}={\stackrel{-}{\mu }}_{ik}-{p}_{ik}{\stackrel{-}{I}}_{i} $

(11)

$ {p}_{ik}=\frac{\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{\mu }_{jk}{I}_{j}-\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{\mu }_{jk}\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{I}_{j}}{\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{I}_{j}^{2}-\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{I}_{j}\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{I}_{j}+\varepsilon } $

(12)

其中：$ {\stackrel{-}{\mu }}_{ik}=\sum\limits _{j\in {w}_{i}}{\mu }_{jk}/\left|{n}_{w}\right| $；$ {\stackrel{-}{I}}_{i}=\sum\limits _{j\in {w}_{i}}{I}_{j}/\left|{n}_{w}\right| $。

根据方差公式$ {s}_{i}^{2}=\sum\limits _{j\in {w}_{i}}({x}_{j}-{\stackrel{-}{x}}_{i})({x}_{j}-{\stackrel{-}{x}}_{i})/\left|{n}_{w}\right| $，$ {\stackrel{-}{x}}_{i}=\sum\limits _{j\in {w}_{i}}{x}_{j}/\left|{n}_{w}\right| $，可将式（12）化简如下：

$ {p}_{ik}=\left(\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{\mu }_{jk}{I}_{j}-{\stackrel{-}{\mu }}_{ik}{\stackrel{-}{I}}_{i}\right)/({s}_{i}^{2}+\varepsilon ) $

(13)

2.2 利用隶属度引导像素滤波的引导滤波

隶属度引导像素滤波的目的是得到新的像素值，将模糊聚类得到的隶属度传入利用隶属度引导像素滤波的代价函数，得到引导滤波的线性系数$ {p}_{1i} $和$ {q}_{1i} $，从而得到新的像素值。隶属度引导像素滤波的目标函数如下：

$ {E}_{1}(\boldsymbol{U}, {\boldsymbol{p}}_{1}, {\boldsymbol{q}}_{1})=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}\sum\limits _{j\in {w}_{i}}\left[\right({p}_{1i}{\mu }_{jk}+{q}_{1i}-{I}_{j}{)}^{2}+\varepsilon {p}_{1i}^{2}] $

(14)

其中：$ {p}_{1i} $和$ {q}_{1i} $是当像素中心为$ i $时引导滤波的线性系数，1是为了和前面的$ {p}_{ik} $、$ {q}_{ik} $作区分。

该目标函数分别对$ {p}_{1i} $和$ {q}_{1i} $求偏导并令其等于0，于是获得：

$ {q}_{1i}={\stackrel{-}{I}}_{i}-{p}_{1i}{\stackrel{-}{\mu }}_{ik} $

(15)

$ {p}_{1i}=\frac{\sum\limits _{k=1}^{c}\left(\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{I}_{j}\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{\mu }_{jk}\right)}{\sum\limits _{k=1}^{c}\left(\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{\mu }_{jk}^{2}-{\left(\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{\mu }_{jk}\right)}^{2}\right)+\varepsilon } $

(16)

令$ {s}_{1i}^{2}=\sum\limits _{k=1}^{c}\sum\limits _{j\in {w}_{i}}^{}({\mu }_{jk}-{\stackrel{-}{\mu }}_{ik})({\mu }_{jk}-{\stackrel{-}{\mu }}_{ik})/\left|{n}_{w}\right| $，$ {\stackrel{-}{\mu }}_{ik}=\sum\limits _{j\in {w}_{i}}{\mu }_{jk}/\left|{n}_{w}\right| $，于是可将式（16）化简如下：

$ {p}_{1i}=\left(\frac{1}{\left|{n}_{w}\right|}\sum\limits _{j\in {w}_{i}}{\mu }_{jk}{I}_{j}-{\stackrel{-}{\mu }}_{ik}{\stackrel{-}{I}}_{i}\right)/({s}_{1i}^{2}+\varepsilon ) $

(17)

2.3 KWFLICM算法与交替引导滤波的结合 2.3.1 算法流程

_GFKWFLICM_GF算法流程如图 1所示。原始加噪图像被当作引导图像，先利用KWFLICM算法对原始加噪图像进行聚类，将得到的隶属度分别传入像素引导隶属度滤波和隶属度引导像素滤波模块，分别得到新的隶属度和新的像素值，并通过新的隶属度对聚类中心进行更新，再将更新后的隶属度和像素值传入KWFLICM算法，循环上述过程，直到算法收敛。

2.3.2 优化过程

_GFKWFLICM_GF算法的迭代过程即优化以下多目标函数的过程：

$ \begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}{J}_{m}(\boldsymbol{U}, \boldsymbol{V})\mathrm{ }=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}{\mu }_{ik}^{m}\left[\right(1-K({x}_{i}, {v}_{k}))+{G}_{ik}^{\text{'}}]\\ E(\boldsymbol{U}, \boldsymbol{p}, \boldsymbol{q})\mathrm{ }=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}\sum\limits _{j\in {w}_{i}}\left[\right({p}_{ik}{I}_{j}+{q}_{ik}-{\mu }_{jk}{)}^{2}+\varepsilon {p}_{ik}^{2}]\\ L\left({\boldsymbol{\mu }}^{\text{'}}\right)\mathrm{ }=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}\sum\limits _{j\in {w}_{i}}({p}_{jk}{x}_{i}+{q}_{jk}-{\mu }_{ik}^{\text{'}}{)}^{2}\\ {E}_{1}(\boldsymbol{U}, {\boldsymbol{p}}_{1}, {\boldsymbol{q}}_{1})\mathrm{ }=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}\sum\limits _{j\in {w}_{i}}\left[\right({p}_{1i}{\mu }_{jk}+{q}_{1i}-{I}_{j}{)}^{2}+\varepsilon {p}_{1i}^{2}]\\ {L}_{1}\left({\boldsymbol{x}}^{\text{'}}\right)\mathrm{ }=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}\sum\limits _{j\in {w}_{i}}({p}_{1j}{\mu }_{ik}+{q}_{1j}-{x}_{i}^{\text{'}}{)}^{2}\end{array} $

s.t.$ 0\le {\mu }_{ik}\le 1 $，$ i=\mathrm{1, 2}, \cdots , n $，$ k=\mathrm{1, 2}，\cdots , c $

$ \sum\limits _{k=1}^{c}{\mu }_{ik}=1 \text{，} i=\mathrm{1, 2}, \cdots , n $

$ 0\le \sum\limits _{i=1}^{n}{\mu }_{ik}\le n \text{，} k=\mathrm{1, 2}, \cdots , c $

(18)

其中：$ {\mu }_{ik}^{\text{'}} $和$ {x}_{i}^{\text{'}} $分别为利用交替引导滤波修正后的隶属度和像素值。

在此多目标模型中，$ {J}_{m}(\boldsymbol{U}, \boldsymbol{V}) $是KWFLICM算法的目标函数，通过初始化隶属度和聚类中心可得到算法隶属度和聚类中心的迭代更新公式，将隶属度分别传入像素引导隶属度滤波和隶属度引导像素滤波的代价函数$ E(\boldsymbol{U}, \boldsymbol{p}, \boldsymbol{q}) $和$ {E}_{1}(\boldsymbol{U}, {\boldsymbol{p}}_{1}, {\boldsymbol{q}}_{1}) $，求得像素引导隶属度滤波的线性系数$ {p}_{ik} $、$ {q}_{ik} $和隶属度引导像素滤波的线性系数$ {p}_{1i} $、$ {q}_{1i} $，分别将线性系数$ {p}_{ik} $、$ {q}_{ik} $传入隶属度回归模型$ L\left({\boldsymbol{\mu }}^{\text{'}}\right) $，$ {p}_{1i} $、$ {q}_{1i} $传入像素值回归模型$ {L}_{1}\left({\boldsymbol{x}}^{\text{'}}\right) $，利用最小二乘法可得到新的隶属度值$ {\mu }_{ik}^{\text{'}} $和新的像素值$ {x}_{i}^{\text{'}} $，利用新的隶属度值$ {\mu }_{ik}^{\text{'}} $对聚类中心进行更新，最终将新的隶属度和像素值传入KWFLICM算法，循环此过程直到算法收敛。

多目标函数的优化过程具体如下：

1）在KWFLICM算法中，邻域位置信息和灰度信息被引入模糊因子$ {G}_{ik}^{\text{'}} $，利用最小二乘法得到隶属度$ {\mu }_{ik} $和聚类中心$ {v}_{k} $的更新公式如式（8）和式（9）所示。

2）KWFLICM算法得到的隶属度被传入像素引导隶属度滤波的目标函数$ E(\boldsymbol{U}, \boldsymbol{p}, \boldsymbol{q})=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}\sum\limits _{j\in {w}_{i}}\left[\right({p}_{ik}{I}_{j}+{q}_{ik}-{\mu }_{jk}{)}^{2}+\mathrm{\varepsilon }{p}_{ik}^{2}] $，目标函数分别对$ {p}_{ik} $和$ {q}_{ik} $求偏导，通过最小二乘法得到线性系数$ {p}_{ik} $和$ {q}_{ik} $的计算公式如式（13）和式（11）所示。

3）利用隶属度回归目标函数$ L\left({\boldsymbol{\mu }}^{\text{'}}\right)=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}\sum\limits _{j\in {w}_{i}}({p}_{jk}{x}_{i}+{q}_{jk}-{\mu }_{ik}^{\text{'}}{)}^{2} $对$ {\mu }_{ik}^{\text{'}} $求偏导并令其等于0，得到：

$ {\mu }_{ik}^{\text{'}}=\sum\limits _{j\in {w}_{i}}({p}_{jk}{x}_{i}+{q}_{jk})={\stackrel{-}{p}}_{ik}{x}_{i}+{\stackrel{-}{q}}_{ik} $

(19)

其中：$ {\stackrel{-}{p}}_{ik}=\sum\limits _{j\in {w}_{i}}{p}_{jk}/\left|{n}_{w}\right| $；$ {\stackrel{-}{q}}_{ik}=\sum\limits _{j\in {w}_{i}}{q}_{jk}/\left|{n}_{w}\right| $。

4）KWFLICM算法得到的隶属度被传入隶属度引导像素滤波的目标函数$ E1(\boldsymbol{U}, {\boldsymbol{p}}_{1}, {\boldsymbol{q}}_{1})=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}\sum\limits _{j\in {w}_{i}}\left[\right({p}_{1i}{\mu }_{jk}+{q}_{1i}-{I}_{j}{)}^{2}+\varepsilon {p}_{1i}^{2}] $，目标函数分别对$ {p}_{1i} $和$ {q}_{1i} $求偏导，通过最小二乘法得到线性系数$ {p}_{1i} $和$ {q}_{1i} $的计算公式如式（17）和式（15）所示。

5）像素值回归目标函数$ {L}_{1}\left({\boldsymbol{x}}^{\text{'}}\right)=\sum\limits _{i=1}^{n}\sum\limits _{k=1}^{c}\sum\limits _{j\in {w}_{i}}({p}_{1j}{\mu }_{ik}+ $ $ {q}_{1j}-{x}_{i}^{\text{'}}{)}^{2} $对$ {x}_{i}^{\text{'}} $求偏导并令其等于0，得到：

$ {x}_{i}^{\text{'}}=\sum\limits _{k=1}^{c}\sum\limits _{j\in {w}_{i}}({p}_{1j}{\mu }_{ik}+{q}_{1j})=\sum\limits _{k=1}^{c}{\stackrel{-}{p}}_{1i}{\mu }_{ik}+{\stackrel{-}{q}}_{1i} $

(20)

其中：$ {\stackrel{-}{p}}_{1i}=\sum\limits _{j\in {w}_{i}}{p}_{1j}/\left|{n}_{w}\right| $；$ {\stackrel{-}{q}}_{1i}=\sum\limits _{j\in {w}_{i}}{q}_{1j}/\left|{n}_{w}\right| $。

2.3.3 算法步骤

_GFKWFLICM_GF算法步骤具体如下：

1）输入原始噪声图$ \boldsymbol{X} $，初始化隶属度$ {\mu }_{ik}^{\left(t\right)} $和聚类中心$ {v}_{k}^{\left(t\right)} $，迭代计数器$ t=0 $，最大迭代次数$ {t}_{\mathrm{m}\mathrm{a}\mathrm{x}} $和迭代终止阈值$ \sigma $。

2）$ t=t+1 $。

3）根据初始化的隶属度和聚类中心，得到KWFLICM算法隶属度的更新公式$ {\mu }_{ik}^{\left(t\right)} $如式（8）所示。

4）$ \boldsymbol{X} $是原始噪声图，$ \boldsymbol{I} $是引导图，令$ \boldsymbol{I}=\boldsymbol{X} $，将隶属度$ {\mu }_{ik}^{\left(t\right)} $传入像素引导隶属度滤波的目标函数，根据式（13）和式（11）可得线性系数$ {p}_{ik} $和$ {q}_{ik} $，然后根据式（19）可得到新的隶属度$ {\mu }_{ik}^{\text{'}} $。

5）将隶属度$ {\mu }_{ik}^{\left(t\right)} $传入隶属度引导像素滤波的目标函数，根据式（17）和式（15）可得线性系数$ {p}_{1i} $和$ {q}_{1i} $，然后根据式（20）可得新的像素值$ {x}_{i}^{\text{'}} $。

6）通过新的隶属度，利用式（9）对聚类中心进行更新。

7）将新的像素值$ {x}_{i}^{\text{'}} $和新的隶属度$ {\mu }_{ik}^{\text{'}} $重新传入KWFLICM算法。

8）当$ t > {t}_{\mathrm{m}\mathrm{a}\mathrm{x}} $或$ \left|\right|{v}_{k}^{\left(t\right)}-{v}_{k}^{(t-1)}\left|\right| < \sigma $时，算法迭代结束，否则，转入步骤2。

2.4 收敛性分析

本文算法的本质是通过引导滤波利用邻域像素值对KWFLICM算法的隶属度进行修正，同时利用像素隶属度对像素值进行引导滤波，引导滤波修正后的隶属度和像素值分别为$ {\mu }_{ik}^{\text{'}}={\stackrel{-}{p}}_{ik}{x}_{i}+{\stackrel{-}{q}}_{ik} $和$ {x}_{i}^{\text{'}}= $ $ \sum\limits _{k=1}^{c}{\stackrel{-}{p}}_{1i}{\mu }_{ik}+{\stackrel{-}{q}}_{1i} $，其中，$ {\stackrel{-}{p}}_{ik} $和$ {\stackrel{-}{q}}_{ik} $为引导滤波利用像素引导隶属度滤波的邻域线性系数均值，$ {\stackrel{-}{p}}_{1i} $和$ {\stackrel{-}{q}}_{1i} $为引导滤波利用隶属度引导像素滤波的邻域线性系数均值，以此来提高算法的鲁棒性。文献[30]利用邻域像素隶属度对当前像素聚类所得到的隶属度进行修正，再调整聚类中心，以此来提高FCM分割算法的鲁棒性。文献[31]提出条件空间约束的模糊聚类算法并用于MRI图像的分割，其本质仍是利用邻域像素隶属度信息对当前像素的隶属度进行约束。文献[32]提出基于组合隶属度的快速模糊聚类算法，通过构造组合隶属度函数，对迭代中的隶属度直接滤波处理，提高了算法的实用性。本文算法的基本思想与文献[30-32]相同，且通过实验分析也证实了本文算法是收敛的。

3 实验结果

MRI图像的质量对临床医学的诊断有着重要意义，是由高斯噪声污染的实部图像和虚部图像经取模运算而得到的模图像，因此最终得到的图像是受莱斯噪声污染的^[33]。MRI图像分布与信号有关，当图像的信噪比较大时莱斯分布趋向于高斯分布，当图像的信噪比较小时莱斯分布呈现瑞利分布，这使得去除莱斯噪声更加困难。因此，本文对受莱斯干扰的合成图、MRI图像及偏置图进行分割测试。

采用误分率（Misclassified Error，ME）、精确度（Accuracy，ACC）、峰值信噪比（Peak Signal to Noise Ratio，PSNR）、划分系数（Partition Coefficient，PC）、划分熵（Partition Entropy，PE）和Jaccard相似（Jaccard Similarity，JS）系数作为算法分割的评价指标。将像素引导隶属度滤波和隶属度引导像素滤波的交替引导滤波与KWFLICM算法相结合。为了验证本文算法的优越性，对基于交替引导滤波的核模糊聚类算法和传统核模糊聚类算法的分割结果进行比较。对合成图、MRI图像和偏置图分别添加不同强度的莱斯噪声进行测试。

3.1 莱斯噪声干扰的合成图、MRI图、偏置图测试

对合成图添加均方差为45的莱斯噪声，采用KFCM、KFCM_S1、KFCM_S2、KFLICM、KWFLICM算法及它们应用交替引导滤波进行改进的算法对图像进行分割，所得分割结果和评价指标分别如图 2及表 1所示，其中最优指标值用加粗字体标示。

对2幅灰度MRI图像分别添加均方差为40和55的莱斯噪声，采用KFCM、KFCM_S1、KFCM_S2、KFLICM、KWFLICM算法及它们采用交替引导滤波进行改进的算法对图像进行分割，所得分割结果和评价指标分别如图 3、图 4及表 2所示，其中最优指标值用加粗字体标示。

对2幅灰度偏置图（灰度设置为20%和0%）分别添加均方差为40和40的莱斯噪声，采用KFCM、KFCM_S1、KFCM_S2、KFLICM、KWFLICM算法及它们采用交替引导滤波进行改进的算法对图像进行分割，所得分割结果和评价指标分别如图 5、图 6及表 3所示，其中最优指标值用加粗字体标示。

3.2 实验结果分析

从图 2~图 6可以看出：对于受莱斯噪声干扰的图像，KFCM算法忽略了像素的空间信息，导致分割结果不理想，分割结果中含有大量噪声点；KFCM_S1、KFCM_S2、KFLICM、KWFLICM算法的分割结果中存在的噪声点较少，在一定程度上增强了算法的抗噪性、提高了算法的分割性能，其中KFCM_S2引入邻域中值信息，抑制了原图像中的大部分莱斯噪声，KFLICM算法虽然引入了具有空间邻域信息的局部模糊因子，但分割结果却不理想，KWFLICM算法在KFLICM算法的基础上改进了局部模糊因子，分割结果较好，但仍含有部分噪声，在引入交替引导滤波之后，相比于传统核模糊聚类算法，分割结果中的噪声点明显减少，算法评价指标也明显变好；将KWFLICM算法和交替引导滤波相结合的本文算法获得的图像几乎不含噪声且细节清晰，分割效果较好。

从表 1~表 3可以看出：本文算法在莱斯噪声干扰下的误分率、精确度、峰值信噪比、划分系数、划分熵、Jaccard相似系数等评价指标上表现更优，表明本文算法在莱斯噪声抑制方面具有一定的优势。

4 噪声干扰对算法性能的影响分析

为了验证交替引导滤波的优越性，对MRI图 1分别添加不同强度的莱斯噪声，测试不同算法在莱斯噪声干扰下的误分率、精确度、峰值信噪比、划分系数、划分熵、Jaccard相似系数的曲线变化，如图 7所示。

从图 7可以看出，分割结果的误分率和划分熵会随着噪声强度的增加而增加，精确度、峰值信噪比、划分系数、Jaccard相似系数会随着噪声强度的增加而减小。在核空间内，基于交替引导滤波的模糊聚类算法的噪声变化曲线图均优于传统核模糊聚类算法的噪声变化曲线图。本文算法将引入空间位置信息与邻域灰度信息的KWFLICM算法与交替引导滤波相结合，随着噪声强度的增加，除个别评价指标外，其余各项评价指标均优于其他算法，且噪声变化曲线能保持相对稳定性，大幅提升了算法的性能，无论是从视觉效果还是定量评价指标来看，本文算法在抗莱斯噪声方面具有明显优势。

5 算法复杂度分析

为比较上述算法的运行效率，对不同算法的计算复杂度进行分析。对$ n $个像素且可分为$ c $类的图像进行分割，$ {w}^{2} $代表像素的邻域窗口大小，$ t $代表迭代次数。FCM算法的计算复杂度来自算法的迭代复杂度$ O\left(nct\right) $；FCM_S1或FCM_S2算法的计算复杂度主要由两部分组成：一部分来自每个像素的均值或中值信息的计算复杂度$ O\left(n{w}^{2}\right) $，另一部分为算法迭代复杂度$ O\left(nct\right) $；FLICM算法的复杂度为算法迭代复杂度$ O\left(nct{w}^{2}\right) $；KWFLICM算法在FLICM算法的基础上改进了模糊因子，计算复杂度为$ O(nct{w}^{2}+n{(w+1)}^{2}) $。在引入双重引导滤波后，算法的计算复杂度由引导滤波器的计算复杂度$ O\left(nct\right) $和算法本身的计算复杂度组成。不同算法的计算复杂度如表 4所示。为了直观地展示本文算法的运行效率，分别对受莱斯噪声污染的MRI图 1测试所有算法的时间开销。不同算法的时间开销柱状图如图 8所示。

6 结束语

为实现基于引导滤波的模糊聚类算法在图像分割中的应用并降低图像中噪声的敏感性，本文提出一种基于交替引导滤波的鲁棒模糊聚类算法。该算法将像素引导隶属度和隶属度引导像素的交替引导滤波与KWFLICM算法相结合，有效抑制了图像中的噪声，进一步提升了聚类性能。实验结果表明，本文算法对于莱斯噪声图像的分割效果明显优于传统核模糊聚类算法，证明了其可对噪声环境下的图像进行有效分割。后续将结合并行紧致差分进化方法、Spark平台和GPU系统进一步加快本文聚类算法的运行速度，以满足实时图像处理需求。