基于马尔可夫随机场的非监督声呐图像分割方法

叶秀芬，张元科

(哈尔滨工程大学自动化学院，黑龙江哈尔滨150001)

对水下声呐图像进行目标分割是非常复杂和困难的，它不仅取决于被分割的不同目标区域，还与海底混响噪声、背景区域等有着紧密的联系［1］。对声呐图像分割的目的就是要从复杂的海底混响区域中提取出目标和阴影，并尽量保留图像原始边缘信息，它是图像分析的关键步骤。声呐在民用上可以搜寻失事飞机、船只的残骸等目标;在军用上可以用于探测各类军事目标［2］。因此，如何才能有效地对水下声呐图像进行分割是国内外研究者们研究的热点与难点。

国内外的研究者们已对马尔可夫随机场(Markov random field，MRF)分割方法在声呐图像上的应用进行了深入的研究，取得了重要的研究成果。如文献［3］为此作了开创性的工作，基于Gibbs分布和MRF的一致性建立了关于重建图像及其边缘的联合先验分布模型;文献［4］提出了基于分层的MRF声呐图像分割模型，文献［5］尽管使用马尔可夫随机场模型对声呐图像分割取得了较好的分割结果，但是初始分割都需要根据图像来人工选择窗口的大小，并且算法运算较为复杂，很难达到水下目标分割的自动和实时性的要求;文献［6］中使用了快速的K均值聚类以及FCM算法对MRF初始化的参数进行估计，但是也要先人为地确定聚类的数目，而且算法速度较慢，很难用于实时性的声呐图像分割;文献［7］中提出基于灰度直方图的谱聚类图像分割方法，该方法是一种较为新颖的图像分割方法，但是，该方法在水下声呐图像分割上的效果却不是很理想;文献［8］提出了一种非监督的声呐图像分割算法，也是在假设类数已知的情况下进行的。

综上，目前的研究大都是人工确定MRF模型的参数或者是人为地确定聚类的类别数目，而没有一种完全自动的声呐图像分割模型。针对这些问题，本文提出了一种新的基于MRF的非监督声呐图像的自动分割方法，不仅能够自动地确定MRF模型的初始化参数，而且还能自动地确定声呐图像要分割的类别以及类别数。最后，采用条件迭代算法(iterative conditional estimation，ICE)对声呐图像进行了分割实验，得到了较好的分割结果。

1 声呐图像直方图分析

声呐图像一般是灰度图像，而灰度图像的直方图反映了图像中某种灰度出现的频率，特别是声呐图像目标、背景以及阴影的灰度级有一定的差别。理论上，根据声呐图像的直方图就能够把图像的目标、背景以及阴影分割开来，进而也能够确定声呐图像要分割的类别以及类别数。

1.1 图像的直方图处理与分析

灰度级为［0，L－1］范围的数字图像的直方图是离散函数h(rk)=nk，这里rk是第k级灰度，nk是图像中的灰度级为rk的像素个数。经常以图像中的像素点总数(用n表示)来除以它的每一个值得到归一化的直方图。因此，一个归一化的直方图由P(rk)=nk/n给出，这里k=0，1，…，L－1。简单地说，P(rk)给出了灰度级为rk发生的概率估计值［9］。图1为声呐原始图像。图2为原始声呐图像的直方图。

图1 原始声呐图像Fig.1 Original sonar image

图2 原始声呐图像的直方图Fig.2 Histogram of the original sonar image

从图2声呐图像的直方图以及其他大量的实验结果中可以看出，声呐图像的直方图总体上具有高斯分布的特点［10］，但是直方图的分布是离散化的，这种效果不利于使用本文的自动分割算法。因此，有必要对其进行平滑处理。基于多尺度的高斯金字塔模型通过对图像进行平滑，从而达到平滑离散分布的直方图的目的，有利于后续算法自动地确定声呐图像的分割类别数(即确定图像中含有阴影、目标和背景这3类中的几类)，并根据图像中存在的类别数进一步使用原始图像进行分割。

1.2 图像的高斯金字塔模型

图像处理的金字塔方法是将原始图像分解成不同空间分辨率的子图像，高分辨率(尺度大)的子图像放在下层，低分辨率(尺度小)的图像放在上层，从而形成一个金字塔的形状［11］，如图3所示。

图3 图像处理的金字塔模型Fig.3 Pyramid model of image processing

图像的高斯金字塔模型算法:

1)先对图像fl(i，j)进行高斯卷积，下标l表示金字塔的层数，再对图像进行降采样:

2)先对1)处理后的图像进行上采样，再对图像进行高斯卷积:

算法结束。其中高斯卷积核gσ为

对原始的声呐图像进行高斯金字塔处理后的图像如图4所示。图5为高斯金字塔处理后的声呐图像直方图。可以看出，经过平滑处理后的声呐图像与图2相比直方图也得到了平滑，可以很明显地看出背景区域服从高斯分布;其中，背景左边谷值的左侧主要是阴影区域的像素点，背景右边谷值的右侧主要是目标的像素点。如果不存在目标和阴影的话，则较高和较低像素级的直方图的像素点会很少;假设图像中存在目标或阴影，则其像素点在直方图中所占的比率会较多;根据这一特点，本文后面将提出一种声呐图像自动分类的算法，用于确定图像中是否存在目标以及阴影，进而，也就能自动地确定声呐图像分割的类别数。

图4 原始声呐图像的高斯金字塔处理Fig.4 Gaussian pyramid of original sonar image

图5 高斯金字塔处理后的图像直方图Fig.5 Image histogram after Gaussian pyramid

2 理论与算法

2.1 图像分割中的Markov随机场模型

最大后验概率(MAP)是图像处理中最常用的最优化准则，也是MRF建模中最常用的最优化准则。MRF模型与MAP准则结合在一起就称作MAP－MRF体系［12］。

假定观测到的图像数据为F，图像上所有像素点的集合记为S。图像的分割问题即要求解的问题满足最大后验概率准则，对每个像素的分类标号(标号场)，记为ω。这样由Bayes后验概率准则:

式中:P(F)为观测数据的先验分布，当数据给定后为常数，所以不参与计算过程，可以不予以考虑;P(ω)是标号场的先验联合Gibbs分布，即满足马尔可夫性。

假设C表示S所有的集簇，c表示C中的元素，U2(ω)为能量函数，Vc(ω)是与集簇相关的势函数，那么

P(F|ω)是似然概率，在很多情况下，假定为各个位置的像素是独立同分布的，即满足

当假定每个P(Fs|ωs)是高斯分布时，每个类的类参数都是由2个参数唯一确定该分布，即为λ和σ。将 Bayes后验概率准则:ω)取对数，得到的目标函数为lnP(ω)+lnP(F|ω)，即要求的是使该表达式最大的时候ω的估计ϖ，将似然函数和先验Gibbs分布的表达式大代入，这时可以形成MRF-MAP下的目标函数的最优解问题:

根据声呐图像的性质，选取一阶的邻域系统，且势函数为ISING模型的势函数，一阶邻域系统如图6。

图6 平面上的一阶邻域系统Fig.6 First-order neighborhood system of plane

势函数:

式中:β为耦合系数。

对每个像素点选取P(Fs|ωs)服从高斯分布，那么取对数后其函数形式为

取lnP(ω)+lnP(F|ω)为目标函数，得到的分割结果为

其中，U2(ω)为能量函数，且

2.2 MRF模型的ICE迭代算法

条件迭代算法是典型的确定松弛算法，算法流程如下:

1)通过训练样本得到所需似然函数P(F|ω)参数集合，即对应不同分类λ情况有

初始化势函数中的耦合系数β，经验值是取区间［0.5，1］之间的数值。

2)依据似然概率即P(F|ω)最大化的准则选取初始的标记场ω0，即对每一个像素点s取，遍历整个图像得到整个图像的初始分割ϖ0;

3)根据目标函数计算当前分割结果:取k为当前的迭代次数，每一个像素点s取:

遍历整个图像得到标记场ϖk;

4)判断收敛条件:以每次迭代过程中全局能量的变化量为收敛条件，计算当前全局能量的值:

如果Δ≥Ek－Ek－1认为全局能量变化很小，标记场ϖk为最后的分割结果，式中Δ为常数，即预先设定的能量改变量的阈值，得到了ICM算法的分割结果;否则取k=k+1，转到步骤2)。

2.3 声呐图像自动分类的模型与算法

本文提出一种能够自动确定声呐图像分类及分类个数的模型如下:

式中:p(rk)为图像经过高斯金字塔处理后模型的归一化统计直方图。y1为图像是否含有阴影区的判别函数，y1∈{0，1};y2为图像是否含有背景区的判别函数，y2∈{0，1};y3为图像是否含有目标区的判别函数，y3∈{0，1};n为图像分类的个数，此处声呐图像n∈ {1，2，3}。

