程一峰, 劉增力(昆明理工大學(xué) 信息工程與自動(dòng)化學(xué)院, 云南 昆明 650500)

1 引 言

近年來(lái),稀疏表示在圖像去噪中的應(yīng)用越來(lái)越廣泛,成為前沿研究課題。Mallat等[1]提出超完備字典,基于超完備字典的稀疏表示是圖像稀疏去噪的重要研究方向,而字典的設(shè)計(jì)則是該方向上重要的研究環(huán)節(jié)。全局字典以及自適應(yīng)字典是字典的兩類(lèi)重要分支,最近幾年,超完備字典通過(guò)學(xué)習(xí)和訓(xùn)練而獲得的方法得到很大的進(jìn)展。Aharon M,Elad M等[2]提出了K-SVD(K-singular value decomposition)算法,達(dá)到了良好的去噪效果。然而傳統(tǒng)的K-SVD算法在一些方面仍然存在不足,具體表現(xiàn)如下:字典在更新過(guò)程中會(huì)出現(xiàn)每一列無(wú)法全部更新或部分列失效的問(wèn)題,原因是初始化字典的選擇不恰當(dāng),最終結(jié)果會(huì)導(dǎo)致信號(hào)利用率低;在噪聲比較強(qiáng)的時(shí)候,字典中含有與噪聲匹配的大量噪聲原子,導(dǎo)致字典在學(xué)習(xí)過(guò)程中無(wú)法有效地將噪聲去除;K-SVD算法運(yùn)算時(shí)間較長(zhǎng),體現(xiàn)在圖像訓(xùn)練過(guò)程中會(huì)占用大量的時(shí)間;傳統(tǒng)最速下降算法鋸齒現(xiàn)象嚴(yán)重。

本文針對(duì)傳統(tǒng)的K-SVD算法的幾點(diǎn)不足,對(duì)其進(jìn)行了改進(jìn):首先,通過(guò)稀疏貝葉斯學(xué)習(xí)(sparse Bayesian learning,SBL)[3]對(duì)圖像進(jìn)行預(yù)處理,提高圖像信號(hào)的利用率;然后,對(duì)字典進(jìn)行優(yōu)化,并結(jié)合Bartlett檢驗(yàn)法[4]將噪聲原子裁剪,提高去噪后圖像質(zhì)量;第三,針對(duì)運(yùn)算時(shí)間較長(zhǎng)的問(wèn)題,研究凸優(yōu)化算法中的最速下降法,將改進(jìn)的最速下降算法與正交匹配追蹤(orthogonal matching pursuit,OMP)相結(jié)合;第四,將最速下降算法進(jìn)行一維搜索,有效地抑制了鋸齒現(xiàn)象。實(shí)驗(yàn)結(jié)果表明,通過(guò)本文的圖像去噪方法,圖像精度和峰值信噪比(peak signal to noise ratio,PSNR)值都得到了提高。

2 圖像的稀疏表示模型

2.1 K-均值聚類(lèi)算法

通過(guò)K-均值(K-means)聚類(lèi)算法將圖像塊進(jìn)行聚類(lèi)處理,將具有相似集合形態(tài)的圖像塊聚集在一起,再通過(guò)K-SVD算法對(duì)每一個(gè)相似塊訓(xùn)練出各自的優(yōu)化字典,最終得到的字典更具有針對(duì)性。

K-均值聚類(lèi)算法的主要方法步驟如下:首先從全體n個(gè)數(shù)據(jù)中選擇k個(gè)聚類(lèi)中心。根據(jù)剩下的n-k個(gè)對(duì)象與k個(gè)聚類(lèi)中心的相似度,將每個(gè)對(duì)象分配到k個(gè)聚類(lèi)中心中與其最為相似的聚類(lèi)中去。如果類(lèi)中心有變化則對(duì)其進(jìn)行反復(fù)聚類(lèi)計(jì)算,直到類(lèi)中心不再發(fā)生變化以及平方誤差函數(shù)收斂為止。平方誤差函數(shù)的定義可以表示為[5,6]:

(1)

式中:xji為第j類(lèi)第i個(gè)樣本;kj為第j類(lèi)的均值或稱(chēng)聚類(lèi)中心;nj為第j類(lèi)樣本個(gè)數(shù);E為數(shù)據(jù)對(duì)象及其相應(yīng)的聚類(lèi)中心的均方差總和。

