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

葉秀芬，張元科

(哈爾濱工程大學(xué)自動化學(xué)院，黑龍江哈爾濱150001)

對水下聲吶圖像進(jìn)行目標(biāo)分割是非常復(fù)雜和困難的，它不僅取決于被分割的不同目標(biāo)區(qū)域，還與海底混響噪聲、背景區(qū)域等有著緊密的聯(lián)系［1］。對聲吶圖像分割的目的就是要從復(fù)雜的海底混響區(qū)域中提取出目標(biāo)和陰影，并盡量保留圖像原始邊緣信息，它是圖像分析的關(guān)鍵步驟。聲吶在民用上可以搜尋失事飛機、船只的殘骸等目標(biāo);在軍用上可以用于探測各類軍事目標(biāo)［2］。因此，如何才能有效地對水下聲吶圖像進(jìn)行分割是國內(nèi)外研究者們研究的熱點與難點。

國內(nèi)外的研究者們已對馬爾可夫隨機場(Markov random field，MRF)分割方法在聲吶圖像上的應(yīng)用進(jìn)行了深入的研究，取得了重要的研究成果。如文獻(xiàn)［3］為此作了開創(chuàng)性的工作，基于Gibbs分布和MRF的一致性建立了關(guān)于重建圖像及其邊緣的聯(lián)合先驗分布模型;文獻(xiàn)［4］提出了基于分層的MRF聲吶圖像分割模型，文獻(xiàn)［5］盡管使用馬爾可夫隨機場模型對聲吶圖像分割取得了較好的分割結(jié)果，但是初始分割都需要根據(jù)圖像來人工選擇窗口的大小，并且算法運算較為復(fù)雜，很難達(dá)到水下目標(biāo)分割的自動和實時性的要求;文獻(xiàn)［6］中使用了快速的K均值聚類以及FCM算法對MRF初始化的參數(shù)進(jìn)行估計，但是也要先人為地確定聚類的數(shù)目，而且算法速度較慢，很難用于實時性的聲吶圖像分割;文獻(xiàn)［7］中提出基于灰度直方圖的譜聚類圖像分割方法，該方法是一種較為新穎的圖像分割方法，但是，該方法在水下聲吶圖像分割上的效果卻不是很理想;文獻(xiàn)［8］提出了一種非監(jiān)督的聲吶圖像分割算法，也是在假設(shè)類數(shù)已知的情況下進(jìn)行的。

綜上，目前的研究大都是人工確定MRF模型的參數(shù)或者是人為地確定聚類的類別數(shù)目，而沒有一種完全自動的聲吶圖像分割模型。針對這些問題，本文提出了一種新的基于MRF的非監(jiān)督聲吶圖像的自動分割方法，不僅能夠自動地確定MRF模型的初始化參數(shù)，而且還能自動地確定聲吶圖像要分割的類別以及類別數(shù)。最后，采用條件迭代算法(iterative conditional estimation，ICE)對聲吶圖像進(jìn)行了分割實驗，得到了較好的分割結(jié)果。

1 聲吶圖像直方圖分析

聲吶圖像一般是灰度圖像，而灰度圖像的直方圖反映了圖像中某種灰度出現(xiàn)的頻率，特別是聲吶圖像目標(biāo)、背景以及陰影的灰度級有一定的差別。理論上，根據(jù)聲吶圖像的直方圖就能夠把圖像的目標(biāo)、背景以及陰影分割開來，進(jìn)而也能夠確定聲吶圖像要分割的類別以及類別數(shù)。

1.1 圖像的直方圖處理與分析

灰度級為［0，L－1］范圍的數(shù)字圖像的直方圖是離散函數(shù)h(rk)=nk，這里rk是第k級灰度，nk是圖像中的灰度級為rk的像素個數(shù)。經(jīng)常以圖像中的像素點總數(shù)(用n表示)來除以它的每一個值得到歸一化的直方圖。因此，一個歸一化的直方圖由P(rk)=nk/n給出，這里k=0，1，…，L－1。簡單地說，P(rk)給出了灰度級為rk發(fā)生的概率估計值［9］。圖1為聲吶原始圖像。圖2為原始聲吶圖像的直方圖。

圖1 原始聲吶圖像Fig.1 Original sonar image

圖2 原始聲吶圖像的直方圖Fig.2 Histogram of the original sonar image

