郝 巖 馮象初 許建樓②

①(西安電子科技大學(xué)理學(xué)院 西安 710071)

②(河南科技大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院 洛陽 471003)

1 引言

圖像去噪是圖像處理中最基本的問題之一，其目的就是從所得到的帶噪圖像中恢復(fù)出原圖像，即求原圖像的某種最優(yōu)意義下的估計值。通常，一個較好的去噪方法應(yīng)該是在去除噪聲的同時又能較好地保持圖像的原有信息。傳統(tǒng)的圖像去噪方法，如中值濾波、高斯濾波等，主要是將圖像的高頻成分濾除， 因此，這些方法不可避免地會將圖像的一些細(xì)節(jié)特征去掉，從而使得所恢復(fù)出來的圖像在紋理區(qū)域看起來比較模糊。

近年來，變分、偏微分方程方法在數(shù)字圖像處理中得到了廣泛應(yīng)用[1-5]。其中圖像去噪的各項異性擴(kuò)散模型(Perona Malik,PM)[1]以及全變差模型(Rudin Osher Fatemi,ROF)[2]是比較典型的代表。目前，ROF模型因其具有良好的保邊性能而備受青睞，但是該模型在圖像恢復(fù)過程中常常會出現(xiàn)階梯效應(yīng)，因此，繼ROF模型之后，許多方法被提出去處理這個問題。特別地，Lysaker等人[6]提出用高階偏微分方程(Lysaker Lundervold Tai,LLT模型)來克服階梯效應(yīng)，但高階偏微分方程演化圖像的邊緣比ROF模型快， 所以經(jīng)常在圖像的邊緣引起模糊現(xiàn)象。為了在去除噪聲的同時既保持圖像的邊緣，又避免階梯現(xiàn)象，Lysaker等人[7]又提出一個兩步去噪模型(Lysaker Osher Tai,LOT模型)。實驗表明LOT模型能夠彌補(bǔ)ROF模型及LLT模型的不足，然而，該模型的收斂速度比較慢。因此，對LOT模型給出有效的加速算法顯得十分必要。最近，一些改進(jìn)的 LOT模型及算法被相繼提出[8,9]。其中，文獻(xiàn)[9]提出在 LOT模型中用一個輔助變量θ代替法向量n，然后再利用分裂 Bregman方法[10-14]對其兩步分別進(jìn)行數(shù)值求解。這樣在一定程度上降低了原有算法的計算復(fù)雜度，從而實現(xiàn)加速收斂。

本文在研究LOT模型的基礎(chǔ)上，提出了一種新的變分去噪模型。新模型由于耦合了兩個變量，因此，在求解時，先利用交替極小化方法化原模型為兩個簡單的子模型，然后再利用分裂Bregman方法對其進(jìn)行交替求解。由于兩子模型對LOT模型中的兩步從構(gòu)成方式及計算方式上分別進(jìn)行了改進(jìn)，因此，實驗表明，新方法不僅具有較快的收斂速度，而且在去噪過程中其去噪能力及保邊性能均好于LOT模型。

2 LOT模型

從數(shù)學(xué)的角度來講，圖像去噪問題是一個不適定的線性逆問題。其數(shù)學(xué)模型為

其中u0是觀察到的帶噪圖像，u是原圖像，η是加性噪聲。對于該不適定問題，一般可通過正則化方法來求解。關(guān)于式(1)，文獻(xiàn)[7]建議用一個兩步去噪模型(LOT模型)進(jìn)行求解。其基本思想是：首先構(gòu)造出圖像等照度線的法向場，然后再根據(jù)該向量場重構(gòu)圖像。因此， LOT模型的第1步就是解下面優(yōu)化問題

其中第1項為法向量逼近項，第2項為圖像逼近項，β＞0 為懲罰參數(shù)。

LOT模型中由于使用低階的偏微分方程，因此，它能夠避免求解高階偏微分方程時所引起的邊緣模糊缺陷。實驗表明，LOT模型在去噪的同時不僅能較好地保持圖像的邊緣，而且還能減緩階梯效應(yīng)。然而它的收斂速度比較慢。

針對上述缺陷，最近，文獻(xiàn)[9]在對LOT模型改進(jìn)的基礎(chǔ)上提出了一種加速方法。通過用一個輔助變量θ，文獻(xiàn)[9]將LOT模型等價地轉(zhuǎn)換成下面模型

式(4)，式(5)與LOT模型中式(2)，式(3)相比，不但可以減少一個未知量的計算，而且還可以利用一些快速算法進(jìn)行求解。文獻(xiàn)[9]通過利用分裂Bregman方法實現(xiàn)了加速收斂。然而，在去噪效果上，文獻(xiàn)[9]中方法并沒有太大改進(jìn)。

3 本文模型及其算法

