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

0 概述

图像分割是图像处理领域的重要研究方向。常用的图像分割方法主要有基于阈值的分割方法、基于聚类的图像分割方法和基于模型的分割方法。近年来，基于模型的分割方法^[1]成为了研究热点，其中高斯混合模型(Gaussian Mixture Model，GMM)因其建模简单有效而受到广泛关注。而基于GMM图像分割方法通常假设图像中的像素点相互独立，没有考虑图像像素间的空间约束关系，导致其分割准确性欠佳。为解决该问题，研究人员通过对马尔科夫随机场(Markov Random Field，MRF)模型引入图像的空间约束关系，建立的高斯马尔科夫随机场(Gaussian MRF，GMRF)可有效提升图像分割效果，并且在图像分割、图像恢复以及图像匹配等领域得到广泛应用。然而，MRF模型经常使用Potts类型先验知识建立图像的空间约束关系，这种简单的像素点对关系不能准确描述复杂图像的先验知识，且在分割过程中易产生过平滑现象。

为提高MRF模型对复杂图像先验知识的描述能力，需引入更多邻域节点来表达图像的空间一致性约束。文献[2]针对传统GMM难以自动获取类属数和对噪声敏感的问题，提出一种基于可变类空间约束GMM的遥感图像分割方法。利用MRF模型将GMM的像素类属性建模为图像的先验概率，以提高GMM的抗噪性。文献[3]提出一种局部空间自适应MRF模型，通过在带权值的自适应空间GMM中引入像素的局部空间信息，使得图像的噪声或者纹理突变具有一定的鲁棒性，但由于区域模型的平滑作用，将导致分割结果出现边缘带现象。文献[4]提出一种快速鲁棒的改进GMM(Modified GMM，MGMM)图像分割模型，该模型利用欧氏距离度量局部窗内像素间的相似性并引入图像局部空间信息。为保持图像的边缘结构，该模型设定了较小的MRF模型先验参数$ \beta $，但是对于自然图像而言，其全局非同态性使得具有较小先验惩罚的MRF模型对复杂图像强度信息的突变较为敏感，不能得到更平滑的分割结果。文献[5]提出基于图像片权重MRF分割模型，对邻域内不同的图像片根据相似度计算对应权重，同时通过K-L(Kullback-Leibler)距离引入熵的惩罚，从而解决过分割问题。文献[6]提出基于改进高斯混合模型的运动图像目标检测算法，将邻域像素值加权和作为像素重塑后的值，并利用像素空间距离建立滤波窗口的各像素权值，以有效抑制噪声并且更好地保留图像细节。文献[7]提出空间约束非对称GMM，采用K-L距离计算邻域像素先验分布，引入了图像空间信息，并根据像素类内和类间先验知识和后验概率建立图像似然特征的非对称分布模型，有效拟制噪声干扰。

上述模型利用局部空间信息提高了分割模型对噪声或纹理突变干扰的鲁棒性，但由于这些模型没有考虑到图像的边缘信息，因此局部区域能量的最小化易导致分割结果过平滑，使得图像分割结果常出现边缘模糊、边界不清的问题。为解决该问题，研究人员将图像边缘信息融合到MRF模型中，以提升图像分割的边缘效果。文献[8]基于隐马尔科夫随机场(Hidden MRF，HMRF)理论提出了HMRF-EM-edge模型，利用邻域像素相互作用引入空间信息，并使用图像边缘特征来保持分割边缘，但对于图像中边缘细小的区域，分割结果依然不理想。文献[9]提出一种结合形状信息的静态MST区域划分和RHMRF-FCM算法的图像分割方法，并定义一种MST静态同质区域划分准则，通过借助MST引入边界信息，可较好地表达图像不规则边界和形状信息，有效抑制几何噪声。文献[10]提出一种结合MRF能量和活动轮廓模型的方法，利用Garbo纹理特征、DCE-MRI时域特征和灰度特征构成的特征向量与聚类中心向量的距离构建一种模糊速度函数，并将其引入到活动轮廓模型中作为MRF能量模型的边缘检测项。文献[11]针对多尺度马尔科夫随机场(MRMRF)模型非重叠区域在最优化过程中造成的块效应，提出一种具有边缘保持的MRMRF模型图像分割方法，有效抑制MRF能量模型在分割过程中向图像的边缘收敛，但这种Cauchy分布对于图像弱边缘的惩罚较大，不利于保持图像的细节边缘特征。

为有效表达复杂图像局部统计特征并保持图像分割区域的边缘结构，本文提出一种基于边缘约束局部区域MRF(Edge Constrained Local Region MRF，ECLRMRF)的分割模型。该模型先利用图像局部区域内像素的欧氏距离度量像素的相似度，构建局部空间约束高斯混合模型(Local Spatial Constrained Gaussian Mixture Model，LSCGMM)，以提取图像的局部区域似然特征。为约束图像分割区域的边缘，利用Canny边缘检测算子提取图像边缘特征和图像分割区域的边缘结构，并在MRF模型框架下将图像的局部区域似然特征和边缘特征相融合，建立具有边缘约束的局部区域MRF分割能量模型，再采用Gibbs算法优化MRF模型并得到最终分割结果。

1 基于GMRF的图像分割模型

假设$ S=\left\{\mathrm{1, 2}, \cdots , N\right\} $表示一个格点集合，其中N表示集合S中的格点总数。令$ X=\left\{{\mathit{\boldsymbol{x}}}_{s}\right|s\in S\} $表示观测图像，其中$ {\mathit{\boldsymbol{x}}}_{s} $是指图像中第$ s $个像素的强度值。令$ Y=\left\{{y}_{s}\right|s\in S\} $表示分割图像的标签场，$ {y}_{s}\in \left\{\mathrm{1, 2}, \cdots , K\right\} $表示第$ s $个像素的标签，其中$ K $为图像的类别总数，图像中每一个像素$ {\mathit{\boldsymbol{x}}}_{s} $对应唯一的标签$ {y}_{s} $。

