練秋生 王 艷

(燕山大學信息科學與工程學院 秦皇島 066004)

1 引言

圖像壓縮感知(CS)是近幾年圖像處理領域的重大突破，它是由Donoho和Candes等人[1,2]在稀疏表示和優(yōu)化理論的基礎上提出的一種成像理論。其本質是直接將采樣與壓縮相結合，在對圖像進行隨機投影得到少量觀測值后，利用圖像表示的稀疏性先驗知識，從觀測值中恢復原始圖像的過程。圖像稀疏性是指其小波域只有少數(shù)與圖像邊緣紋理相對應的系數(shù)不為零或具有較大值，而與圖像光滑部分相對應的大多數(shù)系數(shù)為零或近似為零。CS理論利用圖像稀疏性先驗知識，所需觀測值遠遠小于Nyquist采樣定理所要求的數(shù)目。目前的CS重構算法包括最優(yōu)化算法和貪婪算法兩大類。最優(yōu)化算法將重構轉化為特定優(yōu)化問題，利用各類優(yōu)化算法對其求解，常用算法有凸優(yōu)化算法、非凸優(yōu)化算法及光滑?0范數(shù)(SL0)算法等。凸優(yōu)化算法用?1范數(shù)等價代替?0范數(shù)，將重構轉化為帶約束的凸優(yōu)化問題，然后求其最優(yōu)解。迭代硬閾值法(IHT)[3]，凸集交替投影法(POCS)[4]，迭代加權法[5]等都屬于此類算法。非凸優(yōu)化算法[6]將?0范數(shù)轉化為?p(p＜1)范數(shù)，將重構轉化為非凸優(yōu)化問題，然后求其最優(yōu)解，而光滑?0范數(shù)法[7]則通過構建一個光滑函數(shù)來逼近?0范數(shù)，然后對光滑函數(shù)進行優(yōu)化實現(xiàn)重構。貪婪算法在每次迭代過程中選取一個或幾個與觀測值余量最大相關的基，尋找一組與觀測值匹配的最稀疏基，從而實現(xiàn)信號的重構。匹配追蹤(MP)、正交匹配追蹤法(OMP)及它們的各種改進方法(如CoSaMP[8]，ROMP[9])均為此類算法。與傳統(tǒng)基于Nyquist采樣的重構相比，CS是一個重大突破。但標準CS重構只利用圖像稀疏性先驗知識，重構自由度大。實際上，圖像的小波系數(shù)除具有稀疏性外，還具有一些結構上的分布特征，如與邊緣紋理相對應的大系數(shù)往往呈樹狀分布，小波系數(shù)幅值在尺度間呈指數(shù)衰減等，將這些先驗知識引入的CS重構中可改善重構效果，減少重構所需的采樣值數(shù)目。Baraniuk在文獻[10]中提出了基于模型的CS重構理論，證明了將描述信號小波系數(shù)分布的結構模型引入到CS重構理論的正確性及有效性，為基于模型的CS圖像重構提供了理論依據(jù)。文獻[11]提出了TMP(Treebased Matching Pursuit)和TOMP(Tree-based Orthogonal Matching Pursuit)算法，將信號小波系數(shù)呈樹狀分布結構用于MP或OMP算法中。這兩種算法先利用MP或OMP算法尋找位于小波樹頂端的大系數(shù)，然后沿著有大系數(shù)的子樹向下搜索它的子系數(shù)，而對于沒有大系數(shù)的子樹不再進行搜索，與MP或OMP相比，此方法極大地縮小了搜索空間，減少了計算量，實現(xiàn)信號的快速重構。文獻[12]提出了BOMP算法(Block-based Orthogonal Matching Pursuit)，將信號小波系數(shù)的塊分布模型引入到OMP算法中。TMP,TOMP和BOMP算法相比與標準的CS重構有了較大改進，但上述算法引入的模型較為簡單，不能精確刻畫信號小波系數(shù)的結構分布特征，且上述算法對于一維信號的重構效果較好，但對于圖像重構，算法計算開銷大，重構效果不理想。圖像小波域的隱馬爾可夫模型(HMT)能較為精確地描述圖像小波系數(shù)的分布特征，它已廣泛應用于圖像去噪[13]和圖像復原[14]中。文獻[15]將該模型引入到一維信號的CS重構中，采用重置權重的?1范數(shù)最優(yōu)化算法，直接將HMT模型參數(shù)引入的重置權重的修正中。但文獻[15]未對算法的合理性做相應的理論推導，缺少理論依據(jù)。本文將HMT模型引入到圖像的CS的重構中，從最大后驗概率(MAP)出發(fā)，經(jīng)過理論推導，將基于HMT的圖像重構轉化為與標準圖像CS重構形式相似的優(yōu)化問題，并提出采用基于貝葉斯優(yōu)化的凸集交替投影法(Projections Onto Convex Sets，POCS)求解該優(yōu)化問題。另外本文將具有更多方向選擇性的雙樹小波域通用HMT(universal HMT，uHMT)模型代替小波域HMT模型，省去了使用HMT模型的計算開銷巨大的訓練過程。為進一步提高圖像重構質量，本文對uHMT模型進行修正，提出了改進的uHMT模型(improved uHMT, iuHMT)，它能更精確刻畫小波系數(shù)的局部分布特性，實驗結果表明了算法的有效性。

