自适应双边滤波超声斑点降噪

摘要：斑点噪声压缩和组织结构保留对超声图像的临床诊断非常重要，因此文章提出一种自适应方差双边滤波超声斑点降噪方法。该方法用4方向梯度计算局部方差的变化量，以此变化量将超声图像归类为斑点噪声和组织结构，并用双边滤波滤除噪声。物理体膜图像和人类肝脏图像的实验结果表明该方法可以有效滤除斑点噪声和保持图像边缘。

Abstract： Speckle suppression and tissue retained are very important to clinical diagnosis of ultrasound images， therefore， a ultrasound speckle algorithm with adaptive variances bilateral filter is puts forward in this paper. Gradient operator with four directions is used to compute the differences of variances and the whole ultrasound image is classified into speckle and tissue by these differences. At last， speckle is filtered out applying bilateral filter. Ultrasound phantom testing and in vivo liver imaging show that the proposed method can effectively preserve edges and reduce speckle.

关键词：超声图像；斑点噪声；降噪；双边滤波；自适应方差

Key words： ultrasound image；speckle；noise reduction；bilateral filter；adaptive variances

中图分类号：TP391 文献标识码：A 文章编号：1006-4311（2016）04-0220-04

0 引言

医学超声图像具有安全、方便、实时、低成本和高效率的特点，因而成为最重要的诊断工具之一。然而医学超声图像等相干图像中广泛存在着一类颗粒状的斑点噪声，其表现为像素灰度值的突变，掩盖了和模糊了图像的细节，而医学超声图像的细节信息对很多临床诊断至关重要，因而斑点应当被抑制。

超声包络信号幅值服从瑞利分布[1-3]，然而从超声机器获取的图像已经被对数压缩，其信号模型可以用Loupas等人[4]提出的加性噪声模型表示。很多斑点抑制方法被提出，其中非线性滤波因其具有边缘保持和均匀区域平滑的特性而得到广泛应用，这些方法包括SRAD模型（Speckle Reducing Anisotropic Diffusion）[5]、NCD模型（nonlinear complex diffusion）[6]、LPND模型（Laplacian pyramid-based nonlinear diffusion）[7]、CTCD模型（contourlet transform based complex diffusion）[8]等。SRAD模型改进了原始各向异性扩散模型的边缘检测和扩散系数，在降噪和边缘定位等性能上有很大提高。NCD模型的基本思想是将实数域的扩散方程推广到复数域，并将虚部成分作为一个二阶的高斯导数来控制扩散的过程。NCD模型的边缘检测算子对噪声更具鲁棒性。LPND模型的基本思想是将原始超声图像映射到拉普拉斯多分辨率空间中，然后对每个尺度的图像进行不同门限下的各向异性扩散，可更有效地滤除噪声和保持图像信息。CTCD模型使用轮廓波作为多分辨工具，因而取得比LPND模型更好的结果。

上述各向异性方法在斑点噪声抑制方面取得很好的效果，但因为需要手动或根据经验设置门限值，难以准确估计像素点受斑点噪声污染程度。为了准确探测斑点噪声，根据医学超声图像在斑点区域的局部统计特征，提出一种自适应方差双边滤波超声斑点抑制方法。

1 基于局部特征变化分析的斑点噪声抑制

1.1 局部斑点统计特征

Dutt和Greenleaf提出了一种广义瑞利密度函数，即K分布[9]。斑点可以通过一个近似的K分布函数滤除，然而这个处理的实现是复杂的。对数压缩的超声图像满足加性噪声模型[4]。当前研究中假设所有的斑点完全发育。用x表示真实信号，n表示0均值乘性噪声，y表示观察到的信号，则斑点噪声模型可以表示为：

假设在均匀斑点区域，真实信号x为常数m，噪声n方差为σn2，则观察信号y的方差为：

由式（2）可知相邻斑点区域的局部方差近似不变，因而可以用局部方差的变化速率探测斑点与边缘。

在进行滤波处理前必须选择一个合适大小的窗口，为了获得可靠的斑点特征，窗口大小至少包含一个斑点分辨单元。医学超声图像中，根据图像系统参[9-10]，单元斑点尺寸可按下式计算：

Slateral/Dl=1.447・z0/fcDSaxial/Da=2.36・fs/B.W （3）

式中Slateral和Saxial分别为横向和轴向的平均斑点尺寸，Dl和Da分别为横向和轴向上每个像素的包络抽样距离，fs是轴向包络抽样频率，fc是中心频率，D是线性传感器在横向上的孔径尺寸，z0为传感器到聚焦带的距离，B.W为包络谱的带宽。当fc=3.5MHz，z0=100mm，D=16mm，Dl=0.5mm，fs=3MHz，B.W=1.5MHz，可以计算出Slateral/Dl=5.17pixels，Saxial/Da=4.72pixels，因而斑点特征的最小统计窗口应取7x7像素。