在基于模型的图像分割方法中，GMM是描述图像像素分配标签过程的理想模型之一。假定分割过程是一个高斯过程，即图像像素的强度值与其对应标签的概率分布可用GMM概率密度函数描述为：

$ P\left({\mathit{\boldsymbol{x}}}_{s}|{y}_{s}=k, \Theta \right)=\sum\limits_{k=1}^{K}\pi \left({y}_{s}=k\right)p\left({x}_{s}|{\theta }_{k}\right) $

(1)

其中，$ \pi \left({y}_{s}=k\right) $表示给一个像素分配为第$ k $类图像标签的先验概率，且满足约束条件$ 0\le \pi \left({y}_{s}=k\right)\le 1 $和$ \sum\limits_{k=1}^{K}\pi ({y}_{s}=k)=1 $；$ k $为分割类别且$ k\in \left\{\mathrm{1, 2}, \cdots , K\right\} $；$ \Theta =\left\{{\theta }_{1}, {\theta }_{2}, \cdots , {\theta }_{K}\right\} $是GMM参数。$ p\left({x}_{s}|{\theta }_{k}\right) $表示GMM第$ k $类高斯分布的概率密度函数，计算公式如下：

$ p\left({x}_{s}|{\theta }_{k}\right)=\frac{1}{{\left(2\mathrm{\pi }\right)}^{1/2}{\sigma }_{k}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{1}{2}{\left({\mathit{\boldsymbol{x}}}_{s}-{\mathit{\boldsymbol{\mu }}}_{k}\right)}^{\mathrm{T}}{\sigma }_{k}^{-1}\left({\mathit{\boldsymbol{x}}}_{s}-{\mathit{\boldsymbol{\mu }}}_{k}\right)\right] $

(2)

其中，$ {\theta }_{k}=\left({\mathit{\boldsymbol{\mu }}}_{k}, {\sigma }_{k}\right) $分别为第$ k $类高斯分布的均值和方差。

传统GMM通常假定图像像素是独立同分布的，即不考虑邻域像素点的空间约束关系时，由式(1)可得基于GMM的对数似然图像分割模型为：

$ L\left[P\left(X|Y, \Theta \right)\right]=\sum\limits_{s=1}^{N}\mathrm{l}\mathrm{n}\sum\limits_{k=1}^{K}\pi \left({y}_{s}=k\right)p\left({\mathit{\boldsymbol{x}}}_{s}|{\theta }_{k}\right) $

(3)

为了在基于GMM图像分割模型中引入邻域像素的相关性，常利用MRF建模图像先验概率$ \pi \left({y}_{s}=k\right) $刻画邻域像素的空间约束关系。假设$ N\left(s\right)=\left\{{N}_{s}|\forall s\in S\right\} $为定义在图像域$ X $上的一个邻域系统，$ {N}_{s} $为第$ s $个节点的邻域像素集合，且$ s\notin {N}_{s} $。根据Hammersley-Clifford定理^[15]，基于MRF描述的图像空间Gibbs先验分布可表示为：

$ \begin{array}{l}\pi \left({y}_{s}=k\right)\approx P\left({y}_{s}=k|{y}_{{N}_{s}}\mathrm{ }, \beta \right)=\\ \frac{1}{Z\left(\beta \right)}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{1}{T}{E}_{s}\left({y}_{s}=k|{y}_{{N}_{s}}\mathrm{ }, \beta \right)\right]\end{array} $

(4)

其中，$ Z\left(\beta \right) $为归一化常数，$ T $为温度参数，$ \beta $为模型先验参数且$ \beta >0 $，$ {E}_{s}\left({y}_{s}=k|{y}_{{N}_{s}}, \beta \right) $是一个Gibbs先验能量项。根据贝叶斯定理并结合式(1)和式(4)，可得到GMRF后验概率密度联合分布为：

$ \begin{array}{l}P\left(Y|X, \beta , \Theta \right)\propto P\left(X|Y, \Theta \right)P\left(Y|\beta \right)\propto \\ \prod _{s=1}^{N}\prod _{k=1}^{K}p\left({\mathit{\boldsymbol{x}}}_{s}|{y}_{s}\mathrm{ }, {\theta }_{k}\right)P\left({y}_{s}=k|{y}_{{N}_{s}}\mathrm{ }, \beta \right)\end{array} $

(5)

其中，$ P\left(Y|X, \beta , \Theta \right) $为标签场$ Y $在给定观察条件$ X, \beta , \Theta $下的后验概率分布，$ P\left(X|Y, \Theta \right) $为观察图像场$ X $关于标签场$ Y $的似然概率分布，$ \Theta =\left\{{\theta }_{1}, {\theta }_{2}, \cdots , {\theta }_{K}\right\} $为GMM参数，$ P\left(Y|\beta \right) $为标签场$ Y $的先验概率分布，$ \beta $为先验参数，$ k $为类别数。

GMRF模型同时考虑了图像强度特征和空间约束先验知识，能够更好地对图像进行分割。但是常规Pairwise MRF模型中的先验参数$ \beta $使空间邻域内像素点趋于取相同的类标签，对于含有丰富特征的图像易产生细节特征的过平滑问题，且不能很好地描述边缘，这是由邻域内每个像素对中心像素的惩罚相同的假设条件不合理造成的^[16]。

2 局部区域一致性MRF的分割模型