3.1 模型建立

考慮到交替迭代方法的諸多優(yōu)點，本文對式(1)提出一個同時包含變量θ和u的耦合去噪模型

式中前2項是正則項，后3項是保真項，λ,β1,β2和μ是正則參數(shù)。

由于式(6)中耦合了兩個變量，因此首先采用交替迭代策略將其轉(zhuǎn)化為下面兩個簡單的子模型：

顯然，通過式(7)可以構(gòu)建角θ，而通過式(8)能夠重構(gòu)與θ角相適應(yīng)的圖像u。與式(4)，式(5)相比，式(7)，式(8)主要有以下兩個不同點：

第一，模型的構(gòu)成方式不同。式(7)在構(gòu)建角θ時，保真項中不但要求θ逼近θ0，而且還要求(c osθ,sinθ)逼近法向量因此，式(7)比式(4)多出一個L2法向量逼近項，并且該逼近項還充分用到了上層重構(gòu)出的圖像nu。隨著迭代的交替進(jìn)行，包含nu的這個逼近項不斷地對θ角進(jìn)行校正，從而使其得到優(yōu)化。然而式(4)在構(gòu)建θ角時，保真項中僅僅要求θ逼近θ0，并沒有利用到圖像u。此外，在第2步重構(gòu)圖像時，式(8)比式(5)不僅多了一個全變差(Total Variation,TV)正則項而且還將原來的法向量逼近項換成了現(xiàn)在的L2法向量逼近項，并且此逼近項還能隨θn+1不斷地進(jìn)行更新，最終使得重構(gòu)圖的法向量場逼近較優(yōu)的向量場(cos?θ,sin ?θ)，這一點式(5)是無法比擬的。

第二，模型的計算方式不同。式(4)，式(5)在計算方式上是先計算式(4)，然后再根據(jù)式(4)結(jié)果計算式(5)，當(dāng)達(dá)到停止標(biāo)準(zhǔn)時便結(jié)束迭代，并輸出圖像u*，在計算上這兩步之間并不相互制約。然而,本文的式(7)，式(8)在計算上卻相互制約，相互影響，交替計算。隨著交替迭代的進(jìn)行，重構(gòu)出的圖像u能夠促使角θ進(jìn)一步得到優(yōu)化，反過來，優(yōu)化后的θ角也使式(8)中的重構(gòu)圖像u向更優(yōu)的方向演化，直到獲得一個較好的圖像為止。

3.2 算法提出

分裂Bregman方法是由文獻(xiàn)[10]提出的一種高效的迭代方法，常用于求解帶有L1項的優(yōu)化問題。由于該方法具有編程簡單，數(shù)值求解過程比較穩(wěn)定，在計算過程中保持正則化參數(shù)為一個常數(shù)，占有內(nèi)存小且具有較快的計算速度和收斂速度等優(yōu)勢，因此，這種方法已經(jīng)被廣泛地運(yùn)用到圖像處理中。本小節(jié)將采用該方法對式(7)，式(8)進(jìn)行數(shù)值求解。

首先，在第 1個子模型式(7)中引入輔助變量ρx,ρy來處理不可微項然后再引入變量bx,by來更新迭代過程。最終，求解式(7)的迭代公式如下：

同樣，在第 2個子模型式(8)中引入輔助變量ηx,ηy,cx,cy，用類似于上面的方法，可以推得式(8)的迭代公式為

其中

λ2為懲罰參數(shù)。結(jié)合式(9)和式(10)，給出本文算法。

新模型的分裂Bregman算法：

步驟 4 停止迭代：對ε＞0，當(dāng)＜ε時，則停止迭代，并輸出復(fù)原圖像否則，令，轉(zhuǎn)步驟2。

注：對于式(9)，式(10)中的第1個式子，本文都是運(yùn)用Gauss-Seidel迭代求得。其具體的迭代公式分別為

此外，類似于文獻(xiàn)[13]中的結(jié)果，新算法是收斂的。

4 數(shù)值實驗

為了驗證新算法的有效性，本節(jié)分別對兩幅大小為256×256的帶噪圖像進(jìn)行實驗，并將新算法分別和文獻(xiàn)[7,9]中方法進(jìn)行比較。本文將從視覺效果與定量指標(biāo)兩個方面對圖像質(zhì)量進(jìn)行比較。采用的定量指標(biāo)有峰值信噪比(Peak Signal to Noise Ratio,PSNR)，均方誤差(Mean Squared Error,MSE)和結(jié)構(gòu)相似性(Structural SIMilarity,SSIM)。其定義分別如下：