1.2 局部特征变化分析

假设任意斑点噪声区域的噪声方差都为常数且斑点完全发育，由超声图像灰度变化易知，当像素点位于斑点区域时，σ2变化平缓；当像素点位于边缘时，σ2变化剧烈。上述分析说明，当滤波窗口沿像素点滑动时，σ2的值变化缓慢，像素点极大可能位于斑点噪声区域或均匀区域，应对图像进行平滑处理；σ2的值变化剧烈，像素点位于边缘或边缘附近，应对图像不做或少做平滑处理。

超声图像局部特征的变化速率可以用梯度来度量。令Im为原始超声图像，I为由Im计算局部方差所得的特征图像，I与Im具有相同大小尺寸。i和j分别表示图像第一维与第二维索引，则I（i，j）表示超声图像Im中索引为（i，j）的像素点的局部特征值σ2 （i，j），由偏微分方程定义易得水平方向、垂直方向、45度方向、135度方向梯度幅值Ilateral（i，j）、Iaxial（i，j）、I45（i，j）、I135（i，j），其描述如下：

Ilateral（i，j）=│I（i，j+1）-I（i，j-1）│Iaxial（i，j）=│I（i+1，j）-I（i-1，j）│I45（i，j）=│I（i-1，j+1）-I（i+1，j-1）│I135（i，j）=│I（i-1，j-1）-I（i+1，j+1）│（4）

上式中符号|・|表示取绝对值。由（4）式可得点I（i，j）的梯度近似值P（i，j），其描述如下：

P（i，j）=Ilateral（i，j）+Iaxial（i，j）+I45（i，j）+I135（i，j） （5）

P（i，j）的值可以表示超声图像中点Im（i，j）的局部特征值变化速率。P（i，j）取值越大，表示σ2变化越剧烈，像素点Im（i，j）越可能为边缘；P（i，j）取值越小，表示σ2变化越小，像素点Im（i，j）越可能为斑点。根据P（i，j）的变化特征，相似度函数描述如下：

b为用户自定义参数，不同的b值满足不同的应用。Sij∈[0，1]，大的相似值（接近1）表示像素点类似于斑点，反之，小的相似值表示像素点更可能为边缘。

1.3 双边滤波

双边滤波器[11]是一种可以在平滑图像的同时保留锐利和明显边缘的非线性滤波器。它同时考虑了像素点之间的亮度距离和空间距离，可以将图像分解成不同尺度而不产生光晕，因而广泛应用于计算摄影应用。双边滤波器表现出比线性空间不变滤波器更好的边缘保持和噪声去除效果，其输入输出关系可以用下式表示：

其中D是以c为中心的局部像素集，q是D的元素。||q-c||表示局部像素q与中心像素c间的欧几里得距离，|xq-xc|表示q与c间的亮度距离，xq是像素q的灰度值。Gσr和Gσz分别是标准差为σr，σz的高斯函数。

双边滤波器具有3个参数：窗口直径、σz和σr。其中，双边滤波器窗口尺寸应当和斑点特征统计窗口大小一致。斑点噪声的灰度变化较大，为了获取足够的平滑，需要将σz设置为一个较大值以排除欧氏距离的影响。因而影响超声图像降噪效果的参数为局部窗口直径和σr。σr可以根据相似度函数自适应，其表达式为：

σr=254Sij+1 （9）

1.4 算法具体步骤

算法具体实施步骤总结如下：

①用r×r窗口计算原始超声图像Im的局部方差得到特征图像I。

②用式（4）、式（5）、式（6）计算图像I得到相似度图像。

③用式（9）计算步骤②中的结果得到Im中每一点的自适应方差σr，用直径r，方差σr的双边滤波过滤图形Im，得到滤波结果。

④用滤波结果替代Im并重复上述步骤。当重复一定次数或取得满意的效果后停止重复，得到最终降噪结果。

2 实验结果

信噪比SNR[12-13]（signal-to-noise ratio）和对比噪声比CNR[14]（contrast-to-noise ratio）是评价超声图像质量的两个主要标准。SNR和CNR描述如下：

其中μt和σt为分别为感兴趣区域的均值和方差，μb和σb分别为和感兴趣区域大小相等的背景区域的均值和方差。图像质量指标的增长率用upSNR和upCNR度量，其描述为：

upSNR=100・（SNRresult/SNRoriginal-1）（12）

upCNR=100・（CNRresult/CNRoriginal-1）（13）

为了验证本文方法的性能，实验中将本文方法与SRAD、DPAD[15]（detail preserving anisotropic diffusion）进行比较。在SRAD方法中，步长λ设置为0.5，迭代次数为100。在DPAD方法中，步长λ设置为0.8，迭代次数为200。本文方法中，斑点特征统计窗口和双边滤波器窗口大小都为9 × 9，迭代次数为15，b=12，σz=50。从图像系统参数可知，统计窗口最小为7x7并且小于一个合理值。