2 稀疏變換域的uHMT模型及其改進

2.1 小波域HMT模型

圖像的小波系數(shù)分布具有以下特征：(1)小波系數(shù)的非高斯性分布：即只有極少數(shù)小波系數(shù)具有較為顯著的較大值，而絕大多數(shù)小波系數(shù)的幅值較小，近似為零。(2)小波系數(shù)狀態(tài)的尺度持續(xù)性：小波系數(shù)“大”或“小”的狀態(tài)具有在尺度間傳遞的特性。圖像小波域HMT模型[14]能有效描述小波系數(shù)的上述特性。圖像小波系數(shù)呈樹狀圖分布，其HMT模型為每個小波系數(shù)wi設置對應的隱狀態(tài)變量Si。HMT模型采用混合高斯模型來模擬小波系數(shù)分布的非高斯性：如果已知第i個結點隱狀態(tài)變量的概率分布，則小波系數(shù)wi的概率密度與其他小波系數(shù)及隱狀態(tài)變量無關，即式中pSi(m)=p(Si=m)為隱狀態(tài)變量概率分布，它們滿足為高斯概率密度函數(shù)，m為小波系數(shù)隱狀態(tài)變量的一個實現(xiàn)，隱狀態(tài)變量在1,…,M中取值，本文中M=2。m=1的高斯成分具有較小的方差，刻畫與圖像光滑部分對應的小波系數(shù)分布特征； m=2的高斯成分具有較大方差，刻畫與圖像邊緣紋理對應的大系數(shù)的分布特征。

小波系數(shù)狀態(tài)沿尺度的持續(xù)性可采用隱狀態(tài)變量沿尺度的狀態(tài)轉移概率表征：任一結點小波系數(shù)的隱狀態(tài)僅依賴于其父結點ρ(i)的隱狀態(tài)，這種依賴關系由條件概率=p(m| S=r)來表Si| Sρ(i)ρ(i)征，它構成了此結點狀態(tài)轉移矩陣。上述參數(shù)、、μi,m、加上根結點處隱狀態(tài)變量概率分布pS1(m)構成了HMT模型參數(shù)，記為θ。文獻[16]假設同一級小波系數(shù)具有相同的模型參數(shù)，因此參數(shù)個數(shù)近似為4J(J為小波分解的級數(shù))，它們可以通過EM算法訓練得到。已知小波系數(shù)w和模型參數(shù)θ，可通過Upward-Downward算法獲得小波系數(shù)的后驗概率p(Si=m|w, θ)。在HMT模型下，由于小波變換的正交性，小波3個子帶相互獨立，小波系數(shù)的聯(lián)合概率分布可表示為

2.2 小波及雙樹小波域uHMT模型