其中P為原始無噪圖像，P'為恢復(fù)后的圖像，M和N為圖像尺寸大小，i,j為圖像像素下標(biāo)，σP',σP,μP',μP分別為P',P的標(biāo)準(zhǔn)差和均值，σP'P為P'和P之間的協(xié)方差。SSIM可以用來測量兩圖像之間的結(jié)構(gòu)相似性，SSIM值越大，說明兩圖像的結(jié)構(gòu)越相似。

圖1 Lena 圖實驗結(jié)果

實驗 1 選取細(xì)節(jié)豐富的結(jié)構(gòu)圖“Lena”作為測試圖像。實驗開始先對“Lena”圖加入均值為0、方差為10的高斯白噪聲。試驗發(fā)現(xiàn)，新算法對參數(shù)β1不是很敏感，對參數(shù)λ,β2和μ較為敏感。大量的數(shù)值仿真結(jié)果表明，參數(shù)β1在0～1之間取值時效果較優(yōu)。本實驗中取參數(shù)λ=0 .5,β1=0 .1,β2=0.1,μ=0 .001。此外，為便于比較，取停止標(biāo)準(zhǔn)中參數(shù)ε=0 .01。圖1給出了各種方法的去噪結(jié)果，其中圖1(a)為含噪圖，圖1(b)為文獻(xiàn)[7]方法(時間步長Δt=0.1)去噪后的結(jié)果，圖1(c)為文獻(xiàn)[9]方法(α=0 .15,β=0 .23)去噪后的結(jié)果，圖1(d)為本文方法的處理結(jié)果，圖1(e)，圖1(f)，圖1(g)分別為圖1(b)，圖1(c)，圖1(d)的局部放大。比較3種去噪結(jié)果，可以看出，圖1(b)和圖1(c)中仍存在少量的塊效應(yīng)，如 Lena的額頭、臉頰和下巴，而圖1(d)看起來比較光滑、自然，并且對圖像的邊緣處理得也較好，如Lena的鼻梁比圖1(b)和圖1(c)都得到較好地保持，從放大的圖像中可以很清楚地看到這一點。 圖2給出了3種方法的峰值信噪比和均方誤差隨重構(gòu)圖時的迭代次數(shù)的變化曲線圖。從圖2(a)可以看出，本文方法在重構(gòu)圖像時僅需很少的迭代就可以達(dá)到穩(wěn)定，并且其峰值信噪比的最大值高于其它兩種方法，說明本文方法比其它兩種方法具有較好的去噪效果。除此之外，本文方法還具有較快的收斂速度，這從圖2(b)可以清楚地看到。

圖2 實驗1中峰值信噪比和均方誤差隨重構(gòu)圖時的迭代次數(shù)的變化曲線

實驗 2 選取比較難處理的“Cameraman”圖像作為測試目標(biāo)，并對其加入均值為0、方差為15的高斯白噪聲?！癈ameraman”圖像中不但含有像素跳躍區(qū)域(相機(jī)支架)、圖像漸變區(qū)域(天空)，而且還含有圖像震蕩區(qū)域，或者說是紋理(如草坪)。試驗中取參數(shù)λ=0 .8,β1=0 .05，參數(shù)β2,μ的取值和實驗1的相同。圖3給出了文獻(xiàn)[7]、文獻(xiàn)[9]方法和本文方法的去噪結(jié)果，其中圖3(a)為噪聲圖，圖3(b)為文獻(xiàn)[7]方法(時間步長Δt=0.1)的去噪結(jié)果，圖3(c)為文獻(xiàn)[9]方法(α=0 .15,β=0 .14)的去噪結(jié)果，圖3(d)為本文方法的處理結(jié)果，圖3(e)，圖3(f)，圖3(g)分別為圖3(b)，圖3(c)，圖3(d)與圖3(a)的差圖。從3種差圖可以明顯地看到，圖3(e)中含有較多的邊，而圖3(f)與圖3(e)相比，雖然含有較少的邊，但也含有較少的噪聲，這說明文獻(xiàn)[9]中方法在保邊方面要優(yōu)于文獻(xiàn)[7]中方法，但其去噪性能略弱于文獻(xiàn)[7]。與前兩種差圖相比較，本文方法相應(yīng)的差圖中不僅含有較少的邊，而且還含有較多的噪聲，這說明了本文方法在去噪和邊緣保持方面均優(yōu)于其它兩種方法。表1中的實驗數(shù)據(jù)也從客觀上表明了本文方法的優(yōu)越性。

圖3 Cameraman 圖實驗結(jié)果

表1列出了兩實驗中各種方法峰值信噪比，均方誤差和結(jié)構(gòu)相似性的比較結(jié)果。從表中數(shù)據(jù)可以看出，文獻(xiàn)[7,9]幾乎有著同樣的復(fù)原結(jié)果，而本文方法的PSNR和SSIM均有明顯地提高。表2給出了各種方法達(dá)到停止標(biāo)準(zhǔn)時所需的執(zhí)行時間和兩步中分別所需的迭代次數(shù)。由各種方法的迭代次數(shù)可以看出，本文方法通過有效使用交替迭代以及分裂Bregman方法，使其計算復(fù)雜度大大降低，導(dǎo)致它的計算時間較其它兩種方法要少得多，從而使其達(dá)到加速收斂的目的。

