0 引言

目前，很多科学研究文献已经提出了多种图像分割方法。其中，图像聚类方法主要受到数据分组初始安排规则影响。近些年来的主要科学研究成果均集中在图形理论方法研究、基于均值漂移(Mean Shift，MS)的相关图像分割算法和Rate-distortion理论方法^[1]。

针对有限混合模型(Finite Mixture Model, FMM)的概率密度函数(Probability Density Function, PDF)的像素点相关属性(如：强度、纹理等属性)建模到群体数据上是自然形式，因为它会根据组件自动提供产生分组混合结构。此外，针对聚类性能来说，FMM概率密度函数是米制单位^[2]。基于FMM的概率密度函数建模方法已经成功应用在生物信息学^[3]、图像检索^[4]等相关领域。FMM模型参数能够通过极大似然估计融合期望最大化算法求得^[5]。

图像先验属性针对强化空间平滑程度操作有着重要意义，而强化空间平滑度是图像处理应用的关键问题^[6]。图像处理应用包括图像恢复、图像去噪、图像分割、图像优化等各种问题。在常见概率密度函数框架中，图像平滑操作是针对图像先验特征的具体方法。

本文提出的模型与现有方法有区别。首先，文本、场景和目标分类均为有监督学习问题且相关分割方法也是无监督性质;同时，在现有科研成果来看，针对隐狄利克雷分布(Latent Dirichlet Allocation, LDA)参数估计通常是由变分推理或者简化Logistic模型来完成。本文提出的改进模型优点是在E-step步骤中可表示在密闭形式上，也明确说明本文模型假设改进空间约束贝叶斯网络模型的概率向量具体比例。通过期望最大化(Expectation Maximization, EM)算法的推理将执行LDA参数的多项式方程，因此，密闭形式M-step步骤所得参数能够满足概率约束所需条件。

1 空间变化有限混合模型

假设xⁱ=(x₁ⁱ, x₂ⁱ, …, x_Lⁱ)^T具有向量特征(例如：强度特征、文理特征)，xⁱ代表第i个像素点具体位置，i=1, 2, …, N，将d维向量图像建模成独立同分布随机向量。空间变化有限混合模型(Spatially Variant Finite Mixture Model, SVFMM)^[7]提供改进FMM方法^[5]的像素点标记。它假设存在K个分量的每个具有θ_j密度的函数来定义参数向量混合模型。

像素点i的特征概率向量是πⁱ=(π₁ⁱ, π₂ⁱ, …, π_Kⁱ)^T。其中，K是成员数量。定义Π={(π¹)^T, (π²)^T, …, (π^N)^T}作为概率矢量集合和Θ={θ₁, θ₂, …, θ_K}作为成员参数。参数π_jⁱ(j=1, 2, …, K)是空间变化有限混合模型对于每个像素点表示第 i个像素点概率，属于第j类必须满足约束条件如下：

$\sum\limits_{j=1}^{K}{\pi _{j}^{i}=1};0\le \pi _{j}^{i}\le 1,i=1,2,\cdots ,N,j=1,2,\cdots ,K$

(1)

假定像素点位置xⁱ的标准有限混合模型概率密度函数可以表示为式(2) 形式:

$f\left( {{\mathit{\boldsymbol{x}}^i}|\mathit{\Pi };\mathit{\Theta }} \right) = \sum\limits_{j = 1}^K {\pi _j^i\phi \left( {{\mathit{\boldsymbol{x}}^i}|{\theta _j}} \right)} $

(2)

φ(xⁱ|θ_j)的高斯分布参数是θ_j = {μ_j , Σ_j}，其中μ_j=(μ_{j, 1}, μ_{j, 2}, …, μ_{j, L})是平均矢量，Σ_j是L维高斯分布的协方差矩阵，Π可定义为随机变量和 Θ参数。

空间变化有限混合模型使用先验密度分布p(Π)的随机变量Π，因此，X表示该组像素点的特征向量{xⁱ}，当i=1, 2，…, N时，本文假设像素点是独立统计，应用贝叶斯模型条件，将所得后验概率密度函数由式(3) 计算：