为有效保持图像分割的边缘，并提高基于GMRF的图像分割方法在分割过程中对噪声、纹理等突变信号的鲁棒性，本文提出一种局部区域一致性MRF分割模型。在GMRF模型中，图像的标签过程缺乏像素间的空间信息，易受噪声或图像非同态性的干扰。为解决上述问题，本文提出模型在图像观察场X中利用K-L距离引入邻域内像素空间关系，建立局部空间约束高斯混合模型(Local Spatial Constrained GMM，LSCGMM)。在标签场Y中，针对常规MRF的点对交互结构不能有效表达复杂图像的局部先验知识问题，本文提出一种基于区域的部分二阶Potts先验模型，并对局部区域内标签取不同值时进行不一致性惩罚。在MRF框架下，将图像的局部空间约束高斯似然特征和局部区域一致性先验知识相结合来建立图像的局部区域一致性MRF分割模型。

2.1 局部空间约束高斯混合模型

常规GMM模型的独立性假设忽略了图像的空间先验关系，为解决此问题，建立图像局部区域内像素与标签带空间约束的似然关系并提出一种LSCGMM模型。令$ {\mathit{\boldsymbol{x}}}_{s} $为图像局部区域内中心像素点，定义其邻域$ {N}_{s} $内的像素点为$ {\mathit{\boldsymbol{x}}}_{r} $，即$ {\mathit{\boldsymbol{x}}}_{r}= $ $ \left\{{\mathit{\boldsymbol{x}}}_{r}\right.| $ $ r\left.\in {N}_{s}\right\} $。在局部区域内，倾向于取相同标签的相邻像素强度值相近，反之则相邻像素强度值相差较大。基于此假设，在距离度量方式中引入局部区域内像素的空间约束关系，获得中心像素$ {x}_{s} $与邻域像素$ {x}_{r} $的空间相似性权重，并将该权重引入到GMM中，得到LSCGMM似然模型。

$ P\left(X\right|Y, \Theta )=\prod _{s=1}^{N}\left[p\left({\mathit{\boldsymbol{x}}}_{s}|{y}_{s}, \theta \right)\prod _{r\in {N}_{s}}p{\left({\mathit{\boldsymbol{x}}}_{r}|{y}_{r}, \theta \right)}^{\frac{w\left({y}_{r}\right)}{{w}_{r}}}\right] $

(6)

其中，$ {w}_{r}=\sum\limits_{r\in {N}_{s}}w\left({y}_{r}\right) $为一个归一化函数，对于区域内所有像素，当$ w\left({y}_{r}\right)=1 $时，$ {w}_{r} $为区域内像素的个数。包含空间信息的权重函数$ w\left({y}_{r}\right) $可表示为：

$ w\left({y}_{r}\right)=\frac{1}{\sqrt{2\mathrm{\pi }}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\left|{\mathit{\boldsymbol{x}}}_{r}-{\mathit{\boldsymbol{x}}}_{s}\right|}^{2}}{2{\delta }^{2}}\right) $

(7)

其中，$ {\delta }^{2} $是一个高斯核函数的方差。显然，在邻域系统$ {N}_{s} $中，中心像素点$ {\mathit{\boldsymbol{x}}}_{s} $与其邻域像素点$ {\mathit{\boldsymbol{x}}}_{r} $之间距离越大，则相似度越低、赋予像素权重越小。通过引入邻域系统像素间相似性，在分割过程中引入像素的空间相关性信息，能够将像素类的分布根据邻域空间信息自适应调整，约束图像区域内具有相似像素强度值的标签一致性。从式(6)、式(7)可知，LSCGMM模型利用像素间的距离描述图像的局部空间关系，通过空间相似性权重约束邻域内像素的标签一致性，使得同质区域像素具有相同类标签，避免了噪声对分割结果的干扰。

2.2 部分二阶区域Potts先验模型

常规Pairwise MRF模型的简单点对结构限制了MRF模型的表达能力，使其无法表达复杂自然图像的区域和全局先验知识^[17]，而高阶MRF(Higher-order MRF，HMRF)虽然具有较强的表达能力，能够表达更加复杂的图像先验知识和统计信息，但是由于高阶势能项学习和推理较高的复杂度限制了HMRF模型在图像处理领域中的进一步应用。为解决此问题，本文提出一种部分二阶区域Potts先验模型。该模型首先对图像进行区域划分，建立大小为$ w\times w $的局部区域集合$ W $，假定第s个局部区域$ {w}_{s} $的标签场$ {Y}_{{w}_{s}} $是一个MRF分布，局部区域内标签场$ {Y}_{{w}_{s}} $的先验使用部分二阶区域Potts先验模型描述，根据MRF与Gibbs分布的等价性，则建立的部分二阶区域Potts先验模型可描述为：

$ P\left(Y|\beta \right)=\prod _{s}P\left({y}_{s}|{y}_{{w}_{s}}, \beta \right)=\prod _{s}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\beta {G}_{d}\sum\limits_{s, r\in {w}_{s}}\left[1-\delta \left({y}_{s}, {y}_{r}\right)\right]\right\}/Z\left(\beta \right) $

(8)

其中，$ Z\left(\beta \right) $是归一化常量，$ \beta $为Potts模型先验参数，控制局部区域分割的平滑度。$ {G}_{d} $为局部区域内中心节点在水平、垂直、对角和反对角等8个方向上的惩罚因子，对于水平和垂直方向的邻域节点，令$ {G}_{d}=1 $，对于对角和反对角方向的邻域节点，由于其距离中心节点较远，因此令$ {G}_{d}=1/\sqrt{2} $。$ \delta \left({y}_{s}, {y}_{r}\right) $为delta函数，可表示为：