從圖2聲吶圖像的直方圖以及其他大量的實驗結(jié)果中可以看出，聲吶圖像的直方圖總體上具有高斯分布的特點［10］，但是直方圖的分布是離散化的，這種效果不利于使用本文的自動分割算法。因此，有必要對其進(jìn)行平滑處理?；诙喑叨鹊母咚菇鹱炙Ｐ屯ㄟ^對圖像進(jìn)行平滑，從而達(dá)到平滑離散分布的直方圖的目的，有利于后續(xù)算法自動地確定聲吶圖像的分割類別數(shù)(即確定圖像中含有陰影、目標(biāo)和背景這3類中的幾類)，并根據(jù)圖像中存在的類別數(shù)進(jìn)一步使用原始圖像進(jìn)行分割。

1.2 圖像的高斯金字塔模型

圖像處理的金字塔方法是將原始圖像分解成不同空間分辨率的子圖像，高分辨率(尺度大)的子圖像放在下層，低分辨率(尺度小)的圖像放在上層，從而形成一個金字塔的形狀［11］，如圖3所示。

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

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

1)先對圖像fl(i，j)進(jìn)行高斯卷積，下標(biāo)l表示金字塔的層數(shù)，再對圖像進(jìn)行降采樣:

2)先對1)處理后的圖像進(jìn)行上采樣，再對圖像進(jìn)行高斯卷積:

算法結(jié)束。其中高斯卷積核gσ為

對原始的聲吶圖像進(jìn)行高斯金字塔處理后的圖像如圖4所示。圖5為高斯金字塔處理后的聲吶圖像直方圖?？梢钥闯?，經(jīng)過平滑處理后的聲吶圖像與圖2相比直方圖也得到了平滑，可以很明顯地看出背景區(qū)域服從高斯分布;其中，背景左邊谷值的左側(cè)主要是陰影區(qū)域的像素點，背景右邊谷值的右側(cè)主要是目標(biāo)的像素點。如果不存在目標(biāo)和陰影的話，則較高和較低像素級的直方圖的像素點會很少;假設(shè)圖像中存在目標(biāo)或陰影，則其像素點在直方圖中所占的比率會較多;根據(jù)這一特點，本文后面將提出一種聲吶圖像自動分類的算法，用于確定圖像中是否存在目標(biāo)以及陰影，進(jìn)而，也就能自動地確定聲吶圖像分割的類別數(shù)。

圖4 原始聲吶圖像的高斯金字塔處理Fig.4 Gaussian pyramid of original sonar image

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

2 理論與算法

2.1 圖像分割中的Markov隨機場模型

最大后驗概率(MAP)是圖像處理中最常用的最優(yōu)化準(zhǔn)則，也是MRF建模中最常用的最優(yōu)化準(zhǔn)則。MRF模型與MAP準(zhǔn)則結(jié)合在一起就稱作MAP－MRF體系［12］。

假定觀測到的圖像數(shù)據(jù)為F，圖像上所有像素點的集合記為S。圖像的分割問題即要求解的問題滿足最大后驗概率準(zhǔn)則，對每個像素的分類標(biāo)號(標(biāo)號場)，記為ω。這樣由Bayes后驗概率準(zhǔn)則:

式中:P(F)為觀測數(shù)據(jù)的先驗分布，當(dāng)數(shù)據(jù)給定后為常數(shù)，所以不參與計算過程，可以不予以考慮;P(ω)是標(biāo)號場的先驗聯(lián)合Gibbs分布，即滿足馬爾可夫性。

假設(shè)C表示S所有的集簇，c表示C中的元素，U2(ω)為能量函數(shù)，Vc(ω)是與集簇相關(guān)的勢函數(shù)，那么

P(F|ω)是似然概率，在很多情況下，假定為各個位置的像素是獨立同分布的，即滿足

當(dāng)假定每個P(Fs|ωs)是高斯分布時，每個類的類參數(shù)都是由2個參數(shù)唯一確定該分布，即為λ和σ。將 Bayes后驗概率準(zhǔn)則:ω)取對數(shù)，得到的目標(biāo)函數(shù)為lnP(ω)+lnP(F|ω)，即要求的是使該表達(dá)式最大的時候ω的估計?，將似然函數(shù)和先驗Gibbs分布的表達(dá)式大代入，這時可以形成MRF-MAP下的目標(biāo)函數(shù)的最優(yōu)解問題:

根據(jù)聲吶圖像的性質(zhì)，選取一階的鄰域系統(tǒng)，且勢函數(shù)為ISING模型的勢函數(shù)，一階鄰域系統(tǒng)如圖6。

圖6 平面上的一階鄰域系統(tǒng)Fig.6 First-order neighborhood system of plane

勢函數(shù):

式中:β為耦合系數(shù)。