$p\left( \mathit{\Pi }|X;\mathit{\Theta } \right)\propto \prod\limits_{i=1}^{N}{p\left( \mathit{\Pi } \right)f\left( {{\mathit{\boldsymbol{x}}}^{i}}|\mathit{\Pi },\mathit{\Theta } \right)}$

(3)

根据密度对数函数，可推导出式(4) :

$L\left( \mathit{\Pi }|X;\mathit{\Theta } \right)=\sum\limits_{i=1}^{N}{\left\{ \lg f\left( {{\mathit{\boldsymbol{x}}}^{i}}|\mathit{\Pi },\mathit{\Theta } \right)+\lg \left( p\left( \mathit{\Pi } \right) \right) \right\}}$

(4)

p(Π)是Gauss-Markov随机域^[8]表达的典型实例，具体形式如式(5) 所示:

$p\left( \mathit{\Pi } \right)\propto \prod\limits_{j=1}^{K}{\beta _{j}^{-N}\exp \left[ -\frac{1}{2}\frac{\sum\limits_{i=1}^{N}{\sum\limits_{m\in {{N}_{i}}}{{{\left( \pi _{j}^{i}-\pi _{j}^{m} \right)}^{2}}}}}{\beta _{j}^{2}} \right]}$

(5)

因为参数β_j可获取集群数据的空间平滑性和执行不同平滑度，所以在每个集群数据上都找到j以便适应数据模型。

SVFMM空间变化模型具体如图 1所示。在标准FMM模型中，对于给定的像素点特征向量x取决于离散隐变量z来表示混合部分，最终生成特征向量x。如果z_j=1，像素点 x则属于j类。在这种情况下，对于给定类混合比例属于该类别像素点百分比。在SVFMM情况下，每个像素点i都有自身固定混合比例πⁱ，称为像素点标签概率。这些上下文混合比例均由平滑先验知识来执行空间约束。

最大后验概率(Maximum A Posteriori, MAP)估计模型参数估计的EM算法^[9]需要在E-step迭代步长t计算隐变量的条件期望值，具体如式(6) 所示：

$z_{j}^{i\left( t \right)}=\frac{\pi _{j}^{i\left( t \right)}\phi \left( {{\mathit{\boldsymbol{x}}}^{i}}|\theta _{j}^{\left( t \right)} \right)}{\sum\limits_{p=1}^{K}{\pi _{j}^{i\left( t \right)}\phi \left( {{\mathit{\boldsymbol{x}}}^{i}}|\theta _{p}^{\left( t \right)} \right)}}$

(6)

在M-step步骤中，考虑到完整数据是隐变量的线性似然对数，那么完整数据的最大化对数似然估计模型参数如式(7) 所示:

$\begin{align} & Q\left( \mathit{\Pi },\mathit{\Theta }|{{\mathit{\Pi }}^{\left( t \right)}},{{\mathit{\Theta }}^{\left( t \right)}} \right)= \\ & \sum\limits_{i=1}^{N}{\sum\limits_{j=1}^{K}{z_{j}^{i\left( t \right)}\left\{ \lg \left( \pi _{j}^{i} \right)+\lg \left( \phi \left( {{x}^{i}}|{{\theta }^{j}} \right) \right) \right\}+\lg \left( p\left( \mathit{\Pi } \right) \right)}} \\ \end{align}$

(7)

式(7) 中的函数Q(·)能够针对每个参数执行独立最大化操作，并提供以下在t+1步骤的混合模型参数。