$ \delta \left({y}_{s}, {y}_{r}\right)=\left\{\begin{array}{cc}1, {y}_{s}\ne {y}_{r}& \\ 0, {y}_{s}={y}_{r}& \end{array}\right. $

(9)

由式(8)可知，本文所提先验模型在局部区域$ {w}_{s} $内中心节点标签$ {y}_{s} $与邻域节点标签$ {y}_{r} $不相等时，这个不一致惩罚为$ \beta {G}_{d}\sum\limits_{s, r\in {w}_{s}}\left[1-\delta \left({y}_{s}, {y}_{r}\right)\right] $。由此可知，该不一致性惩罚值随局部区域内不一致标签数量的增加而呈线性增长，即区域内较多像素点的类标签趋于一致时，具有较小的惩罚值，且避免了常规的区域同一性惩罚造成的过平滑现象。仅当中心节点标签$ {y}_{s} $与邻域节点标签$ {y}_{r} $全部相等时，对局部区域内标签的不一致性惩罚为0。与传统基于四邻域结构的Pairwise MRF模型相比，部分二阶区域结构通过对八邻域结构引入更多的图像边缘线过程特征，有助于描述更为复杂的图像边缘特征。虽然本文所提先验模型引入的较大局部区域包含了图像更多的局部先验知识，但是并未显著增加本文所提先验模型的计算复杂度。

2.3 基于局部区域一致性约束MRF分割模型

在常规的MRF分割框架下，根据贝叶斯定理，基于LSCGMM似然特征和部分二阶区域Potts先验构造图像的MRF分割能量模型为：

$E\left( {Y|X,B,\Theta ,\beta } \right) = \underbrace {E\left( {X|Y,\Theta } \right)}_{{E_{{\rm{Data}}}}\left( {X,Y} \right)} + \underbrace {E\left( {Y|\beta } \right)}_{{E_{{\rm{Smooth}}}}\left( Y \right)}$

(10)

其中，令$ {E}_{\mathrm{D}\mathrm{a}\mathrm{t}\mathrm{a}}\left(X, Y\right)=-\mathrm{l}\mathrm{n}P\left(X\right|Y, \Theta ) $表示图像观察场与标签场的似然特征能量项，根据提出的LSCGMM似然模型式(6)，则$ {E}_{\mathrm{D}\mathrm{a}\mathrm{t}\mathrm{a}}\left(X, Y\right) $可定义为：

$\begin{array}{l} {E_{{\rm{Data}}}}\left( {X,Y} \right) = - {\rm{ln}}P(X|Y,\Theta ) = \\ \sum\limits_{s \in S} - \left[ {{\rm{ln}}\sqrt {2\pi \sigma _k^2} + \frac{1}{2}{{\left( {{\mathit{\boldsymbol{x}}_s} - {\mathit{\boldsymbol{\mu }}_k}} \right)}^{\rm{T}}}\sigma _k^{ - 1}\left( {{\mathit{\boldsymbol{x}}_s} - {\mathit{\boldsymbol{\mu }}_k}} \right)} \right.\\ \left. { + \sum\limits_{r \in {N_s}} {\frac{{w\left( {{y_r}} \right)}}{{2{w_r}}}} {{\left( {{\mathit{\boldsymbol{x}}_r} - {\mathit{\boldsymbol{\mu }}_k}} \right)}^{\rm{T}}}\sigma _k^{ - 1}\left( {{\mathit{\boldsymbol{x}}_r} - {\mathit{\boldsymbol{\mu }}_k}} \right)} \right] \end{array}$

(11)

令$ {E}_{\mathrm{S}\mathrm{m}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{h}}\left(Y, \beta \right)=-\mathrm{l}\mathrm{n}P\left(Y\right|\beta ) $表示图像标签场的空间先验能量约束项，根据所提部分二阶区域Potts先验模型式(8)，则$ {E}_{\mathrm{S}\mathrm{m}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{h}}\left(Y, \beta \right) $可定义为：

$ {E}_{\mathrm{S}\mathrm{m}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{h}}\left(Y\right)=-\mathrm{l}\mathrm{n}\prod _{s}P\left({y}_{s}\right|{y}_{{w}_{s}}, \beta )=\sum\limits_{s}\beta {G}_{d}\sum\limits_{r\in {w}_{s}}\left[1-\delta \left({y}_{s}, {y}_{r}\right)\right] $

(12)

根据式(12)可知，在局部区域$ {w}_{s} $内，$ {E}_{\mathrm{D}\mathrm{a}\mathrm{t}\mathrm{a}}\left(X, Y\right) $利用局部像素空间相关性信息，可有效约束更相似的邻域像素取相同的标签值，对于图像的噪声或者纹理突变信息的干扰具有鲁棒性。$ {E}_{\mathrm{S}\mathrm{m}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{h}}\left(Y, \beta \right) $对于局部区域内标签的不一致惩罚的大小随着邻接标签取值不一致个数的增大而呈线性增长趋势，这种类别不一致性惩罚避免了局部区域内异质像素相同的惩罚值造成的过平滑现象，且能够有效保证图像局部区域一致性约束。因此，本文提出的局部区域一致性MRF分割模型不仅能对噪声图像或纹理信号的干扰具有鲁棒性，而且可有效提升图像分割结果的边缘效果，从而得到更平滑的分割结果。

3 基于边缘约束局部区域MRF的分割模型

虽然本文建立的局部区域一致性约束MRF模型对噪声或纹理突变信息干扰具有较强的鲁棒性，但是由于图像分割对象的边缘呈高度尖峰态分布，而局部区域一致性模型常约束区域内像素倾向于取相同的标签，因此本文所提模型难以描述图像强烈突变的高维复杂性边缘特征，易导致边缘区域的过平滑，使得对含有丰富细节特征的图像产生较大的误分割。图像边缘是图像的基本特征之一，通常携带有图像的大部分信息，这种边缘信息能够有效约束图像边缘分割的准确性。为解决局部区域一致性约束MRF模型对图像边缘的过平滑问题，本文提出利用图像边缘特征建立分割对象的一致性边缘约束，并融合局部区域一致性模型表示的区域特征，建立一种ECLRMRF分割模型。

3.1 边缘能量约束先验模型

为建立ECLRMRF模型，基于Canny边缘检测算子提取图像的边缘特征，提出一种边缘能量约束先验模型，并融合局部区域一致性MRF图像分割模型，建立所提ECLRMRF模型的边缘能量项，且在分割过程中对图像目标分割的边界进行一致性边缘约束。本文所提边缘约束先验模型使用Canny边缘检测算子提取图像的边缘特征，得到图像的边缘信息$ B=\left\{{b}_{1}, {b}_{2}, \cdots , {b}_{N}\right\} $，令$ {b}_{s} $表示边缘图像中第$ s $个节点，则图像在该点的边缘特征$ {b}_{s} $可表示为：

$ {b}_{s}=\left\{\begin{array}{cc}0, s\notin \mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}& \\ 1, s\in \mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}& \end{array}\right. $

