多尺度形態(tài)梯度和標(biāo)記分水嶺的時頻譜圖分割

郭海濤，徐 雷，趙紅葉，田原嫄，焦圣喜

(1.內(nèi)蒙古大學(xué)電子信息工程學(xué)院，內(nèi)蒙古 呼和浩特 010021；2.東北電力大學(xué)電氣工程學(xué)院，吉林 吉林 132012；3.東北電力大學(xué)自動化工程學(xué)院，吉林 吉林 132012；4.東北電力大學(xué)機(jī)械工程學(xué)院，吉林 吉林 132012)

?

多尺度形態(tài)梯度和標(biāo)記分水嶺的時頻譜圖分割

郭海濤1，2，徐 雷3，趙紅葉3，田原嫄4，焦圣喜3

(1.內(nèi)蒙古大學(xué)電子信息工程學(xué)院，內(nèi)蒙古 呼和浩特 010021；2.東北電力大學(xué)電氣工程學(xué)院，吉林 吉林 132012；3.東北電力大學(xué)自動化工程學(xué)院，吉林 吉林 132012；4.東北電力大學(xué)機(jī)械工程學(xué)院，吉林 吉林 132012)

時頻譜圖干擾強，目標(biāo)之間、目標(biāo)與干擾之間有重疊，其分割是重要而困難的問題.提出一種基于圖像熵定義的時頻譜圖多尺度形態(tài)梯度圖像融合方法，將該方法與標(biāo)記分水嶺分割結(jié)合形成一種基于多尺度形態(tài)梯度和標(biāo)記分水嶺的時頻譜圖分割方法.實驗結(jié)果表明，與基于單尺度形態(tài)梯度和標(biāo)記分水嶺的分割方法相比，新方法實用性更強;與Otsu法相比，新方法分割更準(zhǔn)確.

時頻譜圖;多尺度梯度圖像;圖像熵;水嶺算法;圖像分割

眾所周知，利用時頻分析方法可以獲得時間信號(包括噪聲)的時頻分布.將該時頻分布作為圖像，即為信號的時頻譜圖.時頻譜圖的處理有助于對信號進(jìn)行時頻分析，而信號的時頻分析有廣闊的應(yīng)用前景.信號的種類(應(yīng)用的領(lǐng)域)不同，時頻分析的方法也可以不同，因此時頻譜圖的種類多種多樣，但其處理本質(zhì)上屬于圖像處理范疇.時頻譜圖分割的目的是分離出時頻曲線.顯然，時頻譜圖的分割是時頻譜圖低層處理的核心與關(guān)鍵.

時頻譜圖分割是重要而困難的問題，有關(guān)其分割研究的報道甚少.不同領(lǐng)域的學(xué)者以不同的應(yīng)用為背景研究這個問題.文獻(xiàn)[1]利用Otsu閾值化方法分割柴油機(jī)振動信號時頻譜圖，文獻(xiàn)[2]利用Otsu閾值化方法結(jié)合區(qū)域分裂與增長實現(xiàn)心音信號時頻譜圖的分割.文獻(xiàn)[1-2]中的方法適用于直方圖為理想雙峰的時頻譜圖的分割，而對于直方圖為非理想雙峰的時頻譜圖的分割效果較差.文獻(xiàn)[3]利用腐蝕、膨脹、開閉等數(shù)學(xué)形態(tài)學(xué)運算對二值化后的心音信號時頻譜圖進(jìn)行分割，其方法不適用于灰度時頻譜圖.文獻(xiàn)[4]利用BP神經(jīng)網(wǎng)絡(luò)分割時頻譜圖，該方法存在以下原理上缺欠:(1)需要大量的學(xué)習(xí)樣本;(2)樣本的選擇無理論指導(dǎo);(3)方法的性能與多種因素有關(guān)，存在不確定性.文獻(xiàn)[5]利用分水嶺算法分割海豚音等自然信號的Capon時頻譜圖，文獻(xiàn)[6]利用分水嶺算法分割語音信號時頻譜圖，文獻(xiàn)[7]應(yīng)用文獻(xiàn)[6]的方法分割步態(tài)周期時頻譜圖.文獻(xiàn)[5-7]中的方法容易得到連續(xù)的單像素寬邊界，但其梯度算子容易受到干擾和量化誤差的影響而產(chǎn)生局部極小值，從而導(dǎo)致過分割.為了解決過分割，文獻(xiàn)[8-9]將分水嶺算法與區(qū)域合并、區(qū)域增長結(jié)合使用，取得較好的效果.這些方法均是基于分水嶺算法的，筆者因此受到啟示，考慮利用分水嶺算法分割時頻譜圖，但使用的分水嶺方法是基于多尺度形態(tài)梯度的，這一點與前述方法不同.