$\left\{ \begin{align} & \mu _{j,l}^{\left( t+1 \right)}=\left( \sum\limits_{i=1}^{N}{z_{j}^{i\left( t \right)}x_{l}^{i\left( t \right)}} \right)/\left( \sum\limits_{i=1}^{N}{z_{j}^{i\left( t \right)}} \right) \\ & \sigma _{j,l}^{\left( t+1 \right)}=\sqrt{\left( \sum\limits_{i=1}^{N}{z_{j}^{i\left( t \right)}{{\left[ x_{l}^{i\left( t \right)}-\mu _{j,l}^{\left( t+1 \right)} \right]}^{2}}} \right)/\left( \sum\limits_{i=1}^{N}{z_{j}^{i\left( t \right)}} \right)} \\ \end{align} \right.$

(8)

通过设定$\frac{\partial Q}{\partial \pi _{j}^{i}}=0$来生成式(9) ，针对π_jⁱ计算得到概率π_i^j。

$\left| {{N}_{i}} \right|{{\left( \pi _{j}^{i} \right)}^{2}}-\left[ \sum\limits_{m\in {{N}_{i}}}{\pi _{j}^{m}} \right]\left( \pi _{j}^{i} \right)-z_{j}^{i}\beta _{j}^{2}=0$

(9)

其中：|N_i|是在邻近第i个像素点的像素点数量。

这样通过验证式(9) 始终是非负解为π_i^j，然而，SVFMM模型缺点是没有明确考虑πⁱ，这是概率向量(π_jⁱ≥0，$\sum\limits_{j=1}^{K}{\pi _{j}^{i}=1,}i=1,2,\cdots ,N$)的空间平滑条件。为了达到目标，在M-step步骤中引入补偿计算来执行非负解π_i^j，以满足这些约束条件变量。文献[7]的梯度投影算法和文献[9]的二次规划算法都类似于此方法。此方法可提高式(7) 和模型性能。

2 改进空间约束贝叶斯网络模型

为了克服SVFMM算法局限性，本章提出改进空间约束贝叶斯网络模型。本改进空间约束贝叶斯网络模型依据LDA，基于上层分布混合结构来提出图像分割上下文混合模型Π。LDA是多项式形式，概率向量πⁱ的相关参数是由LDA分布产生^[5]。相似先验知识已经在改进空间的上下文约束条件中提出，其中LDA参数估计通过迭代梯度下降来进行方法优化^[7]。此外，空间平滑通过隐Dirichlet分布参数的方程形式以密闭形式计算实施具有真正非负解。

假设生成图像模型为了产生第i个像素点以求达到第j个组件目的，实现LDA过程，因此，式(6) 的隐变量zⁱ具有第j个分量。在这种情况下，后验概率混合模型结构可表述为式(10) 形式：

$\begin{matrix} \left\{ \begin{align} & \pi _{j}^{i}=p\left( z_{j}^{i}=1|{{\mathit{\boldsymbol{x}}}^{i}} \right)=1 \\ & \pi _{m}^{i}=p\left( z_{m}^{i}=1|{{\mathit{\boldsymbol{x}}}^{i}} \right)=0 \\ \end{align} \right.; & m=1,\cdots ,K,m\ne j \\ \end{matrix}$

(10)

考虑LDA实现过程(M=1) 与Γ(x+1) =xΓ(x)，式(10) 中标签的第i个像素点概率可改写成为式(11) ：

$\pi _{j}^{i}=p\left( z_{j}^{i}=1|{{\mathit{\boldsymbol{x}}}^{i}} \right)=\alpha _{j}^{i}/\sum\limits_{l=1}^{K}{\alpha _{l}^{i}},j=1,2,\cdots ,K$

(11)

新型模型有可能通过引入LDA分布参数A达到空间变化目标。假设高斯-马尔可夫随机域来估计闭合形式具体参数^[8]。

$p\left( A \right)\propto \prod\limits_{j=1}^{K}{\beta _{j}^{-N}\exp \left[ -\frac{1}{2}\frac{\sum\limits_{i=1}^{N}{\sum\limits_{m\in {{N}_{i}}}{{{\left( \alpha _{j}^{i}-\alpha _{j}^{m} \right)}^{2}}}}}{\beta _{j}^{2}} \right]}$

(12)

已有先验知识主要特征是为了提供更好的先期适应数据，强制每个数据簇中拥有不同程度的平滑程度。

图 2代表这种分层办法提出的图形模型。本文参考LDA分布的空间变化有限混合模型(LDA-SVFMM)模型，LDA-SVFMM产生图像模型工作原理如下：首先产生样本ξ(概率向量)使用LDA分布相关参数，从而获得多项式分布参数ξ。隐变量z表示观察点类的x变量，它是用ξ参数进行多项式计算所得结果。LDA分布空间参数的约束平滑条件α需要根据标准化SVFMM算法执行。

3 MAP-EMM方法

本文模型主要思想是在最大后验概率(Maximum A Prior, MAP)算法上应用EM方法。应用式(12) 中的Gauss-Markov参数A产生下面所提MAP函数，以最大化EM算法的M-step目标。关于参数A，本文给出定义如下：

$\begin{align} & Q\left( A,\mathit{\Theta }|{{A}^{\left( t \right)}},{{\mathit{\Theta }}^{\left( t \right)}} \right)=\sum\limits_{i=1}^{N}{\sum\limits_{j=1}^{K}{\left\{ z_{j}^{i}\left[ \lg \left( \alpha _{j}^{i}/{{A}^{i}} \right)+ \right. \right.}} \\ & \lg \left( \phi \left( {{x}^{i}}|{{\theta }^{j}} \right) \right)-\lg \left( \beta _{j}^{N} \right)-\frac{1}{2}\sum\limits_{m\in {{N}_{i}}}{{{\left( \alpha _{j}^{i}-\alpha _{j}^{m} \right)}^{2}}}/\beta _{j}^{2} \\ \end{align}$

(13)

其中：${{A}^{i}}=\sum\limits_{j=1}^{K}{\alpha _{j}^{i}}$。为了计算模型参数α_j^i(t+1) ，本文在M-step中最大化式(13) 从而达到提高α_jⁱ参数的目的。计算α_jⁱ 的偏导数并将结果设置为0。设定邻近像素点数量|N_i|=8和$\frac{\partial Q}{\partial \alpha _{j}^{i}}=0$，本文给出多项式方程从而提高α_j^i(t+1) 参数以求达到i=1, 2，…, N和j=1, 2，…, K的要求。

$\begin{align} & {{\left( \alpha _{j}^{i} \right)}^{3}}+\left[ A_{-j}^{i}-\frac{\sum\limits_{m\in {{N}_{i}}}{\alpha _{j}^{m}}}{\left| {{N}_{i}} \right|} \right]{{\left( \alpha _{j}^{i} \right)}^{2}}+\frac{A_{-j}^{i}\sum\limits_{m\in {{N}_{i}}}{\alpha _{j}^{m}}}{\left| {{N}_{i}} \right|}\left( \alpha _{j}^{i} \right)- \\ & \frac{z_{j}^{i}{{A}^{i}}-j\beta _{j}^{2}}{2\left| {{N}_{i}} \right|}=0 \\ \end{align}$

(14)

最后，通过$\frac{\partial Q}{\partial \beta _{j}^{2}}=0$来解决β_j²参数的时间步长为(t+1) 步时的相关解决方案。

$\beta _{j}^{2}=\frac{1}{N}\sum\limits_{i=1}^{N}{\sum\limits_{m\in {{N}_{i}}}{{{\left( \alpha _{j}^{i}-\alpha _{j}^{m} \right)}^{2}}}},j=1,\cdots ,K$

(15)

本文MAP-EMM算法的具体过程如下所述：

步骤1 初始化算法参数和LDA分布向量参数αⁱ。

步骤2 应用式(13) 计算MAP函数：

步骤2.1 (E-Step)：应用式(13) 的后验概率z_j^i(t)，计算第i个像素点是否属于第j类图像。

步骤2.2 (M-Step)：

步骤2.2.1 应用式(8) 计算改进空间约束贝叶斯网络模型参数；

步骤2.2.2 应用式(14) 所得非负解针对LDA分布参数进行替换；

步骤2.2.3 通过像素点标签替换像素点已有概率参数；

步骤2.2.4 应用式(15) 计算像素点标签参数。

步骤3 直至MAP函数收缩至无穷小范围，算法自动达到结束条件。

4 实验结果与分析

本文实验环境为Matlab 2014b，图像大小是256×256。本文针对MAP-EMM改进算法的迭代次数和初始化随机生成条件，确保提供对数似然函数的最大值。这里考虑EM算法结束条件准则是收敛定义成在式(4) 的对数似然变化的百分比在两个连续迭代次数之间小于0.001％。

为了验证MAP-EMM改进算法提出的先验LDA的必要性，针对SVFMM^[7]提出在不同复杂程度下的相同分段图像分割的具体实例。经过实验证明，在人工图像和自然图像中，高斯噪声对于MAP-EMM改进算法的鲁棒性影响很小。

本文针对MAP-EMM改进算法在Berkeley图像数据库的300张图像中进行图像分割^[6]。本文MAP-EMM改进算法与JSEG^[10]、MAP-ML(MM)^[3]、CTM^[5]进行效果比较。

应用标准滤波器进行特征描述，在最近文献中的纹理描述主要方法包括Blobworld特征^[11]和MRF特征^[12]。Blobworld特征通过生成六维颜色矢量和纹理信息数据的过程关键特点是正确估计纹理规模。MRF特征只应用于PCA进行降维操作，其次是像素点终止窗口的向量化操作。

图 3(a)~(d)和图 4(a)~(d)展示本文方法和其他三种方法分割300幅Berkeley图像数据库部分图像的结果比较。从中可得，本文方法分割结果中，噪声区域比较少，边界保持较好；其次，从本文方法和MM分割结果的比较可见，两种方法都采用了图像优化方法，分割结果的边界保持较好(如图 3中MM的分割结果第1~3及第5行)，但图像的标签数量对分割结果影响较大，尽管MM在分割过程中能够根据能力的变换自动调整标签数，但容易造成过分割或欠分割。

对比其他方法，JSTG能够得到较为同质的区域，但已得到过分割结果，并且不能很好地区分视觉差异不明显区域。CTM采用基于超像素的区域合并策略对图像进行分割，从图 3和图 4可见，其分割结果边界不光滑、错位，采用的最小描述长度准则并不能较好地适应Berkeley数据库300幅图像，造成过分割或欠分割现象。

为了更好评价各比较方法的分割性能，采用四个常用评价指标函数：PRI(Probabilistic Rand Index)、VoI(Variation of Information)、GCE(Global Consistency Error)和BDE(Boundary Displacement Error)对分割结果进行评价。其中:PRI是统计机器分割和多个人工分割之间标签一致的像素对的个数占整个像素对个数的比率；VoI则把机器分割和人工分割之间的距离定义为在给定人工分割的条件下机器分割的平均条件熵，它能够测量机器分割中不能被人工分割所解释的随机性的量；GCE测量一个分割可被看作为另外一个分割的程度；BDE则是测量两个分割结果中边界像素的平均位移误差。量化结果中PRI值越大，VoI、GCE和BED值越小，则机器分割结果与人工分割结果越接近。

表 1给出了4种方法分割300幅图像的结果在评价指标上的量化分析。可见，本文方法在PRI、VoI、GCE指标上优于其他3种方法，在BDE指标上仅次于CTM，相对于JSEG、MM和CTM，本文方法的分割结果更加接近于人工分割结果。

5 结语

本文根据贝叶斯理论和空间分层建模约束混合模型提出MAP-EMM改进算法。本文模型应用LDA概率密度模型和高斯—马尔可夫定理的随机域参数混合过程来实现参数平滑。本文方法根据空间信息先验平滑变换操作，在待处理像素点的上下文混合结构中引入LDA符合多项式分布，从而用来替换传统期望最大化算法中映射操作。本文仿真结果对比图显示空间变化的混合模型性能方面比采用投影步骤的混合结构约束方法有很大改进。