基于馬爾可夫隨機場的非監督聲吶圖像分割方法

葉秀芬，張元科

(哈爾濱工程大學自動化學院，黑龍江哈爾濱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.