(13)

其中，$ {b}_{s}=1 $表明第$ s $个节点是边缘节点，$ {b}_{s}=0 $表明第$ s $个节点不是边缘节点。

本文定义边缘先验约束能量项为：

$ E\left(B|\gamma \right)=\gamma \sum\limits_{s}\sum\limits_{r\in {N}_{s}}{\phi }_{sr}\left({b}_{s}, {b}_{r}\right) $

(14)

其中，$ \gamma $表示边缘先验参数，其大小衡量边缘约束的重要性。$ {N}_{s} $表示第$ s $个节点的邻域节点集合，二阶边缘势函数$ {\phi }_{sr}\left({x}_{s}, {x}_{r}\right) $可表示为：

$ {\phi }_{sr}\left({b}_{s}, {b}_{r}\right)=\left|{b}_{s}-{b}_{r}\right| $

(15)

式(15)表明，对于邻接节点对$ \left(s, r\right) $，如果2个邻域节点均是边缘上的点或者不是边缘上的点，则边缘约束项的惩罚为0，以保证分割结果不过于平滑，有利于保持图像的边缘约束，反之则给该点增加一个边缘惩罚$ \gamma $。因此，边缘特征的引入增加了对图像目标边界的有效约束，在能量模型优化过程中，通过边缘能量项的最小化约束图像分割结果向边缘收敛，能够有效提高图像分割边缘的准确度。

3.2 边缘约束局部区域MRF模型

基于MRF框架，在能量函数中融合LSCGMM模型表达的似然特征、部分二阶区域Potts模型建立的区域先验知识和边缘能量先验约束模型描述的图像边缘特征，建立基于边缘约束的局部区域MRF(ECLRMRF)模型可表示为：

$\begin{array}{*{20}{l}} {E\left( {Y|X,B,\Theta ,\beta ,\gamma } \right) = \underbrace {E\left( {X|Y,\Theta } \right)}_{{E_{{\rm{Data}}}}\left( {X,Y} \right)} + \underbrace {E\left( {Y|\beta } \right)}_{{E_{{\rm{Smooth}}}}\left( Y \right)} + \underbrace {E\left( {B|\gamma } \right)}_{{E_{{\rm{Edge}}}}\left( B \right)} = }\\ {\underbrace {\sum\limits_{s \in S} {\left[ {{\rm{ln}}\sqrt {2\pi \sigma _k^2} + \frac{1}{2}{{\left( {{\mathit{\boldsymbol{x}}_s} - {\mathit{\boldsymbol{\mu }}_k}} \right)}^{\rm{T}}}\sigma _k^{ - 1}\left( {{\mathit{\boldsymbol{x}}_s} - {\mathit{\boldsymbol{\mu }}_k}} \right) + \sum\limits_{r \in {N_s}} {\frac{{w\left( {{y_r}} \right)}}{{2{w_r}}}} {{\left( {{\mathit{\boldsymbol{x}}_r} - {\mathit{\boldsymbol{\mu }}_k}} \right)}^{\rm{T}}}\sigma _k^{ - 1}\left( {{\mathit{\boldsymbol{x}}_r} - {\mathit{\boldsymbol{\mu }}_k}} \right)} \right]} }_{{E_{{\rm{Data}}}}\left( {X,Y} \right)} + }\\ {\underbrace {\sum\limits_s {\sum\limits_{r \in {N_s}} \beta } {G_d}\left[ {1 - \delta \left( {{y_s},{y_r}} \right)} \right]}_{{E_{{\rm{Smooth}}}}\left( Y \right)} + \underbrace {\gamma \sum\limits_s {\sum\limits_{r \in {N_s}} {{\phi _{sr}}} } \left( {{b_s},{b_r}} \right)}_{{E_{{\rm{Edge}}}}\left( B \right)}} \end{array}$

(16)

本文提出的模型由以下3个部分组成：局部空间约束高斯混合似然能量项$ {E}_{\mathrm{D}\mathrm{a}\mathrm{t}\mathrm{a}}\left(X, Y\right) $、部分二阶区域Potts先验项$ {E}_{\mathrm{S}\mathrm{m}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{h}}\left(Y\right) $和边缘能量约束项$ {E}_{\mathrm{E}\mathrm{d}\mathrm{g}\mathrm{e}}\left(B\right) $。似然能量项$ {E}_{\mathrm{D}\mathrm{a}\mathrm{t}\mathrm{a}}\left(X, Y\right) $基于相似性测量的方式引入了融合邻域空间信息的高斯混合似然特征，保持区域内像素与分配标签的似然一致性。先验项$ {E}_{\mathrm{S}\mathrm{m}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{h}}\left(Y\right) $根据局部区域内和中心标签值不一致数量呈正比的线性惩罚建立局部区域先验的一致性平滑约束，并且引入更多的先验空间约束，有效提高分割过程的鲁棒性。边缘能量约束项$ {E}_{\mathrm{E}\mathrm{d}\mathrm{g}\mathrm{e}}\left(B\right) $通过引入图像边缘特征，增加对目标分割边界的约束，提高图像分割边缘的准确度。