考慮到噪聲和過多偽極小值點的存在可導(dǎo)致分水嶺算法時頻譜圖嚴(yán)重的過分割，因此先對時頻譜圖進(jìn)行濾波.這里利用文獻(xiàn)[10]中的形態(tài)學(xué)開閉重構(gòu)濾波器.分水嶺算法圖像分割的性能在很大程度上依賴于待分割圖像的梯度圖像[11]，因此梯度圖像的形成是分水嶺算法圖像分割中的重要環(huán)節(jié).考慮到形態(tài)梯度圖像使圖像灰度級階躍更為急劇、圖像邊緣更加突出，多尺度梯度圖像能夠更好地保留圖像邊緣細(xì)節(jié)，因此在時頻譜圖分割中利用多尺度形態(tài)梯度圖像.為了形成時頻譜圖的多尺度形態(tài)梯度圖像，筆者提出一種基于圖像熵定義的時頻譜圖多尺度形態(tài)梯度圖像融合方法.上述時頻譜圖的濾波和梯度圖像的形成統(tǒng)稱為時頻譜圖的預(yù)處理.

1 開閉重構(gòu)濾波器

利用形態(tài)學(xué)開閉重構(gòu)濾波器[10]對時頻譜圖進(jìn)行平滑，抑制圖像中的噪聲及偽極小值點.形態(tài)學(xué)開閉重構(gòu)濾波器的步驟如圖1所示.首先對圖像進(jìn)行測地腐蝕得到標(biāo)記圖像，然后以原始圖像作為掩膜圖像對標(biāo)記圖像進(jìn)行形態(tài)學(xué)開重構(gòu).在此基礎(chǔ)上，對開重構(gòu)后的圖像進(jìn)行測地膨脹，將膨脹得到的結(jié)果作為掩膜圖像，對開重構(gòu)后的圖像進(jìn)行形態(tài)學(xué)閉重構(gòu).

圖1 開閉重構(gòu)濾波器

圖1中相關(guān)操作的定義如下:若f為掩膜圖像，b為結(jié)構(gòu)元素，r為標(biāo)記圖像，則測地膨脹和測地腐蝕運算可分別表示為

dr(f)=min{d(f)，r}，er(f)=max{e(f)，r}，

其中d(f)和e(f)分別表示利用結(jié)構(gòu)元素b對圖像f作膨脹和腐蝕運算.相應(yīng)的膨脹重構(gòu)和腐蝕重構(gòu)分別表示為

膨脹重構(gòu)和腐蝕重構(gòu)均為迭代運算，當(dāng)運算收斂時迭代終止.開重構(gòu)與閉重構(gòu)分別表示為

2 多尺度形態(tài)梯度圖像融合方法

基于分水嶺變換的圖像分割方法，其性能在很大程度上依賴于形成梯度圖像的梯度算法[11].傳統(tǒng)的梯度算子，如高斯濾波器一階偏微分和形態(tài)學(xué)梯度算子存在嚴(yán)重的缺陷，它們會受噪聲和量化誤差的影響，在均勻一致的區(qū)域內(nèi)部產(chǎn)生過多的局部最小值，從而導(dǎo)致過分割現(xiàn)象.針對時頻譜圖干擾強的特點，筆者提出一種基于圖像熵定義的時頻譜圖多尺度形態(tài)梯度圖像融合方法.

文獻(xiàn)[12-13]定義單尺度形態(tài)梯度算子為

Grad(f)=(f⊕b)-(fΘb)，

(1)

