基于自適應(yīng)雙邊濾波的太赫茲圖像去噪算法

郭俊文，宋貴才，李興廣

（1.長春理工大學(xué) 理學(xué)院，長春 130022；2.長春理工大學(xué) 電子信息工程學(xué)院，長春 130022）

太赫茲時域光譜技術(shù)由于具有良好的時間分辨率和寬光譜覆蓋范圍，已成為太赫茲光譜、成像和介質(zhì)無損檢測的主要手段［1-3］。由于激光器的抖動和介質(zhì)的不均勻?qū)е路瓷鋻呙鑸D像中背景灰度分布不均勻同時容易受到噪聲的污染，導(dǎo)致圖像邊緣分辨率較差、圖像質(zhì)量不高［4-6］。在提高圖像質(zhì)量方面有大量的文獻(xiàn)［7-9］，其研究主要集中在信噪比的提高，而在有效去除噪聲的同時保留更多的圖像細(xì)節(jié)信息則是太赫茲圖像濾波算法的關(guān)鍵。

傳統(tǒng)的雙邊濾波（BF）［10］，人為修改參數(shù)達(dá)到濾波效果，但人為設(shè)置參數(shù)往往取值不合理，以致出現(xiàn)不能完全濾除噪聲或者圖像模糊的現(xiàn)象，即使參數(shù)設(shè)置合理，在一幅圖像中的不同區(qū)域也會存在降噪和保留邊緣特征之間的矛盾，文獻(xiàn)［11］提出了雙邊濾波算法在太赫茲共焦掃描圖像的應(yīng)用，雖然取得了良好的效果，但在設(shè)置相關(guān)參數(shù)時仍取經(jīng)驗(yàn)值，無法實(shí)現(xiàn)自適應(yīng)太赫茲圖像的降噪。文獻(xiàn)［12］針對太赫茲圖像的分辨率差，噪聲嚴(yán)重的問題，提出了一種自適應(yīng)流行高維濾波的降噪算法，但在降噪過程中圖像中的殘留噪聲和背景信息也被增強(qiáng)，導(dǎo)致圖像的模糊。文獻(xiàn)［13］P.V.Sudeep中提出了CT圖像中的自適應(yīng)雙邊濾波，通過極大似然噪聲評估的方法對多幀CT圖像進(jìn)行噪聲評估，合理選擇濾波參數(shù)。

本文針對太赫茲反射圖像中的高斯白噪聲，較好的平衡了降噪和細(xì)節(jié)信息間的關(guān)系問題。對傳統(tǒng)雙邊濾波進(jìn)行分析確定影響降噪性能的主要參數(shù)指標(biāo)，明確雙邊濾波的濾波效果主要取決于亮度標(biāo)準(zhǔn)差，建立亮度標(biāo)準(zhǔn)差和圖像噪聲方差之間的線性關(guān)系。通過主成分分析（PCA）的方法迭代估計太赫茲圖像中高斯白噪聲的噪聲方差，由估計的噪聲方差確定亮度標(biāo)準(zhǔn)差進(jìn)而對太赫茲噪聲圖像進(jìn)行降噪處理。

1 自適應(yīng)雙邊濾波算法

傳統(tǒng)的雙邊濾波是Tomasi和Manduchi于1998年提出的非線性濾波，雙邊濾波的優(yōu)勢在于平滑噪聲的同時能夠很好地保護(hù)圖像的邊緣細(xì)節(jié)［10］。雙邊濾波定義如下：

式中，(m)為像素點(diǎn)m濾波后的值；I(n)為輸入的像素值；N(m)為m附近的方形領(lǐng)域；Z是歸一化因子；n是鄰域中像素點(diǎn)的位置；wD是空間函數(shù)；wR是亮度函數(shù)，定義如下：

式中，σr為亮度標(biāo)準(zhǔn)差，σd為空間標(biāo)準(zhǔn)差，Im和In分別為m和n處的像素值。

由文獻(xiàn)［13］可知，推導(dǎo)自適應(yīng)雙邊濾波器的第一步是雙邊濾波器轉(zhuǎn)換為自適應(yīng)雙邊濾波器時參數(shù)的自適應(yīng)選擇，雖然兩個參數(shù)決定了濾波器的性能，但是空間標(biāo)準(zhǔn)差對噪聲不敏感所以需要自適應(yīng)設(shè)置亮度標(biāo)準(zhǔn)差減少噪聲，定義σr如下：

式中，?2表示從輸入圖像估計的噪聲方差，τ為固定常數(shù)，由文獻(xiàn)［13］可知σd在為1.8時，信噪比隨著σr的增加而增大，意味著圖像的去噪能力在增強(qiáng)。因此亮度標(biāo)準(zhǔn)差可作為對圖像處理算法效果優(yōu)劣的評價標(biāo)準(zhǔn)，其結(jié)果主要由?2的取值決定。

