張墨華,彭建華
(1.國家數(shù)字交換系統(tǒng)工程技術(shù)研究中心,鄭州 450002;2.河南財經(jīng)政法大學(xué) 計算機(jī)與信息工程學(xué)院,鄭州 450002)
圖像去噪在數(shù)學(xué)上屬于不適定問題,需要使用圖像先驗(yàn)進(jìn)行正則化,從而生成有意義的解[1]。自然圖像的一些自身屬性可以作為有用的圖像先驗(yàn)用于圖像去噪,如稀疏性、多尺度自相似與非局部自相似性等。文獻(xiàn)[2-4]采用圖像的稀疏性作為先驗(yàn)信息,應(yīng)用小波與曲波等稀疏變換,或是通過學(xué)習(xí)圖像塊得到完備字典以及圖像塊的稀疏表示。文獻(xiàn)[5]采用圖像像素在空間域和強(qiáng)度域中表現(xiàn)出的相似性作為先驗(yàn)信息,通過相似塊進(jìn)行協(xié)同處理。文獻(xiàn)[6]利用同一幅圖像中相同尺度或不同尺度的相似子塊進(jìn)行去噪。文獻(xiàn)[7-8]將圖像的非局部自相似性作為先驗(yàn)信息,對相隔較遠(yuǎn)的相似塊進(jìn)行非局部處理。文獻(xiàn)[9-11]則是使用多種圖像先驗(yàn)組合的方法進(jìn)行圖像去噪。
除了使用預(yù)定義的圖像先驗(yàn),文獻(xiàn)[12]還提出了從自然圖像中學(xué)習(xí)先驗(yàn)的方法。生成式圖像先驗(yàn)學(xué)習(xí)方法通常從一組干凈的外部圖像中學(xué)習(xí)先驗(yàn)?zāi)P?并將其應(yīng)用于給定的退化圖像[13-14],或者從給定的退化圖像中學(xué)習(xí)先驗(yàn)[15]。近年來,判別式圖像先驗(yàn)學(xué)習(xí)方法也開始得到廣泛應(yīng)用。該方法根據(jù)清晰-噪聲圖像對學(xué)習(xí)去噪模型,可進(jìn)一步細(xì)分為基于深度學(xué)習(xí)的方法[16]、基于隨機(jī)場的方法[17]與基于反應(yīng)擴(kuò)散的方法[18]。
在上述方法中,文獻(xiàn)[12]提出的使用高斯混合模型(Gaussian Mixture Model,GMM)先驗(yàn)進(jìn)行自然圖像塊建模的EPLL方法取得了較大成功。這一方法的主要思想是從干凈圖像塊中學(xué)習(xí)先驗(yàn),并將其用于求解圖像復(fù)原問題。但EPLL方法的混合分量數(shù)目固定(文獻(xiàn)[12]中為200),學(xué)習(xí)訓(xùn)練完畢后分量不易進(jìn)行擴(kuò)展。針對這個問題,文獻(xiàn)[19]提出非參貝葉斯方法,為結(jié)構(gòu)化數(shù)據(jù)(如文本文檔、時間序列和圖像)的無監(jiān)督建模提供了靈活框架,其模型的復(fù)雜度可以根據(jù)數(shù)據(jù)規(guī)模增減進(jìn)行自適應(yīng)變化。狄利克雷過程(Dirichlet Process,DP)是最流行的非參貝葉斯方法,其采樣結(jié)果是一個分布,而不是變量,因此DP混合模型可以看作無窮多個混合分量的混合模型[20]。作為分層貝葉斯框架下的非參先驗(yàn),DP可以根據(jù)具體參數(shù)值對該模型生成的數(shù)據(jù)進(jìn)行劃分,具有自然解釋特性,分量的數(shù)量是隨機(jī)的,可以隨著觀測到的新數(shù)據(jù)增長。但DP的學(xué)習(xí)是一個復(fù)雜的最優(yōu)化問題,難以精確求解,只能進(jìn)行近似求解,常用的方法之一是變分推理。變分推理中基于傳統(tǒng)均值場的坐標(biāo)上升算法按照先后順序?qū)?shù)據(jù)集中所有的局部變分參數(shù)與全局變分參數(shù)進(jìn)行更新,因此在大數(shù)據(jù)集上效率較低,同時由于更新緩慢,容易陷入局部最優(yōu)解。
針對上述問題,本文構(gòu)建一種基于狄利克雷過程的可擴(kuò)展高斯混合模型,從清晰圖像數(shù)據(jù)庫中學(xué)習(xí)外部通用先驗(yàn),使模型復(fù)雜度根據(jù)數(shù)據(jù)規(guī)模自適應(yīng)變化。同時,為改善模型的推理過程,提出一種基于批次更新方式的可擴(kuò)展變分算法,將數(shù)據(jù)集分為若干批次,在每個批次訪問中對圖像的局部變分參數(shù)與全局變分參數(shù)進(jìn)行更新,從而求解圖像去噪中所有隱變量的變分后驗(yàn)分布,實(shí)現(xiàn)通用先驗(yàn)學(xué)習(xí)。
定義1(狄利克雷過程) 令G0是集合Θ上的分布,α是正值實(shí)數(shù)。對于Θ上的有限劃分A1,A2,…,Ar,則向量G(A1,A2,…,Ar)是隨機(jī)的。如果有:
G(A1,A2,…,Ar)~Ddir(αG0(A1),αG0(A2),…,αG0(Ar)
(1)
則稱G是基分布G0和集中參數(shù)α上的狄利克雷過程(DP)所生成的分布,記作:G~DDP(α,G0)。
假定從G中獨(dú)立生成N個隨機(jī)變量ηn,則有:
ηn~G,n∈{1,2,…,N}
(2)
聯(lián)合分布{η1,η2,…,ηn}服從波利亞壇子模型(Polya’s urn scheme)[21]。
定理1令G~DDP(α,G0),對于任何集合A,有E[G(A)]=G0(A)
Var[G(A)]=Go(A)(1-Go(A))/(α+1)
從定理1可知,隨著α值增大,方差逐步減小。這樣DP將更多質(zhì)量集中于均值附近。α值越小,G就更為離散。證明可以參見文獻(xiàn)[22]。
定義2(狄利克雷過程混合模型) 如果ηn是第n個觀測變量分布F的參數(shù),即有F(x|ηn),則DPMM可以看作為無窮多個混合分量的混合模型,即有:
G|{α,G0}~DDP(α,G0)
ηn|G~G
xn|ηn~F(x|ηn)
(3)
在DP混合模型中,DP作為分層貝葉斯框架下的非參先驗(yàn),可以根據(jù)具體參數(shù)值對該模型生成的數(shù)據(jù)進(jìn)行劃分,因此,DP混合模型擁有自然解釋特性,其中分量的數(shù)量(劃分中原子的數(shù)目)是隨機(jī)的,并且能夠隨著觀測到的新數(shù)據(jù)進(jìn)行增長。
文獻(xiàn)[23]根據(jù)折棍子構(gòu)造的過程,對DP給出了更為清晰的特征表述。考慮2個獨(dú)立隨機(jī)變量無限集合,vi~Bbeta(1,α),ηi~G0,i={1,2,…},它們的折棍子過程表示如下:
(4)
由式(4)可見,G是離散的。G的支持集包括可數(shù)無限原子集合,這些原子集合從G0獨(dú)立生成?;旌媳壤耰(v)通過將單元長度的“棍子”連續(xù)折斷成無限數(shù)目的分段生成,每分段的大小與棍子剩余部分成正比,從Bbeta(1,α)分布獨(dú)立生成。
在DP混合中,向量βi(v)是一個包含比例的無限向量,{η1,η2,…}是表示混合分量的原子。令zn表示數(shù)據(jù)點(diǎn)xn所關(guān)聯(lián)的混合分量的分派變量。數(shù)據(jù)生成過程如下:
1)生成vi|α~Bbeta(1,α),i={1,2,…}。
2)生成ηi|G0~G0,i={1,2,…}。
3)對于第n個數(shù)據(jù)點(diǎn):
(1)生成zn|{v1,v2,…}~Mmult(β(v))。
(2)生成xn|zn,ηn~F(xn|ηzn)。
在圖像去噪任務(wù)中,每個圖像具有其特定的聚類頻度wi={wi1,wi2,…,wiK,…},通過施加分層DP先驗(yàn),聚類頻度可以看作一個有限的狄利克雷分布:
[wi1,wi2,…,wiK,wi>K]~
Ddir(αβ1,αβ2,…,αβK,αβ>K)
(5)
從式(5)可知,wi的均值為β,方差由集中參數(shù)α所決定,下標(biāo)>K表示除了前K個分量外其余分量累加之和。本文約定觀測數(shù)據(jù)從指數(shù)家族分布生成,DP的基分布是對應(yīng)的共軛先驗(yàn)。
使用基于DP的高斯混合模型對圖像的生成過程進(jìn)行建模。自然圖像可以看作圖像塊組合,每個塊包含D個像素,并且有其對齊基點(diǎn)b,b的取值在D中,即:b~Ccat(1/D,…,1/D)。一幅圖像由所有對齊基點(diǎn)下各個塊重疊生成,每個對齊基點(diǎn)下,可能有部分塊處于部分觀測的狀態(tài),即部分參與最終圖像的生成,部分不參與。
參照文獻(xiàn)[12]的方法,將每個圖像塊去平均化,對齊基點(diǎn)下每個圖像塊gibn看成由高斯混合分布生成,即每個塊表示為均值為0、精度矩陣為S的高斯混合,如下式所示:
(6)
其中,i表示第i幅圖像(共N幅),b表示對齊基點(diǎn)(共D個),n表示圖像i在對齊基點(diǎn)b中的第n個塊(共Nib個)。精度矩陣S的先驗(yàn)為Wishart分布,滿足:
S~Wwis(μ,V)
(7)
圖像塊均值mibn滿足高斯分布:
mibn~Nnorm(τ,σ)
(8)
每個圖像塊可以指派為高斯混合模型中某個分量(聚類)中的一個,用zibn滿足分類分布,表示為:
zibn~Ccat(w1,w2,…,wk)
(9)
給定通過對齊基點(diǎn)b生成的、塊均值為mibn的塊gibn,圖像xi的采樣可以通過式(10)生成:
(10)
在圖像去噪任務(wù)中,觀測圖像yi是干凈圖像xi的退化版本,是由向xi中加入高斯白噪聲后所得。其生成定義如下:
yi~Nnorm(yi|xi,σ2I)
(11)
其中,σ2為噪聲方差。
整個圖像生成過程如圖1(a)所示,其變分參數(shù)的推理過程如圖1(b)所示,第3節(jié)將給出詳細(xì)過程。
圖1 圖像生成有向圖模型及推理過程Fig.1 Directed graph model of image generation andits inference process
本文通過變分分布的參數(shù)推理實(shí)現(xiàn)圖像先驗(yàn)?zāi)P偷膶W(xué)習(xí)。給定干凈觀測圖像xi,需要學(xué)習(xí)模型中的各種隱變量,如全局圖像變量與局部圖像變量。全局圖像變量包括折棍子過程的折斷比例vk,與每個聚類的精度矩陣Sk;局部圖像變量包括聚類頻度wi、圖像塊gibn、塊均值mibn與塊聚類指派zibn。直接計算后驗(yàn)分布難度較大,因此常用變分推理進(jìn)行近似求解。變分推理是一種近似似然和后驗(yàn)的確定性方法[24],采用優(yōu)化方法嘗試尋找一個合適變分分布q,它與真實(shí)的后驗(yàn)在KL散度上無限接近。
定義待學(xué)習(xí)的隱變量集合W={g,S,wi,mi,zi},觀測圖像變量x={x1,x2,…,xN},超參數(shù)θ={α,μ,V}。隱變量的后驗(yàn)分布與其近似分布q的KL散度為:
D(q(W)‖p(W|x,θ))=Eq[lbq(W)]-
Eq[lbp(W,x|θ)]+lbp(x|θ)
(12)
式(12)中KL散度的最小化求解可以轉(zhuǎn)換為對數(shù)邊緣似然下界(ELBO)的最大化求解[25]。如下所示:
lbp(x|θ)≥Eq[lbp(W,x|θ)]-
Eq[lbq(W)]?L
(13)
將L定義為待優(yōu)化求解的目標(biāo)函數(shù),通過參數(shù)的學(xué)習(xí)過程,盡量使L值最大化,從而盡可能減小變分近似后驗(yàn)q與真實(shí)后驗(yàn)的差值。采用均值場方法[25],假定后驗(yàn)完全獨(dú)立,q可以分解為指數(shù)家族密度函數(shù)的乘積。
(14)
其中,隨機(jī)度量的隨機(jī)性通過無限集合{v1,v2,…}和{S1,S2,…}體現(xiàn),而隨機(jī)性的生成過程則通過經(jīng)典的截斷折棍子過程表示[21]。假設(shè)固定截斷位置為T,則混合比例wt(v)在t>T時等于0。T是變分參數(shù),可以自由設(shè)置,并非先驗(yàn)?zāi)P偷囊徊糠帧?/p>
下面具體給出所有變量的變分后驗(yàn)分布,這些分布參數(shù)稱為變分參數(shù),是要進(jìn)行優(yōu)化求解的對象。
在局部變量因子中,圖像塊g、圖像塊均值m、塊分量指派z、圖像聚類頻度w、圖像x的變分后驗(yàn)分布分別為:
(15)
在全局變量因子中,折斷比例變量v的變分分布為:
(16)
另一個全局變量因子精度矩陣S的變分分布為:
(17)
對于自然圖像,所有的對齊基點(diǎn)具有相似的值,設(shè)計對齊基點(diǎn)B的變分后驗(yàn)為
q(B)=Ccat(1/D,…,1/D)
(18)
該設(shè)計簡化了更新操作,同時仍然可以避免不重疊塊所造成的偽影情況。
圖像通用先驗(yàn)的學(xué)習(xí)是基于外部干凈圖像數(shù)據(jù)庫進(jìn)行的,此時圖像x是可觀測變量。根據(jù)式(12),結(jié)合3.1節(jié)所述各變分分布,目標(biāo)函數(shù)L可以分解為3個部分:
L?Limage+Lentropy+Lweight
其中:
Limage?Eq[lbp(x|g,m,b)+
(19)
(20)
(21)
(22)
(23)
其中,期望Eq[β]通過式(4)進(jìn)行計算,Cik是對于圖像i中所有圖像塊使用聚類k的情況的統(tǒng)計量,定義為:
(24)
(25)
(26)
從上式中抽取出兩個統(tǒng)計量Ck和Fk,定義為:
(27)
(28)
這兩個統(tǒng)計量是訓(xùn)練圖像集中所有分派給聚類k的數(shù)據(jù)統(tǒng)計量,將在3.3節(jié)中詳述。
(29)
(30)
(31)
(32)
其中,Oibn表示在對齊基點(diǎn)b下第n個塊中可觀測的像素的數(shù)目。
(33)
αEq[βk]Tk
(34)
其中,Tk定義為:
(35)
(36)
傳統(tǒng)的坐標(biāo)上升算法先對數(shù)據(jù)集中所有局部變分參數(shù)進(jìn)行更新,再對全局變分參數(shù)進(jìn)行更新,因此在大數(shù)據(jù)集上效率較低。同時,由于全局參數(shù)在第一次參與局部參數(shù)更新時是初始化值,計算意義不大。此外,傳統(tǒng)算法更新較為緩慢,容易陷入局部最優(yōu)解。針對這些問題,本文利用式(27)、式(28)、式(35)中所得3個統(tǒng)計量的可累加性,采用批次更新方式完成坐標(biāo)上升算法。將數(shù)據(jù)集分成若干批次,對所有批次完成訪問稱為一輪。在每輪批次訪問中,首先對批次中每個圖像的局部變分參數(shù)進(jìn)行更新,然后進(jìn)行全局變分參數(shù)的更新。
在每個批次中都對式(27)、式(28)、式(35)所得的3個統(tǒng)計量的值進(jìn)行記錄(表示為Ubt),然后累積到全局統(tǒng)計量Uo中,其定義如下:
其中,B表示總批次數(shù)。在訪問完每個批次bt之后,執(zhí)行增量更新來累積統(tǒng)計量,從而反映新批次的匯總情況,并刪除該批次bt之前的值。
(37)
其中,ISbt是圖像數(shù)據(jù)集IS的批次數(shù)據(jù)集。綜上,通用先驗(yàn)訓(xùn)練過程如算法1所示。
算法1通用先驗(yàn)訓(xùn)練的變分推理算法
輸入圖像數(shù)據(jù)集IS
1.Repeat
2.For each batch ISbtin IS do
12.End for
13.Until ELBO收斂
3.4.1 新增機(jī)制
在通用先驗(yàn)訓(xùn)練過程中,通過分量新增機(jī)制可以增加新的有用分量,幫助規(guī)避局部最優(yōu)解。但即使這些分量得到整個圖像數(shù)據(jù)庫的支持,每個批次也無法提供缺失分量所需要的足夠樣本,因此,整個新增過程需要兩輪數(shù)據(jù)訪問。第一輪采集目標(biāo)數(shù)據(jù)樣本,第二輪生成新的分量,并且使用擴(kuò)展的模型對每個批次進(jìn)行更新。具體步驟如下:
2)在第2輪數(shù)據(jù)訪問之前生成新的分量。通過運(yùn)行有限次數(shù)的變分推理,將數(shù)據(jù)集IS′擬合為包括K′(設(shè)置為10)個混合分量的模型。從而將整個模型擴(kuò)展為K+K′個分量。期間并沒有對這些新分量所產(chǎn)生的ELBO的變化進(jìn)行評估,依賴隨后的歸并操作來移除不需要的分量。
3)在第2輪中訪問每個批次,并且對擴(kuò)展后的K+K′個混合模型執(zhí)行局部和全局參數(shù)更新。通過對目標(biāo)數(shù)據(jù)集X′分析得到統(tǒng)計量U′,并入全局統(tǒng)計量Uo中。數(shù)據(jù)集X′有兩種分派結(jié)果:分派給最初的分量(最有可能是k′),或者分派給全新的分量。第2輪結(jié)束后,從Uo中去掉U′,使得Uo和全局參數(shù)與數(shù)據(jù)集IS一致。
一次新增多個分量有助于規(guī)避局部最優(yōu)解。盡管新增分量可能會因?yàn)樵黾硬槐匾姆至慷鴮?dǎo)致ELBO小幅下降,但隨后歸并操作會拒絕較差的新增分量。
3.4.2 歸并機(jī)制
歸并操作有助于優(yōu)化全局?jǐn)?shù)據(jù)目標(biāo)。由于推理的代價隨分量數(shù)目K的增加而增加,因此為了保持較小的K,設(shè)計歸并操作時可以將兩個分量歸并為一個分量。
具體來說,通過隨機(jī)選擇分量ka和kb,比較現(xiàn)有模型q和候選模型q′的ELBO值。由于L(q′)的計算中除了Lentropy項,其他都是U統(tǒng)計量的線性函數(shù),因此可以提前對所有可能的歸并對計算Lentropya,b項,每個批次最多需要保存K(K-1)/2個標(biāo)量。這種預(yù)計算會加快后續(xù)的歸并操作,有助于改善模型質(zhì)量。通過在每輪數(shù)據(jù)訪問中執(zhí)行一次新增操作和若干次歸并操作,幾輪訪問后可以提升先驗(yàn)?zāi)P偷馁|(zhì)量,使其結(jié)構(gòu)更加緊湊。
(38)
(39)
在完成{ν,S}的更新后,固定全局變分因子,對所有局部變分變量參數(shù)進(jìn)行多輪更新,直至最大化變分目標(biāo)L′:
L′?Eq[lb(y,x,IS,.)-lbq(x,.)]
(40)
綜上,算法2給出圖像去噪算法。
算法2圖像去噪算法
輸入噪聲圖像y,通用先驗(yàn)?zāi)P?/p>
1.Repeat
7.Until ELBO收斂
8.Return
本文從BSDS數(shù)據(jù)集[28]中挑選的200張圖像中均勻采樣出200萬個圖像塊作為訓(xùn)練集,用來學(xué)習(xí)通用外部先驗(yàn)知識。實(shí)驗(yàn)測試用圖選取部分經(jīng)典合成圖像及BSDS圖像集中部分圖像,如圖2所示。
圖2 部分測試用圖Fig.2 Part of images for testing
利用最大似然估計來計算圖像塊均值m的參數(shù)τ和σ,類似地,使用計算塊的經(jīng)驗(yàn)協(xié)方差矩陣來估計精度矩陣S中的超參數(shù)n和V,并利用變分學(xué)習(xí)算法使模型能根據(jù)觀測數(shù)據(jù)進(jìn)行自適應(yīng)調(diào)整。
采用峰值信噪比(PSNR)作為算法的客觀度量,PSNR計算過程如下:
(41)
圖像塊大小在模型中扮演重要的角色。過大會導(dǎo)致邊緣模糊,丟失紋理細(xì)節(jié),過小將產(chǎn)生鋸齒效果。本文通過實(shí)驗(yàn)對比,選擇8×8的塊大小,能較好地適配不同的噪聲水平。
使用標(biāo)準(zhǔn)差分別為30、40、50、75的加性高斯噪聲對原始圖像進(jìn)行污染。將本文方法的去噪結(jié)果與BM3D[4]、LSSC[15]、EPLL[12]、PGPD[10]、NL-Bayes[8]、EPPGIC[11]等優(yōu)秀去噪算法進(jìn)行對比,結(jié)果如表1所示。表中數(shù)值為每個噪聲標(biāo)準(zhǔn)差下所有測試圖像的平均PSNR值,加粗?jǐn)?shù)據(jù)為最高值。
表1 去噪性能對比Table 1 Comparison of denoising performance dB
對于所選的經(jīng)典合成測試圖片,許多優(yōu)秀算法經(jīng)過調(diào)優(yōu)取得了出色結(jié)果。而本文方法在效果上優(yōu)于大部分算法,并且在更大的BSDS測試集合中取得了更好的結(jié)果,這顯示了本文模型對于大圖片集上非參學(xué)習(xí)的應(yīng)用價值。
EPLL可以看作是本文模型的簡化,其模型的分量數(shù)目是固定的。對于所有的噪聲級別和數(shù)據(jù)集,本文模型內(nèi)外部先驗(yàn)結(jié)合的方法均優(yōu)于EPLL,在性能上有所提升,表明貝葉斯非參學(xué)習(xí)方法能夠得到合適的模型聚類分量數(shù)目。
圖3為各方法對噪聲標(biāo)準(zhǔn)差為30時的Barbara圖像的去噪結(jié)果對比。本文方法在PSNR值上優(yōu)于LSSC、EPLL 、NL-Bayes、PGPD方法,相比BM3D和EPPGIC方法在PSNR值上略弱,但是在局部恢復(fù)細(xì)節(jié)方面有著突出的表現(xiàn),如圖3(e)所示,原始Barbara圖中腿部的兩處黑點(diǎn)區(qū)域,在上述兩種方法的復(fù)原圖中被過渡平滑了,而本文方法在恢復(fù)紋理的同時,上述兩處黑點(diǎn)仍保持清晰可見。EPLL方法盡管也復(fù)原了黑點(diǎn),但是腿部條紋恢復(fù)不太理想。
圖3 去噪視覺效果對比1Fig.3 Comparison 1 of visual effects after denoising
在更高噪聲級別下,各種方法的復(fù)原都出現(xiàn)了一些偽像。如圖4所示,BM3D、PDGD背景有大量的波紋,EPPGIC的背景有較多混雜的偽像,NL-Bayes方法在蛇身周邊有較多偽像,EPLL方法背景平滑度較弱,本文方法在背景有少量偽影,但無論在PSNR還是視覺觀感上都具有一定優(yōu)勢。圖4噪聲標(biāo)準(zhǔn)差為75。
圖4 去噪視覺效果對比2Fig.4 Comparison 2 of visual effects after denoising
本文構(gòu)建一種基于狄利克雷過程的可擴(kuò)展高斯混合模型用于圖像去噪。該模型從干凈圖像數(shù)據(jù)庫中學(xué)習(xí)通用先驗(yàn),借助聚類分量新增及歸并機(jī)制和模型中統(tǒng)計量的可累加性,使模型復(fù)雜性可以自適應(yīng)于訓(xùn)練觀測圖像的變化。實(shí)驗(yàn)結(jié)果表明,相比傳統(tǒng)去噪模型,該模型能取得更高的峰值信噪比,復(fù)原效果更佳。后續(xù)將從退化圖像中學(xué)習(xí)內(nèi)部先驗(yàn),以此捕獲圖像內(nèi)部自相似性,實(shí)現(xiàn)內(nèi)外部先驗(yàn)結(jié)合,使模型可用于圖像填充、去模糊與超分辨等任務(wù)。