2.2 稀疏貝葉斯學(xué)習(xí)

在K-SVD算法中,字典在更新過(guò)程中經(jīng)常會(huì)出現(xiàn)每一列不能全部更新或部分列失效的現(xiàn)象,為了解決這個(gè)問(wèn)題,可以通過(guò)貝葉斯學(xué)習(xí)[7]對(duì)信號(hào)進(jìn)行預(yù)處理,通過(guò)貝葉斯概率模型進(jìn)行迭代直到信號(hào)的稀疏表示。

1) 計(jì)算βij的協(xié)方差矩陣與βij:

βij協(xié)方差矩陣為

(2)

式中:βij為稀疏向量;Γ為k×k階對(duì)角矩陣,γ1,γ2,…,γk為對(duì)角元素;D為過(guò)完備字典矩陣;σ為噪聲標(biāo)準(zhǔn)差。

(3)

式中:Xij為圖像子塊。

2) 更新噪聲:

(4)

式中:Rij為提取子圖的矩陣;X為含噪圖像塊。

(5)

2.3 字典噪聲原子裁剪

傳統(tǒng)K-SVD的方法在訓(xùn)練含噪圖像之后,得到的字典中會(huì)出現(xiàn)與噪聲匹配的噪聲原子,因此對(duì)于字典進(jìn)行優(yōu)化的關(guān)鍵步驟就是刪除這些噪聲原子。然而關(guān)鍵問(wèn)題是:如何對(duì)噪聲原子進(jìn)行檢測(cè)。文獻(xiàn)[8]中提出利用假設(shè)檢驗(yàn)來(lái)檢測(cè)噪聲原子。噪聲原子采用對(duì)角、反對(duì)角、水平、垂直4個(gè)掃描方向進(jìn)行對(duì)比掃描,更能提高相關(guān)性檢測(cè)的能力。

設(shè)μi為特征向量的方差,在噪聲原子4個(gè)方向上用Bartlett檢驗(yàn)法測(cè)試特征向量的方差μi是否相等,或者至少有2個(gè)原子的方差不等,Bartlett檢驗(yàn)統(tǒng)計(jì)量H的定義為:

(6)

式中:n為原子元素?cái)?shù)量。

對(duì)原子庫(kù)中的原子進(jìn)行判斷,若某個(gè)原子滿足H<χ2(β;3),則認(rèn)為該原子為噪聲原子,并且進(jìn)行刪除。χ2(β;3)表示自由度為3的卡方分布上的β百分點(diǎn)位置。與通過(guò)挖掘原子庫(kù)的冗余性來(lái)減小原子數(shù)目的傳統(tǒng)字典優(yōu)化方法對(duì)比可以發(fā)現(xiàn),該方法不會(huì)對(duì)信息原子產(chǎn)生影響,只會(huì)從字典中刪除噪聲原子,從而提高了逼近質(zhì)量以及去噪效果。

2.4 改進(jìn)的最速下降正交匹配追蹤

傳統(tǒng)的最速下降算法雖然擁有很多的優(yōu)點(diǎn)和可取之處,但是該算法依然存在缺點(diǎn)。在迭代過(guò)程中迭代點(diǎn)呈現(xiàn)了相應(yīng)的“之”字型,此現(xiàn)象為鋸齒現(xiàn)象,對(duì)收斂速度會(huì)產(chǎn)生影響,但是鋸齒現(xiàn)象又是無(wú)法避免的。因此可以通過(guò)抑制鋸齒現(xiàn)象從而可以提高收斂速度。

稀疏分解圖像的目標(biāo)是找到一個(gè)矩陣α并且滿足:

(7)

式中:||α||0為α的非零元個(gè)數(shù)。

目前,稀疏分解方法可分為基追蹤(BP),匹配追蹤(MP)和正交匹配追蹤(OMP)。OMP是最好的匹配追蹤算法并且具有最適用的基礎(chǔ)表達(dá)功能。式(7)可以轉(zhuǎn)換為無(wú)約束最優(yōu)化問(wèn)題,即:

(8)

式中:λ為正則化參數(shù)。

盡管OMP算法比其他方法快,但其仍以損失一定的精度為代價(jià)。考慮到凸優(yōu)化算法具有較高的精度但收斂速度較慢的特點(diǎn),本文將改進(jìn)的最速下降算法與OMP算法相結(jié)合,即改進(jìn)的最速下降OMP。利用改進(jìn)的最速下降算法去求解式(8)。

