基于自适应双边滤波的太赫兹图像去噪算法

郭俊文，宋贵才，李兴广

（1.长春理工大学 理学院，长春 130022；2.长春理工大学 电子信息工程学院，长春 130022）

太赫兹时域光谱技术由于具有良好的时间分辨率和宽光谱覆盖范围，已成为太赫兹光谱、成像和介质无损检测的主要手段［1-3］。由于激光器的抖动和介质的不均匀导致反射扫描图像中背景灰度分布不均匀同时容易受到噪声的污染，导致图像边缘分辨率较差、图像质量不高［4-6］。在提高图像质量方面有大量的文献［7-9］，其研究主要集中在信噪比的提高，而在有效去除噪声的同时保留更多的图像细节信息则是太赫兹图像滤波算法的关键。

传统的双边滤波（BF）［10］，人为修改参数达到滤波效果，但人为设置参数往往取值不合理，以致出现不能完全滤除噪声或者图像模糊的现象，即使参数设置合理，在一幅图像中的不同区域也会存在降噪和保留边缘特征之间的矛盾，文献［11］提出了双边滤波算法在太赫兹共焦扫描图像的应用，虽然取得了良好的效果，但在设置相关参数时仍取经验值，无法实现自适应太赫兹图像的降噪。文献［12］针对太赫兹图像的分辨率差，噪声严重的问题，提出了一种自适应流行高维滤波的降噪算法，但在降噪过程中图像中的残留噪声和背景信息也被增强，导致图像的模糊。文献［13］P.V.Sudeep中提出了CT图像中的自适应双边滤波，通过极大似然噪声评估的方法对多帧CT图像进行噪声评估，合理选择滤波参数。

本文针对太赫兹反射图像中的高斯白噪声，较好的平衡了降噪和细节信息间的关系问题。对传统双边滤波进行分析确定影响降噪性能的主要参数指标，明确双边滤波的滤波效果主要取决于亮度标准差，建立亮度标准差和图像噪声方差之间的线性关系。通过主成分分析（PCA）的方法迭代估计太赫兹图像中高斯白噪声的噪声方差，由估计的噪声方差确定亮度标准差进而对太赫兹噪声图像进行降噪处理。

1 自适应双边滤波算法

传统的双边滤波是Tomasi和Manduchi于1998年提出的非线性滤波，双边滤波的优势在于平滑噪声的同时能够很好地保护图像的边缘细节［10］。双边滤波定义如下：

式中，(m)为像素点m滤波后的值；I(n)为输入的像素值；N(m)为m附近的方形领域；Z是归一化因子；n是邻域中像素点的位置；wD是空间函数；wR是亮度函数，定义如下：

式中，σr为亮度标准差，σd为空间标准差，Im和In分别为m和n处的像素值。

由文献［13］可知，推导自适应双边滤波器的第一步是双边滤波器转换为自适应双边滤波器时参数的自适应选择，虽然两个参数决定了滤波器的性能，但是空间标准差对噪声不敏感所以需要自适应设置亮度标准差减少噪声，定义σr如下：

式中，̂2表示从输入图像估计的噪声方差，τ为固定常数，由文献［13］可知σd在为1.8时，信噪比随着σr的增加而增大，意味着图像的去噪能力在增强。因此亮度标准差可作为对图像处理算法效果优劣的评价标准，其结果主要由̂2的取值决定。

1.1 亮度标准差自适应计算

由式（5）可知，噪声方差决定了自适应灰度标准差，噪声方差估计算法流程图如图1所示，假设图像和噪声是不相关的，x代表大小为S1×S2的无噪图像，y=x+n是均值为零的与信号无关的高斯白噪声的图像。对图像提取子集块，子集块大小为M1×M2，共包含的子集块，对子集块重新排列成有M=M1M2元素的向量，分别被xi,ni,yi,i=1...N认为是随机向量的实现，n是加性高斯白噪声且服从标准正态分布同时X和Y相互独立N～ΝM(0,σ2I)行间距，cov(X,Ν)=0。利用文献［14］中计算噪声方差最大值，同时最大值也被用作迭代过程中的初始值，从而细化了噪声方差的估计。噪声方差可以估计为图像块协方差矩阵的最小特征值。因此，在获取下一代噪声方差时采用主成分分析（PCA）的方法计算样本协方差矩阵的最小特征值。

图1 噪声方差估计算法完整流程图

1.2 主成分分析方法

PCA（Principal Component Analysis）是一种可用于噪声方差估计的方法［15，16］，证明了噪声方差可以被估计为图像块协方差矩阵的最小特征值。本文首先对图像提取子集块，集块大小为M1×M2。SX和SY分别为X和Y的样本协方差矩阵，有以下性质：

式中，s2表示样本方差，表示特征值，表示特征向量。随机向量X的值只存在子空间VM-m，这意味着X分量之间存在线性相关性，即：

当考虑主成分时，协方差cov(X,Ν)=0可知ΣY=ΣX+ΣN，ΣX，ΣN和ΣY分别为X，N和Y的样本协方差矩阵。此外，ΣN=σ2I和ΣX的最小特征值为零。因此，ΣY的最小特征值等于σ2，即噪声方差可以估计为图像块协方差矩阵的最小特征值：

因此可知收敛于σ2，那么噪声方差能够被估计。在此基础上可以利用式（9）计算差的期望值：

进而在迭代计算噪声方差过程中需要满足以下条件：

其中，根据实验过程，对算法的部分参数取值为：M1=M2=5，C0=2.9，m=7，p=0.06，最大迭代次数imax=10，N代表对M图像块提取满足式（10）的有效图像块。最终，采用迭代运算来计算图像噪声方差，具体步骤如下：

（1）以图像的每个像素点为左顶点，可划分出多个5×5的图像块；