其中⊕和Θ分別為膨脹和腐蝕運算，其性能取決于結(jié)構(gòu)元素b的大小.小尺度的結(jié)構(gòu)元素去噪能力弱，但能有效地檢測梯度邊緣細(xì)節(jié);大尺度的結(jié)構(gòu)元素去噪能力強，但梯度圖像邊緣細(xì)節(jié)有所丟失.為了更好地利用小尺度結(jié)構(gòu)元素和大尺度結(jié)構(gòu)元素的各自優(yōu)點，近年來學(xué)者們提出平均加權(quán)的多尺度形態(tài)梯度算子[11].這種簡單的取均值進(jìn)行多尺度形態(tài)梯度圖像融合的方法不能體現(xiàn)出不同尺度梯度圖像的差別，因此不能很好地反映不同尺度下梯度圖像的綜合信息.結(jié)合圖像熵的定義，筆者提出一種新的多尺度形態(tài)梯度算子.該算子的實現(xiàn)步驟是:分別求取單尺度形態(tài)梯度圖像，利用圖像熵的定義計算各尺度形態(tài)梯度圖像的加權(quán)值，然后通過加權(quán)求和將單尺度形態(tài)梯度圖像融合為多尺度形態(tài)梯度圖像.多尺度形態(tài)梯度圖像表示為

(2)

利用(2)式即可融合得到時頻譜圖多尺度形態(tài)梯度圖像.與文獻(xiàn)[11]不同，這里選擇圓盤形結(jié)構(gòu)元素bj.取3個不同的尺度，即取序號j =1，2，3，相應(yīng)的圓盤形結(jié)構(gòu)元素bj的半徑分別為1，2，3.對應(yīng)于這3個尺度下的形態(tài)梯度圖像的熵值和加權(quán)值分別為E1，E2，E3和a1，a2，a3.

3 標(biāo)記分水嶺分割方法

3.1 傳統(tǒng)分水嶺分割

分水嶺算法是基于拓?fù)淅碚摰臄?shù)學(xué)形態(tài)學(xué)分割方法，其基本思想是將圖像看作測地學(xué)上的拓?fù)涞孛?分水嶺算法有很多種，具有代表性的是基于浸沒模擬(Immersion Simulation)[13]的過程.該方法將灰度圖像看成地形學(xué)上的地貌，是一種自下而上形成區(qū)域的方法.文獻(xiàn)[13]給出的基于浸沒模擬的算法級定義是較為嚴(yán)格的遞歸定義，具有一定的代表性.該算法包括2個部分:第一部分為排序，第二部分為泛洪.算法可簡述如下(其中步驟(ⅰ)為排序，步驟(ⅱ)—(ⅵ)為泛洪):

(ⅰ)首先計算圖像中各點的梯度，然后掃描整幅圖像得到各梯度的概率密度.各像素點在排序數(shù)組中的位置由梯度分布的累積概率與該像素點的梯度值計算得到.計算出所有像素點的排序位置并將其存入排序數(shù)組.在排序后的數(shù)組中，梯度值越低的點存放的位置越靠前.

(ⅱ)像素點按梯度值從低到高的順序處理，相同梯度值的點作為一個梯度層級.

(ⅲ)處理一個梯度層級hcur(當(dāng)前層)，將該層中所有鄰域已被標(biāo)識的點加入到一個先進(jìn)先出的隊列中去.

(ⅳ)若先進(jìn)先出隊列非空，則彈出隊列的首元素作為當(dāng)前處理像素.順序處理當(dāng)前像素所有高度為hcur的相鄰點.若鄰點已被標(biāo)識，則根據(jù)該鄰點標(biāo)識刷新當(dāng)前像素點的標(biāo)識.若鄰點尚未標(biāo)識，則將該鄰點加入到先進(jìn)先出隊列中去.循環(huán)執(zhí)行本步驟直至隊列空為止.

(ⅴ)再一次掃描當(dāng)前梯度級的像素點，檢查是否仍有未標(biāo)識點.此時的未標(biāo)識點意味著一個新的極小區(qū).因此，若發(fā)現(xiàn)未標(biāo)識點，則將當(dāng)前區(qū)域標(biāo)識值加1，并將該值賦為未標(biāo)識點的標(biāo)識值.然后從該點出發(fā)執(zhí)行與步驟(ⅳ)相同的泛洪步驟，標(biāo)識該極小區(qū)的所有像素點.

