基于小波框架的稀疏正則化方法及其在圖像復(fù)原中的應(yīng)用

袁存林 宋義壯

( 山東師范大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，250358，濟(jì)南 )

1 引 言

實(shí)際獲取的圖像總是受到模糊和噪聲的干擾，這種干擾降低了圖像的質(zhì)量. 因此，直接從成像設(shè)備中得到的圖像總是為退化圖像. 在一些諸如醫(yī)學(xué)成像的領(lǐng)域中，這種退化圖像會(huì)影響醫(yī)生的診斷，甚至可能會(huì)導(dǎo)致誤診而引起醫(yī)療事故.因此，從含有模糊和噪聲的退化圖像中復(fù)原原始圖像有著重要的應(yīng)用價(jià)值.

記u∈Rm1×m2為原始圖像，觀測(cè)的低質(zhì)量圖像f∈Rm1×m2由圖像u受模糊和噪聲的干擾而形成.在數(shù)學(xué)上，圖像復(fù)原可表示為以下問(wèn)題的反問(wèn)題：

f=Au+η,

(1)

其中η為高斯白噪聲，A為與模糊算子對(duì)應(yīng)的矩陣.

本文旨在從觀測(cè)圖像f重建高分辨率圖像u，該問(wèn)題可通過(guò)求解問(wèn)題(1)的最小二乘問(wèn)題

(2)

來(lái)求解.問(wèn)題(2)通過(guò)變分法可以轉(zhuǎn)化為

ATAu=ATf.

(3)

而矩陣ATA是病態(tài)的，其最小特征值在零附近，在受噪聲干擾的情況下，直接求解(3)式會(huì)放大噪聲.

(4)

其中，第一項(xiàng)為保真項(xiàng)，第二項(xiàng)為正則化項(xiàng)，權(quán)重參數(shù)λ>0.這是一個(gè)N-P難問(wèn)題.