最终根据MAP准则，本文所提ECLRMRF模型的最优估计标签$ {Y}^{\mathrm{*}} $为：

$ {Y}^{\mathrm{*}}=\underset{Y}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}E\left(Y|X, B, \Theta , \beta , \gamma \right)=\underset{Y}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left[E\left(X|Y, \Theta \right)+E\left(Y|\beta \right)+E\left(B|\gamma \right)\right] $

(17)

3.3 基于ECLRMRF模型的Gibbs分割算法

对本文所提ECLRMRF模型求解最优估计标签$ {Y}^{\mathrm{*}} $时，首先需要估计ECLRMRF模型的参数$ \left(\Theta , \beta , \gamma \right) $。部分二阶区域Potts先验模型参数$ \beta $通常采用人工试错法进行估计，选择分割结果最好的值作为分割图像的最优参数值，选择的原则是对噪声强或者纹理较复杂的图像选择较大的$ \beta $值，反之则选择较小的$ \beta $值。边缘先验约束参数$ \gamma $由人工给定。局部空间约束GMM参数$ {\Theta }^{\mathrm{*}}=\left\{{\mu }_{1}, {\mu }_{2}, \cdots , {\mu }_{K};{\sigma }_{1}, {\sigma }_{2}, \cdots , {\sigma }_{K}\right\} $，$ {\mathit{\boldsymbol{\mu }}}_{k}, {\sigma }_{k} $分别是基于LSCGMM的似然模型的第k类分布的均值和方差。由于LSCGMM是一个含有未知量$ Y $的参数模型，对于该模型，期望最大化(Expectation Maximization，EM)算法是一种有效的估计算法。EM算法由期望步(E步)和最大化步(M步)组成。E步根据在前一次迭代估计的参数$ {\Theta }^{\left(t-1\right)} $计算隐含变量$ Y $的期望，M步是对似然函数用最大似然法求解新参数$ {\Theta }^{\left(t\right)} $的估计值。EM算法求解LSCGMM模型参数的迭代公式为：

$ {\mathit{\boldsymbol{\mu }}}_{k}^{t}=\frac{\sum\limits_{s\in S}\sum\limits_{r\in {N}_{s}}w\left({y}_{r}\right){\mathit{\boldsymbol{x}}}_{r}P\left[{\mathit{\boldsymbol{x}}}_{r}|{y}_{r}, {\mathit{\boldsymbol{\mu }}}_{k}^{t-1}, {\sigma }_{k}^{t-1}\right]}{\sum\limits_{s\in S}\sum\limits_{r\in {N}_{s}}w\left({y}_{r}\right)P\left[{\mathit{\boldsymbol{x}}}_{r}|{y}_{r}, {\mathit{\boldsymbol{\mu }}}_{k}^{t-1}, {\sigma }_{k}^{t-1}\right]} $

(18)

$ {\sigma }_{k}^{t}=\frac{\sum\limits_{s\in S}\sum\limits_{r\in {N}_{s}}w\left({y}_{r}\right)\left[{\mathit{\boldsymbol{x}}}_{r}-{\mathit{\boldsymbol{\mu }}}_{k}^{t-1}\right]{\left[{\mathit{\boldsymbol{x}}}_{r}-{\mathit{\boldsymbol{\mu }}}_{k}^{t-1}\right]}^{\mathrm{T}}P\left[{\mathit{\boldsymbol{x}}}_{r}|{y}_{r}, {\mathit{\boldsymbol{\mu }}}_{k}^{t-1}, {\sigma }_{k}^{t-1}\right]}{\sum\limits_{s\in S}\sum\limits_{r\in {N}_{s}}w\left({y}_{r}\right)P\left[{\mathit{\boldsymbol{x}}}_{r}|{y}_{r}, {\mathit{\boldsymbol{\mu }}}_{k}^{t-1}, {\sigma }_{k}^{t-1}\right]} $

(19)

本文提出基于ECLRMRF分割模型的Gibbs采样算法步骤为：

步骤 1  给定输入图像$ X $，对待分割图像进行参数初始化：人工给定分类数$ K $，部分二阶区域Potts先验模型参数$ \beta $和边缘先验参数$ \gamma $，Gibbs采样算法的温度参数初值$ {T}_{0}=4 $℃，局部区域大小给定为$ 3\times 3 $。

步骤 2  LSCGMM的参数$ \left({\mu }_{k}^{\left(0\right)}, {\sigma }_{k}^{\left(0\right)}\right) $初始化由K-means算法估计。

步骤 3  采用EM算法根据式(18)和式(19)估计具有LSCGMM参数$ {\Theta }^{\mathrm{*}}={\left\{{\mu }_{K}, {\sigma }_{K}\right\}}_{K=\mathrm{1, 2}, \cdots ，k} $，并根据式(8)计算权重$ w\left({y}_{n}\right) $。

步骤 4  由步骤1给定的参数$ \beta $、步骤3估计的LSCGMM参数$ {\Theta }^{\mathrm{*}}={\left\{{\mu }_{K}, {\sigma }_{K}\right\}}_{K=\mathrm{1, 2}, \cdots ，k} $以及权重$ w\left({y}_{n}\right) $，按照最大似然法根据式(20)对图像进行初始化分割。