(ⅵ)返回步驟(ⅲ)處理下一梯度層級，直至將所有梯度層級都處理完為止.

3.2 標(biāo)記分水嶺分割

3.2.1 標(biāo)記提取 為了改善傳統(tǒng)分水嶺過分割現(xiàn)象，采用標(biāo)記的方法分別提取圖像的內(nèi)部標(biāo)記和外部標(biāo)記，輔助修改多尺度形態(tài)梯度圖像，抑制無意義的局部極小值.

(1)內(nèi)部標(biāo)記提取:對多尺度形態(tài)梯度圖像用Extend-H-maxima變換[15]計算局部最大值，將局部最大值圖像轉(zhuǎn)化為二值圖像，從而提取內(nèi)部標(biāo)記.

(2)外部標(biāo)記提取:對多尺度形態(tài)梯度圖像用歐氏距離變換后提取外部標(biāo)記[15].

3.2.2 圖像分割 采用數(shù)學(xué)形態(tài)學(xué)中強制最小技術(shù)[16]修改多尺度形態(tài)梯度圖像.具體步驟是:將內(nèi)部標(biāo)記和外部標(biāo)記強制作為多尺度形態(tài)梯度圖像的局部最小值，強制修改多尺度形態(tài)梯度圖像，對修正后的梯度圖像進(jìn)行分水嶺分割得到最終分割結(jié)果.

4 時頻譜圖分割步驟

時頻譜圖分割具體步驟為:

(ⅰ)利用形態(tài)重構(gòu)濾波器去除時頻譜圖中的噪聲.結(jié)合時頻譜圖噪聲特點，這里選用半徑為1的圓盤形結(jié)構(gòu)元素進(jìn)行開閉重構(gòu)濾波，濾波中令r=f-1.

(ⅱ)利用(1)式求圓盤形結(jié)構(gòu)元素半徑分別為1，2，3時的單尺度形態(tài)梯度圖像，并利用(2)式對3個單尺度形態(tài)梯度圖像融合得到多尺度形態(tài)梯度圖像.

(ⅲ)對步驟(ⅱ)中得到的多尺度形態(tài)梯度圖像進(jìn)行Extend-H-maxima變換[15](參數(shù)h= 10).

(ⅳ)對步驟(ⅱ)中得到的多尺度形態(tài)梯度圖像進(jìn)行歐氏距離變換，提取外部標(biāo)記.

(ⅴ)對步驟(ⅲ)，(ⅳ)得到的內(nèi)部標(biāo)記及外部標(biāo)記先進(jìn)行閉操作，后膨脹操作.

(ⅵ)對步驟(ⅴ)中得到的標(biāo)記圖像利用Matlab中bwareaopen函數(shù)去除周邊產(chǎn)生的孤立像素.

(ⅶ)采用形態(tài)學(xué)極小值標(biāo)定技術(shù)修正多尺度形態(tài)梯度圖像.

(ⅷ)對修正后的梯度圖像進(jìn)行分水嶺算法分割.

5 實驗結(jié)果與討論

圖2a是海豚音的Capon時頻譜圖[17]，水平方向表示時間，垂直方向表示頻率.譜圖中大體呈現(xiàn)為水平方向的間斷的曲線是時頻曲線，即為目標(biāo)，其他部分均為背景.從圖2a可以看出，干擾非常嚴(yán)重，存在大量的呈現(xiàn)出不完全紋理現(xiàn)象的強干擾(大體垂直方向的不規(guī)則條紋).干擾的強度與目標(biāo)的強度相當(dāng)，且與目標(biāo)有重疊.對于這樣的時頻譜圖，其分割是一個很大的挑戰(zhàn).

按時頻譜圖分割的步驟(ⅰ)得到去濾波后的時頻譜圖(圖2b).按步驟(ⅱ)得到時頻譜圖的多尺度形態(tài)梯度圖像(圖2c).按步驟(ⅲ)—(ⅷ)得到新方法分割結(jié)果(圖2d).將圖2d中分割得到的目標(biāo)置為白色，其他部分置為黑色，得到圖2e.

