小波域Wiener濾波和Perona-Malik融合去噪的新算法*

尹 芳，付自如，于曉洋

1.哈爾濱理工大學 計算機科學與技術學院，哈爾濱 150080

2.哈爾濱理工大學 儀器科學與技術博士后科研流動站，哈爾濱 150080

1 引言

去噪是圖像預處理的一個重要課題。由于小波變換具有多分辨率和局部時頻分化的特性，從而能夠做到高頻處時間細化，低頻處頻率細分，并且小波變換還能夠直接處理一些圖像特征，便于靈活地搭建各種數學模型，進而得到廣泛的應用。維納濾波是一種基于最小均方誤差準則，對噪聲敏感的最佳線性濾波器。傳統的小波域維納濾波去噪就是兩種算法融合的一種應用，但是該算法的去噪效果不理想，并且存在Gibbs偽影。針對此問題，許多改進的去噪算法被研究者提出。文獻[1]提出一種基于小波變換的圖像去噪算法，其思想是利用脈沖耦合神經網絡（pulse coupled neural network，PCNN）的同步脈沖特性對圖像小波系數進行局部加窗、加權處理，結合維納濾波進行局部去噪，并且對每個像素的連接系數β進行修正，以保護圖像邊緣和細節信息，提高圖像質量，但是PCNN模型中參數都需要根據不同的圖像由經驗而定。文獻[2]提出一種基于維納濾波空間域噪聲方差自適應的圖像濾波算法，其思想是利用Canny算子把圖像分成平滑區和邊緣區，計算各平滑區域的噪聲方差，并將求得平均值作為平滑區域的噪聲估計值，進行維納濾波去噪處理。該算法只對平滑區域進行去噪，忽略邊緣區域的噪聲影響，因此僅對邊緣細節少的圖像有較好的去噪效果。文獻[3]提出一種模糊濾波器和維納濾波器結合去噪的方案，其思想是首先由小波變換提取圖像的4個子帶系數矩陣，然后利用維納濾波器對逼近子帶的小波系數進行去噪，水平、垂直、對角方向的子帶系數則由模糊濾波器處理。該算法雖然圖像質量得到提升，但是小波系數不連續，圖像的低頻區仍存在Gibbs偽影。文獻[4]將含噪圖像進行小波分解，利用改進的維納濾波窗口對高頻小波分量去噪，使用綜合的偏微分方程模型對處理后的圖像進行二次去噪，去除Gibbs偽影，但該算法計算復雜，且偏微分方程迭代次數過多，造成圖像部分呈現階梯效應。文獻[5]提出基于雙樹復小波變換的維納濾波收縮去噪技術，其思想是把雙樹復小波變換的強方向、多尺度平移不變性同維納濾波的局部強處理能力融合進行去噪。該算法雖然可避免圖像產生Gibbs偽影，使得圖像質量得到提高，但是由于雙樹復小波變換去噪會得到大量冗余信息，造成該算法的時間性能差。文獻[6]利用小波在單尺度上的平移不變性，采用軟硬閾值收縮小波系數后重構新的信號，其思想是通過周期性的變換-去噪-變換方式獲得大量小波系數，使其平均值作為最終的小波系數。雖然小波基之間的關聯性得到減弱，但是該算法僅對一維信號的去噪效果有明顯作用，而對二維多尺度圖像的去噪效果就不明顯。文獻[7]在小波域采用局部維納濾波器平滑該區域的塊效應，并根據領域方差大小自適應調整像素的處理方式以達到保護細節的目的。雖然塊效應與偽Gibbs效應都是因小波系數不連續造成圖像重構產生的偽影，但在Gibbs偽影抑制方面得不到理想的效果。上述去噪算法各有千秋和不足，在圖像去噪領域提供了新的思路。

