基于序列二次規劃優化閾值的NSCT高斯噪聲圖像濾波方法

楊 晨，董希旺，李青東，任 章

(北京航空航天大學 自動化科學與電氣工程學院，飛行器控制一體化國防科技重點實驗室，北京 100191)

0 引言

圖像在獲取和傳播的過程中通常會出現被各種噪聲污染的情況。為了得到高質量的圖像，圖像去噪作為圖像處理的重要部分受到了人們的廣泛關注。圖像去噪的關鍵在于如何在去除噪聲的同時盡可能地保留圖像的原始信息。小波變換(Wavelet Transformation,WT)是常用的圖像去噪算法之一，具有很好的處理點狀奇異性的能力[1]。然而，對于二維圖像等高維對象，小波變換缺乏良好的方向性和各向異性。研究已經證明小波變換更適合處理各向同性的對象，但對于處理數字圖像中邊緣和線特征等各向異性對象時表現不佳[2]。這也是許多基于小波變換的方法在圖像壓縮和去噪等應用中會引入模糊的原因[3]。近年來的研究中，基于多尺度變換的圖像去噪方法受到了人們的廣泛關注，許多基于Ridgelet、Curvelet和Contour-let變換的方法都表現出了很好的應用價值。

2002年，Donoho等提出了一種新的二維圖像表示方法——Contourlet變換[4],該方法能夠使用輪廓片段得到圖像多分辨率、多方向性的表示。然而，由于基礎的Contourlet變換在奇異點附近缺乏移不變性，因而引入了偽吉布斯效應[5]。針對這一問題，2006年Cunha等提出了非下采樣Contourlet變換(Nonsubsampled Contourlet Transform,NSCT)[6]。該方法實現了完全移不變性、多方向多尺度的圖像表示，通過非下采樣金字塔分解(Nonsubsampled Pyramids, NSP)和非下采樣方向濾波器(Nonsubsampled Directional Filter Banks, NSDFB)將圖像分解為不同的子帶。經過分解之后，噪聲部分相對于源圖像將具有較小的Contourlet系數[7]。基于Contourlet變換的圖像去噪的關鍵就在于如何選擇合適的閾值和閾值去噪函數對噪聲進行分離。

本文提出了一種基于NSCT的圖像去除高斯白噪聲的新方法。由于大多數圖像的紋理特征都是非對稱的，為每個方向子帶選取各自的閾值會比使用全局閾值的方法得到更好的去噪效果[8]。本文利用交叉驗證準則作為目標函數，基于序列二次規劃(Sequential Quadratic Programming,SQP)優化算法得到各個方向子帶的去噪閾值。在閾值確定之后，為了避免硬閾值和軟閾值去噪函數各自的缺點，本文采用非線性閾值函數對圖像的Contou-rlet系數進行處理，進而通過Contoutlet逆變換得到去噪之后的圖像。

1 非下采樣Contourlet變換

經典的Contourlet變換首先通過拉普拉斯金字塔(Laplacian Pyramid,LP)對圖像進行分解以捕獲點奇異，接著由方向濾波器(Directional Filter Banks,DFB)將分布在同一個方向的奇異點合成為一個系數[4]。Contourlet變換允許每層金字塔變換有不同的分解方向數目，因此相對小波變換有更好的靈活性,只要選擇合適的去噪閾值就能獲得比小波變換更好的圖像去噪效果。然而，由于金字塔變換和方向濾波器的下采樣和上采樣過程，經典的Contourlet變換缺乏移不變性而在奇異點附近導致了偽吉布斯效應。

為了克服這個缺點，基于Contourlet變換提出了非下采樣Contourlet變換。NSCT的結構分為非下采樣金字塔分解和非下采樣方向濾波器組分解兩部分。與Contourlet不同的是，NSCT在圖像的分解和重構過程中沒有對NSP及NSDFB分解后的信號進行分析濾波后的降采樣以及綜合濾波前的上采樣。而是對相應的濾波器進行上采樣，再對信號進行分析濾波和綜合濾波。NSCT分解的結構如圖1所示，其中NSP為非下采樣金字塔分解，NSDFB為非下采樣方向濾波器組，分別實現了圖像的多尺度、多方向性分解。

圖1 NSCT分解的結構示意圖Fig.1 The structure of NSCT decomposition

1.1 非下采樣金字塔分解