表1 各種方法的恢復(fù)圖像定量指標(biāo)比較

表2 各種方法達(dá)到停止標(biāo)準(zhǔn)所需的時間和兩步中分別所需的迭代次數(shù)

5 結(jié)論

本文在研究LOT模型的基礎(chǔ)上，提出了一種新的變分去噪模型。新模型由于耦合了兩個變量，因此，在求解時，先利用交替極小化方法化原模型為兩個簡單的子模型，然后再利用分裂Bregman方法對其進(jìn)行交替求解。由于兩子模型對LOT模型中的兩步從構(gòu)成方式及計算方式上分別進(jìn)行了改進(jìn)，因此，數(shù)值實驗表明，新方法不僅能夠有效地提高恢復(fù)圖像的質(zhì)量，而且與文獻(xiàn)[7]中方法相比收斂速度有明顯改進(jìn)。

[1]Perona P and Malik J.Scale space and edge detection using anisotropic diffusion[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1990,12(7)：629-639.

[2]Rudin L,Osher S,and Fatemi E.Nonlinear total variation based noise removal algorithms[J].Physica D,1992,60(1-4)：259-268.

[3]張偉斌,馮象初,王衛(wèi)衛(wèi).圖像恢復(fù)的小波域加速 Landweber迭代閾值方法[J].電子與信息學(xué)報,2011,33(2)：342-346.Zhang Wei-bin,Feng Xiang-chu,and Wang Wei-wei.Accelerated Landweber iterative thresholding algorithm in wavelet domain for image restoration[J].Journal ofElectronics＆Information Technology,2011,33(2)：342-346.

[4]武曉玥,郭寶龍,李雷達(dá).一種新的結(jié)合非下采樣 Contourlet與自適應(yīng)全變差的圖像去噪方法[J].電子與信息學(xué)報,2010,32(2)：360-365.Wu Xiao-yue,Guo Bao-long,and Li Lei-da.A new image denoising method conbining the nonsubsampled contourlet transform and adaptive total variation[J].Journal of Electronics＆Information Technology,2010,32(2)：360-365.

[5]張軍,韋志輝.SAR圖像去噪的分?jǐn)?shù)階多尺度變分PDE模型及自適應(yīng)算法[J].電子與信息學(xué)報,2010,32(7)：1654-1659.Zhang Jun and Wei Zhi-hui.Fractional-order multiscale variation PDE model and adaptive algorithm for SAR image denoising[J].Journal of Electronics＆Information Technology,2010,32(7)：1654-1659.

[6]Lysaker M,Lundervold A,and Tai X C.Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time[J].IEEE Transactions on Image Processing,2003,12(12)：1579-1590.

[7]Lysaker M,Osher S,and Tai X C.Noise removal using smoothed normals and surface fitting[J].IEEE Transactions on Image Processing,2004,13(10)：1345-1357.

[8]Hahn J,Tai X C,Borok S,et al..Orientation-matching minimization for image denoising and inpainting[J].International Journal of Computer Vision,2011,92(3)：308-324.

[9]Yang Yu-fei,Pang Zhi-feng,Shi Bao-li,et al..Split Bregman method for the modified LOT model in image denoising[J].Applied Mathematics and Computation,2011,217(12)：5392-5403.

[10]Goldstein T and Osher S.The split Bregman method for L1 regularized problems[J].SIAM Journal on Imaging Sciences,2009,2(2)：323-343.

[11]Goldstein T,Bresson X,and Osher S.Geometric applications of the split Bregman method：segmentation and surface reconstruction[J].Journal of Scientific Computing,2010,45(1-3)：272-293.

[12]Wu Chun-lin and Tai X C.Augmented Lagrangian method,dual methods and split-Bregman iterations for ROF,vectorial TV and higher order models[J].SIAM Journal on Imaging Science,2010,3(3)：300-339.

[13]Cai J F,Osher S,and Shen Z W.Split Bregman methods and frame based image restoration[J].Multiscale Modeling and Simulation,2009,8(2)：337-369.

[14]孫玉寶,費(fèi)選,韋志輝,等.稀疏性正則化的圖像泊松恢復(fù)模型及分裂 Bregman迭代算法[J].自動化學(xué)報,2010,36(11)：1512-1519.Sun Yu-bao,Fei Xuan,Wei Zhi-hui,et al..Image restoration model under poisson noise using sparse representations and split Bregman iteration algorithm[J].Acta Automatica Sinica,2010,36(11)：1512-1519.