图1为物理体膜超声图像，图2为人类肝脏超声图像。每一个图（a）都表示原始超声图像，每一个图（b）都为SRAD降噪结果，每一个图（c）都为DPAD降噪结果，每一个图（d）都为本文方法降噪结果，每一个图中的黑色像素围绕方块为背景噪声区域，每一个图中的白色像素围绕方块为感兴趣区域。SNR用黑色方块计算，CNR用白色和黑色两个方块共同计算。SNR和CNR被用作实验性能标准。物理体膜超声图像和人类肝脏超声图像的性能标准分别展示于表1和表2中。

数量上，所有降噪结果的CNR值和SNR值都是原始超声图像的数倍，比较而言，本文方法的数值最大，性能优于其他方法。视觉上，四种方法都可以有效的抑制斑点噪声和保留图像边缘，然而本文方法的噪声抑制程度和边缘保留程度更好。综合数量标准和视觉效果，本文方法优于SRAD和DPAD。

3 结束语

本文首先通过分析得出斑点噪声的局部方差可以作为其统计特征，再通过统计特征的变化速率探测像素点是否为斑点，最后通过自适应方差双边滤波滤除噪声。和文章中其他方法相比，本文的斑点探测方法更为简单，斑点探测准确率更高，噪声滤除效果更好，边缘保留程度更高，并且不需要人工设置阈值，有效避免了准确估计阈值的难题。实验结果表明，本文算法比文章中其他算法取得了更好的噪声抑制和边缘保留效果。

参考文献：

[1]Aysal T C， Barner K E. Rayleigh-maximum-likelihood filtering for speckle reduction of ultrasound images[J]. Medical Imaging， IEEE Transactions on， 2007， 26（5）： 712-727.

[2]Krissian K， Westin C F， Kikinis R， et al. Oriented speckle reducing anisotropic diffusion[J]. Image Processing， IEEE Transactions on， 2007， 16（5）： 1412-1424.

[3]Martin-Fernandez M， Alberola-Lopez C， Ruiz-Alzola J， et al. Sequential anisotropic Wiener filtering applied to 3D MRI data[J]. Magnetic resonance imaging， 2007， 25（2）： 278-292.

[4]Chen Yan， Yin Ruming， Flynn P， et al. Aggressive region growing for speckle reduction in ultrasound images[J]. Pattern Recognition Letters， 2003， 24（4）： 677-691.

[5]Yu Yongjian， Acton S T. Speckle reducing anisotropic diffusion[J]. Image Processing， IEEE Transactions on， 2002， 11（11）： 1260-1270.

[6]Gilboa G， Sochen N， Zeevi Y Y. Image enhancement and denoising by complex diffusion processes[J]. Pattern Analysis and Machine Intelligence， IEEE Transactions on， 2004， 26（8）： 1020-1036.

[7]Zhang Fan， Yoo Y M， Mong K L， et al. Nonlinear diffusion in Laplacian pyramid domain for ultrasonic speckle reduction[J]. Medical Imaging， IEEE Transactions on， 2007， 26（2）： 200-211.

[8]Yu Jinhua and Wang Yuanyuan. Contourlet transform based complex diffusion for ultrasonic speckle reduction[C]. Proceedings of the Sixth IASTED International Conference on Biomedical Engineering. ACTA Press， 2008： 425-430.

[9]Dutt V， Greenleaf J F. Adaptive speckle reduction filter for log-compressed B-scan images[J]. Medical Imaging， IEEE Transactions on， 1996， 15（6）： 802-813.

[10]Koo J I， Park S B. Speckle reduction with edge preservation in medical ultrasonic images using a homogeneous region growing mean filter （HRGMF）[J]. Ultrasonic imaging， 1991， 13（3）： 211-237.

[11]Bae S， Paris S， Durand F. Two-scale tone management for photographic look[C]. ACM Transactions on Graphics （TOG）. ACM， 2006， 25（3）： 637-645.

[12]Behar V， Adam D， Friedman Z. A new method of spatial compounding imaging[J]. Ultrasonics， 2003， 41（5）： 377-384.

[13]Shao Dangguo， Liu P and Liu D C. Local characteristic matching based adaptive iterative bilateral filter to ultrasound speckle reduction[J]. ICIC express letters. An international journal of research and surveys. Part B， Applications， 2012， 3（1）： 77-82.

[14]Li Xiaoying， Liu D C. Ultrasound speckle reduction based on image segmentation and diffused region growing[C]. Proceedings of the 11th Joint Conference on Information Sciences，2008.

[15]Aja-Fernández S， Alberola-López C. On the estimation of the coefficient of variation for anisotropic diffusion speckle filtering[J]. Image Processing， IEEE Transactions on， 2006， 15（9）： 2694-2701.