作為比較，這里給出結(jié)構(gòu)元素半徑分別為1，2，3時單尺度形態(tài)梯度標(biāo)記分水嶺分割后的結(jié)果(圖2f—h);另外，還給出最常用的Otsu法分割后的結(jié)果(圖2i)，閾值為98.在上述的分割結(jié)果2e—i中僅僅是閾值化的結(jié)果，沒有去除孤立區(qū)，但這并不影響考察分割的結(jié)果.總的來看，各種方法分割出的時頻曲線均有一定的失真，沒有分割出完整的時頻曲線;這是因為時頻譜圖中干擾過于嚴(yán)重，干擾與時頻曲線嚴(yán)重重疊.相比較而言，Otsu法分割出的時頻曲線嚴(yán)重變寬，失真較大，這是因為時頻譜圖的直方圖不是理想的雙峰型，Otsu法不適用造成的.結(jié)構(gòu)元素半徑為2和3時單尺度形態(tài)梯度標(biāo)記分水嶺分割出的時頻曲線殘缺嚴(yán)重，丟失過半.在上述這些方法中，新方法和結(jié)構(gòu)元素半徑為1時單尺度形態(tài)梯度標(biāo)記分水嶺方法的分割效果相對較好.這兩者相比，在分割出的時頻曲線寬度失真方面，新方法失真較嚴(yán)重;在分割出的時頻曲線細(xì)節(jié)保持方面，新方法較好;在方法的實用性方面，由于合適的單一尺度的選擇相對比較困難，盲目性大，而多尺度的選擇相對容易，因此新方法更有優(yōu)勢.綜合來看，新方法是一個可取的折中.

a 海豚音的Capon時頻譜圖

b 濾波后時頻譜圖

c 多尺度形態(tài)學(xué)梯度圖像

d 新方法分割

e 圖2d結(jié)果重新賦值

f 結(jié)構(gòu)元素半徑為1時的分水嶺分割

g 結(jié)構(gòu)元素半徑為2時的分水嶺分割

h 結(jié)構(gòu)元素半徑為3時的分水嶺分割

i Otsu法分割

6 結(jié)語

提出多尺度形態(tài)梯度圖像的一種融合方法，并將該方法與標(biāo)記分水嶺分割結(jié)合形成一種基于多尺度形態(tài)梯度和標(biāo)記分水嶺的時頻譜圖分割方法.新方法思路是:利用圖像熵定義對時頻譜圖3個不同尺度下的形態(tài)梯度圖像進(jìn)行融合，獲得時頻譜圖的多尺度形態(tài)梯度圖像，然后利用標(biāo)記分水嶺方法分割時頻譜圖的多尺度形態(tài)梯度圖像.

新方法較最常用的Otsu法更適合于時頻譜圖分割，較基于單尺度形態(tài)梯度和標(biāo)記分水嶺的時頻譜圖分割方法更具實用性.由于時頻譜圖中干擾過于嚴(yán)重，干擾與時頻曲線嚴(yán)重重疊，因此新方法分割出的時頻曲線同樣存在一定程度的失真.綜合考慮方法的失真程度和實用性，新方法是一個較好的折中選擇.

基于多尺度形態(tài)梯度和標(biāo)記分水嶺的時頻譜圖分割方法的分割效果還有待提高.認(rèn)為以下2個問題需要進(jìn)一步深入研究:(1)形成時頻譜圖的多尺度形態(tài)梯度圖像時，如何根據(jù)時頻譜圖自身的特點選取不同尺度和不同形狀的結(jié)構(gòu)元素;(2)標(biāo)記分水嶺分割時，如何更有效地進(jìn)行時頻譜圖多尺度形態(tài)梯度圖像的標(biāo)記.這2個問題的解決可望大幅度提高時頻譜圖的分割效果.

[1] 蔡艷平，李艾華，王 濤，等.基于時頻譜圖與圖像分割的柴油機(jī)故障診斷[J].內(nèi)燃機(jī)學(xué)報，2011，29(2):181-186.

[2] GAVROVSKA A M，PASKAS M P，DUJKOVIC D，et al.Region-Based Phonocardiogram Event Segmentation in Spectrogram Image[C].IEEE:The 10th Symposium on Neural Network Applications in Electrical Engineering，2010:69-72.

[3] GAVROVSKA A M，PASKAM P，RELJIN I S.Determination of Morphologically Characteristic PCG Segments from Spectrogram Image[J].Telfor Journal，2010，2(2):4-77.