使用l1或lq(0

(5)

小波框架變換W一般可表示為以下形式：

(6)

圖像u在小波框架下分解，得到系數(shù)Wαu，其中α∈Λ={(α1,α2)∈Z2|0≤α1,α2≤r}.當(dāng)α=0時(shí)，W0,0u為稠密系數(shù)，它刻畫了圖像的全局特征，即每個(gè)像素取為局部像素的加權(quán)平均；當(dāng)α≠0時(shí)，由于Wαu中大部分元素為零或集中在零附近，稱Wαu為稀疏系數(shù)，它刻畫了圖像的局部特征，即圖像的奇點(diǎn)或細(xì)節(jié)特征，比如邊緣和隱式邊緣.系數(shù)大部分為零或集中在零附近的現(xiàn)象稱為稀疏性.模型(5)中參數(shù)γ的選取需要反映系數(shù)Wu的稀疏性，在本文中，當(dāng)α=(0, 0)時(shí)，Wαu為刻畫圖像全局特征的稠密系數(shù)，取γ[α]=0;當(dāng)α≠(0, 0)時(shí)，Wαu為刻畫圖像局部特征的稀疏系數(shù)，取γ[α]=1.

下面將使用交替方向乘子法(Alternating Direction Method of Multipliers，ADMM)[7]求解模型(5).

2 模型求解

對(duì)于小波框架模型(5)，取d=Wu，則模型(5)可轉(zhuǎn)化為具有約束條件的極小化問(wèn)題：

(7)

問(wèn)題(7)的增廣拉格朗日函數(shù)為

(8)

其中y稱為拉格朗日乘子，ρ>0為懲罰參數(shù).對(duì)(8)式適當(dāng)變形有

(9)

極小化問(wèn)題(9)的ADMM型算法(算法1)如下：

(10)

(11)

yn+1=yn+ρ0(Wun-dn),

(12)

ρn+1=ρn·C.

(13)

其中，初始參數(shù)ρ0>0，常數(shù)C>1，y0,d-1為給定初值. 在上述算法1中，第1步與第2步相當(dāng)于交替更新變量u,d求Lρn(u,d,yn)的最小值解，第3步使用梯度上升法不斷更新拉格朗日乘子y，第4步使參數(shù)ρ趨于無(wú)窮以保證算法的穩(wěn)定性，這一點(diǎn)將在數(shù)值實(shí)驗(yàn)中體現(xiàn).

對(duì)(10)式右端求Fréchet導(dǎo)數(shù)，并令左端為零，可得

(14)

由d=Wu及(6)式中W的定義，可知dα=Wαu，其中α∈Λ.故而(11)式可寫為以下形式：

(15)

當(dāng)α=(0, 0)時(shí)，由于γ[α]=0，可得

(16)

當(dāng)α≠(0,0)時(shí)，γ[α]=1，對(duì)(15)式使用半閾值公式[5]求解可得

(17)

(18)

接下來(lái)將給出圖像恢復(fù)領(lǐng)域常用的基于小波框架的l1正則化模型[6]，用以衡量小波框架下l1/2正則化模型(5)對(duì)圖像的恢復(fù)效果.小波框架下l1正則化模型如下：

(19)

求解模型(19)的ADMM算法(算法2)如下：

(20)

(21)

yn+1=yn+ρ0(Wun-dn).

(22)

(20)式可通過(guò)求Fréchet導(dǎo)數(shù)求解，(21)式可使用軟閾值公式求解. 在接下來(lái)的數(shù)值實(shí)驗(yàn)中，將對(duì)比算法1和算法2對(duì)圖像的恢復(fù)效果，以驗(yàn)證小波框架下l1/2正則化方法的可行性.

3 數(shù)值實(shí)驗(yàn)

首先介紹一下數(shù)值實(shí)驗(yàn)中使用的小波框架和衡量圖像恢復(fù)效果的標(biāo)準(zhǔn). 小波框架使用B樣條緊小波框架[6]，在一維情況下，有濾波器

二維情況下的濾波器qα可通過(guò)一維情況下的濾波器做張量積得到qα[k1,k2]=qα1[k1]qα2[k2]. 小波框架變換Wα作用于圖像u相當(dāng)于濾波器qα與圖像u做卷積：Wαu=qα[-·]*u，*為循環(huán)卷積. 圖像的恢復(fù)效果通過(guò)峰值信噪比(Peak Signal to Noise Ratio，PSNR)衡量，峰值信噪比PSNR定義如下：

PSNR=10×log10((max(u))2/MSE),

接下來(lái)將探究l1/2、l1兩種正則化方法對(duì)圖像的恢復(fù)效果. 四幅測(cè)試圖像分別為Shepp-Logan、Cameraman、Lenna和Fingerprint(圖1)，圖像大小均為256×256，對(duì)其施加模糊核大小為7×7且標(biāo)準(zhǔn)差為4的高斯模糊以及40 dB的高斯白噪聲.

圖1 測(cè)試所用圖像

由表1可以看出，四幅圖像使用l1/2正則化方法所耗用的時(shí)間均低于使用l1正則化方法. 在恢復(fù)效果上，F(xiàn)ingerprint圖像使用l1正則化方法的恢復(fù)效果略好一些，Shepp-Logan、Cameraman和Lenna三幅圖像使用l1/2正則化方法的恢復(fù)效果較好，且Shepp-Logan圖像使用l1/2正則化方法有著較大的優(yōu)勢(shì), 這種現(xiàn)象可能與圖像的稀疏程度有關(guān). 文獻(xiàn)[8]使用圖像的梯度特征衡量其稀疏程度，考慮到小波框架下濾波器可以沿不同方向捕獲圖像的各類型局部特征，在本文中，使用稀疏度S(u)衡量圖像u的稀疏程度，S(u)定義如下：

表1 小波框架下l1/2、l1正則化方法的恢復(fù)效果(PSNR值)及耗用時(shí)間

其中m1×m2為圖像像素的總個(gè)數(shù)，指標(biāo)集K={(i,j)∈Z2|0≤i,j≤r,i2+j2≠0}，(r+1)2-1為集合K中元素的個(gè)數(shù)，Γ(Wαu)為Wαu中絕對(duì)值小于或等于1的元素總個(gè)數(shù)(圖像像素值在0到255之間).通過(guò)上述定義可知，稀疏度S(u)衡量了圖像u在變換域中系數(shù)的稀疏程度，S(u)越大，說(shuō)明變換域中的系數(shù)越稀疏.為衡量?jī)煞N正則化方法間的相對(duì)恢復(fù)效果，定義某圖像l1/2的離差為使用l1/2正則化方法恢復(fù)的PSNR值與使用兩種正則化方法得到的PSNR均值之間的差，類似可定義l1的離差. 離差越大，表示該正則化方法的相對(duì)恢復(fù)效果越好.