1.1 亮度標(biāo)準(zhǔn)差自適應(yīng)計算

由式（5）可知，噪聲方差決定了自適應(yīng)灰度標(biāo)準(zhǔn)差，噪聲方差估計算法流程圖如圖1所示，假設(shè)圖像和噪聲是不相關(guān)的，x代表大小為S1×S2的無噪圖像，y=x+n是均值為零的與信號無關(guān)的高斯白噪聲的圖像。對圖像提取子集塊，子集塊大小為M1×M2，共包含的子集塊，對子集塊重新排列成有M=M1M2元素的向量，分別被xi,ni,yi,i=1...N認(rèn)為是隨機(jī)向量的實(shí)現(xiàn)，n是加性高斯白噪聲且服從標(biāo)準(zhǔn)正態(tài)分布同時X和Y相互獨(dú)立N～ΝM(0,σ2I)行間距，cov(X,Ν)=0。利用文獻(xiàn)［14］中計算噪聲方差最大值，同時最大值也被用作迭代過程中的初始值，從而細(xì)化了噪聲方差的估計。噪聲方差可以估計為圖像塊協(xié)方差矩陣的最小特征值。因此，在獲取下一代噪聲方差時采用主成分分析（PCA）的方法計算樣本協(xié)方差矩陣的最小特征值。

圖1 噪聲方差估計算法完整流程圖

1.2 主成分分析方法

PCA（Principal Component Analysis）是一種可用于噪聲方差估計的方法［15，16］，證明了噪聲方差可以被估計為圖像塊協(xié)方差矩陣的最小特征值。本文首先對圖像提取子集塊，集塊大小為M1×M2。SX和SY分別為X和Y的樣本協(xié)方差矩陣，有以下性質(zhì)：

式中，s2表示樣本方差，表示特征值，表示特征向量。隨機(jī)向量X的值只存在子空間VM-m，這意味著X分量之間存在線性相關(guān)性，即：

當(dāng)考慮主成分時，協(xié)方差cov(X,Ν)=0可知ΣY=ΣX+ΣN，ΣX，ΣN和ΣY分別為X，N和Y的樣本協(xié)方差矩陣。此外，ΣN=σ2I和ΣX的最小特征值為零。因此，ΣY的最小特征值等于σ2，即噪聲方差可以估計為圖像塊協(xié)方差矩陣的最小特征值：

因此可知收斂于σ2，那么噪聲方差能夠被估計。在此基礎(chǔ)上可以利用式（9）計算差的期望值：

進(jìn)而在迭代計算噪聲方差過程中需要滿足以下條件：

其中，根據(jù)實(shí)驗(yàn)過程，對算法的部分參數(shù)取值為：M1=M2=5，C0=2.9，m=7，p=0.06，最大迭代次數(shù)imax=10，N代表對M圖像塊提取滿足式（10）的有效圖像塊。最終，采用迭代運(yùn)算來計算圖像噪聲方差，具體步驟如下：

（1）以圖像的每個像素點(diǎn)為左頂點(diǎn)，可劃分出多個5×5的圖像塊；

（2）利用初始噪聲圖像，計算其噪聲方差的上邊界

（3）估計下一代噪聲方差滿足同時滿足輸出否則繼續(xù)步驟（3）；

（5）停止運(yùn)算，輸出

（6）令，轉(zhuǎn)至步驟（2）。

實(shí)驗(yàn)表明，經(jīng)過3次迭代運(yùn)算即可估計出噪聲方差，由式（5）圖像的噪聲方差和亮度標(biāo)準(zhǔn)差為線性關(guān)系，因此可以自適應(yīng)亮度標(biāo)準(zhǔn)差。

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

本文所有的降噪方法均在CPU為Core i7，3.30GHZ，64GRAM，64位的WIN7操作系統(tǒng)，Matlab R2016a中進(jìn)行測試。本實(shí)驗(yàn)中，首先對仿真Lena圖像和太赫茲圖像添加不同等級的高斯白噪聲，對比實(shí)際噪聲方差和估計的噪聲方差，驗(yàn)證噪聲方差估計的準(zhǔn)確性，由下圖2給出具體的實(shí)驗(yàn)數(shù)據(jù)。

圖2 噪聲方差估計結(jié)果對比圖，lena和thz image分別添加不同等級的噪聲方差，估計的噪聲方差分布在實(shí)際噪聲方差附近