上述HMT模型在使用中需經(jīng)過大量的訓練才能獲得具體的參數(shù)值，計算量較大。文獻[17]提出一種小波域通用HMT(universal HMT，uHMT)模型，它根據(jù)圖像小波系數(shù)的幅值沿尺度呈現(xiàn)指數(shù)衰減特性得到HMT模型的各級高斯成分方差有以下關系：式中j表示小波變換級數(shù)，j=1,2,…,J 。小波系數(shù)的“大”或“小”的狀態(tài)持續(xù)性沿尺度呈指數(shù)增強的特性，其狀態(tài)轉移矩陣可表示如下：

根據(jù)上述特性將HMT模型參數(shù)由原來的4J個減少為9個，即文獻[17]提出的小波域uHMT模型參數(shù)為

該模型參數(shù)對于大多數(shù)自然圖像都是適用的，因此用uHMT代替HMT能省去計算量較大的訓練步驟，有效減少了計算開銷。

普通的正交或雙正交2維小波變換的方向選擇性差(只有水平、垂直和對角線3個方向)，并且不具有平移不變性，它不能實現(xiàn)圖像幾何信息(如邊緣和紋理)的最優(yōu)表示。為提高圖像壓縮感知系統(tǒng)的性能，本文應用具有六方向選擇性和平移不變性的雙樹小波[18]來實現(xiàn)圖像稀疏表示。雙樹小波有實小波和復小波兩種，為減少冗余度和計算量，本文應用冗余度為2的雙樹實小波。雙樹小波的uHMT模型參數(shù)與小波uHMT模型參數(shù)類似，經(jīng)過對多幅圖像的訓練獲得雙樹小波域uHMT模型參數(shù)為

2.3 改進的uHMT(iuHMT)模型

上述uHMT模型將同級所有小波系數(shù)的uHMT模型的高斯成分簡化為具有相同的方差。然而由于小波系數(shù)的分布具有局部性，不同子帶不同位置的高斯成分方差并不相同，且成分方差與小波系數(shù)的局部能量相關?；谶@一特性本文提出了局部自適應性的改進uHMT模型(iuHMT,improved uHMT)，將高斯成分的方差與小波系數(shù)的局部能量相關聯(lián)。設w為圖像小波系數(shù)或雙樹小波系數(shù)，圖像的噪聲為2σ，則某點處小波系數(shù)wi的局部能量表示為：，其中Ni是以i點為中心的鄰域，M為鄰域大小，本文取M=3×3。wi兩個高斯成分的方差可分別表示為：則iuHMT模型的參數(shù)可改為：θ=。其中參數(shù)kS和kL分別表示兩個高斯成份的方差與該點小波系數(shù)局部能量的比例關系。通過對多幅圖像的統(tǒng)計結果，kL為0.1，kS分別為0.2和0.25時小波與雙樹小波域的iuHMT性能最好。本文中小波域的iuHMT模型的參數(shù)取值為

雙樹小波域的iuHMT模型的參數(shù)取值為

iuHMT的模型參數(shù)比uHMT少兩個，并且由于高斯成分的方差具有局部性和自適應性，它能更精確描述小波或雙樹小波系數(shù)的分布特征。

3 基于雙樹小波uHMT模型的壓縮感知圖像重構

3.1 圖像壓縮感知

設信號x的稀疏變換為x=Ψw，Ψ為信號的稀疏基，w為信號的稀疏系數(shù)。采用投影矩陣Φ對信號進行投影，得到信號的觀測值y=Φx。理論證明只要Θ=ΦΨ滿足限制等容性[1,2](Restricted Isometry Property, RIP)，利用觀測值y，即可重構信號x。信號x的重構可表示為求解以下優(yōu)化問題：

式中η為觀測噪聲。式(3)表示在滿足觀測值條件下，尋找x的非零變換系數(shù)最少的解，即x的最稀疏解。由于上述優(yōu)化問題是一個典型的NP-hard問題，通常采用?1范數(shù)等價代替?0范數(shù)[1,2]，即優(yōu)化問題式(3)等價為

求解上式的典型算法包括迭代硬閾值法(IHT)和迭代軟閾值法(IST)。

3.2 重構模型的推導

引入HMT模型后，根據(jù)最大后驗概率(MAP)準則，壓縮感知圖像重構可表示為

