基于序列二次規(guī)劃優(yōu)化閾值的NSCT高斯噪聲圖像濾波方法

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

(北京航空航天大學(xué) 自動(dòng)化科學(xué)與電氣工程學(xué)院，飛行器控制一體化國(guó)防科技重點(diǎn)實(shí)驗(yàn)室，北京 100191)

0 引言

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

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

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

1 非下采樣Contourlet變換

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

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

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

1.1 非下采樣金字塔分解

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

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

1.2 非下采樣方向?yàn)V波器組

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

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

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

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

y=x+η

(1)

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

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

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

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

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

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

(2)

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

(3)

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

(4)

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

(5)

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

xk+1=xk+αkdk

(6)

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

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

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

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

T=argmin(GCV(δ))

(10)

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

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

(11)

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

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

(13)

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

(14)

2.4 本文圖像去噪方法流程

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

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

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

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

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

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

3 實(shí)驗(yàn)結(jié)果與分析

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

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

(15)

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

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

(16)

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

(17)

(18)

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

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

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

圖4所示為不同去噪方法對(duì)于噪聲標(biāo)準(zhǔn)差40的Barbara圖像的去噪結(jié)果對(duì)比。圖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。通過(guò)對(duì)比可以發(fā)現(xiàn)，本文圖像去噪方法相對(duì)于另外兩種去噪方法有著更好的去噪效果。

(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的圖像高斯白噪聲去除方法。首先通過(guò)SQP優(yōu)化算法對(duì)GCV目標(biāo)函數(shù)進(jìn)行尋優(yōu)，從而為每個(gè)NSCT圖像子帶選取合適的去噪閾值，進(jìn)而使用非線性閾值函數(shù)對(duì)分解后的NSCT系數(shù)進(jìn)行處理。該方法不需要圖像噪聲方差等先驗(yàn)信息，能夠僅利用觀測(cè)噪聲圖像實(shí)現(xiàn)閾值的選取。實(shí)驗(yà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]. 計(jì)算機(jī)應(yīng)用研究, 2009, 26(9):3542-3544.