在去噪實(shí)驗(yàn)中，大多數(shù)情況下較高的噪聲方差估計精度會決定更高的去噪質(zhì)量。因此本文算法在自適應(yīng)參數(shù)上采用主成分分析的方法對圖像的整體噪聲進(jìn)行估計，從圖2可以看出估計的噪聲方差和實(shí)際噪聲方差之間接近。

為了驗(yàn)證本文算法的改進(jìn)之處，本文以傳統(tǒng)雙邊濾波、三維塊匹配濾波和引導(dǎo)濾波為比較對象。選用Lena和太赫茲樹葉兩幅大小為256x256的灰度圖像進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)步驟為：首先為圖像加入不同水平的高斯白噪聲得到噪聲圖像；利用PCA噪聲方差估計的方法來確定噪聲方差；其次根據(jù)噪聲方差計算噪聲標(biāo)準(zhǔn)差，利用整幅圖像的噪聲方差估算出亮度標(biāo)準(zhǔn)差；最后根據(jù)亮度標(biāo)準(zhǔn)差對整幅圖像進(jìn)行濾波。將實(shí)驗(yàn)結(jié)果與三種參照算法作比較，為了客觀評價算法性能，分別對仿真圖像和實(shí)驗(yàn)圖像進(jìn)行實(shí)驗(yàn)采用峰值信噪比（Peak Signal to Noise Ratio，PSNR）和結(jié)構(gòu)相似度（Structural Similarity Index Measurement，SSIM）（圖像的輪廓和細(xì)節(jié)）對圖像質(zhì)量進(jìn)行評價。實(shí)驗(yàn)結(jié)果如下圖所示。

式中，MSE表示處理后圖像I(i,j)與原圖像I0(i,j)之間的均方誤差，M和N表示圖像的大小。在一般情況下，PSNR值越高表示濾波算法的去噪效果越好。

式中，μf和分別為f的均值和f的方差；σf和σfg分別為圖像f的標(biāo)準(zhǔn)差以及圖像的協(xié)方差；C1和C2為小的常數(shù)，避免分母為零出現(xiàn)不穩(wěn)定現(xiàn)象。

圖3（a）是原始圖像，圖3（b）是添加方差為15的高斯白噪聲的Lena圖像，從圖中可以看出，四種算法在一定程度上去除噪聲。圖3（e）是雙邊濾波后的圖像，傳統(tǒng)雙邊濾波由于參數(shù)固定，處理后的圖像噪聲分布密集，效果不理想。圖3（d）是引導(dǎo)濾波后的圖像，引導(dǎo)濾波在設(shè)置參數(shù)時選擇了最優(yōu)值，比雙邊濾波可以更有效的去除噪聲，同時也出現(xiàn)了過平滑部分圖像細(xì)節(jié)被平滑。圖3（c）是三維塊匹配濾波后的圖像，三維塊匹配通過在三維空間進(jìn)行濾波處理，再將結(jié)果反變換融合到二維，該算法去噪效果顯著但是時間復(fù)雜度較高。圖3（f）是本文算法濾波后的圖像，本文算法在去除噪聲的同時沒有出現(xiàn)過度平滑，也盡可能的保留了圖像的細(xì)節(jié)信息。

圖3 均值為0，方差為15噪聲圖像濾波效果

圖4 太赫茲圖像濾波效果

圖4（a）是太赫茲反射掃描圖像，圖4（b-e）分別為三維塊濾波、引導(dǎo)濾波、雙邊濾波和本文算法去除噪聲后的圖像。從圖像上可以看出，圖4（c）和圖4（d）為經(jīng)引導(dǎo)濾波和雙邊濾波處理，出現(xiàn)了邊緣平滑問題，同時仍然殘留噪聲，圖像表現(xiàn)較為模糊，圖4（b）和圖4（e）分別為經(jīng)三維塊匹配和本文算法處理后的圖像，在達(dá)到相同處理效果上，本文算法在結(jié)構(gòu)相似度上要優(yōu)于三維塊匹配算法。因此，本文算法在時間復(fù)雜度、峰值信噪比和結(jié)構(gòu)相似度均優(yōu)參照算法，同時圖像細(xì)節(jié)得到更好的保留，為圖像的檢測和識別提供了高質(zhì)量的太赫茲圖像。

由表1可知，通過與參照算法的對比實(shí)驗(yàn)，可以得出以下結(jié)論：提出的自適應(yīng)雙邊濾波算法在整體的濾波效果上有明顯的優(yōu)勢，在峰值信噪比和結(jié)構(gòu)相似度方面有明顯提高，這說明既保留了圖像的細(xì)節(jié)信息也達(dá)到了濾波效果。該算法也只是針對與太赫茲圖像中的高斯噪聲分布，需要改進(jìn)的地方是對于乘性噪聲和泊松噪聲具有一定的局限性，需要對噪聲模型做進(jìn)一步改進(jìn)。

