郭俊文，宋貴才，李興廣

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

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

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

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

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

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

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

式中，σr為亮度標準差，σd為空間標準差，Im和In分別為m和n處的像素值。

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

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

1.1 亮度標準差自適應(yīng)計算

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

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

1.2 主成分分析方法

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

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

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

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

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

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

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

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

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

（5）停止運算，輸出

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

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

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

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

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

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

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

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

式中，μf和分別為f的均值和f的方差；σf和σfg分別為圖像f的標準差以及圖像的協(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)了過平滑部分圖像細節(jié)被平滑。圖3（c）是三維塊匹配濾波后的圖像，三維塊匹配通過在三維空間進行濾波處理，再將結(jié)果反變換融合到二維，該算法去噪效果顯著但是時間復(fù)雜度較高。圖3（f）是本文算法濾波后的圖像，本文算法在去除噪聲的同時沒有出現(xiàn)過度平滑，也盡可能的保留了圖像的細節(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)三維塊匹配和本文算法處理后的圖像，在達到相同處理效果上，本文算法在結(jié)構(gòu)相似度上要優(yōu)于三維塊匹配算法。因此，本文算法在時間復(fù)雜度、峰值信噪比和結(jié)構(gòu)相似度均優(yōu)參照算法，同時圖像細節(jié)得到更好的保留，為圖像的檢測和識別提供了高質(zhì)量的太赫茲圖像。

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

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

3 結(jié)論

綜上所述，本文提出了基于PCA的自適應(yīng)雙邊濾波算法，依據(jù)圖像塊的高斯噪聲分布，對不同等級的圖像噪聲方差算法自行進行準確的估計，避免噪聲方差參數(shù)的不平滑，最優(yōu)的選取了空間標準差和亮度標準差，滿足了去噪后圖像細節(jié)特征信息保留的情況下，圖像的峰值信噪比相比于雙邊濾波提高1.5963dB，圖像的結(jié)構(gòu)相似度相比于雙邊濾波提高0.06?？梢越鉀Q在太赫茲反射掃描成像過程中，傳統(tǒng)的去噪處理算法帶來的細節(jié)信息丟失問題。實驗結(jié)果表明，本文算法對太赫茲反射掃描圖像有良好的降噪效果，結(jié)果很好地保留了圖像的邊緣細節(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］.長春理工大學學報：自然科學版，2017，40（3）：25-31.

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

［6］徐利民.高分辨太赫茲圖像處理［D］.西安：中國科學院大學，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］.激光與光電子學進展，2015，52（12）：93-98.

［12］朱書進，李躍華.基于自適應(yīng)流形高維濾波的太赫茲圖像濾波算法［J］.太赫茲科學與電子信息學報，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.