通常將式(8)看成是一種多元函數(shù)f(x)的無(wú)約束優(yōu)化問(wèn)題,其尋優(yōu)表達(dá)式為[9,10]:

αk+1=αk+tkpk

(9)

式中:tk為優(yōu)化迭代的步長(zhǎng);pk為第k次迭代搜索方向。

通常認(rèn)為在一個(gè)點(diǎn)上最快的下降方向是這一點(diǎn)的負(fù)梯度方向。于是可以定義:

pk=-f(αk)

(10)

每次優(yōu)化迭代的步長(zhǎng)tk可以通過(guò)式(11)求得,即:

(11)

如果tk-1f(X(k+1))Tp(k-1)<0,則轉(zhuǎn)入式(12),否則轉(zhuǎn)入式(14)。

令p′k=X(k)-X′,從X(k)出發(fā),沿著p′k搜索方向進(jìn)行一維搜索,求出步長(zhǎng)t′k,使得

(12)

若f(X(k)+t′kp′k)

X(k+1)=X(k)+tkpk

(13)

置k=k+1,轉(zhuǎn)到式(10)。

該算法通過(guò)在p=X(k+1)-X(k-1)方向上進(jìn)行一維搜索,得到一個(gè)比原始算法更好的迭代點(diǎn)。

(14)

式中:G=DTD;b=DTX。

對(duì)式(14)求導(dǎo)可得:

f(Xk)=Gα-b=-DTr

(15)

pk=-f(Xk)=DTr

(16)

式中：r為殘差矩陣。

(17)

結(jié)合改進(jìn)之后的最速下降算法在凸優(yōu)化算法中的優(yōu)點(diǎn)及貪婪迭代算法中OMP的優(yōu)點(diǎn)。改進(jìn)的最速下降OMP可以概括如下:

OMP的分解方法的目的是解決

(18)

步驟1:初始化。殘差r0=X,αx是X(X是含噪圖像塊)的稀疏系數(shù),Dx是X的字典,Dx∈D,ε是誤差控制參數(shù),原子下標(biāo)索引集合I=φ。

步驟2:迭代L次,更新原子下標(biāo)索引集I=I∪{i}。

1)選擇原子i和相應(yīng)的稀疏系數(shù)αx,通過(guò)集合I使得目標(biāo)函數(shù)最小化

(19)

2)更新I,I=I∪{i}

3)更新矩陣的最優(yōu)稀疏系數(shù)αx

(20)

4)更新殘差向量

rk+1=rk-Dxtkpk

(21)

步驟3:更新迭代殘差,當(dāng)||rk+1||2<ε停止迭代,其它情況則返回步驟2。

3 改進(jìn)的K-SVD圖像去噪算法

本文利用稀疏貝葉斯學(xué)習(xí)對(duì)含噪圖像進(jìn)行預(yù)處理,提高其信號(hào)利用率。并且結(jié)合改進(jìn)的最速下降正交匹配追蹤以及噪聲原子裁剪對(duì)字典進(jìn)行更新。具體算法步驟如下:

步驟1:初始化階。首先利用K-means聚類(lèi)算法對(duì)大小為N的圖像Y進(jìn)行聚類(lèi),將圖像分為n類(lèi)(n?N),并令有k個(gè)超參數(shù)的向量為γ=[γ1,γ2,…,γk]T,得到圖像子塊Xij=RijX,i=1,2,…,k;j=1,2,…,n;Rij為圖像子塊矩陣。最后,初始化字典D(0)∈Rn×k為過(guò)完備DCT基,并令J=1。

步驟2:稀疏貝葉斯學(xué)習(xí)。當(dāng)J≤Jmax時(shí),針對(duì)每一類(lèi)圖像塊Xij用稀疏貝葉斯學(xué)習(xí)不斷更新噪聲大小,并且求得稀疏系數(shù)向量αij。

(22)

1)定義稀疏系數(shù)不為零的圖像塊集合為ωl={(i,j)|αij(l)≠0}。

步驟5:字典優(yōu)化。進(jìn)行噪聲監(jiān)測(cè)并刪除噪聲原子。

步驟6:圖像去噪。在去噪過(guò)程中存在著圖像塊的重疊現(xiàn)象,因此應(yīng)該對(duì)圖像塊的重疊部分做平均,去噪之后的圖像的近似解為:

選取PSNR較大的圖像作為最終去噪輸出圖像。

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

實(shí)驗(yàn)是將本文算法與經(jīng)典圖像去噪算法:小波閾值去噪,基于非自適應(yīng)冗余DCT字典去噪,基于傳統(tǒng)K-SVD方法去噪在MATLAB上進(jìn)行仿真對(duì)比。實(shí)驗(yàn)結(jié)果見(jiàn)圖1、圖2和圖3。

圖1 Boat原圖與加噪圖像

圖2 不同字典對(duì)比

從仿真結(jié)果可以看出,本文方法對(duì)于含噪圖像去噪效果明顯,能夠很好地保留圖像的邊緣紋理信息。與傳統(tǒng)去噪方法進(jìn)行對(duì)比之后,可以發(fā)現(xiàn),本文算法的去噪效果優(yōu)于傳統(tǒng)去噪,圖像精度提高。

對(duì)Barbara,Flatland和Boat 3種圖像,選用不同去噪方法的峰值信噪比值比較結(jié)果見(jiàn)表1。

圖3 不同方法的去噪對(duì)比

表1 不同方法的峰值信噪比值比較 dB

通過(guò)表1可以發(fā)現(xiàn),不論是稀疏性較好的Barbara圖亦或是稀疏性較弱的Flatland和Boat圖,利用本文方法,峰值信噪比都明顯比其他方法提高不少,相對(duì)于傳統(tǒng)的K-SVD,本文方法提高了大概1 dB。因此,本文方法相對(duì)于其他去噪方法,綜合去噪能力顯著提高。

5 結(jié) 論

首先通過(guò)稀疏貝葉斯學(xué)習(xí)對(duì)含噪圖像進(jìn)行預(yù)處理, 提高信號(hào)利用率。接著將改進(jìn)之后的OMP算法,以及原子裁剪方法相結(jié)合,構(gòu)成了改進(jìn)之后的K-SVD算法對(duì)圖像進(jìn)行去噪。通過(guò)實(shí)驗(yàn)分析可以看出,相對(duì)于傳統(tǒng)去噪算法,本文方法能夠更好地提高圖像的精度,圖像邊緣和紋理部分保留完好,峰值信噪比顯著提高,并且去噪效率也有所提高。

[參考文獻(xiàn)]

[1] Mallat S G, ZHANG Z F. Matching pursuit with time-frequency dictionary[J].IEEETransactionsonSignalProcessing, 1993, 41(12):3397-3415.

[2] Aharon M, Elad M, Bruckstrein A M. TheK-SVD: an algorithm for designing of over-complete dictionaries for sparse representation[J].IEEETransactiononImageProcessing, 2006, 54(11): 4311-4322.

[3] Wipf D P, Rao B D. Sparse Bayesian learning for basisselection[J].SignalProcessingIEEETranson, 2004, 52,(8):2153-2164.

[4] Bartlett M S. Properties of Sufficiency and Statistical Tests[J].ProceedingsoftheRoyalSocietyofLondon.SeriesA,MathematicalandPhysicalSciences,1937, 160(901):268-282.

[5] 蔣帥.K-均值聚類(lèi)算法研究[D]. 西安:陜西師范大學(xué), 2010.

[6] 田源,王洪濤. 基于量子核聚類(lèi)算法的圖像邊緣特征提取研究[J]. 計(jì)量學(xué)報(bào)，2016,37(6):582-586.

[7] ZHOU F, LI L. Modified sparse Bayesian dictionary learningK-SVD algorithm for video image sparse representation[J].JournalofComputationalInformationSystems, 2015, 11(12):4567-4580.

[8] Lee J, Kim Y H, Nam J H. Adaptive noise reduction algorithms based on statistical hypotheses tests[J].IEEETransactionsonConsumerElectronics, 2008, 54(3):1406-1414.

[9] 龍建輝, 胡錫炎, 張磊. 最速下降法解二次矩陣方程[J]. 湖南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2008, 54(3):85-88.

[10] 鄭慧峰,喻桑桑,王月兵,等. 基于經(jīng)驗(yàn)?zāi)B(tài)分解和奇異值分解的振動(dòng)聲調(diào)制信號(hào)分析方法研究[J]. 計(jì)量學(xué)報(bào)，2016,37(4):398-401.

[11] 周燦梅. 基于壓縮感知的信號(hào)重建算法研究[D].北京:北京交通大學(xué), 2010.