面向圖像先驗建模的可擴展高斯混合模型

張墨華,彭建華

(1.國家數字交換系統工程技術研究中心,鄭州 450002;2.河南財經政法大學 計算機與信息工程學院,鄭州 450002)

0 概述

圖像去噪在數學上屬于不適定問題,需要使用圖像先驗進行正則化,從而生成有意義的解[1]。自然圖像的一些自身屬性可以作為有用的圖像先驗用于圖像去噪,如稀疏性、多尺度自相似與非局部自相似性等。文獻[2-4]采用圖像的稀疏性作為先驗信息,應用小波與曲波等稀疏變換,或是通過學習圖像塊得到完備字典以及圖像塊的稀疏表示。文獻[5]采用圖像像素在空間域和強度域中表現出的相似性作為先驗信息,通過相似塊進行協同處理。文獻[6]利用同一幅圖像中相同尺度或不同尺度的相似子塊進行去噪。文獻[7-8]將圖像的非局部自相似性作為先驗信息,對相隔較遠的相似塊進行非局部處理。文獻[9-11]則是使用多種圖像先驗組合的方法進行圖像去噪。

除了使用預定義的圖像先驗,文獻[12]還提出了從自然圖像中學習先驗的方法。生成式圖像先驗學習方法通常從一組干凈的外部圖像中學習先驗模型,并將其應用于給定的退化圖像[13-14],或者從給定的退化圖像中學習先驗[15]。近年來,判別式圖像先驗學習方法也開始得到廣泛應用。該方法根據清晰-噪聲圖像對學習去噪模型,可進一步細分為基于深度學習的方法[16]、基于隨機場的方法[17]與基于反應擴散的方法[18]。

在上述方法中,文獻[12]提出的使用高斯混合模型(Gaussian Mixture Model,GMM)先驗進行自然圖像塊建模的EPLL方法取得了較大成功。這一方法的主要思想是從干凈圖像塊中學習先驗,并將其用于求解圖像復原問題。但EPLL方法的混合分量數目固定(文獻[12]中為200),學習訓練完畢后分量不易進行擴展。針對這個問題,文獻[19]提出非參貝葉斯方法,為結構化數據(如文本文檔、時間序列和圖像)的無監督建模提供了靈活框架,其模型的復雜度可以根據數據規模增減進行自適應變化。狄利克雷過程(Dirichlet Process,DP)是最流行的非參貝葉斯方法,其采樣結果是一個分布,而不是變量,因此DP混合模型可以看作無窮多個混合分量的混合模型[20]。作為分層貝葉斯框架下的非參先驗,DP可以根據具體參數值對該模型生成的數據進行劃分,具有自然解釋特性,分量的數量是隨機的,可以隨著觀測到的新數據增長。但DP的學習是一個復雜的最優化問題,難以精確求解,只能進行近似求解,常用的方法之一是變分推理。變分推理中基于傳統均值場的坐標上升算法按照先后順序對數據集中所有的局部變分參數與全局變分參數進行更新,因此在大數據集上效率較低,同時由于更新緩慢,容易陷入局部最優解。

針對上述問題，本文構建一種基于狄利克雷過程的可擴展高斯混合模型，從清晰圖像數據庫中學習外部通用先驗，使模型復雜度根據數據規模自適應變化。同時，為改善模型的推理過程，提出一種基于批次更新方式的可擴展變分算法，將數據集分為若干批次，在每個批次訪問中對圖像的局部變分參數與全局變分參數進行更新，從而求解圖像去噪中所有隱變量的變分后驗分布，實現通用先驗學習。

1 理論定義

定義1(狄利克雷過程) 令G0是集合Θ上的分布,α是正值實數。對于Θ上的有限劃分A1,A2,…,Ar,則向量G(A1,A2,…,Ar)是隨機的。如果有:

G(A1,A2,…,Ar)～Ddir(αG0(A1),αG0(A2),…,αG0(Ar)

(1)

則稱G是基分布G0和集中參數α上的狄利克雷過程(DP)所生成的分布,記作:G～DDP(α,G0)。

假定從G中獨立生成N個隨機變量ηn,則有:

ηn～G,n∈{1,2,…,N}

(2)

聯合分布{η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將更多質量集中于均值附近。α值越小,G就更為離散。證明可以參見文獻[22]。

定義2(狄利克雷過程混合模型) 如果ηn是第n個觀測變量分布F的參數,即有F(x|ηn),則DPMM可以看作為無窮多個混合分量的混合模型,即有:

G|{α,G0}～DDP(α,G0)

ηn|G～G

xn|ηn～F(x|ηn)

(3)

在DP混合模型中,DP作為分層貝葉斯框架下的非參先驗,可以根據具體參數值對該模型生成的數據進行劃分,因此,DP混合模型擁有自然解釋特性,其中分量的數量(劃分中原子的數目)是隨機的,并且能夠隨著觀測到的新數據進行增長。

文獻[23]根據折棍子構造的過程,對DP給出了更為清晰的特征表述。考慮2個獨立隨機變量無限集合,vi～Bbeta(1,α),ηi～G0,i={1,2,…},它們的折棍子過程表示如下:

(4)

由式(4)可見,G是離散的。G的支持集包括可數無限原子集合,這些原子集合從G0獨立生成。混合比例βi(v)通過將單元長度的“棍子”連續折斷成無限數目的分段生成,每分段的大小與棍子剩余部分成正比,從Bbeta(1,α)分布獨立生成。

在DP混合中,向量βi(v)是一個包含比例的無限向量,{η1,η2,…}是表示混合分量的原子。令zn表示數據點xn所關聯的混合分量的分派變量。數據生成過程如下:

1)生成vi|α～Bbeta(1,α),i={1,2,…}。

2)生成ηi|G0～G0,i={1,2,…}。

3)對于第n個數據點:

(1)生成zn|{v1,v2,…}～Mmult(β(v))。

(2)生成xn|zn,ηn～F(xn|ηzn)。

在圖像去噪任務中,每個圖像具有其特定的聚類頻度wi={wi1,wi2,…,wiK,…},通過施加分層DP先驗,聚類頻度可以看作一個有限的狄利克雷分布:

[wi1,wi2,…,wiK,wi>K]～

Ddir(αβ1,αβ2,…,αβK,αβ>K)

(5)

從式(5)可知,wi的均值為β,方差由集中參數α所決定,下標>K表示除了前K個分量外其余分量累加之和。本文約定觀測數據從指數家族分布生成,DP的基分布是對應的共軛先驗。

2 圖像生成模型

使用基于DP的高斯混合模型對圖像的生成過程進行建模。自然圖像可以看作圖像塊組合,每個塊包含D個像素,并且有其對齊基點b,b的取值在D中,即:b～Ccat(1/D,…,1/D)。一幅圖像由所有對齊基點下各個塊重疊生成,每個對齊基點下,可能有部分塊處于部分觀測的狀態,即部分參與最終圖像的生成,部分不參與。

參照文獻[12]的方法,將每個圖像塊去平均化,對齊基點下每個圖像塊gibn看成由高斯混合分布生成,即每個塊表示為均值為0、精度矩陣為S的高斯混合,如下式所示:

(6)

其中,i表示第i幅圖像(共N幅),b表示對齊基點(共D個),n表示圖像i在對齊基點b中的第n個塊(共Nib個)。精度矩陣S的先驗為Wishart分布,滿足:

S～Wwis(μ,V)

(7)

圖像塊均值mibn滿足高斯分布:

mibn～Nnorm(τ,σ)

(8)

每個圖像塊可以指派為高斯混合模型中某個分量(聚類)中的一個,用zibn滿足分類分布,表示為:

zibn～Ccat(w1,w2,…,wk)

(9)

給定通過對齊基點b生成的、塊均值為mibn的塊gibn,圖像xi的采樣可以通過式(10)生成:

(10)

在圖像去噪任務中,觀測圖像yi是干凈圖像xi的退化版本,是由向xi中加入高斯白噪聲后所得。其生成定義如下:

yi～Nnorm(yi|xi,σ2I)

(11)

其中,σ2為噪聲方差。

整個圖像生成過程如圖1(a)所示,其變分參數的推理過程如圖1(b)所示,第3節將給出詳細過程。

圖1 圖像生成有向圖模型及推理過程Fig.1 Directed graph model of image generation andits inference process

3 圖像先驗學習

3.1 變分分布

本文通過變分分布的參數推理實現圖像先驗模型的學習。給定干凈觀測圖像xi,需要學習模型中的各種隱變量,如全局圖像變量與局部圖像變量。全局圖像變量包括折棍子過程的折斷比例vk,與每個聚類的精度矩陣Sk;局部圖像變量包括聚類頻度wi、圖像塊gibn、塊均值mibn與塊聚類指派zibn。直接計算后驗分布難度較大,因此常用變分推理進行近似求解。變分推理是一種近似似然和后驗的確定性方法[24],采用優化方法嘗試尋找一個合適變分分布q,它與真實的后驗在KL散度上無限接近。

定義待學習的隱變量集合W={g,S,wi,mi,zi},觀測圖像變量x={x1,x2,…,xN},超參數θ={α,μ,V}。隱變量的后驗分布與其近似分布q的KL散度為:

D(q(W)‖p(W|x,θ))=Eq[lbq(W)]-

Eq[lbp(W,x|θ)]+lbp(x|θ)

(12)

式(12)中KL散度的最小化求解可以轉換為對數邊緣似然下界(ELBO)的最大化求解[25]。如下所示:

lbp(x|θ)≥Eq[lbp(W,x|θ)]-

Eq[lbq(W)]?L

(13)

將L定義為待優化求解的目標函數,通過參數的學習過程,盡量使L值最大化,從而盡可能減小變分近似后驗q與真實后驗的差值。采用均值場方法[25],假定后驗完全獨立,q可以分解為指數家族密度函數的乘積。

(14)

其中,隨機度量的隨機性通過無限集合{v1,v2,…}和{S1,S2,…}體現,而隨機性的生成過程則通過經典的截斷折棍子過程表示[21]。假設固定截斷位置為T,則混合比例wt(v)在t>T時等于0。T是變分參數,可以自由設置,并非先驗模型的一部分。

下面具體給出所有變量的變分后驗分布,這些分布參數稱為變分參數,是要進行優化求解的對象。

在局部變量因子中,圖像塊g、圖像塊均值m、塊分量指派z、圖像聚類頻度w、圖像x的變分后驗分布分別為:

(15)

在全局變量因子中,折斷比例變量v的變分分布為:

(16)

另一個全局變量因子精度矩陣S的變分分布為:

(17)

對于自然圖像,所有的對齊基點具有相似的值,設計對齊基點B的變分后驗為

q(B)=Ccat(1/D,…,1/D)

(18)

該設計簡化了更新操作,同時仍然可以避免不重疊塊所造成的偽影情況。

3.2 變分參數的更新

圖像通用先驗的學習是基于外部干凈圖像數據庫進行的,此時圖像x是可觀測變量。根據式(12),結合3.1節所述各變分分布,目標函數L可以分解為3個部分:

L?Limage+Lentropy+Lweight

其中:

Limage?Eq[lbp(x|g,m,b)+

(19)

(20)

(21)

(22)

(23)

其中,期望Eq[β]通過式(4)進行計算,Cik是對于圖像i中所有圖像塊使用聚類k的情況的統計量,定義為:

(24)

(25)

(26)

從上式中抽取出兩個統計量Ck和Fk,定義為:

(27)

(28)

這兩個統計量是訓練圖像集中所有分派給聚類k的數據統計量,將在3.3節中詳述。

(29)

(30)

(31)

(32)

其中,Oibn表示在對齊基點b下第n個塊中可觀測的像素的數目。

(33)

αEq[βk]Tk

(34)

其中,Tk定義為:

(35)

(36)

3.3 通用先驗學習算法

傳統的坐標上升算法先對數據集中所有局部變分參數進行更新,再對全局變分參數進行更新,因此在大數據集上效率較低。同時,由于全局參數在第一次參與局部參數更新時是初始化值,計算意義不大。此外,傳統算法更新較為緩慢,容易陷入局部最優解。針對這些問題,本文利用式(27)、式(28)、式(35)中所得3個統計量的可累加性,采用批次更新方式完成坐標上升算法。將數據集分成若干批次,對所有批次完成訪問稱為一輪。在每輪批次訪問中,首先對批次中每個圖像的局部變分參數進行更新,然后進行全局變分參數的更新。

在每個批次中都對式(27)、式(28)、式(35)所得的3個統計量的值進行記錄(表示為Ubt),然后累積到全局統計量Uo中,其定義如下:

其中,B表示總批次數。在訪問完每個批次bt之后,執行增量更新來累積統計量,從而反映新批次的匯總情況,并刪除該批次bt之前的值。

(37)

其中,ISbt是圖像數據集IS的批次數據集。綜上,通用先驗訓練過程如算法1所示。

算法1通用先驗訓練的變分推理算法

輸入圖像數據集IS

1.Repeat

2.For each batch ISbtin IS do

12.End for

13.Until ELBO收斂

3.4 聚類分量的新增與歸并機制

3.4.1 新增機制

在通用先驗訓練過程中,通過分量新增機制可以增加新的有用分量,幫助規避局部最優解。但即使這些分量得到整個圖像數據庫的支持,每個批次也無法提供缺失分量所需要的足夠樣本,因此,整個新增過程需要兩輪數據訪問。第一輪采集目標數據樣本,第二輪生成新的分量,并且使用擴展的模型對每個批次進行更新。具體步驟如下:

2)在第2輪數據訪問之前生成新的分量。通過運行有限次數的變分推理,將數據集IS′擬合為包括K′(設置為10)個混合分量的模型。從而將整個模型擴展為K+K′個分量。期間并沒有對這些新分量所產生的ELBO的變化進行評估,依賴隨后的歸并操作來移除不需要的分量。

3)在第2輪中訪問每個批次,并且對擴展后的K+K′個混合模型執行局部和全局參數更新。通過對目標數據集X′分析得到統計量U′,并入全局統計量Uo中。數據集X′有兩種分派結果:分派給最初的分量(最有可能是k′),或者分派給全新的分量。第2輪結束后,從Uo中去掉U′,使得Uo和全局參數與數據集IS一致。

一次新增多個分量有助于規避局部最優解。盡管新增分量可能會因為增加不必要的分量而導致ELBO小幅下降,但隨后歸并操作會拒絕較差的新增分量。

3.4.2 歸并機制

歸并操作有助于優化全局數據目標。由于推理的代價隨分量數目K的增加而增加,因此為了保持較小的K,設計歸并操作時可以將兩個分量歸并為一個分量。

具體來說,通過隨機選擇分量ka和kb,比較現有模型q和候選模型q′的ELBO值。由于L(q′)的計算中除了Lentropy項,其他都是U統計量的線性函數,因此可以提前對所有可能的歸并對計算Lentropya,b項,每個批次最多需要保存K(K-1)/2個標量。這種預計算會加快后續的歸并操作,有助于改善模型質量。通過在每輪數據訪問中執行一次新增操作和若干次歸并操作,幾輪訪問后可以提升先驗模型的質量,使其結構更加緊湊。

3.5 圖像去噪算法

(38)

(39)

在完成{ν,S}的更新后,固定全局變分因子,對所有局部變分變量參數進行多輪更新,直至最大化變分目標L′:

L′?Eq[lb(y,x,IS,.)-lbq(x,.)]

(40)