近年來偏微分方程去噪算法得到迅速發展，它通常對采樣信號進行局部熱擴散，經過多次迭代運算后逼近得到去噪的圖像，從而該算法具有良好地光滑信號的作用，對Gibbs偽影有很好的抑制能力，但若圖像去噪迭代次數過多，則會產生階梯效應。本文將小波域維納濾波作為圖像去噪的迭代基，迭代過程中加入自適應的噪聲權系數η以保護圖像信息不遭受過多損失，并修正去噪量的大小。由于閾值收縮去噪的圖像去除小值小波系數后，圖像在重構時易在邊緣等奇異處引起震蕩從而產生Gibbs偽影，最后使用Perona-Malik[8]算法作為本文算法的后處理階段去除偽Gibbs效應。

2 偏微分方程去噪算法

Perona-Malik[8]是由Perona和Malik基于偏微分方程各向同性擴散提出來的各向異性擴散模型，其迭代表達式為：

其中，?y為梯度，值為；div表示擴散算子；g是與梯度成反比的函數，稱為擴散系數。g在0處等于1，而在無窮處等于0，由于圖像邊緣的梯度很大，此時g很小，則擴散較弱；圖像同質區域的梯度較小，此時g很大，則擴散較強。該算法可以實現在同質區域圖像被平滑，而在邊緣區域圖像細節保持不變，但其在同質區域處又有可能把噪聲當作邊緣而造成階梯偽影，降低了圖像的主觀質量。因此偏微分方程對閾值收縮去噪產生的Gibbs偽影有很好的抑制作用，但若去噪中迭代次數過多則很大可能產生階梯效應。

3 小波域維納濾波去噪算法

圖1去噪方程可表示為y=y0-h，其中：y0為原含噪圖像，h表示方差為?2的零均值高斯噪聲。令W為小波變換中的變換矩陣，由小波正交性變換可得Y=Y0-H，其中Y=Wy，Y0=Wy0，H=Wh，H表示小波域上方差為?2的零均值高斯噪聲。由小波變換原理可知，H與Y0是互為不相關的，因此小波域上的維納濾波系數的形式可簡化為標量關系，如式（2）：

Fig.1 Denoising process of wavelet domain Wiener filter圖1 小波域維納濾波去噪的過程

E可利用及其領域值估計得到，不失一般性，可以利用一個(2R+1)×(2R+1)方形窗中的平均值求得：

即將qi,j作為的近似估計值，把qi,j代入式（3）可得：

維納濾波的系數可由式（5）代入到式（2）中得到：

4 融合去噪新方案

4.1 小波域維納濾波和Perona-Malik融合算法

本文將對含噪圖像y0做連續的迭代，小波域維納濾波被視為一種擴散迭代基（即在每個迭代階段采用不同的小波系數，以使圖像得到不同尺度的處理），迭代中噪聲權系數η起到修正?2大小及保護圖像信息的作用，最后迭代得到的圖像yn再由Perona-Malik[8]算法進行二次去噪，平滑偽Gibbs效應。本文將此算法稱為小波域維納濾波和Perona-Malik融合算法（wavelet-domain Wiener withPerona-Malikfusion，WWPM），其算法流程如圖2所示。

圖2中A算法為迭代的小波域維納濾波，ηi(i=0,1,…)為迭代階段的噪聲權系數，它是由η權值生成器產生的。式（7）是本文算法給出的圖像去噪方程：

式中，η表示噪聲權系數。將圖像關系式（7）變換到小波域上的迭代方程為：

Fig.2 Process of this paper algorithm圖2 本文算法流程圖

x表示(i,j)像素點處的第m層小波分解，這里記為x(m:i,j)；k表示第k次迭代。根據小波變換的正交性原理可知，在各次迭代中方差為的零均值高斯噪聲經小波變換到小波域里仍是方差為的零均值高斯噪聲；因小波變換具有去相關性，可假設信號成分(Yk)i,j(x)(k=0，1,…)相互獨立。

小波域上的維納濾波系數的形式簡化為標量關系如下：