如图5高斯金字塔处理后的图像直方图所示，定义像素峰值rpv为图像统计直方图取得极大值时的灰度级;左边像素谷值rvvl为图像统计直方图从像素峰值rpv开始向左统计直方图和变化较小时的灰度级;同理，右边像素谷值rvvr为图像统计直方图从像素峰值rpv开始向右统计直方图和变化较小时的灰度级。并定义判别某类别是否存在的函数:

式中:θ为判断阈值。

自动分类模型算法的实现:

1)设定迭代条件iter=0.01，计数个数count=0，迭代步数l=0;

2)计算像素峰值:

并计算灰度级为rpv时的概率:pl=p(rpv);

3)计算灰度区间［rpv－l，rpv+l］的概率统计和:pl+1=p(［rpv－l，rpv+l］);若|pl+1－pl|＜iter，则count++;count大于给定一整数m，则转4)，否则，l=l+1，并转3);

4)计算左右边像素谷值:

若p(［rvvl，rvvr］)＞1－iter，则令y1=0，y2=1，y3=0，n=1，算法终止;

5)计算图像的判别函数yi及类别个数n:

在式(1)中，假定背景区域像素点占整幅图像的比率超过80%。

2.4 局部能量极化的原理与算法

令(x，y)为某一图像中像素的坐标，令Sxy表示某一确定大小的邻域(子图像)，其中心在(x，y)。则在Sxy中像素的平均值、能量和方差能以下面的式子计算［9］:

局部能量极化即把图像分成很多大小相同的区域，再分别计算源图像中各个区域的均值、能量和方差。然后根据能量函数的大小，以及图像分类的个数和判别函数来得到MRF分割模型的初始化参数。

其算法如下:

1)根据前面的算法得到图像的高斯金字塔模型gl(i，j)、左右边像素谷值rvvl和rvvr以及判别函数yi;

2)分别以rvvl和rvvr为阈值得到金字塔模型的阈值化函数，即阴影模型gshadow、目标模型gtarget和背景模型gback;

3)根据上面得到的各个模型和判别函数yi可以估计MRF模型初始化参数，其中所选取的区域的长L和宽W;

4)然后根据区域能量公式计算各个区域的能量函数，并将能量最小时的区域均值μ1和均方差σ1作为MRF模型的阴影区域参数;将能量取得极大值时的区域均值μ2和均方差σ2作为MRF模型的目标区域参数;而选取一块既非目标也非阴影的区域的均值μ3和均方差σ3作为MRF模型的背景区域参数;即可得到MRF模型的初始化参数均值μλ，均方差σλ，其中λ∈{1，2，3}为不同的分割区域。

最后，用得到的初始化参数初始化MRF模型，并利用ICE对声呐图像进行分割，将分割后的模型分别与阴影模型gshadow和目标模型gtarget进行简单的与运算与及高斯低通滤波处理，即可以得到最终的声呐图像分割结果。

3 实验结果与分析

通过对图1原始声呐图像和图2原始声呐图像的直方图以及大量实验数据的分析可以看出，声呐图像含有比较复杂的海底混响噪声，因此，对声呐图像的分割是比较困难的。本文在对声呐图像进行大量实验的基础上，得出声呐图像经过适当的预处理之后也是有规律可循的结论，在此基础上提出了能够自动确定分类个数以及自动确定MRF模型初始化参数的分割方法，不仅分割效果较好，而且时间复杂度也较小，有利于实时声呐图像数据处理。

根据本文提出的算法，在实验运行环境为Celeron(R)2.93GHz CPU、1.00GB 内存、Windows XP 操作系统下的VS2008环境中进行了分割实验。实现了对图2既有目标又有阴影的声呐图像以及图7较为复杂的声呐图像的分割实验，验证了本文算法的有效性。其中图2像素尺寸大小为160×160，图7像素尺寸大小为100×100。

该实验结果与其他文献算法的分割结果相比较，文献［6］利用FCM快速聚类算法对MRF模型的参数进行初始化，但仍旧需要首先指定聚类的个数，文献［7］同样也需要指定分类的类别数，而且两者的算法的时间复杂度较大，不利于实时性的声呐图像数据处理。图8为不同算法对图2原始声呐图像的分割效果。

图7 较为复杂的声呐图像Fig.7 More sophisticated sonar image