綜上,算法2給出圖像去噪算法。

算法2圖像去噪算法

輸入噪聲圖像y,通用先驗模型

1.Repeat

7.Until ELBO收斂

8.Return

4 實驗結果與分析

本文從BSDS數據集[28]中挑選的200張圖像中均勻采樣出200萬個圖像塊作為訓練集,用來學習通用外部先驗知識。實驗測試用圖選取部分經典合成圖像及BSDS圖像集中部分圖像,如圖2所示。

圖2 部分測試用圖Fig.2 Part of images for testing

利用最大似然估計來計算圖像塊均值m的參數τ和σ,類似地,使用計算塊的經驗協方差矩陣來估計精度矩陣S中的超參數n和V,并利用變分學習算法使模型能根據觀測數據進行自適應調整。

采用峰值信噪比(PSNR)作為算法的客觀度量,PSNR計算過程如下:

(41)

圖像塊大小在模型中扮演重要的角色。過大會導致邊緣模糊,丟失紋理細節,過小將產生鋸齒效果。本文通過實驗對比,選擇8×8的塊大小,能較好地適配不同的噪聲水平。

使用標準差分別為30、40、50、75的加性高斯噪聲對原始圖像進行污染。將本文方法的去噪結果與BM3D[4]、LSSC[15]、EPLL[12]、PGPD[10]、NL-Bayes[8]、EPPGIC[11]等優秀去噪算法進行對比,結果如表1所示。表中數值為每個噪聲標準差下所有測試圖像的平均PSNR值,加粗數據為最高值。

表1 去噪性能對比Table 1 Comparison of denoising performance dB

對于所選的經典合成測試圖片,許多優秀算法經過調優取得了出色結果。而本文方法在效果上優于大部分算法,并且在更大的BSDS測試集合中取得了更好的結果,這顯示了本文模型對于大圖片集上非參學習的應用價值。

EPLL可以看作是本文模型的簡化,其模型的分量數目是固定的。對于所有的噪聲級別和數據集,本文模型內外部先驗結合的方法均優于EPLL,在性能上有所提升,表明貝葉斯非參學習方法能夠得到合適的模型聚類分量數目。

圖3為各方法對噪聲標準差為30時的Barbara圖像的去噪結果對比。本文方法在PSNR值上優于LSSC、EPLL 、NL-Bayes、PGPD方法,相比BM3D和EPPGIC方法在PSNR值上略弱,但是在局部恢復細節方面有著突出的表現,如圖3(e)所示,原始Barbara圖中腿部的兩處黑點區域,在上述兩種方法的復原圖中被過渡平滑了,而本文方法在恢復紋理的同時,上述兩處黑點仍保持清晰可見。EPLL方法盡管也復原了黑點,但是腿部條紋恢復不太理想。

圖3 去噪視覺效果對比1Fig.3 Comparison 1 of visual effects after denoising

在更高噪聲級別下,各種方法的復原都出現了一些偽像。如圖4所示,BM3D、PDGD背景有大量的波紋,EPPGIC的背景有較多混雜的偽像,NL-Bayes方法在蛇身周邊有較多偽像,EPLL方法背景平滑度較弱,本文方法在背景有少量偽影,但無論在PSNR還是視覺觀感上都具有一定優勢。圖4噪聲標準差為75。

圖4 去噪視覺效果對比2Fig.4 Comparison 2 of visual effects after denoising

5 結束語

本文構建一種基于狄利克雷過程的可擴展高斯混合模型用于圖像去噪。該模型從干凈圖像數據庫中學習通用先驗,借助聚類分量新增及歸并機制和模型中統計量的可累加性,使模型復雜性可以自適應于訓練觀測圖像的變化。實驗結果表明,相比傳統去噪模型,該模型能取得更高的峰值信噪比,復原效果更佳。后續將從退化圖像中學習內部先驗,以此捕獲圖像內部自相似性,實現內外部先驗結合,使模型可用于圖像填充、去模糊與超分辨等任務。