（2）利用初始噪声图像，计算其噪声方差的上边界

（3）估计下一代噪声方差满足同时满足输出否则继续步骤（3）；

（5）停止运算，输出

（6）令，转至步骤（2）。

实验表明，经过3次迭代运算即可估计出噪声方差，由式（5）图像的噪声方差和亮度标准差为线性关系，因此可以自适应亮度标准差。

2 实验结果与分析

本文所有的降噪方法均在CPU为Core i7，3.30GHZ，64GRAM，64位的WIN7操作系统，Matlab R2016a中进行测试。本实验中，首先对仿真Lena图像和太赫兹图像添加不同等级的高斯白噪声，对比实际噪声方差和估计的噪声方差，验证噪声方差估计的准确性，由下图2给出具体的实验数据。

图2 噪声方差估计结果对比图，lena和thz image分别添加不同等级的噪声方差，估计的噪声方差分布在实际噪声方差附近

在去噪实验中，大多数情况下较高的噪声方差估计精度会决定更高的去噪质量。因此本文算法在自适应参数上采用主成分分析的方法对图像的整体噪声进行估计，从图2可以看出估计的噪声方差和实际噪声方差之间接近。

为了验证本文算法的改进之处，本文以传统双边滤波、三维块匹配滤波和引导滤波为比较对象。选用Lena和太赫兹树叶两幅大小为256x256的灰度图像进行实验。实验步骤为：首先为图像加入不同水平的高斯白噪声得到噪声图像；利用PCA噪声方差估计的方法来确定噪声方差；其次根据噪声方差计算噪声标准差，利用整幅图像的噪声方差估算出亮度标准差；最后根据亮度标准差对整幅图像进行滤波。将实验结果与三种参照算法作比较，为了客观评价算法性能，分别对仿真图像和实验图像进行实验采用峰值信噪比（Peak Signal to Noise Ratio，PSNR）和结构相似度（Structural Similarity Index Measurement，SSIM）（图像的轮廓和细节）对图像质量进行评价。实验结果如下图所示。

式中，MSE表示处理后图像I(i,j)与原图像I0(i,j)之间的均方误差，M和N表示图像的大小。在一般情况下，PSNR值越高表示滤波算法的去噪效果越好。

式中，μf和分别为f的均值和f的方差；σf和σfg分别为图像f的标准差以及图像的协方差；C1和C2为小的常数，避免分母为零出现不稳定现象。

图3（a）是原始图像，图3（b）是添加方差为15的高斯白噪声的Lena图像，从图中可以看出，四种算法在一定程度上去除噪声。图3（e）是双边滤波后的图像，传统双边滤波由于参数固定，处理后的图像噪声分布密集，效果不理想。图3（d）是引导滤波后的图像，引导滤波在设置参数时选择了最优值，比双边滤波可以更有效的去除噪声，同时也出现了过平滑部分图像细节被平滑。图3（c）是三维块匹配滤波后的图像，三维块匹配通过在三维空间进行滤波处理，再将结果反变换融合到二维，该算法去噪效果显著但是时间复杂度较高。图3（f）是本文算法滤波后的图像，本文算法在去除噪声的同时没有出现过度平滑，也尽可能的保留了图像的细节信息。

图3 均值为0，方差为15噪声图像滤波效果

图4 太赫兹图像滤波效果

图4（a）是太赫兹反射扫描图像，图4（b-e）分别为三维块滤波、引导滤波、双边滤波和本文算法去除噪声后的图像。从图像上可以看出，图4（c）和图4（d）为经引导滤波和双边滤波处理，出现了边缘平滑问题，同时仍然残留噪声，图像表现较为模糊，图4（b）和图4（e）分别为经三维块匹配和本文算法处理后的图像，在达到相同处理效果上，本文算法在结构相似度上要优于三维块匹配算法。因此，本文算法在时间复杂度、峰值信噪比和结构相似度均优参照算法，同时图像细节得到更好的保留，为图像的检测和识别提供了高质量的太赫兹图像。

由表1可知，通过与参照算法的对比实验，可以得出以下结论：提出的自适应双边滤波算法在整体的滤波效果上有明显的优势，在峰值信噪比和结构相似度方面有明显提高，这说明既保留了图像的细节信息也达到了滤波效果。该算法也只是针对与太赫兹图像中的高斯噪声分布，需要改进的地方是对于乘性噪声和泊松噪声具有一定的局限性，需要对噪声模型做进一步改进。

表1 太赫兹图像滤波结果PSNR和SSIM对比

3 结论

综上所述，本文提出了基于PCA的自适应双边滤波算法，依据图像块的高斯噪声分布，对不同等级的图像噪声方差算法自行进行准确的估计，避免噪声方差参数的不平滑，最优的选取了空间标准差和亮度标准差，满足了去噪后图像细节特征信息保留的情况下，图像的峰值信噪比相比于双边滤波提高1.5963dB，图像的结构相似度相比于双边滤波提高0.06。可以解决在太赫兹反射扫描成像过程中，传统的去噪处理算法带来的细节信息丢失问题。实验结果表明，本文算法对太赫兹反射扫描图像有良好的降噪效果，结果很好地保留了图像的边缘细节特征，峰值信噪比也得到了提高。

［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］乔晓利，任姣姣，张丹丹，等.太赫兹波谱多特征参数成像方法［J］.长春理工大学学报：自然科学版，2017，40（3）：25-31.

［5］陈思宏，韩旭，李丽娟，等.基于太赫兹时域光谱技术的塑料检测［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］杨永发，李琦.双边滤波算法的太赫兹共焦扫描图像去噪应用［J］.激光与光电子学进展，2015，52（12）：93-98.

［12］朱书进，李跃华.基于自适应流形高维滤波的太赫兹图像滤波算法［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.