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

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

0 引言

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

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

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

1 非下采樣Contourlet變換

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

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

圖1 NSCT分解的結(jié)構(gòu)示意圖Fig.1 The structure of NSCT decomposition

1.1 非下采樣金字塔分解

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

圖2 三級NSP分解結(jié)構(gòu)示意圖Fig.2 The three-stage NSP decomposition

1.2 非下采樣方向濾波器組

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

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

2 基于SQP優(yōu)化閾值的NSCT圖像去噪

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

y=x+η

(1)

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

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

2.1 序列二次規(guī)劃算法

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

SQP算法主要思想是利用原優(yōu)化問題的信息構(gòu)造一個二次規(guī)劃子問題，通過在迭代過程中求解該子問題對當(dāng)前解進(jìn)行修正來優(yōu)化性能指標(biāo)[11]。二次規(guī)劃子問題的約束由原規(guī)劃問題的約束線性化得到，子問題的目標(biāo)函數(shù)是拉格朗日方程的二次近似。

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

(2)

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

(3)

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

(4)

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

(5)

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

xk+1=xk+αkdk

(6)

步長αk基于可行性優(yōu)先準(zhǔn)則確定，持續(xù)迭代過程即可使得目標(biāo)函數(shù)不斷減小。即使算法選取的初值不滿足約束，通過迭代優(yōu)化后依然可以得到滿足約束的解，如果評價函數(shù)或解的更新梯度小于設(shè)定的閾值則結(jié)束迭代。

2.2 廣義交叉驗證準(zhǔn)則

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

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

(7)

(8)

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

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

(9)

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

T=argmin(GCV(δ))

(10)

2.3 非線性閾值函數(shù)

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

(11)

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

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

(13)

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

(14)

2.4 本文圖像去噪方法流程

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

1)設(shè)定合適的NSP與NSDFB分解層數(shù)，對噪聲圖像進(jìn)行多尺度NSCT分解；

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

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

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

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

3 實驗結(jié)果與分析

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

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

(15)

其中，M代表圖像濾波后某區(qū)域所有像素的平均值，SV代表所有像素的標(biāo)準(zhǔn)差。FI值越高，表示濾波器的平滑作用越強。

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

(16)

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

(17)

(18)

表1 不同去噪方法得到的PSNR結(jié)果

表2 不同去噪方法得到的EPI結(jié)果

表1列出了三種去噪算法對于疊加不同標(biāo)準(zhǔn)差高斯噪聲圖像得到的PSNR結(jié)果，表2所示為三種去噪算法得到的圖像邊緣保持指數(shù)?？梢钥闯霰疚奶岢龅乃惴ㄔ谒袦y試圖像中都得到了最好的PSNR結(jié)果，并且保留了最多的原始圖像邊緣信息。

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

(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 結(jié)論

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

參考文獻(xiàn)

[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準(zhǔn)則SAR圖像去噪[J]. 計算機應(yīng)用研究, 2009, 26(9):3542-3544.