式中f(w| y)為小波系數(shù)的條件概率密度，y=ΦΨw +η為圖像x的觀測值，η為高斯白噪聲：η～N(μ, σ2)。根據(jù)貝葉斯公式式(5)等價為

式中f(w)和f(y)分別為w和y的概率密度函數(shù)，f(y)為常數(shù)，f(w)即為式(2)表示的HMT模型的聯(lián)合概率密度函數(shù)。由于η=y?ΦΨw 為高斯分布，因此有

對式(6)右端取自然對數(shù)得

將式(2)，式(7)代入式(8)并去掉常數(shù)項得

將式(9)的右邊乘?1有

根據(jù)拉格朗日乘子法，式(12)可以寫成下列形式：

式(13)與式(4)類似，均為在滿足觀測值條件下尋找滿足要求的最優(yōu)解，關鍵區(qū)別為式(4)是只利用了信號的稀疏性先驗知識，而式(13)則包含了HMT模型蘊含的稀疏性和樹結構分布先驗知識。

3.3 重構算法的實現(xiàn)

式(13)為凸優(yōu)化問題，本文采用基于貝葉斯優(yōu)化的凸集交替投影(POCS)算法求解，優(yōu)化問題涉及以下兩個凸集：

POCS算法交替向兩個凸集C1和C2進行投影， C1和C2的交點即為式(13)的最優(yōu)解：

POCS算法迭代步驟如下：

(1)將第t次迭代的結果wt向C1投影：

(2)利用梯度下降法將β向C2投影：

梯度下降法需在每次迭代中在下降方向上尋找最優(yōu)步長，即進行一次線性優(yōu)化搜索，計算開銷較大。本文采用貝葉斯優(yōu)化[17]代替梯度下降法。令λ=diag(λi)，其中λi為

根據(jù)式(15)得到t+1w的第i個分量為

將λi的值代入并化簡得

式(17)即為文獻[17]給出的HMT模型的貝葉斯優(yōu)化結果。貝葉斯優(yōu)化方法根據(jù)模型參數(shù)，對系數(shù)進行統(tǒng)計優(yōu)化 ，避免了每次迭代過程對步長的優(yōu)化搜索，減小了計算量。

基于雙樹小波iuHMT模型的壓縮感知重構算法具體步驟包括：

(1)初始化：設置最小誤差Emin和最大迭代次數(shù)tmax, w1=ΨTΦTy ，t=1。

(2)利用式(14)將wt投影到C1得到β。式中Φ表示PDCT(Permuted Discrete Cosine Transform)投影，Ψ和ΨT表示雙樹小波變換及其逆變換。

(3)利用穩(wěn)健中值算子估計β所含噪聲的方差σ2。

(4)利用雙樹小波域iuHMT模型參數(shù)及Upward-Downward算法估計β的隱狀態(tài)變量后驗概率pSi(m)。

(5)將β投影到C2，即利用式(17)獲得wt+1。

(6)xt+1=Ψwt+1，若||xt+1?xt||＜Emin或t＞tmax則轉移到下一步，否則令t=t+1轉到(2)繼續(xù)執(zhí)行。

(7)輸出重構圖像xt+1。

4 實驗結果

為驗證算法有效性，本文將基于uHMT模型和iuHMT模型的圖像壓縮感知重構與標準的圖像壓縮感知重構算法進行比較。本文實驗選取barbara,boat, couple, hill, lena, man, fingerprint 7幅512×512的標準灰度圖像進行重構。實驗中分別用DWTIHT, DWTuHMT, DWTiuHMT, TDWTIHT,DTDWTuHMT, DTDWTiuHMT表示小波域的IHT重構、小波域uHMT模型重構、小波域iuHMT模型重構、雙樹小波域IHT重構、雙樹小波域uHMT模型重構及雙樹小波域iuHMT模型的重構。其中小波域的uHMT及iuHMT重構算法將上部分算法實現(xiàn)中步驟(2)的雙樹小波變換改為Daubechies4正交小波變換，步驟(4)中的參數(shù)相應改為小波域的模型參數(shù)值。