非下采樣金字塔分解實現了圖像的多尺度分解，保證了NSCT的多尺度特性。NSCT采用二通道移不變非下采樣濾波器組實現NSP分解。下一級的濾波器由對前一級濾波器進行上采樣得到，使得NSCT具有了多尺度特性[6]。圖2所示為一個三層金字塔分解結構示意圖。圖2中，一個三級金字塔能夠得到4層子圖像帶(y0,y1,y2,y3)，包括與源圖像尺度一致的1層低頻圖像和3層高頻圖像。金字塔分解的層數為3，H0和H1分別為不同尺度下的低通和帶通輸出。

圖2 三級NSP分解結構示意圖Fig.2 The three-stage NSP decomposition

1.2 非下采樣方向濾波器組

非下采樣方向濾波器組結構實現了NSCT的多方向性。NSDFB去掉了DFB中的降采樣和上采樣過程，使用扇形方向濾波器組將二維頻域平面分割為多個具有方向性的楔形結構，每個楔形塊包含該方向上的圖像細節特征，從而形成一個由多通道NSDFB組成的樹形結構[9]。NSDFB可以對任一尺度下NSP分解后的子帶圖像進行l級方向分解，從而得到2l個與源圖像具有相同尺寸大小的方向子帶圖像。圖3所示為一個四通道NSDFB方向分解示意圖，由幾組扇形濾波器和棋盤濾波器級聯而成，灰色部分代表濾波器的頻域通帶，實現了圖像的多方向分解。

圖3 四通道NSDFB方向分解示意圖Fig.3 Four-channel NSDFB decomposition

2 基于SQP優化閾值的NSCT圖像去噪

高斯白噪聲圖像去噪問題中一般采用式(1)所示的模型[10]

y=x+η

(1)

式中，y是觀測的含噪聲圖像，x為原始圖像，η為噪聲。圖像去噪的目的就是從含噪圖像y中恢復原始圖像x。

NSCT圖像去噪主要包含3個步驟：噪聲圖像的NSCT分解、NSCT系數的處理和NSCT逆變換得到去噪圖像。其中最為重要的部分在于NSCT系數的處理，關系到能否從含噪圖像中有效地去除噪聲并保留原始圖像信息。閾值去噪是一種簡單有效的系數處理方法，在圖像去噪問題中有著廣泛的應用。閾值去噪問題的關鍵在于如何選擇一個合適的閾值，較小的閾值能夠保留更多的圖像細節但同時一部分噪聲也會保留下來；較大的閾值能夠有效地去除噪聲但可能導致源圖像信息的丟失。為了解決這個問題，本文采用序列二次規劃算法為Contourlet域的圖像去噪選擇合適的閾值。

2.1 序列二次規劃算法

序列二次規劃方法是一種求解有約束的非線性優化問題的有效方法，起源于1963年Wilson提出的牛頓拉格朗日方法，后來經過Han與Powelld等的修改完善得到。序列二次規劃法從誕生至今經過了國內外學者非常深入的研究，是一種十分成熟的優化算法，對大多數等式和不等式約束優化問題都能取得較好的求解效果，非常適合解決中小型非線性優化問題，在許多領域有著非常廣泛的應用。

SQP算法主要思想是利用原優化問題的信息構造一個二次規劃子問題，通過在迭代過程中求解該子問題對當前解進行修正來優化性能指標[11]。二次規劃子問題的約束由原規劃問題的約束線性化得到，子問題的目標函數是拉格朗日方程的二次近似。

非線性規劃問題的拉格朗日函數如式(2)所示

(2)

式中，λi是Lagrange乘子，gi為等式與不等式約束條件。在每次迭代中，計算拉格朗日函數的Hessian矩陣Hk的近似。這樣，原問題就轉化為一個如式(3)～式(5)所示的二次規劃子問題。

(3)

s.t.ceqi(xk)+ceqi(xk)Tdk=0i∈E

(4)

ci(xk)+ci(xk)Tdk≥0i∈I

(5)

其中，xk是第k次迭代的解，dk表示第k次迭代時的搜索方向。之后通過求解二次規劃子問題根據式(6)調整當前的解。

xk+1=xk+αkdk

(6)

步長αk基于可行性優先準則確定，持續迭代過程即可使得目標函數不斷減小。即使算法選取的初值不滿足約束，通過迭代優化后依然可以得到滿足約束的解，如果評價函數或解的更新梯度小于設定的閾值則結束迭代。

2.2 廣義交叉驗證準則

本文利用廣義交叉驗證(General Cross Validation,GCV)準則建立序列二次規劃算法的目標函數模型。GCV準則提供了一種僅通過觀測到的含噪圖像獲取去噪閾值的方法，能夠在缺乏噪聲方差等先驗信息的條件下選擇合適的去噪閾值。

