基于SURE無偏估計的圖像自適應稀疏收縮去噪

沙正虎，余 劍，崔 琛，2

1.解放軍電子工程學院 信息工程系，合肥 230037

2.安徽省電子制約技術重點實驗室，合肥 230037

圖像在攝取或傳輸過程中容易受到噪聲的污染，從而影響圖像信息的提取，因此，去噪是圖像處理技術的重要組成部分。圖像去噪既要盡可能去除噪聲分量，又要盡可能保留圖像的細節成分，如紋理和邊緣，而各種圖像去噪其實就是在去除噪聲和保留細節之間進行權衡。為了克服這個問題，Donoho和Johnstone提出了小波收縮法[1]，利用收縮函數對圖像小波域系數有選擇性地進行取舍，以達到去噪的目的。自此以后，無論是在正交基變換（如小波等），還是緊框架變換（如脊波、曲波等），收縮都成為最常用的去噪算法。目前，收縮主要的研究方向集中在兩點：一是收縮曲線的選擇，常用的有硬閾值和軟閾值兩種；另一點則是閾值的選取問題。對于閾值的選取，目前主要有Visu通用閾值[2]、MiniMaxi閾值[3]、SURE閾值[4]（Stein Unbiased Risk Estimator，SURE）等。在上述閾值中，Visu通用閾值（，σ為標準方差，N為信號長度）計算簡單，但其趨向于“過扼殺”變換系數，從而會導致過多的高頻信息流失；MiniMaxi閾值，由于基于悲觀決策的思想，即最大均方誤差最小化，所以也會過扼殺現象；SURE閾值估計是一種基于Stein無偏風險估計準則的閾值選擇，該準則是均方差準則的無偏估計，趨近于理想閾值。

近幾年，稀疏表示作為信號處理領域一個新的工具，越來越引起人們的重視。由于圖像本身是稀疏的，所以利用稀疏表示來處理圖像具有極大的優越性。傳統的小波變換（Wavelet）能夠表示圖像的點狀奇異特征，但對高維信息表示能力不足；Donoho等學者提出的多尺度幾何變換[5]，主要有 Rideglet、Curvelet、Contourlet、Bandlet等變換，這些變換充分考慮了圖像的某些幾何特性，能比小波變換更好地捕捉圖像的邊緣等奇異特性，但每種變換只能對應某些特征表示是稀疏的；圖像的稀疏表示，則是用超完備字典取代傳統中的正交基或緊框架來表示圖像，而字典的這種冗余性恰恰能夠更好地刻畫圖像的各種奇異特性，且可以同時表示多種奇異特性。

針對以上問題，本文研究了稀疏框架下收縮去噪問題，提出一種基于SURE無偏估計的自適應閾值選擇算法（Adaptive SURE，AdpSURE）。本文算法首先推導稀疏框架下收縮去噪的SURE無偏估計目標函數，然后利用黃金分割法搜索全局最小值點，對應的閾值即為自適應選擇的閾值。仿真結果驗證了本文算法的優越性。

1 圖像的稀疏表示

1.1 稀疏表示

自1991年Mallat和Zhang[6]提出匹配追蹤算法（Matching Pursuit，MP）以來，圖像的稀疏表示理論得到極大的發展，并取得了豐碩的成果。其基本思想是超完備字典取代傳統中的正交基，字典的選擇應盡可能地包含被表達圖像所含有的特殊結構。其構成可以沒有任何限制，字典中的元素被稱為原子。

設 y∈?N為N維的含噪圖像序列，A∈?M×N的過完備字典，考慮圖像為加性高斯白噪聲，圖像稀疏表示可以用求解x的l0范數表示：

式中，ε是誤差因子，x是原圖像的稀疏表示。現已證明，精確求解式（1）是一個NP難題[7]。幸運的是，近些年研究人員提出了多種有效的近似求解算法，主要可以概括為貪婪追蹤[6]和凸松弛[8]兩類算法，常見的有MP和BP（Basis Pursuit）等。其中MP算法計算簡單，并且可以獲得較高的重構精度，因此在圖像處理領域得到廣泛應用。下面就簡要描述MP算法的原理。

1.2 匹配追蹤算法