表1為DWTuHMT、DWTiuHMT和DWTIHT算法在不同采樣率下重構的峰值信噪比(PSNR)比較。表中平均值一欄為各算法在各個采樣率下7幅圖像重構PSNR的平均值，黑體數(shù)字表示較高的PSNR。由表1可以看出在采樣率為20%，30%和40%時DWTuHMT算法重構的PSNR平均值要比DWTIHT算法的PSNR平均值分別高出1.39 dB,2.12 dB和2.50 dB，DWTuHMT算法重構效果明顯優(yōu)越于DWTIHT算法；在采樣率10%，DWTuHMT相比DWTIHT高出0.25 dB。在4種采樣率下，DWTiuHMT重構的PSNR平均值較DWTIHT分別高5.31 dB，5.34 dB，5.24 dB，4.93 dB，較DWTuHMT分別高5.06 dB，3.95 dB，3.12 dB，2.43 dB，充分體現(xiàn)了基于iuHMT模型重構的優(yōu)越性。表2 為DTDWTuHMT, DTDWTiuHMT和DTDWTIHT算法對7幅圖像重構PSNR平均值的比較。從表2可以看出在采樣率為20%、30%、40%時，DTDWTuHMT算法重構PSNR平均值比DTDWTIHT算法分別高出1.03 dB，1.51 dB，1.74 dB。在4種采樣率下，DTDWTiuHMT重構的PSNR平均值比DTDWTIHT分別高出了1.30 dB，2.46 dB，2.39 dB，2.20 dB之多，較DTDWTuHMT也分別高出了1.13 dB，1.43 dB，0.88 dB，0.46 dB。基于DTDWTiuHMT模型的重構圖像PSNR的平均值比DTDWTuHMT模型高0.97 dB，同樣體現(xiàn)了基于iuHMT模型重構算法的優(yōu)越性。將表1中DWTuHMT、DWTiuHMT算法對各圖像的重構PSNR與表2中DTDWTuHMT, DTDWTiuHMT算法對各圖像的重構PSNR進行比較可以看出，DTDWTuHMT在4種采樣率下重構PSNR平均值比DWTuHMT分別高出4.82 dB，4.17 dB，3.32 dB，2.47 dB；DTDWTiuHMT的重構PSNR平均值比DWTiuHMT分別高出了0.89 dB，1.65 dB，1.06 dB，1.50 dB，充分表明雙樹小波比普通小波更適合于壓縮感知。

表1 小波域uHMT、iuHMT重構與IHT重構PSNR的比較

表2 雙樹小波域的iuHMT、uHMT重構與IHT重構PSNR的比較

圖1為采樣率為20%時各種算法對barbara重構結果的局部圖像。圖1(a)為小波域IHT算法重構的圖像，該圖像紋理信息丟失嚴重，且含有大量噪聲；圖1(b)為小波域的uHMT模型重構圖像，與圖1(a)相比，該圖像恢復了更多的紋理信息，重構效果有明顯改善；圖1(c)為小波域的iuHMT模型重構圖像，與圖1(a)，1(b)相比，紋理信息進一步增多，且所含噪聲進一步降低；圖1(d)為雙樹小波域的IHT重構圖圖像，它恢復了較多的紋理信息，但圖像過度平滑，造成模糊失真；圖1 (e)為雙樹小波域uHMT算法的重構圖像，包含了更豐富的紋理信息，圖像更清晰，有效改善了圖1(d)的模糊失真；圖1(f)為雙樹小波域iuHMT算法的重構效果圖，與其余5幅圖像相比，它的紋理和邊緣信息最清晰，視覺效果最好。

5 結論

標準的壓縮感知重構只是利用了小波系數(shù)稀疏性的先驗知識，而未利用小波系數(shù)分布結構特征的先驗知識。本文提出一種基于雙樹小波域改進uHMT模型的圖像壓縮感知重構算法，將有效描述圖像小波系數(shù)分布的HMT模型引入的圖像的壓縮感知重構。從最大后驗概率出發(fā)，經(jīng)過理論推導，將基于HMT的重構轉化為型如標準壓縮感知重構的優(yōu)化問題，并提出采用基于貝葉斯優(yōu)化的POCS算法對優(yōu)化問題求解。本文還提出了基于雙樹小波域的uHMT模型，與小波相比，雙樹小波能更加有效地描述圖像的幾何信息，且采用uHMT模型省去了HMT模型所需的計算量巨大的訓練過程，提高了重構質量和重構速度。另外本文對uHMT進行了改進，提出參數(shù)個數(shù)更少，更能精確描述雙樹小波變換系數(shù)分布特征的iuHMT模型，進一步提高了圖像重構質量。