在對式(1)中模型進行非下采樣Contourlet變換之后，可以得到

(7)

(8)

其中，N為圖像子帶中Contourlet系數的總數目，N0代表進行閾值變換之后被置0的Contourlet系數的數目，ω、ωδ分別代表閾值變換前后的Contourlet系數矩陣。

圖像的均方差(Mean Square Error, MSE)定義為

(9)

GCV函數可以看作圖像均方差的估計，可以反映去噪后的圖像質量。Jason等已經證明，當N→∞時通過最小化GCV函數得到的閾值即等于使MSE(δ)最小的去噪閾值[12]。利用這個原理，就能夠在噪聲方差等先驗信息未知的條件下，通過最小化函數GCV(δ)為每個圖像子帶選擇合適的去噪閾值T

T=argmin(GCV(δ))

(10)

2.3 非線性閾值函數

在此前的研究中，硬閾值函數和軟閾值函數是最常用的閾值去噪函數[13]。硬閾值去噪函數如式(11)所示。

(11)

軟閾值去噪函數如式(12)所示：

可以看出，在硬閾值函數去噪中，小于閾值的系數被直接置0，而大于閾值的系數保留不變，這種不連續性會使得去噪圖像有較大的方差；另一方面，在軟閾值函數去噪中，大于閾值的系數整體減去了一個閾值的大小，使得去噪圖像與源圖像有較大的偏差[14]。與軟閾值和硬閾值函數相比，非線性閾值函數是連續的，經過非線性閾值函數處理后的均方誤差曲線比較光滑。本文將非線性閾值函數應用于NSCT圖像去噪問題中，非線性閾值函數如式(13)所示：

(13)

式中，Tm,n為圖像子帶Sm,n的閾值，權重因子αi,j定義為

(14)

2.4 本文圖像去噪方法流程

經過NSCT分解之后，在分解的每個子帶中以GCV函數作為優化目標函數，使用SQP優化算法為每個圖像子帶選擇合適的去噪閾值。在所有閾值確定之后，使用式(13)的非線性閾值函數處理圖像的NSCT系數。本文提出的基于SQP優化閾值的NSCT圖像去噪方法流程如下：

1)設定合適的NSP與NSDFB分解層數，對噪聲圖像進行多尺度NSCT分解；

2)設置SQP算法初始參數，計算每個圖像子帶中最大的Contourlet系數wmax，并將[0,wmax]作為閾值尋優的搜索區間；

3)基于GCV準則，使用SQP優化算法通過最小化函數GCV(δ)為每個圖像子帶選擇合適的去噪閾值；

4)根據以上確定的閾值，使用式(13)的非線性閾值函數處理圖像的NSCT系數；

5)對閾值處理后的系數進行NSCT反變換得到去噪圖像。

3 實驗結果與分析

基于本文提出的圖像去噪方法對一組疊加不同程度高斯白噪聲的標準測試圖像(Lena, Barbaba, Peppers)進行了實驗，這3組圖片包含較多的圖像細節信息，可以對去噪效果進行對比。噪聲標準差分別設為20、40、60，可以模擬一般圖像獲取和傳播的過程中產生的不同程度的高斯噪聲污染。

使用平滑指數(FI)評價圖像經過濾波器后的噪聲平滑能力，其計算公式如下

(15)

其中，M代表圖像濾波后某區域所有像素的平均值，SV代表所有像素的標準差。FI值越高，表示濾波器的平滑作用越強。

經過實驗表明隨著分解層的增加，平滑指數逐漸增大[15]，同時考慮到Contourlet系數數目越多GCV準則也越適用，實驗中NSP和NSDFB的分解層數設為[4 8 16]。實驗結果與標準CT和NSCT去噪方法進行了對比，使用峰值信噪比(Peak Signal to Noise Ratio, PSNR)對圖像去噪效果進行評價

(16)

使用邊緣保持指數(Edge Preservation Index,EPI)評價去噪后圖像對原始圖像水平和垂直方向邊緣的保持能力：

(17)

(18)

表1 不同去噪方法得到的PSNR結果

表2 不同去噪方法得到的EPI結果

表1列出了三種去噪算法對于疊加不同標準差高斯噪聲圖像得到的PSNR結果，表2所示為三種去噪算法得到的圖像邊緣保持指數。可以看出本文提出的算法在所有測試圖像中都得到了最好的PSNR結果，并且保留了最多的原始圖像邊緣信息。