表2給出了各圖像的稀疏度及l(fā)1/2、l1的離差，并通過(guò)圖2展示了l1/2、l1的離差隨著圖像稀疏度的變化趨勢(shì). 通過(guò)表2和圖2可以發(fā)現(xiàn)，圖像的稀疏度和l1/2的離差呈現(xiàn)出正相關(guān)性，即隨著稀疏度的增加，l1/2正則化方法的恢復(fù)效果(PSNR值)逐漸超越l1的恢復(fù)效果，并且優(yōu)勢(shì)越來(lái)越大.l1的離差和圖像的稀疏度呈現(xiàn)出負(fù)相關(guān)性，對(duì)于稀疏度低的圖像，l1正則化方法具有一定的優(yōu)勢(shì). 從而可以根據(jù)圖像稀疏度這一先驗(yàn)信息，選擇合適的正則化方法進(jìn)行圖像恢復(fù)：對(duì)于稀疏度低的圖像，選擇l1正則化方法更合適；對(duì)于稀疏度高的圖像，選擇l1/2正則化方法更合適.

表2 各圖像的稀疏度及在不同正則化方法下的離差

圖2 稀疏度與離差關(guān)系圖

圖3給出了Shepp-Logan和Lenna圖像在兩種正則化方法下的恢復(fù)圖像. 為了進(jìn)一步對(duì)比小細(xì)節(jié)的恢復(fù)效果，本文對(duì)Shepp-Logan和Lenna圖像的局部區(qū)域進(jìn)行放大(圖3第2排和第4排). 可以看出，對(duì)于Shepp-Logan圖像的細(xì)節(jié)特征，l1/2正則化方法能更好地兼顧降噪和去模糊，局部圖案的去模糊效果更好一些，圖案以外區(qū)域降噪效果更明顯；對(duì)于Lenna圖像，l1/2正則化方法對(duì)其臉部輪廓的恢復(fù)效果更好. 鑒于l1/2正則化方法在恢復(fù)圖像細(xì)節(jié)特征方面的表現(xiàn)，對(duì)于小細(xì)節(jié)要求較高的圖像恢復(fù)工作，l1/2正則化方法是一個(gè)可行的選擇. 對(duì)于實(shí)驗(yàn)中使用的數(shù)據(jù)，C取為2，初值y0,d-1為零，求解l1/2、l1正則化模型的兩種算法的迭代次數(shù)分別為20與100. 不同圖像使用的參數(shù)λ,ρ0在表3中給出.

圖3 恢復(fù)圖像對(duì)比圖

表3 不同圖像使用的初始參數(shù)λ,ρ0

最后，為了驗(yàn)證求解l1/2正則化模型的算法1的穩(wěn)定性，定義誤差Ek+1=‖uk-u‖2/‖u‖2，uk由算法1在第k+1次迭代時(shí)生成，k≥0，u為真實(shí)圖像.圖4展示了誤差隨迭代次數(shù)的變化趨勢(shì)，可以看出，隨著迭代次數(shù)的增加，誤差不斷減小，并逐漸達(dá)到穩(wěn)定狀態(tài)，從而證實(shí)了算法1的穩(wěn)定性.

圖4 誤差隨迭代次數(shù)的變化趨勢(shì)

由實(shí)驗(yàn)可知，在小波框架下，用以求解l1/2正則化模型的算法是穩(wěn)定的.同時(shí)，l1/2和l1兩種正則化方法之間的相對(duì)恢復(fù)效果和圖像在小波變換域中系數(shù)的稀疏程度呈現(xiàn)出明顯的相關(guān)性：圖像在小波變換域中的系數(shù)越稀疏，l1/2正則化方法的相對(duì)恢復(fù)效果越好.對(duì)于圖像中受干擾的細(xì)節(jié)特征，l1/2正則化方法相比于l1正則化方法能夠更好地兼顧降噪和去模糊，從而更好地恢復(fù)圖像的細(xì)節(jié)特征.