[4] SHAFI I，AHMAD J，SHAH S I，et al.Time Frequency Image Analysis Using Neural Networks[C].IEEE:IMACS Multiconference on Computational Engineering in Systems Applications，2006:315-320.

[5] LEPRETTRE B，MARTIN N.Extraction of Pertinent Subsets from Time-Frequency Representations for Detection and Recognition Purposes[J].Signal Processing，2002，82(2):229-238.

[6] STEINBERG R，O’SHAUGHNESSY D.Segmentation of a Speech Spectrogram Using Mathematical Morphology[C].IEEE:IEEE International Conference on Acoustics，Speech and Signal Processing，2008:1 637 -1 640.

[7] YUWONO M，SU S W，MOULTON B D，et al.Gait Cycle Spectrogram Analysis Using a Torso-Attached Inertial Sensor[C].IEEE:Annual International Conference of the IEEE Engineering in Medicine and Biology Society，2012:6 539-6 542.

[9] ZHANG Zhenggang，XU Kailiang，TA Dean，et al.Joint Spectrogram Segmentation and Ridge-Extraction Method for Separating Multimodal Guided Waves in Long Bones[J].Science China-Physics，Mechanics and Astronomy，2013，56(7):1 317-1 323.

[10] MUKHOPADHYAY S，CHANDA B.An Edge Preserving Noise Smoothing Technique Using Multiscale Morphology[J].Signal Processing，2002，82:527-544.

[11] WANG D.A Multiscale Gradient Algorithm for Image Segmentation Using Watersheds[J].Pattern Recognition，1997，30(12):2 043-2 052.

[12] SALEMBIER P.Morphological Multiscale Segmentation for Image Coding[J].Signal Processing，1994，38(3):359-386.

[13] VINCENT L，SOILLE P.Watersheds in Digital Spaces:An Efficient Algorithm Based on Immersion Simulations[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，1991，13(6):583-598.

[14] PUN T，KAPUR J N，WESZKA J S.Entropic Thresholding:A New Approach[J].Computer Vision，Graphics and Image Processing，1981，16(3):210-239.

[15] WU Qiang，MERCHANT F A，CASTLEMAN K R.Microscope Image Processing[M].Waltham，Massachusetts:Academic Press Elsevier，2008:135-142.

[16] SOILLE P.Morphological Image Analysis:Principle and Applications[M].New York:Springer-Verlag，2003.

[17] LEPRETTRE B，MARTIN N.Extraction of Pertinent Subsets from Time-Frequency Representations for Detection and Recognition Purposes[J].Signal Processing，2002，82(2):229-238.

(責(zé)任編輯 向陽潔)

Time-Frequency Spectrogram Segmentation Using the Multi-Scale Morphological Gradient and the Marked Watershed Algorithm

GUO Haitao1，2，XU Lei3，ZHAO Hongye3，TIAN Yuanyuan4，JIAO Shengxi3

(1.College of Electronic Information Engineering，Inner Mongolia University，Hohhot 010021，China;2.Electrical Engineering College，Northeast Dianli University，Jilin 132012，Jilin China;3.School of Automation Engineering，Northeast Dianli University，Jilin 132012，Jilin China;4.School of Mechanical Engineering，Northeast Dianli University，Jilin 132012，Jilin China)

It is important and difficult to segment a time-frequency spectrogram due to its strong interference as well as serious overlap among targets and between targets and interference.An image entropy based fusion method for multi-scale morphological gradient image of a time-frequency spectrogram is presented.By combining that method with the marked watershed algorithm，a method for segmenting time-frequency spectrogram based on the multi-scale morphological gradient image and the marked watershed algorithm，is obtained.The experiment results show that it is more practical than the segmentation method which is based on single scale morphological gradient and marker watershed，and more accurate than Otsu method.

time-frequency spectrogram;multi-scale morphological gradient image;image entropy;watershed algorithm;image segmentation

1007-2985(2015)04-0012-06

郭海濤(1965—)，男，黑龍江安達(dá)人，內(nèi)蒙古大學(xué)電子信息工程學(xué)院教授，博士，主要從事圖像處理、模式識別等研究.

TP391.4

A

10.3969/j.issn.1007-2985.2015.04.004