從而維納濾波各系數可由式（12）代入式（9）中得：

將ei,j帶入式（10）可重構去噪圖像，如此迭代n次得到的最后重構圖像極有可能產生偽Gibbs效應。由于Gibbs偽影是一種明顯的振動不光滑信號，可由Perona-Malik[8]算法進行二次去噪，去除Gibbs偽影。

4.2 噪聲權系數的產生

由式（6）可知，信號方差的估算對局部維納濾波性能有很大的影響作用。然而，由于噪聲的影響，這種估算必然不精確。從整體來看，圖像的高斯白噪聲均值雖然為0，但是決不能由此得出局部窗內高斯白噪聲的均值也為0。而式（5）是假定局部窗內噪聲均值為0的，這顯然不能很好地符合事實，于是在式（6）中，局部維納濾波中引入修正因子是很有必要的。

噪聲權系數η的確定是根據本次得到的重構圖像領域的內子帶的邊緣標準差最佳估計值DMAP與前一次迭代的最佳估計值Dopt進行比較得到的。參數η能夠起到保護細節和對噪聲方差?2修正的作用，η越小，保護細節的能力就越強，但是抑噪能力就越弱，反之η越大，抑噪能力就越強，但是保護細節的能力就越弱。由式（7）和（8）可知，若使yk+1圖像信息更加真實，Hk和η值的確定就顯得尤為重要。小波域中噪聲權系數定義公式為：

其中，η0表示η的初始給定值；ηv表示控制η增大速度的參數；t表示η的優化次數。

η值的確定思想：分別計算本次迭代得到的重構圖像領域的內子帶的邊緣標準差估計值的后驗概率和前一次迭代的最佳估計值的后驗概率Popt：

式中，w為真實小波系數；P(w)為小波系數w的先驗分布；P(D(t))為第t次優化后的當前邊緣標準差的概率，若P(D?(t)/w)＞ηPopt，則t.加1，否則就保持不變，由此可以優化計算出參數η。

Fig.3 Optimization of?MAP圖3 ?MAP的優化

該算法的優化示意圖如圖3所示，領域內子帶的邊緣標準差估計值是利用上述思想求得η的關鍵。求解值的估計思想為：大量自然圖像的小波系數具有非高斯特性，各細節子帶的小波系數都基本在0附近對稱分布，這種分布可近似由拉布拉斯描述[9-11]，因此本文將子帶內小波系數視為不同邊緣標準差的拉普拉斯分布，而邊緣標準差又假設為其強局部相關的隨機變量，從而以當前系數為中心，可利用其鄰域窗口中的小波系數來估計當前系數的邊緣標準差，進而可得有局部特性的MAP估計值。計算局部系數邊緣標準差MAP估計值D?MAP步驟如下。

步驟1設噪聲h（高斯噪聲）概率密度分布為：

步驟2要計算出邊緣標準差的值，首先須對小波系數w做MAP估計，而在給定小波觀測系數y的情況下，可計算出使后驗概率密度Pw|y(w|y)最大時的w估計值：

步驟3每一個含噪觀測子帶的方差估計可采用最大似然估計方法得到：

式中，N表示窗口中小波系數的個數；M(i)表示y的領域窗口。則當前系數的邊緣標準差估計為：

式中，可采用Donoho提出的魯棒性中值估計[12]。已知噪聲方差?h為:

模型邊緣標準差以子帶為單位進行估計，則拉普拉斯分布的邊緣標準差?計算如下：

5 仿真實驗

5.1 實驗重要部分說明