其中，gψ0為字典中使殘差能量最小的原子，即展開信號與原子 gψ0的內積最大。

由于 gψ0與 R1y正交，因此滿足

同樣的方法，可以得到：

其中RMy為M項近似的殘差，并且滿足下式：

Mallat已對匹配追蹤的收斂性進行過分析，結果表明匹配追蹤時收斂的，并且在有限維條件下滿足指數級收斂[6]。

2 圖像稀疏收縮去噪框架

收縮是目前研究最為廣泛的去噪策略，其主要思想是：圖像主要由一些平滑的不重疊區域組成，這些區域是以一些邊緣為界限的。因此，去噪過程中不僅要保留含有大量信息的低頻分量，也要盡量保留含有邊緣信息的高頻分量。而收縮通過對收縮函數和閾值的選擇，能夠較好地保留圖像的紋理和邊緣信息，避免去除噪聲同時造成圖像的模糊。

目前，收縮主要應用于傳統的正交基和緊框架。在稀疏框架下，使用超完備字典取代傳統的正交基或緊框架[9]。因此，稀疏框架下收縮的過程，如圖1所示。

圖1 稀疏框架下收縮去噪框圖

設A為列l2標準化的冗余字典，為了估計式（1）的解，采用如下的收縮公式：

式中，A+表示字典 A的Moore-Penrose逆，λ表示閾值，Sλ(?)表示收縮函數，參量為矢量時輸出結果為矢量。收縮去噪的關鍵在于如何選擇合適的收縮函數和收縮閾值。對于收縮函數，Donoho提出了兩種常用的收縮函數：硬閾值函數和軟閾值函數[1]。但是這兩種函數導數都不連續，而很多情況下需要對收縮函數進行一階甚至高階導數運算，通常采用如下高階可導的收縮函數[11]：

其中，k表示偶常數。k取值越大，對應的曲線越接近硬閾值函數，而當k取20時，收縮曲線實際上已經非常接近硬閾值。

閾值選擇問題是本文研究的一個重點，下面介紹本文提出的一種基于SURE無偏估計的自適應閾值選擇算法。

3 基于SURE無偏估計的自適應閾值選擇

3.1 算法的提出

圖像去噪，實質上是盡可能地使圖像的估計值接近本身的無噪圖像。這也是圖像去噪的難點，因為在實際過程中往往很難得到真實圖像本身。Stein提出的Stein無偏風險估計[12]很好地解決了這個問題。SURE無偏估計推導如下：

上式是關于變量λ的隱函數，去噪的目的是使η(λ)最小化。函數 η(λ)在區間（0，+∞）是一個關于λ的凸函數，存在唯一的全局最小值。下面進行簡要證明。

3.2 η(λ)關于λ的凸函數證明

由于 k 是偶常數，易知 f(λ)、g(λ)在區間（0，+∞）都是凸函數。同時根據凸函數性質易知，兩個非負凸函數的正系數線性組合也是凸函數。故根據公式（18）可得，η(λ)在區間（0，+∞）是關于λ的凸函數，即 η(λ)在區間（0，+∞）存在唯一的全局最小值。由于黃金分割法僅需計算函數值，而且每次迭代只需要計算一次函數值，算法簡單，適用范圍廣，因此可以用黃金分割法對公式（15）進行全局最小值搜索，搜索所得到的閾值即為本文的自適應閾值。

3.3 算法分析

設圖像序列 y長度為N，字典 A大小為M×N，稀疏度為K。本文算法中，稀疏表示階段采用MP算法，其算法復雜度為O(K MN )，而在迭代收縮階段，黃金分割法每步迭代只計算一次目標函數值，在第γ步迭代中，計算一次目標函數計算復雜度為O(( k +1)MN+2(k2+3k)N ) 。設收縮階段總迭代次數為γo，則本文算法復雜度為O(KMN+γo((k+1)MN+2(k2+3k)N))。經過大量的仿真表明，γo的數值在14左右，最大不超過20，因此本文算法的復雜度在可接受范圍之內。

4 實驗仿真和結果分析

為驗證本文自適應閾值選擇算法的有效性，對超完備字典下Visu閾值、AdpSURE閾值去噪性能進行實驗對比。圖像采用疊加了不同強度高斯噪聲的大小為512×512的Lena、Barbara圖；字典由冗余Haar小波字典構成，字典冗余度為7；收縮函數如式（9）所示，取參數k為20，采用峰值信噪比作為判定標準。