$ \begin{array}{l}Y=\underset{Y}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}E\left(Y|X, {\Theta }^{\mathrm{*}}, \beta \right)=\\ \underset{Y}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left[E\left(X|Y, {\Theta }^{\mathrm{*}}\right)+E\left(Y|\beta \right)\right]\end{array} $

(20)

步骤 5  根据ECLRMRF模型式(16)，采用Gibbs采样算法进行迭代估计给像素$ {\mathit{x}}_{s} $分配一个标签$ {y}_{s} $时的接受率，Gibbs采样算法接受这个分配的标签的概率如式(21)所示：

$ P\left({y}_{s}=k\right)=\frac{\mathrm{e}\mathrm{x}\mathrm{p}\left[-E\left({y}_{s}=k|{x}_{{w}_{s}}, {b}_{s}, \Theta , \beta , \gamma \right)\right]}{\sum\limits_{\zeta \in \left\{\mathrm{1, 2}, \cdots , K\right\}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-E\left({y}_{s}=\zeta |{x}_{{w}_{s}}, {b}_{s}, \Theta , \beta , \gamma \right)\right]} $

(21)

其中，$ E\left({y}_{s}=k|{x}_{{w}_{s}}, {b}_{s}, \Theta , \beta , \gamma \right) $为像素$ {\mathit{x}}_{s} $分配一个标签为$ k $时边缘约束局部区域一致性能量。

步骤 6  算法迭代运行，如果$ \left|{E}^{\left(t\right)}\left(Y|X, B, \Theta , \beta , \gamma \right)-\right. $ $ \left.{E}^{\left(t-1\right)}\left(Y|X, B, \Theta , \beta , \gamma \right)\right|\le \epsilon $($ \epsilon $是一个很小的值)或达到最大迭代次数，则算法终止，并输出最优分割结果$ {Y}^{\mathrm{*}}=\left\{{y}_{1}, {y}_{2}, \cdots , {y}_{N}\right\} $，否则减小Gibbs采样算法温度$ T=0.95{T}^{\left(t-1\right)} $，并返回步骤5继续迭代。

4 实验结果与分析

为验证本文所提模型的可行性和有效性，以人工合成图像和自然图像作为实验对象，分别采用传统高斯马尔科夫随机场(CGMRF)模型^[20]、HMRF-EM算法^[8]、局部区域高斯马尔科夫随机场(LRGMRF)模型^[21]、RLSI-HMRF模型^[22]以及本文所提ECLRMRF模型进行分割实验。在实验中估计GMM参数时，由于EM算法对随机初值敏感，因而上述模型中的GMM初值都采用K-means算法进行初始化，不同图像Potts先验模型参数$ \beta $和边缘先验约束参数$ \gamma $由试错法给定。根据大量实验验证，参数$ \beta $取值区间(0，6]、参数$ \gamma $取值区间为[0，1.0]可得到较好的分割结果，实验设定$ \gamma =0.5 $。对于噪声弱、较平坦的图像，平滑参数$ \beta $和边缘先验约束参数$ \gamma $可取较小值，而对于噪声强、纹理丰富的图像，可增加$ \beta $和$ \gamma $的取值，分类数$ K $根据图像由人工指定，且邻域窗口大小设置为$ 3\times 3 $。

4.1 人工合成图像分割实验

人工合成图像如图 1(a)所示，图像大小为128像素×128像素、灰度级为256。图像分割为4类，各类均值分别为0、85、170和255。为评估本文所提模型对噪声图像的鲁棒性，人工叠加均值为零、方差为500的高斯噪声，实验中Potts先验模型参数$ \beta $设为1.2，对比模型GMM参数分别用EM算法进行估计。

5种不同模型对人工合成图像的分割图像对比如图 1所示。从图 1可以看出：CGMRF模型对噪声比较敏感且存在较多的误分点，使得分割效果最差；LRGMRF模型由于利用了较多的图像空间信息，因而对噪声具有一定的鲁棒性，但仍然存在一些误分割斑点且较大的区域平滑导致边缘的误分割；HMRF-EM模型受噪声影响出现了误分割现象，且分割结果出现了明显的边缘带；RLSI-HMRF模型由于建立了远距离的空间关系，分割结果更为平滑但缺乏边缘约束，导致分割结果出现边缘带；与上述4种模型相比，在本文所提ECLRMRF模型中，局部区域一致性约束对图像噪声干扰具有更强的鲁棒性，同时边缘特征约束图像分割结果向边缘收敛，有效提高了图像分割边缘的准确度。

为更加准确地表明本文所提模型的有效性，对人工合成噪声图像的5种模型分割结果进行量化评估，并采用分类正确率(Correct Classification Ratio，CCR)^[23]作为量化评估指标，且CCR可定义为：

$ C=\sum\limits_{k=1}^{K}\frac{\left|\mathrm{g}{\mathrm{t}}_{k}\bigcap \mathrm{s}\mathrm{e}{\mathrm{g}}_{k}\right|}{\left|\mathrm{G}\mathrm{T}\right|}\times 100\mathrm{\%} $

(22)

其中，$ \mathrm{g}{\mathrm{t}}_{k} $表示第$ k $个分割的真实结果，$ \mathrm{s}\mathrm{e}{\mathrm{g}}_{k} $表示第$ k $个分割的模型分割结果，$ \mathrm{G}\mathrm{T}={\cup }_{k=1}^{K}\mathrm{g}{\mathrm{t}}_{k} $。$ C\in $ $ \left[\mathrm{0, 1}\right] $，其数值越接近于1，表示分割结果越好。