為驗證本文方法的有效性設計實驗如下：以圖4中的加噪圖像降噪處理為例進行Matlab仿真實驗，實驗采用大小均為512×512像素的圖像。在實驗中的小波域維納濾波去噪階段，設n=5，即圖2中yn為y5，小波域維納濾波共迭代5次，各個迭代階段依次選取的小波基為sym9、sym7、sym5及sym3、sym1，而在每次小波變換中圖像又被分解5層，本文使用F1=[7,7,5,5,3]來表示第一次迭代中的各層正方形卷積窗口，其中7、7、5、5、3依次為第一層、第二層、第三層、第四層、第五層窗口的大小。同樣記法，則第二、三、四、五次迭代各層中的窗口大小可依次記為F2=[7,5,5,3,3]、F3=[5,5,3,3,3]、F4=[5,3,3,3,3]、F5=[3,3,3,3,3]。鑒于η值主要受優化次數t的影響，則η值確定的時間復雜度是O(n2)；考慮到ηv可任意取值的情況，t可能會使得的比較結束之后過大，不易控制，這里取ηv=1.5η0，增強其可操作性，以便降低實驗的復雜性，保證其時間性能。在小波域維納濾波去噪階段結束后，采用Perona-Malik[8]算法對圖像yn進行二次去噪，由經驗驗證該算法只需采取3次迭代。

針對圖4中所有仿真圖像的類型特點及去噪難度進行如下說明：Angle、House等圖像具有平滑區域多，邊緣、細節少的特點，此類圖像的信號易于噪聲區別，去噪過程應該主要進行低通濾波；Lena、Agriculture、Chili等圖像具有平滑區和邊緣區較多，細節、紋理較少的特點，此類圖像的信號與噪聲較易區別，去噪相對容易；Hair、Boat、Racetrack、Barbara、Orangutan等圖像具有豐富的細節、紋理特點，此類圖像去噪比較難，這是因為去噪時細節很容易和噪聲混淆，造成圖像信息損失。

5.2 實驗中η的確定

由4.2節可知，若η初值較大就會失去其優化噪聲量的意義，甚至對圖像信息造成損失，因此實驗中η0的確定思路是：鑒于η0與η呈正相關，則首先給定參數η0一個較小值，然后η0以步長d0=0.2增大。實驗中取得的部分峰值信噪比（peak signal to noise ratio，PSNR）如表1所示，處理后相應的圖像如圖5所示。實驗的其他數據環境如5.1節所述，結合表1與圖5進行如下描述：

Fig.4 Original images and noise images圖4 原始圖與加噪圖

（1）從表1中可以看出，在η0逐漸變大的情況下Angle、Boat、Chili和 Hair的序號為 1～3的圖像的PSNR值是不斷增大的，序號為4～5的圖像的PSNR值開始遞減。

Table 1 PSNR values of images byWWPM in the case ofη0表1 在η0給定值下圖像經WWPM算法得到的PSNR值

（2）從圖5中各個圖像處理后的主觀效果可看出，Angle、Boat、Chili和Hair的1～3序號圖像的噪聲不斷減少并且Gibbs偽影逐漸變淡，圖像細節也不斷地清晰可見，但是序號為4～5的Hair、Boat、Angle圖像，圖像的細節清晰度就明顯不如序號為3的圖像。

上述結果可做如下分析：在η值確認階段，D?MAP估計值的優化次數t大于給定的小波系數觀察值y的個數，從而使得η值偏小，即η0=0.6時，Angle、Boat、Chili和Hair的PSNR值偏小，而當η0=1.4 時，η0的初始值較大，從而η值偏大，使得Angle、Boat、Chili和Hair的PSNR值偏小。由表1可知，η0最佳取值范圍是(0.8,1.2)。

5.3 實驗結果展示

Fig.5 Images of optimal PSNR value by WWPM with different parameters in Table 1圖5 WWPM算法在表1中不同參數下最佳PSNR值時的圖像

