袁存林 宋義壯

( 山東師范大學(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ù)，它刻畫(huà)了圖像的全局特征，即每個(gè)像素取為局部像素的加權(quán)平均；當(dāng)α≠0時(shí)，由于Wαu中大部分元素為零或集中在零附近，稱(chēng)Wαu為稀疏系數(shù)，它刻畫(huà)了圖像的局部特征，即圖像的奇點(diǎn)或細(xì)節(jié)特征，比如邊緣和隱式邊緣.系數(shù)大部分為零或集中在零附近的現(xiàn)象稱(chēng)為稀疏性.模型(5)中參數(shù)γ的選取需要反映系數(shù)Wu的稀疏性，在本文中，當(dāng)α=(0, 0)時(shí)，Wαu為刻畫(huà)圖像全局特征的稠密系數(shù)，取γ[α]=0;當(dāng)α≠(0, 0)時(shí)，Wαu為刻畫(huà)圖像局部特征的稀疏系數(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稱(chēng)為拉格朗日乘子，ρ>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)式可寫(xiě)為以下形式：

(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]使用圖像的梯度特征衡量其稀疏程度，考慮到小波框架下濾波器可以沿不同方向捕獲圖像的各類(lèi)型局部特征，在本文中，使用稀疏度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均值之間的差，類(lèi)似可定義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é)特征.