表 1给出了5种模型在合成图像分割过程中的运行时间、迭代次数以及CCR。从表 1可以看出，在相同的噪声水平下，本文所提ECLRMRF模型的迭代次数低于CGMRF模型、LRGMRF模型和RLSI-HMRF模型，且具有较高的CCR值。这说明在5种对比模型中，本文所提ECLRMRF模型的分割结果准确性最优，并且与CGMRF、LRGMRF、HMRF-EM、RLSI-HMRF模型相比，其分割正确率分别提高了约5.23、1.32、41.92和0.45个百分点。由于每次迭代计算的区域信息量增多，因此和HMRF-EM模型、CGMRF模型、RLSI-HMRF模型相比，本文模型的迭代时间有所增加，但是迭代时间仍小于LRGMRF模型。因此，在分割准确度和算法运行时间方面而言，本文所提模型具有较高的综合性能。

4.2 自然图像分割实验

为进一步验证本文所提ECLRMRF模型的有效性，实验对5幅自然图像进行分割，而实验图像取自Berkeley数据集BSDS500^[24]。如图 2(a)所示，从左到右图像ID分别为80099、135069、253036、28083、55067，分割类别k分别设定为2、2、3、4、6，实验中Potts先验模型参数$ \beta $分别设定为2.1、1.4、2.3、3.2、5.3。从图 2可以看出，CGMRF模型的分割结果较差，这是由于Pairwise MRF模型简单的点对交互结构不能充分描述图像的空间先验特征，因此易受噪声或者纹理突变信号的干扰而出现误分割。例如，第1幅图中的"水面"有较多的分割斑点噪声，第3幅图中的"草地""树"和第4幅图中的"山""草地""树"这样纹理丰富的区域都出现了较为明显的误分割。与CGMRF模型相比，LRGMRF模型可利用更大的局部统计信息包含更多的空间先验特征，因而受噪声或者纹理突变信号的干扰较小，在噪声较弱或纹理平坦区域有较好的分割结果，但是由于该模型采用的局部区域模型的平滑作用，造成图像的细节特征缺失，且在图像边缘处易出现边缘带现象。第1幅图中"水牛"与"水面"之间出现的边缘细节缺失，且"水面"区域出现分割斑点。采用HMRF-EM模型的分割结果在图像纹理丰富区域和特征相近区域效果较差，易出现误分割现象，如第2幅图中"尾巴"，第3幅图中的"天空"和第4幅图中"草地"和"树"的误分割。RLSI-HMRF模型利用鲁棒的Pn Potts模型表达了更复杂的图像局部区域特征，因此图像分割结果背景平滑且轮廓清晰，但在纹理复杂的区域以及相似区域仍然会出现少量斑点以及误分割情况，比如第1幅图中"水面"区域，第2幅图中"尾巴"区域，第3幅图中"草地"和"树"以及第4幅图中"山"和"树"都出现误分割现象。本文所提ECLRMRF模型与其他4种模型相比，不仅对噪声或者纹理突变信息具有较强的鲁棒性，且在纹理细节丰富的区域取得了更为平滑的分割结果，如第1幅图中"水面"、第2幅图的"翅膀"和"尾巴"、第3幅图中的"草地"和"树"以及第4幅图中的"山"均取得了良好的分割结果，而且同时考虑了图像的边缘特征约束，可较好地保持图像的边缘信息，如第5幅图中"山"之间更为清晰的分割边界。

对于自然图像分割结果常用概率随机索引(Probabilistic Rand Index，PRI)指标衡量^[25]，且该指标可定义为：

$ \mathrm{P}\mathrm{R}\mathrm{I}\left({S}_{\mathrm{s}\mathrm{e}\mathrm{g}}, {S}_{\mathrm{g}\mathrm{t}}\right)=\frac{2}{N\left(N-1\right)}\sum\limits_{\underset{s＜r}{s, r}}\left[{c}_{sr}{p}_{sr}+\left(1-{c}_{sr}\right)\left(1-{p}_{sr}\right)\right] $

(23)

其中，$ {S}_{\mathrm{s}\mathrm{e}\mathrm{g}} $表示图像分割结果，$ {S}_{\mathrm{g}\mathrm{t}} $表示标准分割结果。$ {c}_{sr}=1 $表示像素$ s $和$ r $属于同一类别，$ {c}_{sr}=0 $表示像素$ s $和$ r $属于不同类别。$ \mathrm{P}\mathrm{R}\mathrm{I}\in \left[\mathrm{0, 1}\right] $，且其值越大，说明分割准确度越高。

表 2给出了5幅图像在不同模型下分割结果的PRI值。从表 2可以看出，本文所提ECLRMRF模型的PRI值均高于其他模型，这说明本文模型的分割性能优于其他对比模型。

为进一步评价本文所提分割算法的性能，表 3给出了实验图像在5种不同模型下分割过程迭代次数和运行时间的对比。从表 3可以看出，与其他4种模型相比，本文所提ECLRMRF模型具有相对较少的迭代次数和运行时间。

5 结束语

本文提出一种基于边缘约束局部区域MRF(ECLRMRF)模型的图像分割方法。通过对局部空间约束高斯混合模型(LSCGMM)引入邻接像素似然特征的空间约束相关性，提升图像的似然特征模型表达能力。部分二阶区域Potts先验模型利用较大的区域信息显著提升局部区域的先验知识表达能力，使得Canny边缘特征约束分割结果向图像分割边缘收敛。实验结果表明，本文所提模型不仅对图像噪声或者纹理突变信息的干扰具有较强的鲁棒性，而且能够有效保持图像分割结果的边缘信息，验证了该模型的可行性和有效性。Gibbs采样算法在采样过程中由于计算信息增多而使得分割算法耗时较长，因此下一步将在保持分割结果准确度的前提下引入群体智能优化算法，以进一步提高计算效率。