對每個像素點選取P(Fs|ωs)服從高斯分布，那么取對數(shù)后其函數(shù)形式為

取lnP(ω)+lnP(F|ω)為目標(biāo)函數(shù)，得到的分割結(jié)果為

其中，U2(ω)為能量函數(shù)，且

2.2 MRF模型的ICE迭代算法

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

1)通過訓(xùn)練樣本得到所需似然函數(shù)P(F|ω)參數(shù)集合，即對應(yīng)不同分類λ情況有

初始化勢函數(shù)中的耦合系數(shù)β，經(jīng)驗值是取區(qū)間［0.5，1］之間的數(shù)值。

2)依據(jù)似然概率即P(F|ω)最大化的準(zhǔn)則選取初始的標(biāo)記場ω0，即對每一個像素點s取，遍歷整個圖像得到整個圖像的初始分割?0;

3)根據(jù)目標(biāo)函數(shù)計算當(dāng)前分割結(jié)果:取k為當(dāng)前的迭代次數(shù)，每一個像素點s取:

遍歷整個圖像得到標(biāo)記場?k;

4)判斷收斂條件:以每次迭代過程中全局能量的變化量為收斂條件，計算當(dāng)前全局能量的值:

如果Δ≥Ek－Ek－1認(rèn)為全局能量變化很小，標(biāo)記場?k為最后的分割結(jié)果，式中Δ為常數(shù)，即預(yù)先設(shè)定的能量改變量的閾值，得到了ICM算法的分割結(jié)果;否則取k=k+1，轉(zhuǎn)到步驟2)。

2.3 聲吶圖像自動分類的模型與算法

本文提出一種能夠自動確定聲吶圖像分類及分類個數(shù)的模型如下:

式中:p(rk)為圖像經(jīng)過高斯金字塔處理后模型的歸一化統(tǒng)計直方圖。y1為圖像是否含有陰影區(qū)的判別函數(shù)，y1∈{0，1};y2為圖像是否含有背景區(qū)的判別函數(shù)，y2∈{0，1};y3為圖像是否含有目標(biāo)區(qū)的判別函數(shù)，y3∈{0，1};n為圖像分類的個數(shù)，此處聲吶圖像n∈ {1，2，3}。

如圖5高斯金字塔處理后的圖像直方圖所示，定義像素峰值rpv為圖像統(tǒng)計直方圖取得極大值時的灰度級;左邊像素谷值rvvl為圖像統(tǒng)計直方圖從像素峰值rpv開始向左統(tǒng)計直方圖和變化較小時的灰度級;同理，右邊像素谷值rvvr為圖像統(tǒng)計直方圖從像素峰值rpv開始向右統(tǒng)計直方圖和變化較小時的灰度級。并定義判別某類別是否存在的函數(shù):

式中:θ為判斷閾值。

自動分類模型算法的實現(xiàn):

1)設(shè)定迭代條件iter=0.01，計數(shù)個數(shù)count=0，迭代步數(shù)l=0;

2)計算像素峰值:

并計算灰度級為rpv時的概率:pl=p(rpv);

3)計算灰度區(qū)間［rpv－l，rpv+l］的概率統(tǒng)計和:pl+1=p(［rpv－l，rpv+l］);若|pl+1－pl|＜iter，則count++;count大于給定一整數(shù)m，則轉(zhuǎn)4)，否則，l=l+1，并轉(zhuǎn)3);

4)計算左右邊像素谷值:

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

5)計算圖像的判別函數(shù)yi及類別個數(shù)n:

在式(1)中，假定背景區(qū)域像素點占整幅圖像的比率超過80%。

2.4 局部能量極化的原理與算法

令(x，y)為某一圖像中像素的坐標(biāo)，令Sxy表示某一確定大小的鄰域(子圖像)，其中心在(x，y)。則在Sxy中像素的平均值、能量和方差能以下面的式子計算［9］:

局部能量極化即把圖像分成很多大小相同的區(qū)域，再分別計算源圖像中各個區(qū)域的均值、能量和方差。然后根據(jù)能量函數(shù)的大小，以及圖像分類的個數(shù)和判別函數(shù)來得到MRF分割模型的初始化參數(shù)。

其算法如下:

1)根據(jù)前面的算法得到圖像的高斯金字塔模型gl(i，j)、左右邊像素谷值rvvl和rvvr以及判別函數(shù)yi;

2)分別以rvvl和rvvr為閾值得到金字塔模型的閾值化函數(shù)，即陰影模型gshadow、目標(biāo)模型gtarget和背景模型gback;

3)根據(jù)上面得到的各個模型和判別函數(shù)yi可以估計MRF模型初始化參數(shù)，其中所選取的區(qū)域的長L和寬W;