表1 太赫茲圖像濾波結(jié)果PSNR和SSIM對比

3 結(jié)論

綜上所述，本文提出了基于PCA的自適應(yīng)雙邊濾波算法，依據(jù)圖像塊的高斯噪聲分布，對不同等級的圖像噪聲方差算法自行進(jìn)行準(zhǔn)確的估計，避免噪聲方差參數(shù)的不平滑，最優(yōu)的選取了空間標(biāo)準(zhǔn)差和亮度標(biāo)準(zhǔn)差，滿足了去噪后圖像細(xì)節(jié)特征信息保留的情況下，圖像的峰值信噪比相比于雙邊濾波提高1.5963dB，圖像的結(jié)構(gòu)相似度相比于雙邊濾波提高0.06?？梢越鉀Q在太赫茲反射掃描成像過程中，傳統(tǒng)的去噪處理算法帶來的細(xì)節(jié)信息丟失問題。實(shí)驗(yàn)結(jié)果表明，本文算法對太赫茲反射掃描圖像有良好的降噪效果，結(jié)果很好地保留了圖像的邊緣細(xì)節(jié)特征，峰值信噪比也得到了提高。

［1］Shen Y C.Terahertz pulsed spectroscopy and imaging for pharmaceutical applications：a review［J］.International Journal of Pharmaceutics，2011，417（1）：48-60.

［2］Secomartorell C，Lópezdomínguez V，Arauzgarofalo G，etal.Goya'sartwork imaging with Terahertz waves. ［J］.Optics Express， 2013， 21 （15） ：17800-17805.

［3］Redo-Sanchez A，Laman N，Schulkin B，et al.Review ofTerahertzTechnologyReadinessAssessment and Applications［J］.Journal of Infrared Millimeter&Terahertz Waves，2013，34（9）：500-518.

［4］喬曉利，任姣姣，張丹丹，等.太赫茲波譜多特征參數(shù)成像方法［J］.長春理工大學(xué)學(xué)報：自然科學(xué)版，2017，40（3）：25-31.

［5］陳思宏，韓旭，李麗娟，等.基于太赫茲時域光譜技術(shù)的塑料檢測［J］.長春理工大學(xué)學(xué)報：自然科學(xué)版，2017，40（5）：14-16.

［6］徐利民.高分辨太赫茲圖像處理［D］.西安：中國科學(xué)院大學(xué)，2013.

［7］Chen Y，Huang S，Pickwellmacpherson E.Frequency-wavelet domain deconvolution for terahertz reflection imaging and spectroscopy［J］.Optics Express，2010，18（2）：1177.

［8］Galv?o R，Hadjiloucas S，Bowen J，et al.Optimal discrimination and classification of THz spectra in thewaveletdomain［J］.OpticsExpress，2003，11（12）：1462.

［9］Heshmat B，Masnadi-Shirazi M，Lewis R B，et al.Enhanced TerahertzBandwidth and Powerfrom GaAsBi-based Sources［J］.Advanced Optical Materials，2013，1（10）：714-719.

［10］Tomasi C，Manduchi R.Bilateral filtering for gray and colorimages［C］.InternationalConference on Computer Vision.IEEE Computer Society，1998.

［11］楊永發(fā)，李琦.雙邊濾波算法的太赫茲共焦掃描圖像去噪應(yīng)用［J］.激光與光電子學(xué)進(jìn)展，2015，52（12）：93-98.

［12］朱書進(jìn)，李躍華.基于自適應(yīng)流形高維濾波的太赫茲圖像濾波算法［J］.太赫茲科學(xué)與電子信息學(xué)報，2016，14（4）：498-501.

［13］Sudeep P V，Niwas S I，Palanisamy P，et al.Enhancement and bias removal of optical coherence tomography images：An iterative approach with adaptive bilateral filtering［J］.Computers in Biology&Medicine，2016，71（C）：97-107.

［14］Corner B R，Narayanan R M，Reichenbach S E.Noise estimation in remote sensing imagery using data masking［J］.International Journal of Remote Sensing，2003，24（4）：689-702.

［15］Pyatykh S，Hesser J，Zheng L.Image noise level estimation by principalcomponent analysis.［J］.IEEE Transactions on Image Processing A Publication ofthe IEEE SignalProcessing Society，2013，22（2）：687.

［16］I.Jolliffe，“Principal component analysis，”in Encyclopedia of Statistics in Behavioral Science.New York：Springer-Verlag，2002：416-509.