本文與其他參考文獻的算法分別對含有不同零均值方差的高斯噪聲的兩組圖像去噪：含有方差為25 dB的Angle、Hair、House、Lena圖像組和含有方差為20 dB的Chili、Agriculture、Racetrack、Orangutan圖像組。圖6中（a）、（b）、（c）、（d）、（e）為兩組中含噪圖像被不同算法處理后取得表2中最佳峰值信噪比（PSNR）時的圖像組。其中圖組（a）是由小波域維納濾波去噪算法（文獻[13]）處理后的圖像組；圖組（b）是由Perona-Malik[8]算法去噪后的圖像組；圖組（c）是本文算法去噪后的圖像組（參數為η0=0.85,ηv≈1.27，t=28，其他數據環境如5.1節所述）；圖組（d）是 ProbShrink[14]算法處理后的圖像組；圖組（e）是文獻[15]算法去噪后的圖像組。去噪后圖像的主觀效應有以下三點可述：

（1）圖組（b）中Angle、Hair、House、Lena、Barbara、Agriculture、Racetrack、Orangutan圖像有明顯的細節丟失現象；圖組（a）、（d）、（e）中Hair、Racetrack、Barbara圖像在細節保護方面基本相當，且細節明顯比圖組（b）的圖像豐富，沒有明顯丟失現象，但在Orangutan圖像的眼睛和胡須周圍有明顯丟失少許細節的現象，這是由于噪聲與眼睛和胡須周圍的微小細節混淆而造成細節丟失；而圖組（c）全體圖像細節保存得相對較好。

（2）圖組（b）、（c）中全體圖像基本上去除Gibbs偽影；而圖組（a）中Angle、Hair、Lena、Barbara、Orangutan圖像仍有較為明顯的偽Gibbs效應，這是由于小的維納濾波系數被直接歸零造成小波重構后產生偽Gibbs效應；圖組（d）、（e）中Hair、Barbara、Orangutan圖像均有輕微的Gibbs偽影。

（3）根據表2中圖像處理后的PSNR數值對比表明，本文算法在圖像質量提升方面優于另外4種算法，即從主客觀上都能夠說明本文算法具有良好的去噪和抑制偽Gibbs效應的能力。

Fig.6 Images of optimal PSNR value of each method in Table 2圖6 表2中各算法的最佳PSNR值時的圖像

Table 2 Comparison of optimal PSNR and MSE values by difference methods表2 各種算法下的最佳PSNR和MSE值對比

Table 3 Denoising time of difference methods表3 不同算法下的去噪時間

表3是不同算法在同一臺計算機上運行所花費的CPU時間。根據表中Average可以看出，文獻[13]算法（即小波域維納濾波）所花費的時間最少，平均時間為0.214 1 s，但是該算法的去噪效果很差。在其余方法中，Perona-Malik[8]算法花費的平均時間最長，為5.518 8 s，這是因為所測試的圖像中包含細節信號的圖像要多一些，使得在細節系數估計上花費更多的時間。對于本文算法、ProbShrink算法[14]、文獻[15]算法而言，盡管本文算法運行時間比ProbShrink算法[14]、文獻[15]算法多一點，但是其處理的圖像效果卻比較好，這點花費是值得的。

6 結論

本文針對傳統的小波域維納濾波去噪效果不理想，并存在Gibbs偽影的問題，提出以小波域維納濾波為迭代基，迭代過程中又通過噪聲權系數自動修正去噪量，最后使用3次Perona-Malik[8]算法平滑小波重構后產生的Gibbs偽影的方案。實驗結果表明本文算法既能去除高斯噪聲，又能有效地去除因小波重構后產生的偽Gibbs效應，很好地保護了圖像細節和邊緣信息，得到高質量的圖像視覺效果。雖然本文算法在時間性能上不是很好，但是在圖像處理效果方面比小波域維納濾波[13]、Perona-Malik[8]、ProbShrink[14]和文獻[15]算法等優越不少。本文算法對于噪聲是廣義平穩隨機過程的圖像有良好的去噪效果，應用到磁共振取像中能很好地保持圖像的邊緣和紋理細節信息。

