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

葉秀芬，張?jiān)?/p>

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

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

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

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

1 聲吶圖像直方圖分析

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

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

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

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

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

從圖2聲吶圖像的直方圖以及其他大量的實(shí)驗(yàn)結(jié)果中可以看出，聲吶圖像的直方圖總體上具有高斯分布的特點(diǎn)［10］，但是直方圖的分布是離散化的，這種效果不利于使用本文的自動(dòng)分割算法。因此，有必要對(duì)其進(jìn)行平滑處理。基于多尺度的高斯金字塔模型通過(guò)對(duì)圖像進(jìn)行平滑，從而達(dá)到平滑離散分布的直方圖的目的，有利于后續(xù)算法自動(dòng)地確定聲吶圖像的分割類別數(shù)(即確定圖像中含有陰影、目標(biāo)和背景這3類中的幾類)，并根據(jù)圖像中存在的類別數(shù)進(jìn)一步使用原始圖像進(jìn)行分割。

1.2 圖像的高斯金字塔模型

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

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

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

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

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

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

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

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

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

2 理論與算法

2.1 圖像分割中的Markov隨機(jī)場(chǎng)模型

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

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

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

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

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

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

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

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

勢(shì)函數(shù):

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

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

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

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

2.2 MRF模型的ICE迭代算法

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

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

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

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

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

遍歷整個(gè)圖像得到標(biāo)記場(chǎng)?k;

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

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

2.3 聲吶圖像自動(dòng)分類的模型與算法

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

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

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

式中:θ為判斷閾值。

自動(dòng)分類模型算法的實(shí)現(xiàn):

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

2)計(jì)算像素峰值:

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

3)計(jì)算灰度區(qū)間［rpv－l，rpv+l］的概率統(tǒng)計(jì)和: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)計(jì)算左右邊像素谷值:

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

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

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

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

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

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

其算法如下:

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

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

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

4)然后根據(jù)區(qū)域能量公式計(jì)算各個(gè)區(qū)域的能量函數(shù)，并將能量最小時(shí)的區(qū)域均值μ1和均方差σ1作為MRF模型的陰影區(qū)域參數(shù);將能量取得極大值時(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對(duì)聲吶圖像進(jìn)行分割，將分割后的模型分別與陰影模型gshadow和目標(biāo)模型gtarget進(jìn)行簡(jiǎn)單的與運(yùn)算與及高斯低通濾波處理，即可以得到最終的聲吶圖像分割結(jié)果。

3 實(shí)驗(yàn)結(jié)果與分析

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

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

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

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

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

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

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

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

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

4 結(jié)束語(yǔ)

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

［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］陽(yáng)凡林，獨(dú)知行，李家彪，等.基于 MRF場(chǎng)的側(cè)掃聲吶分割方法［J］.海洋學(xué)報(bào)，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］尹芳，陳德運(yùn)，吳銳.改進(jìn)的譜聚類圖像分割方法［J］.計(jì)算機(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é)報(bào)，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.