圖4所示為不同去噪方法對于噪聲標準差40的Barbara圖像的去噪結果對比。圖4(a)和圖4(b)分別是源圖像和含噪圖像，圖4(c)和圖4(d)分別是標準CT和NSCT得到的去噪結果，PSNR值分別為25.150和26.945，圖4(e)為本文方法得到的圖像去噪結果，PSNR值為28.093。通過對比可以發現，本文圖像去噪方法相對于另外兩種去噪方法有著更好的去噪效果。

(a)源圖像 (b)含噪圖像， PSNR=22.203

(c)CT去噪圖像， (d)NSCT去噪圖像， PSNR=25.150 PSNR=26.945

(e)本文方法去噪圖像，PSNR=28.093圖4 不同方法得到的去噪圖像Fig.4 The denoising results of different methods

4 結論

本文提出了一種新的基于NSCT的圖像高斯白噪聲去除方法。首先通過SQP優化算法對GCV目標函數進行尋優，從而為每個NSCT圖像子帶選取合適的去噪閾值，進而使用非線性閾值函數對分解后的NSCT系數進行處理。該方法不需要圖像噪聲方差等先驗信息，能夠僅利用觀測噪聲圖像實現閾值的選取。實驗結果表明本文提出的方法能有效去除圖像的高斯白噪聲，與傳統基于Contourlet變換的圖像去噪方法相比能夠得到更高的峰值信噪比，并較好地保留原始圖像的細節特征。

參考文獻

[1] Wimmer G, Tamaki T, Tischendorf J J, et al. Directional wavelet based features for colonic polyp classification[J]. Medical Image Analysis, 2016, 31:16.

[2] Roux S G, Abry P, Vedel B, et al. Hyperbolic wavelet leaders for anisotropic multifractal texture analysis[C]//IEEE International Conference on Image Processing. IEEE, 2016:3558-3562.

[3] Shao L, Yan R, Li X, et al. From heuristic optimization to dictionary learning: a review and comprehensive comparison of image denoising algorithms[J]. IEEE Transactions on Cybernetics, 2014, 44(7):1001.

[4] Do M N, Vetterli M. The contourlet transform: an efficient directional multiresolution image representation[J]. IEEE Transactions on Image Processing, 2005, 14(12):2091.

[5] Yin M, Liu W, Zhao X, et al. Image denoising using trivariate prior model in nonsubsampled dual-tree complex contourlet transform domain and non-local means filter in spatial domain[J]. Optik -International Journal for Light and Electron Optics, 2013, 124(24):6896-6904.

[6] Cunha A L D, Zhou J, Do M N. The nonsubsampled contourlet transform: theory, design, and applica-tions[J]. IEEE Transactions on Image Processing, 2006, 15(10):3089-3101.

[7] Wang X Y, Yang H Y, Zhang Y, et al. Image denoising using SVM classification in nonsubsampled contourlet transform domain[J]. Information Scie-nces An International Journal, 2013, 246(14):155-176.

[8] Han J, Mirko V D B. Microseismic and seismic denoising via ensemble empirical mode decomposition and adaptive thresholding[J]. Geophysics, 2015, 80(6):KS69-KS80.

[9] Chai Y, Li H, Zhang X. Multifocus image fusion based on features contrast of multiscale products in non-subsampled contourlet transform domain[J]. Optik-International Journal for Light and Electron Optics, 2012, 123(7):569-581.

[10] Luisier F, Blu T, Unser M. Image denoising in mixed poisson-Gaussian noise[M]. IEEE Press, 2011.

[11] Morshed M J, Asgharpour A. Hybrid imperialist competitive-sequential quadratic programming (HIC-SQP) algorithm for solving economic load dispatch with incorporating stochastic wind power: A comparative study on heuristic optimization techniques[J]. Energy Conversion & Management, 2014, 84:30-40.

[12] Jansen M, Bultheel A. Multiple wavelet threshold estimation by generalized cross validation for images with correlated noise[J]. IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society, 1999, 8(7):947-953.

[13] Gnanadurai D, Sadasivam V. An efficient adaptive thresholding technique for wavelet based image denoising[J]. International Journal of Signal Proces-sing, 2006(2).

[14] Chui M, Feng Y, Wang W, et al. Image denoising method with adaptive Bayes threshold in Nonsubsampled Contourlet Domain[J]. Aasri Procedia, 2012, 1(3):512-518.

[15] 楊曉慧, 焦李成. 非下采樣Contourlet域GCV準則SAR圖像去噪[J]. 計算機應用研究, 2009, 26(9):3542-3544.