[1]Gong Xiaolin,Mao Ruiquan.An image de-noising algorithm based on PCNN and local Wiener filter[J].Journal of Beijing University of Posts and Telecommunications,2011,34(5):67-70.

[2]Bai Junqi,Zhao Chunguang,Wang Shoufeng,et al.Adaptive Wiener filtering noise reduction in infrared images[J].Opto-Electronic Engineering,2011,38(11):79-85.

[3]Kethwas A,Jharia B.Image de-noising using fuzzy and Wiener filter in wavelet domain[C]//Proceedings of the 2015 International Conference on Computer and Communication Technologies,Coimbatore,Mar 5-7,2015.Piscataway:IEEE,2015:1-5.

[4]Chen Shen.Research on image denoising algorithm based on wavelet transform and partial differential equation[D].Xi’an:Northwest University,2010.

[5]Naimi H,Adamou-Mitiche A B H,Mitiche L.Medical image denoising using dual tree complex thresholding wavelet transform and Wiener filter[J].Journal of King Saud University:Computer and Information Sciences,2015,27(1):40-45.

[6]Zhong Guosheng,Deng Yingxiang,Ao Liping.Study and application of translation invariant wavelet de-noising for blasting seismic signals[C]//Proceedings of the 2011 International Conference on Multimedia Technology,Hangzhou,Jul 26-28,2011.Piscataway:IEEE,2011:4100-4103.

[7]Nath V K,Hazarika D.Blocking artifacts suppression in wavelet transform domain using local Wiener filtering[C]//Proceedings of the 3rd National Conference on Emerging Trends and Applications in Computer Science,Shillong,Mar 30-31,2012.Piscataway:IEEE,2012:93-97.

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

[9]Tang Chongwu.The key technology of image statistical modeling and analysis of noise[D].Shanghai:Shanghai Jiao Tong University,2015.

[10]Liu J,Moulin P.Information-theoretic analysis of interscale and intrascale dependencies between image wavelet coefficients[J].IEEE Transactions on Image Processing,2001,10(11):1647-1658.

[11]Ha Y M,Yoon J W.An adaptive approximation for Gaussian wavelet kernel[C]//Proceedings of the 18th International Conference on Advanced Communication Technology,Pyeongchang,Jan 31-Feb 3,2016.Piscataway:IEEE,2016:576-580.

[12]Donoho D L,Johnstone I M.Adapting to unknown smoothness via wavelet shrinkage[J].Journal of the American StatisticalAssociation,1995,90(432):1200-1224.

[13]Liu Pengju,Li Hong.Image denoising based on wavelet domain Wiener filtering[J].Computer Simulation,2005,22(9):269-271.

[14]PizuricaA,Philips W.Estimating the probability of the presence of a signal of interest in multiresolution single-and multiband image denoising[J].IEEE Transactions on Image Processing,2006,15(3):654-665.

[15]Hou Yingkun,Liu Mingxia,Yang Deyun.Multilevel block matching transform domain filtering image noise[J].Journal of Computer-Aided Design and Computer Graphics,2014,6(2):225-231.

附中文參考文獻：

[1]宮霄霖,毛瑞全.結合PCNN和局部維納濾波的圖像去噪[J].北京郵電大學學報,2011,34(5):67-70.

[2]白俊奇,趙春光,王壽峰,等.紅外圖像中的自適應維納濾波噪聲抑制技術[J].光電工程,2011,38(11):79-85.

[4]陳燊.基于小波變換與偏微分方程的圖像去噪算法研究[D].西安:西北大學,2010.

[9]唐崇武.圖像統計建模與噪聲分析關鍵技術的研究[D].上海:上海交通大學,2015.

[13]劉鵬舉,李宏.基于小波域維納濾波的圖像降噪技術[J].計算機仿真,2005,22(9):269-271.

[15]侯迎坤,劉明霞,楊德運.多級塊匹配變換域濾波圖像去噪[J].計算機輔助設計與圖形學學報,2014,26(2):225-231.