4)然后根據(jù)區(qū)域能量公式計算各個區(qū)域的能量函數(shù)，并將能量最小時的區(qū)域均值μ1和均方差σ1作為MRF模型的陰影區(qū)域參數(shù);將能量取得極大值時的區(qū)域均值μ2和均方差σ2作為MRF模型的目標(biāo)區(qū)域參數(shù);而選取一塊既非目標(biāo)也非陰影的區(qū)域的均值μ3和均方差σ3作為MRF模型的背景區(qū)域參數(shù);即可得到MRF模型的初始化參數(shù)均值μλ，均方差σλ，其中λ∈{1，2，3}為不同的分割區(qū)域。

最后，用得到的初始化參數(shù)初始化MRF模型，并利用ICE對聲吶圖像進(jìn)行分割，將分割后的模型分別與陰影模型gshadow和目標(biāo)模型gtarget進(jìn)行簡單的與運算與及高斯低通濾波處理，即可以得到最終的聲吶圖像分割結(jié)果。

3 實驗結(jié)果與分析

通過對圖1原始聲吶圖像和圖2原始聲吶圖像的直方圖以及大量實驗數(shù)據(jù)的分析可以看出，聲吶圖像含有比較復(fù)雜的海底混響噪聲，因此，對聲吶圖像的分割是比較困難的。本文在對聲吶圖像進(jìn)行大量實驗的基礎(chǔ)上，得出聲吶圖像經(jīng)過適當(dāng)?shù)念A(yù)處理之后也是有規(guī)律可循的結(jié)論，在此基礎(chǔ)上提出了能夠自動確定分類個數(shù)以及自動確定MRF模型初始化參數(shù)的分割方法，不僅分割效果較好，而且時間復(fù)雜度也較小，有利于實時聲吶圖像數(shù)據(jù)處理。

根據(jù)本文提出的算法，在實驗運行環(huán)境為Celeron(R)2.93GHz CPU、1.00GB 內(nèi)存、Windows XP 操作系統(tǒng)下的VS2008環(huán)境中進(jìn)行了分割實驗。實現(xiàn)了對圖2既有目標(biāo)又有陰影的聲吶圖像以及圖7較為復(fù)雜的聲吶圖像的分割實驗，驗證了本文算法的有效性。其中圖2像素尺寸大小為160×160，圖7像素尺寸大小為100×100。

該實驗結(jié)果與其他文獻(xiàn)算法的分割結(jié)果相比較，文獻(xiàn)［6］利用FCM快速聚類算法對MRF模型的參數(shù)進(jìn)行初始化，但仍舊需要首先指定聚類的個數(shù)，文獻(xiàn)［7］同樣也需要指定分類的類別數(shù)，而且兩者的算法的時間復(fù)雜度較大，不利于實時性的聲吶圖像數(shù)據(jù)處理。圖8為不同算法對圖2原始聲吶圖像的分割效果。

圖7 較為復(fù)雜的聲吶圖像Fig.7 More sophisticated sonar image

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

圖9為不同算法對圖7較為復(fù)雜的聲吶圖像進(jìn)行分割的效果。從圖8和圖9可以看出，本文方法和人工MRF方法對2種不同的聲吶圖像的分割效果較好，而基于譜聚類的文獻(xiàn)［7］的算法對聲吶圖像的分割效果不是很理想。圖10分別為對聲吶圖像圖2和圖7進(jìn)行分割所用時間的性能分析。由表1可以看出，本文所提出的算法用時只比人工確定MRF模型參數(shù)的分割算法多了100 ms左右，所用的時間明顯較少，而且算法復(fù)雜度較低。

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

圖10 不同算法的時間復(fù)雜度分析Fig.10 Time complexity analysis of different algorithms

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

4 結(jié)束語

針對目前聲吶圖像分割方法中存在的缺點與不足，特別是已有的算法模型大多是事先確定需要分割的聲納圖像中所含有的分割類別個數(shù)的問題，本文提出了一種能夠自動確定聲吶圖像分類個數(shù)的模型，并在此基礎(chǔ)上通過一種局部能量極值化的方法來得到MRF模型的初始化參數(shù)，該算法簡單，實現(xiàn)容易，有利于實時性的聲納圖像分割與目標(biāo)識別。

［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場的側(cè)掃聲吶分割方法［J］.海洋學(xué)報，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ìn)的譜聚類圖像分割方法［J］.計算機工程與應(yīng)用，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的非監(jiān)督圖象分割［J］.電子學(xué)報，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.