圖1 采樣率為20%時barbara重構效果比較

[1] Donoho D L. Compressed sensing [J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306.

[2] Candes E J nd Romberg J T. Robust uncertainty principles:Exact signal reconstruction from highly incomplete frequency information [J]. IEEE Transactions on Information Theory, 2006, 52(2): 489-509.

[3] Blumensath T and Davies M E. Iterative thresholding for sparse approximations[J]. Journal of Fourier Analysis and Applications, 2008, 14(5): 629-654.

[4] Bregman L M. The method of successive projection for finding a common point of convex sets [J]. Soviet Math, 1965,6(3): 688-692.

[5] Chartrand R and Yin Wotao. Iteratively reweighted algorithms for compressive sensing[C]. IEEE International Conference on Acoustics, Speech and Signal Processing, Las Vegas, NV, USA, 2008: 3869-3872.

[6] Chartrand R. Exact reconstruction of sparse signals via nonconvex minimization [J]. IEEE Signal Processing Letters,2007, 14(10): 707-710.

[7] Mohimani G H, Babaie-Zadeh M, and Jutten C. A fast sparse approach for overcomplete sparse decomposition based on smoothed ?0norm [J]. IEEE Transactions on Signal Processing, 2009, 57(1): 289-301.

[8] Needell D and Tropp J A. CoSaMP: Iiterative signal recovery from incomplete and inaccurate samples[J]. Applied and Computational Harmonic Analysis, 2008, 26(3): 301-321.

[9] Needell D and Vershynin R. Uniform uncertainty principle and signal recovery via regularized orthogonal matching pursuit[J]. Foundations of Computational Mathematics, 2009,9(3): 317-334.

[10] Baraniuk R G, Cevher Volkan, and Marco T D, et al..Model-Based compressive sensing [J]. IEEE Transactions on Information Theory, 2010, 56(4): 1982-2001.

[11] Duarte M F. Fast reconstruction from random incoherent projections[R]. Rice ECE Department Technical Report Tree,2005.

[12] Yonina C E and Helmut B. Block-sparsity:coherence and efficient recovery[C]. IEEE International Conference on Acoustics, Speech and Signal Processing, Taipei, Taiwan,2009: 2885-2888.

[13] Kivinen J J, Sudderth E B, and Jordan M I. Image denoising with nonparametric Hidden Markov trees[C]. IEEE International Conference on Image Processing, San Antonio,Texas, 2007, 3: 121-124.

[14] 趙書斌，彭思龍. 基于小波域HMT模型的圖像超分辨率重構[J]. 計算機輔助設計與圖形學學報，2003, 15(11): 1347-1352.Zhao Shu-bin and Peng Si-long. Wavelet-domain HMT-based image superresolution[J]. Journal of Computer-aided Design& Computer Graphics, 2003, 15(11): 1347-1352.

[15] Marco F D, Wakin M B, and Baraniuk R G. Wavelet-domain compressive signal reconstruction using a hidden markov tree model[C]. IEEE International Conference on Acoutics,Speech and Signal Processing, Las Vegas, NV, USA, 2008:5137-5140.

[16] Crouse Matthew S, Nowak R D, and Baraniuk R G. Waveletbased statistical signal processing using Hidden Markov Models[J]. IEEE Transactions on Signal Processing, 1996,3(6): 1029-1035.

[17] Romberg J K, Hyeokho C, and Baraniuk R G. Bayesian tree-structured image modeling using wavelet-domain hidden markov models [J]. IEEE Transactions on Image Processing,2001, 10(7): 1056-1068.

[18] Selesnick I W and Baraniuk R G. The dual-tree complex wavelet transform[J]. IEEE Signal Processing Magazine, 2005,22(6): 123-151.