图8 不同算法对图2的分割效果Fig.8 Segmentation results of different algorithms in Fig.2

图9为不同算法对图7较为复杂的声呐图像进行分割的效果。从图8和图9可以看出，本文方法和人工MRF方法对2种不同的声呐图像的分割效果较好，而基于谱聚类的文献［7］的算法对声呐图像的分割效果不是很理想。图10分别为对声呐图像图2和图7进行分割所用时间的性能分析。由表1可以看出，本文所提出的算法用时只比人工确定MRF模型参数的分割算法多了100 ms左右，所用的时间明显较少，而且算法复杂度较低。

图9 不同算法对图7的分割效果Fig.9 Segmentation results of different algorithms in Fig.7

图10 不同算法的时间复杂度分析Fig.10 Time complexity analysis of different algorithms

表1 各种模型方法性能比较表Table 1 Performance comparison of various modeling methods

4 结束语

针对目前声呐图像分割方法中存在的缺点与不足，特别是已有的算法模型大多是事先确定需要分割的声纳图像中所含有的分割类别个数的问题，本文提出了一种能够自动确定声呐图像分类个数的模型，并在此基础上通过一种局部能量极值化的方法来得到MRF模型的初始化参数，该算法简单，实现容易，有利于实时性的声纳图像分割与目标识别。

［1］LANGNER F，KNAUER C，JANS W，et al.Side scan sonar image resolution and automatic object detection，classification and identification［C］//Oceans 2009-Europe.Berlin，Germany，2009:1-8.

［2］LIU G Y，BIAN H Y，SHI H.Sonar image segmentation based on spectral matting using morphological operations［J］.Journal of Jilin University:Engineering and Technology Edition，2012，42(1):228-233.

［3］GEMAN S，GEMAN D.Stochastic relaxation，Gibbs distributions，and the Bayesian restoration of images［J］.IEEE Trans Patten Anal Machine Intel，1984，6:721-741.

［4］MIGNOTTE M，COLLET C，PEREZ P，et al.Sonar image segmentation using an unsupervised hierarchical MRF model［J］.IEEE Transactions on Image Processing，2000，9(7):1216-1231.

［5］阳凡林，独知行，李家彪，等.基于 MRF场的侧扫声呐分割方法［J］.海洋学报，2006，28(4):43-48.YANG Fanlin，DU Ruxing，LI Jiabiao，et al.The segmentation method based on MRF field side scan sonar［J］.Acta Oceanologica Sinica，2006，28(4):43-48.

［6］WANG Xingmei，YE Xiufen，ZHANG Zhehui，et al.A novel automatic segmentation algorithm for sonar imagery［C］//Proceedings of the 2008 IEEE International Conference on Mechatronics and Automation(ICMA 2008).Takamatsu，Japan，2008:336-341.

［7］尹芳，陈德运，吴锐.改进的谱聚类图像分割方法［J］.计算机工程与应用，2011，47(21):185-187.YIN Fang，CHEN Deyun，WU Rui.Improved method of image segmentation using spectral clustering［J］.Computer Engineering and Applications，2011，47(21):185-187.

［8］汪西莉，刘芳，焦李成.基于不完全分层MRF的非监督图象分割［J］.电子学报，2004，32(7):1087-1089.WANG Xili，LIU Fang，JIAO Licheng.Unsupervised image segmentation based on incomplete hierarchical MRF［J］.Acta Electronica Sinica，2004，32(7):1087-1089.

［9］MANDHOUJ I，AMIRI H，MAUSSANG F，et al.Sonar image processing for underwater object detection based on high resolution system［C］//SIDOP 2012:2nd Workshop on Signal and Document Processing.Hammamet，Tunisia，2012:5-10.

［10］YAO K C，MIGNOTTE M，COLLET C，et al.Unsupervised segmentation using a self-organizing map and a noise model estimation in sonar imagery［J］.Pattern Recognition，2000，33(9):1575-1584.

［11］MIELKE M，SCHAFER A，BRUCK R.ASIC implementation of a Gaussian pyramid for use in autonomous mobile robotics［C］//2011 IEEE 54th International Midwest Symposium on Circuits and Systems(MWSCAS). Seoul，Korea，2011:1-4.

［12］JIE F，SHI Y，LI Y，et al.Interactive region-based MRF image segmentation［C］//2011 4th International Congress on Image and Signal Processing(CISP).Shanghai，China，2011:1263-1267.