圖2 兩種閾值選擇對應峰值信噪比曲線的位置圖

圖3 Lena圖像去噪效果局部放大圖

圖4 Barbara圖像去噪效果局部放大圖

實驗分為三個部分：（1）不同噪聲強度下，Visu閾值和AdpSURE閾值去噪結果的峰值信噪比，取50組實驗的平均值；（2）噪聲標準差取20，分別繪出Visu閾值和AdpSURE閾值對應峰值信噪比曲線的位置圖；（3）噪聲標準差取20，分別用Visu閾值和AdpSURE閾值對Lena、Barbara進行收縮去噪的局部放大圖。

實驗結果如圖2～圖4所示。表1給出了不同噪聲強度下兩種閾值選擇去噪的PSNR值，由表可知，不同噪聲強度下，本文的PSNR值均高于Visu的PSNR。圖2給出了兩種閾值選擇對應峰值信噪比曲線的位置，由圖可以看出，Visu閾值存在過扼殺的現象，而本文的非常接近PSNR的極大值點。圖3、圖4給出了對Lena和Barbara圖像的兩種去噪效果的局部對比圖，由圖可以看出，采用Visu圖片整體比較平滑，細節部分比較模糊，而本文較好地克服了這一點，無論是Lena圖像中頭發等邊緣特征，還是Barbara圖像中的紋理信息，恢復得均強于Visu閾值。

表1 不同噪聲強度下兩種閾值收縮去噪的PSNR dB

5 結束語

在圖像稀疏收縮框架下，基于一階可導收縮函數，提出了一種基于SURE無偏估計的自適應閾值選擇算法。算法推導了基于Stein無偏估計的目標優化函數，證明了該函數是關于閾值的凸函數，用黃金分割法搜索全局最小值點，獲得自適應的閾值。本文算法算法的優點體現于如下兩點：（1）算法基于SURE的無偏估計，去噪過程無需真實圖像本身；（2）閾值的選擇接近峰值信噪比曲線的最大值點。實驗結果證實了本文算法的有效性。相比于Visu閾值，本文算法無論是客觀指標還是主觀質量，均顯示出更加優越的性能，但其仍有較大改進空間，如字典的選取等，相關工作正在研究之中。

[1]Donoho D L.De-nosing by soft-thresholding[J].IEEE Trans on Information Theory，1995，41：613-627.

[2]Donoho D L，JohnstoneIM.Idealspatialadaptation by wavelet shrinkage[J].Biometrika，1994，81：425-455.

[3]Jonstone I M，Silverman B W.Wavelet threshold estimators for data with correlated noise[J].Journal of Royal Statistical Society：Series B，1997，59（2）：319-351.

[4]Donoho D L，Johnstone I M.Adapting to unknown smoothness via wavelet shrinkage[J].Amer Stat Assoc，1995，90：1200-1224.

[5]焦李成，譚山.圖像的多尺度幾何分析：回顧和展望[J].電子學報，2003，31（12）：1975-1981.

[6]Mallat S，Zhang Z.Matching pursuit in a time-frequency dictionary[J].IEEE Trans on Signal Process，1993，41（12）：3397-3415.

[7]張春梅，尹忠科，肖明霞.基于冗余字典的信號超完備表示與稀疏分解[J].科技通報，2006，51（6）.

[8]Chen S S，Donoho D L，Saunders M A.Atomic decomposition by basis pursuit[J].SIAM Rev，2001，43（1）：129-159.

[9]Elad M.Why simple shrinkage is still relevant for redundant representations[J].IEEE Transactions on Information Theory，2006，52（12）：5559-5569.

[10]曲天書，戴逸松，王樹勛.基于SURE無偏估計的自適應小波閾值去噪[J].電子學報，2002，30（2）.

[11]Elad M.Sparse and redundant representations：from theory to applications in signal and image processing[M].[S.l.]：Springer，2010：298-300.

[12]Stein C M.Estimation of the mean of a multivariate normal distribution[J].Ann Stat，1981，9（6）：1135-1151.