CN105654434A - 基于统计模型的医学超声图像去噪方法 - Google Patents

Abstract

基于统计模型的医学超声图像去噪方法，包括如下步骤：步骤1)医学超声图像模型的建立；步骤2)对第一步得到的对数变换后的图像进行小波分解，得到四个频域(LL¹、LH¹、HL¹和HH¹)；对低频域LL¹继续进行小波分解，再得到四个频域(LL²、LH²、HL²和HH²)；然后重复这个步骤，直到分解最大层数J；步骤3)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理；步骤4)利用引导滤波器对最后一层的低频部分(LL^J)中的小波系数做滤波处理；步骤5)作小波逆变换处理，得到去噪后的医学超声图像。

Description

基于统计模型的医学超声图像去噪方法

技术领域

本发明应用于医学图像去噪领域。

背景技术

在医学成像领域，超声成像、CT、MRI等成像技术已应用于医学临床诊断中。超声诊断作为医学诊断中的一种诊断技术得以应用是从1972年灰阶超声的问世开始的，特别是近几年来，超声诊断在医学临床诊断中运用广泛，主要是因为：

1)超声波为非电离辐射，对人体无损害，且能够直观地显示人体器官及组织的形态结构；

2)超声波对人体软组织的分辨能力高；

3)能实时地观察人体器官的功能活动；

4)超声成像仪器操作简单，使用方便，价格便宜。

因此超声成像技术相对于其他成像技术更加安全。尤其在观察孕妇体内胎儿成长状况与诊断腹部器官病变等临床应用中，超声成像技术的使用更为重要。

根据美国癌症协会在2013年给出的一份女性患乳腺癌的统计数据，在过去的一年里，美国妇女中共有232340例新的浸润性乳腺癌病例，其中39620人死于乳腺癌。对于人体的乳腺以及腹部的其他器官进行检测的主要技术为超声成像技术，即通常所说的B超图像。因此，为医生提供更加清晰无噪声的医学超声图像具有非常重要的意义。

超声诊断主要是依据超声的良好指向性和与光相似的反射、散射、衰减及多普勒效应等物理特性，并且结合其不同的物理参数，使用不同类型的超声诊断仪器，采用各种扫查方法，将超声发射到体内，并在人体组织中传播，当正常组织和病理组织的声阻抗有一定差异时，它们组成的界面就发生反射和散射，再将此回声信号接收，加以检波等处理，显示为图像。由于各种组织的界面形态、组织器官的运动状况和对超声的吸收程度不同，其回声有一定的共性和某些特性，结合生理、病理解剖知识与临床医学，观察、分析、总结这些不同的反射规律，可对患病的部位、性质或功能障碍程度做出判断。

近年来，超声医学成像技术得到了迅猛的发展，彩色多普勒成像(CDI，ColorDopplerImaging)、组织多普勒成像(TDI，TissueDopplerImaging)、组织谐波成像(THI，TissueHarmonicImaging)、超声造影剂成像等一系列新技术接连出现并应用于临床，带动和促进了超声相关领域的研究和发展。但是，在超声成像设备的临床应用中，超声图像中的斑点噪声严重影响了超声图像的质量。由于成像机制的限制，超声图像无法避免的有质量较差这一主要缺点，特别是由于所要观察的器官或组织在结构上存在其不均匀性，有一些微小的结构超声无法识别。加上声波信号也存在干涉现象，这样在图像的形成上有了其特有的斑点噪声，这是得超声图像的质量大大降低，也使得对图像细节的识别与分析更加困难，在图像上体现为相关的形状各异的斑点噪声，它将掩盖那些灰度差别很小的图像特征。对于临床医生而言，斑点噪声对他们的准确诊断造成了很大的干扰，特别是对于经验不是很丰富的医生造成的影响更大。因此，从临床应用的角度出发，需要研究去除斑点噪声的方法，为医生做出更准确的诊断提供技术支持，降低人工诊断的风险。

为了解决这个难题，人们发展了图像去噪技术。图像去噪的过程是根据已知的降质含噪图像估计原始真实图像，得到原图像某种意义下的最优逼近。图像噪声是一个随机的过程，噪声分量灰度值是一个随机分量，按照其概率密度的统计特征可以分为：高斯噪声，椒盐噪声，泊松噪声，瑞利噪声等。根据实际图像的特点、噪声的统计特征及频谱分布规律，发展了各种不同的医学图像去噪方法。这些方法在进行图像去噪时，注重的是去噪的同时尽可能多的保持图像的基本几何结构边缘、轮廓等。但图像中的精细结构信息，细小边缘和纹理等常会被平滑掉。因此解决图像去噪中的平滑噪声和保持结构信息这一矛盾是图像去噪中经久不衰的课题。医学图像去噪声处理属于图像的预处理阶段，从数字图像处理的技术角度来说属于图像恢复的技术范畴，它的存在有着非常重要的意义，主要表现在：

(1)医学图像由于不同的成像机理，得到的初始图像中都含有大量的不同性质的噪声，这些噪声的存在影响着人们的对医学图像的观察，干扰人们对图像信息的理解。噪声严重时候，图像几乎产生变形，更使得图像失去了存储信息的本质意义。显然对图像进行去噪处理，提高医学超声图像的质量，改善视觉效果，是正确识别图像信息的必要保证。

(2)除了能提高人视觉识别信息的准确性，图像去噪的意义还在于它是对图像作进一步处理的可靠保证。如果对一幅含有噪声的医学图像进行特征提取、配准或者图像融合等处理其结果肯定不能令人满意。特别是对于医学图像处理来说，必须要求每一步有尽可能的准确性和可靠性，方便医生更加准确地针对病灶区域做出判断，降低辅助诊断的风险。所以医学图像去噪是必需的。

在数字图像处理领域，有很多传统的图像去噪方法，它们可能已经被提出以至被应用很久了。在这样的学术背景下依然研究医学图像去噪的意义何在？我想意义依然是有的，在于：

(1)虽然医学图像去噪技术是以一般数字图像处理技术为基础，但是医学图像本身具有自己一些鲜明的特征，这些特征正是医学图像所含有的特殊信息。在对医学图像进行去噪处理的时候必须尽可能地保留这些特征，这就需要我们研究新的算法使得这些算法在保留一般数字图像去噪性能的基础上还能满足医学图像去噪的特殊要求。

(2)在医学图像去噪领域，传统方法呈百花齐放之态，但是这些方法并非十全十美，主要表现在去噪的同时对图像细节的丢失。因此进一步研究新的去噪方法或者完善已有的方法意义依然重大。

(3)不同方法都有着不同的数学理论基础，对图像去噪的效果也表现不同。探求它们的内部机理，寻求相应的关系，研究不同方法之间如何取长补短，以达到更好的去噪效果，也是很有意义的。

(4)研究图像去噪的同时对医学图像其他处理环节性能的提升也有着促进意义。

(5)从实际应用来讲，由于医院资源的局限性，特别是医生每天进行人工诊断病人的数量无法满足社会整个阶层的需求，即面临着病人多医生少情况。因此，各种自动诊断仪器的需求越来越大，自动诊断仪器的出现，一方面可以节约医生资源，另一方面可以方便更多的病人进行诊断。随着当今社会经济的飞跃发展，人们自身健康情况却不容乐观，所以人们对家用型医疗自动诊断仪器的需求也非常大，例如家用超声图像自动诊断仪等。但是超声图像自动诊断仪是建立在高质量图像的基础上，并且自动诊断仪需要对超声图像做后期的智能分析，如特征提取、边缘检测和图像分类识别等。因此，从自动化诊断技术的角度出发，需要研究去除斑点噪声的方法，为图像的后期智能处理提供技术保障，促进超声图像自动化诊断技术的发展，具有不可估量的价值。

医学超声图像中斑点噪声的独特性造就了其对于去噪要求的特殊性，首先斑点噪声在图像中表现出大的“颗粒状”；其次一些病灶细节体现在图像边缘处；为了适用于实时成像系统，去噪的方法的效率要高。所以针对以上要求，提出的去噪方法应该具有以下目标：

(1)具有较强的去除斑点噪声的能力。这也是对于去噪方法的最基本的要求，为用户提供高质量清晰的图像也是建立在较强去噪能力的基础上；

(2)具有边缘保持能力。传统的去噪方法尽管具有较强的去噪能力，但是往往会过多的平滑图像的边缘而使得图像过于模糊，这对超声诊断是非常不利的，特别是在病灶部分，对于医师做出正确判断会造成很大的影响。

(3)不存在梯度变形。此目标是针对双边滤波器存在梯度变形的情况而提出来的。

(4)非迭代。在最新的去噪方法当中，很多迭代的方法应运而生，尽管其具有很强的去噪能力，但是这往往是建立在迭代次数比较多的基础上，相应的运算时间就会比较长，不适用于实时成像系统。

(5)算法复杂度低，快速性。该目标和第4点很类似，就是要提高算法的快速性，以达到成像系统对于实时性的要求。

发明内容

根据医学超声图像对于去噪的要求已经要实现的目标，本发明要克服现有技术的上述缺点，提供一种基于统计模型的医学超声图像去噪方法。

本发明分析当前频率域内和空间域内的去噪方法并结合斑点噪声的模型的特点和医学超声图像的处理需求提出了一种新的去噪方法，即一种基于统计模型的医学超声图像去噪方法。小波变换具有时频分析和多尺度分析等优越性，其已在图像处理领域得到了广泛的应用。在处理加性噪声问题时，小波的去噪效果较好，能够满足一般产品需求。然而，仅仅利用小波变换的去噪方法对医学超声图像中斑点噪声的抑制效果不好。对于引导滤波器，它在处理图像噪声时，一方面具有很强的去噪能力，另一方面能够保持图像边缘细节。虽然双边滤波器在去噪的过程中也能很好的保持边缘信息，但是其效率较低，运行时间将会很长并且存在“梯度反转”现象，难以用于实时系统。随着图像的分辨率越来越大，这在很大的程度上限制了双边滤波的应用空间。因此本发明利用引导滤波器替换掉双边滤波器，极大提高去噪性能和运算效率。因此，本发明将结合小波去噪和引导滤波器的优点。具体思路如下，在传统的小波去噪方法的基础之上，根据小波域内超声图像及斑点噪声的统计特性，改进了小波阈值函数和收缩方法，能够更有效地去除高频部分的斑点噪声。由于医学超声图像在小波域内的低频部分依然存在斑点噪声，因此使用去噪效果好并且效率高的的引导滤波器，在抑制低频域内噪声的同时能够保留低频域内的图像边缘信息。本发明解决斑点噪声问题是通过以下技术方案实现的：

步骤1)医学超声图像模型的建立

如果认为超声成像系统能够对那些影响声波功率的因素做出恰当的动态补偿，则超声成像系统采集的包络信号由两部分组成，一是有意义的体内组织的反射信号，另一部分是噪声信号。其中噪声信号可分为相乘噪声与相加噪声。相乘噪声与超声信号成像的原理有关，主要来源于随机的散射信号。相加噪声认为是系统噪声，如传感器的噪声等。超声成像系统初步得到的包络信号为f^pre，它的一般模型如下

f^pre＝g^pren^pre+w^pre(1)

这里，上标^pre表示系统初步得到的信号。函数g^pre表示无噪声信号，n^pre和w^pre分别表示相乘噪声和相加噪声，式中n^pre是噪声的主要成分。

和相乘噪声n^pre相比，相加噪声w^pre所占比重很小，因此将w^pre忽略后的模型为

f^pre＝g^pren^pre(2)

为了适应超声成像系统显示屏幕的动态显示范围，对超声成像系统采集到的包络信号进行对数压缩处理。此时相乘的式(2)模型将变为相加的模型，如下

log(f^pre)＝log(g^pre)+log(n^pre)(3)

此时，得到的信号log(f^pre)即是通常看到的医学超声图像。

步骤2)对第一步得到的对数变换后的图像进行小波分解，得到四个频域(LL¹、LH¹、HL¹和HH¹)。

LL¹分量是对原始信号LL0的列和行进行小波分解后得到的低频分量，即一级小波分解后近似部分，它包含了原始图像最多的低频信息。

LH¹是一次小波分解后的垂直方向上的高频分量，即它包含了图像水平方向上的近似信息和垂直方向上的边缘等高频信息。

HL¹是一次小波分解后的水平方向上的高频分量，即它包含了图像垂直方向上的近似信息和水平方向上的边缘等高频信息。

HH¹是一次小波分解后对角方向上的高频分量，即它包含了图像水平和垂直方向上的边缘等高频信息。

对低频域LL¹继续进行小波分解，再得到第二层四个频域(LL²、LH²、HL²和HH²)。然后重复这个步骤，直到分解最大层数J。

由于小波变换是线性变换，因此式(3)模型经过二维离散小波变换后得到下面模型：

$W_{l,k}^j\left(log\right(f^{pre}\left)\right)=W_{l,k}^j\left(log\right(g^{pre}\left)\right)+W_{l,k}^j\left(log\right(n^{pre}\left)\right)---\left(4\right)$

其中和分别表示含有噪声图像的小波系数、无噪声图像的小波系数和斑点噪声的小波系数。其中上标j为小波变换的分解层数，下标(l,k)为小波域内的坐标。为了方便表示，将式(4)简化为

$F_{l,k}^j=G_{l,k}^j+N_{l,k}^j---\left(5\right)$

本发明认为经过小波分解后的无噪信号的小波系数服从广义拉普拉斯分布，其概率分布如下

p_G(g)＝C(σ_g,β)exp{-[K(σ_g,β)|g]β},-∞<g<+∞,σ_g>0,β>0(6)

其中

$K({\mathrm{\&sigma;}}_g,\mathrm{\&beta;})={\mathrm{\&sigma;}}_g^{-1}{\&lsqb;\mathrm{\&Gamma;}(3/\mathrm{\&beta;})/\mathrm{\&Gamma;}(1/\mathrm{\&beta;})\&rsqb;}^{1/2}$

C(σ_g,β)＝βK(σ_g,β)/[2Γ(1/β)]

其中K(σ_g，β)是信号能量系数，C(σ_g,β)是归一化因子，是伽马函数。σ_g是无噪信号小波系数的标准差，决定广义拉普拉斯分布概率密度函数的扩散程度；β为形状参数，控制着广义拉普拉斯分布概率密度函数的衰减速度。当β＝1时，式(6)将变为拉普拉斯分布，是广义拉普拉斯分布的特殊模型。

为了更好的描述不同散射信号情况下斑点噪声的特性，斑点噪声的小波系数被认为服从瑞利分布

$p_N\left(n\right)=\frac{n}{{\mathrm{\&sigma;}}_n^2}\exp (-\frac{n^2}{2{\mathrm{\&sigma;}}_n^2})---\left(7\right)$

式中σ_n为小波域内噪声的标准差。

步骤3)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理。

在小波去噪方法中，阈值函数的选择会直接影响到最终的图像去噪结果。当阈值选择较小时，一部分大于该阈值的噪声系数会被当作有用信号保留下来，这就导致去噪后的图像依然存在大量噪声；当阈值选择较大时，会将很多系数很小的有用信息当作噪声而置零，这将使得去噪后的图像变得很平滑，损失很多细节信息。因此选择恰当的小波阈值函数非常重要。

Donoho等人提出了一种典型的阈值选取方法，并且从理论上证明了该阈值与噪声的标准差成正比，改阈值函数又称为统一阈值函数，其公式如下

$T={\mathrm{\&sigma;}}_n\sqrt{2\log M}---\left(8\right)$

其中，M即是对应小波域内小波系数的总体个数，σ_n是噪声的标准差。在这种阈值函数中，阈值T受小波系数的个数影响较大，即当M过大时，较大的阈值可能会平滑掉那些系数较小的有用信息。

在式(8)的基础之上，提出了一种更加适合超声图像的阈值函数，其公式如下

$T={\mathrm{\&alpha;}}_j{\mathrm{\&sigma;}}_n\sqrt{2\log M}---\left(9\right)$

其中，σ_n是噪声的标准差，a_j代表j层的自适应参数。这是种常见的阈值改进的方法，a_j的选取是根据实验决定的，如果选择恰当将会得到更好的效果，试验中选择 $a_j=\frac{1}{ln(1+j)}.$

在小波去噪方法中，首先选定一个给定阈值，然后按照一定的规则对小波系数进行收缩，便完成了对小波系数的去噪。即给定一个阈值，所有绝对值小于这个阈值的系数被当作噪声，然后对其作置零处理；对绝对值大于阈值的小波系数用一定的方法进行缩减，然后得到缩减后的新值。

经典的小波收缩方法有软阈值法和硬阈值法，但是在软阈值法中，较大的小波系数总是被阈值缩减，因此收缩后的信号的数学期望与收缩之前不同，所以处理后的图像相对平滑一些。硬阈值法的缺点是在零值域附近的小波系数被突然置零，导致了小波数据的不连续性，并且这使得信号的方差更大了，这些变换对于图像中的细节影响较大。但是在实际应用中，特别是噪声水平很高时，硬阈值法处理后的图像在不连续点周围会产生震荡，影响图像的去噪效果。

由于经典的阈值收缩方法不能满足对医学超声图像去噪的要求，所以对收缩方法做了改进。

无噪信号的小波系数服从广义拉普拉斯分布，小波域内的斑点噪声部分服从瑞利分布。为了简化计算，本发明选择β＝1，则式(6)变为拉普拉斯分布

$p_G\left(g\right)=\frac{1}{\sqrt{2}{\mathrm{\&sigma;}}_g}\exp (-\frac{\sqrt{2}\left|g\right|}{{\mathrm{\&sigma;}}_g})---\left(10\right)$

为了得到小波域内的信号估计值，使用贝叶斯最大后验估计的方法。在后验概率的计算过程中，使用贝叶斯公式如下

$\begin{array}{c}{p}_{G|F}\left(g|f\right)=\frac{1}{p_F\left(f\right)}{p}_{F|G}\left(f|g\right)\&CenterDot;p_G\left(g\right)=\\ \frac{1}{p_F\left(f\right)}{p}_N\left(f-g\right)\&CenterDot;p_G\left(g\right)\end{array}---\left(11\right)$

将式(7)、式(10)带入上式(11)，得到

$\begin{array}{c}{p}_{G|F}\left(g|f\right)=\frac{1}{p_F\left(f\right)}\&CenterDot;\frac{f-g}{\sqrt{2}{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_g}\&times;\\ \exp \left(-\frac{2\sqrt{2}{\mathrm{\&sigma;}}_N^2\left|g\right|+{\mathrm{\&sigma;}}_g{\left(f-g\right)}^2}{2{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_g}\right)\end{array}---\left(12\right)$

为了得到最大后验概率，将ln(p_G|F(g|f))对g求一次导数的方程置零，最后得到

$\hat{g}=sign\left(f\right)\&CenterDot;max\left(\right|f|-\frac{{\mathrm{\&sigma;}}_n^2+\sqrt{{\mathrm{\&sigma;}}_n^4+2{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_g^2}}{\sqrt{2}{\mathrm{\&sigma;}}_g},0)---\left(13\right)$

为g的估计，并且假设包络信号f和经过对数变换后的无噪信号g同号；这样就得到新的收缩方法

$\hat{g}=\left\{\begin{array}{cc}0& f\&le;T_j\\ sign\left(f\right)\&CenterDot;max\left(\right|f|-\frac{{\mathrm{\&sigma;}}_n^2+\sqrt{{\mathrm{\&sigma;}}_n^4+2{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_g^2}}{\sqrt{2}{\mathrm{\&sigma;}}_g},0)& f>T_j\end{array}\right.---\left(14\right)$

斑点噪声的标准差由小波域内最高频的小波系数(即式5中最高频子带)绝对值的中值计算得到

${\widehat{\mathrm{\&sigma;}}}_N=\frac{median\left(\right|F_{l,k}^{{\mathrm{HH}}^1}\left|\right)}{0.6745}---\left(15\right)$

由式(5)可以得到j层小波系数的方差关系式

${\mathrm{\&sigma;}}_{F,j}^2={\mathrm{\&sigma;}}_{G,j}^2+{\mathrm{\&sigma;}}_N^2---\left(16\right)$

由于和都被建模为零均值的式(6)和式(7)模型，因此小波系数在j层的标准差σ_F,j可由小波系数的标准差计算得到

${\widehat{\mathrm{\&sigma;}}}_{F,j}=\sqrt{\frac{1}{N}\underset{l,k=1}{\overset{m,n}{\&Sigma;}}{\left(F_{l,k}^j\right)}^2}---\left(17\right)$

这里的N＝m×n即是对应小波域内小波系数的总个数，m和n分别是对应小波域内系数的行数和列数。

由式(15)、式(16)和式(17)可以得到

${\widehat{\mathrm{\&sigma;}}}_{G,j}^2=max({\widehat{\mathrm{\&sigma;}}}_{F,j}^2-{\widehat{\mathrm{\&sigma;}}}_N^2,0)---\left(18\right)$

本发明改进的小波收缩函数在曲线图像上表现的更加平滑，尤其当小波系数大于小波阈值的区间范围内。

但是就传统的噪声方差估计方法(式15)，实验表明当图像噪声比较弱的时候，在小波域内最高频分区内，边缘细节的比重相对于噪声的分量要大的多，所以按照传统的噪声方差估计方法的到的σ_n与实际的噪声方差有很大的误差，这直接影响到阈值选取的和收缩算法执行的正确性。另外这种方法是建立在噪声绝对服从加性高斯分布假设的基础上，因而对高斯白噪声会有很好的去噪效果，但是对其他噪声(如脉冲噪声)会产生很大的误差，前面的分析中，医学超声图像中的噪声是服从瑞利分布的，并不是严格意义上的高斯分布。本发明就解决这个问题，提出新的噪声方差的估计方法，以两个Laplacian算子之差作为滤波器对图像进行滤波通过滤波结果的方差估计噪声方差。

步骤4)利用引导滤波器对最后一层的低频部分(LL^J)中的小波系数做滤波处理

一般基于小波的去噪方法，即保留低频域(LL)的小波系数不变，仅对高频域(LH、HL、HH)的小波系数做阈值处理。然而，此方法应用于医学超声图像去噪时表现不佳。经过多次实验，发现低频域内的小波系数依然具有很多斑点噪声，为了更有效地滤除低频域内的斑点噪声，本发明选择引导滤波器对低频域内的小波系数作滤波处理。

引导图像滤波由局部线性模型发展而来，该法基本原理如下式所示

$q_i=\underset{j}{\&Sigma;}{W}_{ij}\left(I\right)p_j---\left(19\right)$

式(8)中，I为引导图像，p为输入图像，q为输出图像，，W_ij为关于引导图像I的函数，i和j为像素点的位置，I由具体问题确定，可以令I＝p。

假设在窗口w_k内，中心点为后k，q为I的线性变换，如式(20)所示

$q_i=a_k{I}_i+b_k,\&ForAll;i\&Element;w_k---\left(20\right)$

在图像滤波中，希望能在达到滤波效果的前提下最小化输入图像和输出图像的差异，减小原始图像细节的损失，故通过最小化p和q的差来确定系数a_k和b_k，即使式(21)最小

$E(a_k,b_k)=\underset{i\&Element;w_k}{\&Sigma;}\&lsqb;{(a_k{I}_i+b_k-p_i)}^2+{\mathrm{\&epsiv;a}}_k^2\&rsqb;---\left(21\right)$

式(21)中，ε是正则化参数，目的是为了防止a_k过大。求解式(21)，得

$a_k=\frac{\frac{1}{\left|w\right|}\underset{i\&Element;w_k}{\&Sigma;}{I}_i{p}_i-{\mathrm{\&mu;}}_k{\stackrel{\&OverBar;}{p}}_k}{{\mathrm{\&sigma;}}_k^2+\mathrm{\&epsiv;}}$

$b_k={\stackrel{\&OverBar;}{p}}_k-a_k{\mathrm{\&mu;}}_k---\left(22\right)$

${\stackrel{\&OverBar;}{p}}_k=\frac{1}{\left|w\right|}\underset{i\&Element;w_k}{\&Sigma;}{p}_i$

式中，μ_k和分别为I在w_k中的均值和方差。|w|为w_k中的像素个数，是输入图像p在w_k中的均值。求出该线性模型后，带入整幅图像，由于每一个像素点会有多个包含该像素点的窗口w_k，所以当在不同窗口w_k计算时，q_i的值会不同。故需要对其进行平均处理

$q_i=\frac{1}{\left|w\right|}\underset{k,i\&Element;w_k}{\&Sigma;}(a_k{I}_i+b_k)=\stackrel{\&OverBar;}{a_i}{I}_i+\stackrel{\&OverBar;}{b_i}---\left(23\right)$

式中，系数a的平均值系数b的平均值

综上所述，核函数W_ij可以定义如下

$W_{ij}=\frac{1}{|w{|}^2}\underset{k,(i,j)\&Element;w_k}{\&Sigma;}(1+\frac{(I_i-{\mathrm{\&mu;}}_k)(I_j-{\mathrm{\&mu;}}_k)}{{\mathrm{\&sigma;}}_k^2+\mathrm{\&epsiv;}})---\left(24\right)$

由以上原理可知，引导滤波器去噪的过程如下：

(1)输入图像p；

(2)输入滤波窗口w_k的尺寸和正则化参数ε；

(3)计算I,p和I*p的均值；

(4)计算(I,p)的协方差；

(5)计算(I*I)的均值并计算I的方差；

(6)计算系数a,b；

(7)分别计算a和b的均值；

(8)得到输出图像q。

但是在选取引导图像I时，由于在实际医学图像中不存在无噪声的图像，也就不存在无噪的引导图像，实验表明当以图像自身作为引导图像时，并不能达到理想的效果。所以为了解决这一问题，本发明利用简单的高斯滤波去对图像进行预处理，将得到的图像作为引导图像，再利用上述引导滤波器去噪的过程，去噪效果大大改善，解决了以图像自身作为引导图像时去噪效果不佳的问题。

步骤5)作小波逆变换处理，得到去噪后的医学超声图像。

如果要得到去噪后的超声包络信号，对第5步得到的超声图像作指数变换即可。与现有技术相比，本发明的新颖性和创造性在于：

本发明一方面对通用的阈值函数和收缩方法进行了改进，对高频域的细小噪声具有很强的去噪能力，在噪声方差估计中提出新的估计方法可以适用于不同水平的斑点噪声；另一方面，由于本发明利用引导滤波器对低频部分进行了过滤处理，因此对于颗粒较大的斑点噪声(存在于低频部分)同样具有很强的抑制能力，同时在引导滤波器的使用上又引入新的引导图像解决了以图像自身作为引导图像时去噪效果不佳的问题。由于本发明的时间的消耗很大部分是由于引入的用于处理低频部分噪声的滤波器，所以在本发明中引导滤波器的引入不仅能很好的提高去噪的性能并且大大提高了处理的效率。同时针对医学超声图像的特点，这种结合的方法不仅能很好的抑制斑点噪声，同时还能够保留图像中病灶边缘等的细节部分，能更好的帮助医师进行病情分析。

本发明的优点是：不仅能很好的抑制斑点噪声，同时还能够保留图像中病灶边缘等的细节部分，能更好的帮助医师进行病情分析。

附图说明

图1是仿真图想的去噪效果的对比图，其中图1a是无噪图像，图1b是噪声图像，图1c是小波软阈值法去噪效果图，图1d是小波-双边滤波器去噪效果图，图1e是本发明方法的法去噪效果图。

图2是临床超声图像的去噪结果对比图，其中图2a是噪声图像，图2b是小波软阈值法去噪效果图，图2c是小波-双边滤波器去噪效果图，图2d是本发明方法的法去噪效果图。

具体实施方式

为使本发明的目的、技术方案和优点更加清晰，下面就对本发明的技术方案作进一步描述。

步骤1)医学超声图像模型的建立

如果认为超声成像系统能够对那些影响声波功率的因素做出恰当的动态补偿，则超声成像系统采集的包络信号由两部分组成，一是有意义的体内组织的反射信号，另一部分是噪声信号。其中噪声信号可分为相乘噪声与相加噪声。相乘噪声与超声信号成像的原理有关，主要来源于随机的散射信号。相加噪声认为是系统噪声，如传感器的噪声等。超声成像系统初步得到的包络信号为f^pre，它的一般模型如下

f^pre＝g^pren^pre+w^pre(1)

这里，上标^pre表示系统初步得到的信号。函数g^pre表示无噪声信号，n^pre和w^pre分别表示相乘噪声和相加噪声，式中n^pre是噪声的主要成分。

和相乘噪声n^pre相比，相加噪声w^pre所占比重很小，因此将w^pre忽略后的模型为

f^pre＝g^pren^pre(2)

为了适应超声成像系统显示屏幕的动态显示范围，对超声成像系统采集到的包络信号进行对数压缩处理。此时相乘的式(2)模型将变为相加的模型，如下

log(f^pre)＝log(g^pre)+log(n^pre)(3)

此时，得到的信号log(f^pre)即是通常看到的医学超声图像。

步骤2)对第一步得到的对数变换后的图像进行小波分解，得到四个频域(LL¹、LH¹、HL¹和HH¹)。

LL¹分量是对原始信号LL0的列和行进行小波分解后得到的低频分量，即一级小波分解后近似部分，它包含了原始图像最多的低频信息。

LH¹是一次小波分解后的垂直方向上的高频分量，即它包含了图像水平方向上的近似信息和垂直方向上的边缘等高频信息。

HL¹是一次小波分解后的水平方向上的高频分量，即它包含了图像垂直方向上的近似信息和水平方向上的边缘等高频信息。

HH¹是一次小波分解后对角方向上的高频分量，即它包含了图像水平和垂直方向上的边缘等高频信息。

对低频域LL¹继续进行小波分解，再得到第二层四个频域(LL²、LH²、HL²和HH²)。然后重复这个步骤，直到分解最大层数J。

由于小波变换是线性变换，因此式(3)模型经过二维离散小波变换后得到下面模型：

$W_{l,k}^j\left(log\right(f^{pre}\left)\right)=W_{l,k}^j\left(log\right(g^{pre}\left)\right)+W_{l,k}^j\left(log\right(n^{pre}\left)\right)---\left(4\right)$

其中和分别表示含有噪声图像的小波系数、无噪声图像的小波系数和斑点噪声的小波系数。其中上标j为小波变换的分解层数，下标(l,k)为小波域内的坐标。为了方便表示，将式(4)简化为

$F_{l,k}^j=G_{l,k}^j+N_{l,k}^j---\left(5\right)$

本发明认为经过小波分解后的无噪信号的小波系数服从广义拉普拉斯分布，其概率分布如下

p_G(g)＝C(σ_g,β)exp{-[K(σ_g,β)|g]^β},-∞<g<+∞,σ_g>0,β>0(6)

其中

$K({\mathrm{\&sigma;}}_g,\mathrm{\&beta;})={\mathrm{\&sigma;}}_g^{-1}{\&lsqb;\mathrm{\&Gamma;}(3/\mathrm{\&beta;})/\mathrm{\&Gamma;}(1/\mathrm{\&beta;})\&rsqb;}^{1/2}$

C(σ_g,β)＝βK(σ_g,β)/[2Γ(1/β)]

其中K(σ_g，β)是信号能量系数，C(σ_g,β)是归一化因子，是伽马函数。σ_g是无噪信号小波系数的标准差，决定广义拉普拉斯分布概率密度函数的扩散程度；β为形状参数，控制着广义拉普拉斯分布概率密度函数的衰减速度。当β＝1时，式(6)将变为拉普拉斯分布，是广义拉普拉斯分布的特殊模型。

为了更好的描述不同散射信号情况下斑点噪声的特性，斑点噪声的小波系数被认为服从瑞利分布

$p_N\left(n\right)=\frac{n}{{\mathrm{\&sigma;}}_n^2}\exp (-\frac{n^2}{2{\mathrm{\&sigma;}}_n^2})---\left(7\right)$

式中σ_n为小波域内噪声的标准差。

步骤3)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理。

在小波去噪方法中，阈值函数的选择会直接影响到最终的图像去噪结果。当阈值选择较小时，一部分大于该阈值的噪声系数会被当作有用信号保留下来，这就导致去噪后的图像依然存在大量噪声；当阈值选择较大时，会将很多系数很小的有用信息当作噪声而置零，这将使得去噪后的图像变得很平滑，损失很多细节信息。因此选择恰当的小波阈值函数非常重要。

Donoho等人提出了一种典型的阈值选取方法，并且从理论上证明了该阈值与噪声的标准差成正比，改阈值函数又称为统一阈值函数，其公式如下

$T={\mathrm{\&sigma;}}_n\sqrt{2\log M}---\left(8\right)$

其中，M即是对应小波域内小波系数的总体个数，σ_n是噪声的标准差。在这种阈值函数中，阈值T受小波系数的个数影响较大，即当M过大时，较大的阈值可能会平滑掉那些系数较小的有用信息。

在式(8)的基础之上，提出了一种更加适合超声图像的阈值函数，其公式如下

$T={\mathrm{\&alpha;}}_j{\mathrm{\&sigma;}}_n\sqrt{2\log M}---\left(9\right)$

其中，σ_n是噪声的标准差，a_j代表j层的自适应参数。这是种常见的阈值改进的方法，a_j的选取是根据实验决定的，如果选择恰当将会得到更好的效果，试验中选择 $a_j=\frac{1}{ln(1+j)}.$

在小波去噪方法中，首先选定一个给定阈值，然后按照一定的规则对小波系数进行收缩，便完成了对小波系数的去噪。即给定一个阈值，所有绝对值小于这个阈值的系数被当作噪声，然后对其作置零处理；对绝对值大于阈值的小波系数用一定的方法进行缩减，然后得到缩减后的新值。

经典的小波收缩方法有软阈值法和硬阈值法，但是在软阈值法中，较大的小波系数总是被阈值缩减，因此收缩后的信号的数学期望与收缩之前不同，所以处理后的图像相对平滑一些。硬阈值法的缺点是在零值域附近的小波系数被突然置零，导致了小波数据的不连续性，并且这使得信号的方差更大了，这些变换对于图像中的细节影响较大。但是在实际应用中，特别是噪声水平很高时，硬阈值法处理后的图像在不连续点周围会产生震荡，影响图像的去噪效果。

由于经典的阈值收缩方法不能满足对医学超声图像去噪的要求，所以对收缩方法做了改进。

无噪信号的小波系数服从广义拉普拉斯分布，小波域内的斑点噪声部分服从瑞利分布。为了简化计算，选择β＝1，则式(6)变为拉普拉斯分布

$p_G\left(g\right)=\frac{1}{\sqrt{2}{\mathrm{\&sigma;}}_g}\exp (-\frac{\sqrt{2}\left|g\right|}{{\mathrm{\&sigma;}}_g})---\left(10\right)$

为了得到小波域内的信号估计值，使用贝叶斯最大后验估计的方法。在后验概率的计算过程中，使用贝叶斯公式如下

$\begin{array}{c}{p}_{G|F}\left(g|f\right)=\frac{1}{p_F\left(f\right)}{p}_{F|G}\left(f|g\right)\&CenterDot;p_G\left(g\right)=\\ \frac{1}{p_F\left(f\right)}{p}_N\left(f-g\right)\&CenterDot;p_G\left(g\right)\end{array}---\left(11\right)$

将式(7)、式(10)带入上式(11)，得到

$\begin{array}{c}{p}_{G|F}\left(g|f\right)=\frac{1}{p_F\left(f\right)}\&CenterDot;\frac{f-g}{\sqrt{2}{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_g}\&times;\\ \exp \left(-\frac{2\sqrt{2}{\mathrm{\&sigma;}}_N^2\left|g\right|+{\mathrm{\&sigma;}}_g{\left(f-g\right)}^2}{2{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_g}\right)\end{array}---\left(12\right)$

为了得到最大后验概率，将ln(p_G|F(g|f))对g求一次导数的方程置零，最后得到

$\hat{g}=sign\left(f\right)\&CenterDot;max\left(\right|f|-\frac{{\mathrm{\&sigma;}}_n^2+\sqrt{{\mathrm{\&sigma;}}_n^4+2{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_g^2}}{\sqrt{2}{\mathrm{\&sigma;}}_g},0)---\left(13\right)$

$\hat{g}$ 为g的估计，并且假设包络信号f和经过对数变换后的无噪信号g同号；这样就得到新的收缩方法

$\hat{g}=\left\{\begin{array}{cc}0& f\&le;T_j\\ sign\left(f\right)\&CenterDot;max\left(\right|f|-\frac{{\mathrm{\&sigma;}}_n^2+\sqrt{{\mathrm{\&sigma;}}_n^4+2{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_g^2}}{\sqrt{2}{\mathrm{\&sigma;}}_g},0)& f>T_j\end{array}\right.---\left(14\right)$

斑点噪声的标准差由小波域内最高频的小波系数(即式5中 $F_{l,k}^j$ 最高频子带 $F_{l,k}^{{\mathrm{HH}}^1}$ )绝对值的中值计算得到

${\widehat{\mathrm{\&sigma;}}}_N=\frac{median\left(\right|F_{l,k}^{{\mathrm{HH}}^1}\left|\right)}{0.6745}---\left(15\right)$

由式(5)可以得到j层小波系数的方差关系式

${\mathrm{\&sigma;}}_{F,j}^2={\mathrm{\&sigma;}}_{G,j}^2+{\mathrm{\&sigma;}}_N^2---\left(16\right)$

由于和都被建模为零均值的式(6)和式(7)模型，因此小波系数在j层的标准差σ_F,j可由小波系数的标准差计算得到

${\widehat{\mathrm{\&sigma;}}}_{F,j}=\sqrt{\frac{1}{N}\underset{l,k=1}{\overset{m,n}{\&Sigma;}}{\left(F_{l,k}^j\right)}^2}---\left(17\right)$

这里的N＝m×n即是对应小波域内小波系数的总个数，m和n分别是对应小波域内系数的行数和列数。

由式(15)、式(16)和式(17)可以得到

${\widehat{\mathrm{\&sigma;}}}_{G,j}^2=max({\widehat{\mathrm{\&sigma;}}}_{F,j}^2-{\widehat{\mathrm{\&sigma;}}}_N^2,0)---\left(18\right)$

本发明改进的小波收缩函数在曲线图像上表现的更加平滑，尤其当小波系数大于小波阈值的区间范围内。

但是就传统的噪声方差估计方法(式15)，实验表明当图像噪声比较弱的时候，在小波域内最高频分区内，边缘细节的比重相对于噪声的分量要大的多，所以按照传统的噪声方差估计方法的到的σ_n与实际的噪声方差有很大的误差，这直接影响到阈值选取的和收缩算法执行的正确性。另外这种方法是建立在噪声绝对服从加性高斯分布假设的基础上，因而对高斯白噪声会有很好的去噪效果，但是对其他噪声(如脉冲噪声)会产生很大的误差，前面的分析中，医学超声图像中的噪声是服从瑞利分布的，并不是严格意义上的高斯分布。本发明就解决这个问题，提出新的噪声方差的估计方法，以两个Laplacian算子之差作为滤波器对图像进行滤波通过滤波结果的方差估计噪声方差。

步骤4)利用引导滤波器对最后一层的低频部分(LL^J)中的小波系数做滤波处理

一般基于小波的去噪方法，即保留低频域(LL)的小波系数不变，仅对高频域(LH、HL、HH)的小波系数做阈值处理。然而，此方法应用于医学超声图像去噪时表现不佳。经过多次实验，发现低频域内的小波系数依然具有很多斑点噪声，为了更有效地滤除低频域内的斑点噪声，本发明选择引导滤波器对低频域内的小波系数作滤波处理。

引导图像滤波由局部线性模型发展而来，该法基本原理如下式所示

$q_i=\underset{j}{\&Sigma;}{W}_{ij}\left(I\right)p_j---\left(19\right)$

式(8)中，I为引导图像，p为输入图像，q为输出图像，W_ij为关于引导图像I的函数，i和j为像素点的位置，I由具体问题确定，可以令I＝p。

假设在窗口w_k内，中心点为后k，q为I的线性变换，如式(20)所示

$q_i=a_k{I}_i+b_k,\&ForAll;i\&Element;w_k---\left(20\right)$

在图像滤波中，希望能在达到滤波效果的前提下最小化输入图像和输出图像的差异，减小原始图像细节的损失，故通过最小化p和q的差来确定系数a_k和b_k，即使式(21)最小

$E(a_k,b_k)=\underset{i\&Element;w_k}{\&Sigma;}\&lsqb;{(a_k{I}_i+b_k-p_i)}^2+{\mathrm{\&epsiv;a}}_k^2\&rsqb;---\left(21\right)$

式(21)中，ε是正则化参数，目的是为了防止a_k过大。求解式(21)，得

$\begin{array}{c}{a}_k=\frac{\frac{1}{\left|w\right|}\underset{i\&Element;w_k}{\&Sigma;}{I}_i{p}_i-{\mathrm{\&mu;}}_k{\stackrel{\&OverBar;}{p}}_k}{{\mathrm{\&sigma;}}_k^2+\mathrm{\&epsiv;}}\\ b_k={\stackrel{\&OverBar;}{p}}_k-a_k{\mathrm{\&mu;}}_k\\ {\stackrel{\&OverBar;}{p}}_k=\frac{1}{\left|w\right|}\underset{i\&Element;w_k}{\&Sigma;}{p}_i\end{array}---\left(22\right)$

式中，μ_k和分别为I在w_k中的均值和方差。|w|为w_k中的像素个数， ${\stackrel{\&OverBar;}{p}}_k$ 是输入图像p在w_k中的均值。求出该线性模型后，带入整幅图像，由于每一个像素点会有多个包含该像素点的窗口w_k，所以当在不同窗口w_k计算时，q_i的值会不同。故需要对其进行平均处理

$q_i=\frac{1}{\left|w\right|}\underset{k,i\&Element;w_k}{\&Sigma;}(a_k{I}_i+b_k)=\stackrel{\&OverBar;}{a_i}{I}_i+\stackrel{\&OverBar;}{b_i}---\left(23\right)$

式中，系数a的平均值系数b的平均值

综上所述，核函数W_ij可以定义如下

$W_{ij}=\frac{1}{|w{|}^2}\underset{k,(i,j)\&Element;w_k}{\&Sigma;}(1+\frac{(I_i-{\mathrm{\&mu;}}_k)(I_j-{\mathrm{\&mu;}}_k)}{{\mathrm{\&sigma;}}_k^2+\mathrm{\&epsiv;}})---\left(24\right)$

由以上原理可知，引导滤波器去噪的过程如下：

(1)输入图像p；

(2)输入滤波窗口w_k的尺寸和正则化参数ε；

(3)计算I,p和I*p的均值；

(4)计算(I,p)的协方差；

(5)计算(I*I)的均值并计算I的方差；

(6)计算系数a,b；

(7)分别计算a和b的均值；

(8)得到输出图像q。

但是在选取引导图像I时，由于在实际医学图像中不存在无噪声的图像，也就不存在无噪的引导图像，实验表明当以图像自身作为引导图像时，并不能达到理想的效果。所以为了解决这一问题，本发明利用简单的高斯滤波去对图像进行预处理，将得到的图像作为引导图像，再利用上述引导滤波器去噪的过程，去噪效果大大改善，解决了以图像自身作为引导图像时去噪效果不佳的问题。

步骤5)作小波逆变换处理，得到去噪后的医学超声图像。

经过阈值收缩处理和引导滤波器处理就可以得到去噪后的小波系数，为了得到去噪后的超声图像，需要对小波系数进行小波逆变换，从而可以得到利于医师分析的去噪后的图像，通过实验也验证了本发明确实可以满足对于医学超声图像去噪的要求。

实验验证

为了客观地评价本发明提出的去噪方法，以峰值信噪比(PSNR)、结构相似度(SSIM)、FoM(Pratt’sFigureofMerit)和运行时间作为图像质量评价标准。峰值信噪比的计算公式如下

$PSNR(X,\hat{X})=10\lg \left(\frac{{255}^2}{MSE}\right)---\left(25\right)$

式中，为信号X的估计值，MSE由下面公式计算得到

$MSE=\frac{1}{MN}\underset{i=1}{\overset{M}{\&Sigma;}}\underset{j=1}{\overset{N}{\&Sigma;}}{(X_{i,j}-{\hat{X}}_{i,j})}^2---\left(26\right)$

这里的M，N分别表示二维信号X的长度与宽度。

结构相似度能够量化两幅图像在结构上的差异，公式定义如下

$SSIM(X,\hat{X})=\frac{(2{\mathrm{\&mu;}}_X{\mathrm{\&mu;}}_{\hat{X}}+c_1)(2{\mathrm{\&sigma;}}_{X,\hat{X}}+c_2)}{({\mathrm{\&mu;}}_X^2+{\mathrm{\&mu;}}_{\hat{X}}^2+c_1)({\mathrm{\&sigma;}}_X^2+{\mathrm{\&sigma;}}_{\hat{X}}^2+c_2)}---\left(27\right)$

式中，μ_X、和分别是参考图像和估计图像的均值和方差。是X和的协方差，c₁和c₂为常量。当c₁和c₂都选择为正数时，SSIM的取值范围为[01]，其中1为最好结果，表示两幅图的结构相同。

FoM能够客观地比较去噪图像的边缘检测质量，公式定义如下

$FoM(X,\hat{X})=\frac{1}{max(N_X,N_{\hat{X}})}\underset{i=1}{\overset{N_X}{\mathrm{\&Sigma;}}}\frac{1}{1+{\mathrm{\&alpha;d}}_i^2}---\left(28\right)$

式中，N_X和分别表示理想的和实际检测到的边缘像素个数。α为常数(通常取α＝1/9)，d_i表示为第i边缘像素点到最近理想边缘像素点的距离。FoM的取值范围为[01]，其中1为最好结果，表示为检测到的图像边缘和理想的图像边缘一致。这里检测边缘像素时使用的是Canny检测算法(高斯滤波器的标准差取值σ＝3)。

当然为了更好的展现本发明的优势，接下来主要是通过对比实验，量化比较各个图像评价标准来使本发明的目的和优点更加清晰。

为了定量估计本发明的去噪的效果，先进行仿真图像的实验。得到如图1和表1中的数据

表1去噪方法的性能比较

用于仿真的图像是通过对无噪图像加入噪声得到的。通过仿真实验定量得出的数据可以看出，如果仅仅利用小波变换的去噪方法在抑制斑点噪声中的效果并不是很好，在低频部分依然存在一定的斑点噪声，所以需要引入滤波器对低频部分做滤波处理，本发明引入的是引导滤波器。通过比较可以发现，本发明的实验结果在直观上也能看出来优于小波软阈值的方法，引入引导滤波器的去噪效果也好于引入双边滤波器，利用评价指标也验证了本发明极大提高了去噪效果。

本方法在对临床超声图像进行试验，选择的是含有病灶的妇女乳腺组织的超声图像，如图2所示。

由于现实中不存在无噪声的超声图像，因此，PSNR等质量指标在这里无法有效地使用，所以就引入另一个指标，无参考图像质量指标(NIQE)。得到表2

表2去噪方法的性能比较

去噪方法	NIQE
		噪声图像	5.1143
小波软阈值	7.1911
		小波‐双边滤波器	8.6848
本发明方法	9.1477

通过比较可以发现本发明在实际应用于医学超声图像中，去噪效果得到了显著提高。在去除噪声的同时能够更多的保留图像边缘信息，从而达到医学超声图像对于去噪的要求。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。

Claims (1)

1.基于统计模型的医学超声图像去噪方法，包括如下步骤：

步骤1)医学超声图像模型的建立；

如果认为超声成像系统能够对那些影响声波功率的因素做出恰当的动态补偿，则超声成像系统采集的包络信号由两部分组成，一是有意义的体内组织的反射信号，另一部分是噪声信号；其中噪声信号可分为相乘噪声与相加噪声；相乘噪声与超声信号成像的原理有关，主要来源于随机的散射信号；相加噪声认为是系统噪声，如传感器的噪声等；超声成像系统初步得到的包络信号为f^pre，它的一般模型如下

f^pre＝g^pren^pre+w^pre(1)

这里，上标^pre表示系统初步得到的信号；函数g^pre表示无噪声信号，n^pre和w^pre分别表示相乘噪声和相加噪声，式中n^pre是噪声的主要成分；

和相乘噪声n^pre相比，相加噪声w^pre所占比重很小，因此将w^pre忽略后的模型为

f^pre＝g^pren^pre(2)

为了适应超声成像系统显示屏幕的动态显示范围，对超声成像系统采集到的包络信号进行对数压缩处理；此时相乘的式(2)模型将变为相加的模型，如下

log(f^pre)＝log(g^pre)+log(n^pre)(3)

此时，得到的信号log(f^pre)即是通常看到的医学超声图像；

步骤2)对第一步得到的对数变换后的图像进行小波分解，得到四个频域(LL¹、LH¹、HL¹和HH¹)；

LL¹分量是对原始信号LL0的列和行进行小波分解后得到的低频分量，即一级小波分解后近似部分，它包含了原始图像最多的低频信息。

LH¹是一次小波分解后的垂直方向上的高频分量，即它包含了图像水平方向上的近似信息和垂直方向上的边缘等高频信息。

HL¹是一次小波分解后的水平方向上的高频分量，即它包含了图像垂直方向上的近似信息和水平方向上的边缘等高频信息。

HH¹是一次小波分解后对角方向上的高频分量，即它包含了图像水平和垂直方向上的边缘等高频信息。

对低频域LL¹继续进行小波分解，再得到第二层四个频域(LL²、LH²、HL²和HH²)；然后重复这个步骤，直到分解最大层数J；

由于小波变换是线性变换，因此式(3)模型经过二维离散小波变换后得到下面模型：

$W_{l,k}^j\left(log\right(f^{pre}\left)\right)=W_{l,k}^j\left(log\right(g^{pre}\left)\right)+W_{l,k}^j\left(log\right(n^{pre}\left)\right)---\left(4\right)$

其中和分别表示含有噪声图像的小波系数、无噪声图像的小波系数和斑点噪声的小波系数；其中上标j为小波变换的分解层数，下标(l,k)为小波域内的坐标；为了方便表示，将式(4)简化为

$F_{l,k}^j=G_{l,k}^j+N_{l,k}^j---\left(5\right)$

经过小波分解后的无噪信号的小波系数服从广义拉普拉斯分布，其概率分布如下

p_G(g)＝C(σ_g,β)exp{-[K(σ_g,β)|g|]^β},-∞<g<+∞,σ_g>0,β>0(6)

其中

$K({\mathrm{\&sigma;}}_g,\mathrm{\&beta;})={\mathrm{\&sigma;}}_g^{-1}{\&lsqb;\mathrm{\&Gamma;}(3/\mathrm{\&beta;})/\mathrm{\&Gamma;}(1/\mathrm{\&beta;})\&rsqb;}^{1/2}$

C(σ_g,β)＝βK(σ_g,β)/[2Γ(1/β)]

其中K(σ_g，β)是信号能量系数，C(σ_g,β)是归一化因子，是伽马函数；σ_g是无噪信号小波系数的标准差，决定广义拉普拉斯分布概率密度函数的扩散程度；β为形状参数，控制着广义拉普拉斯分布概率密度函数的衰减速度；当β＝1时，式(6)将变为拉普拉斯分布，是广义拉普拉斯分布的特殊模型；

为了更好的描述不同散射信号情况下斑点噪声的特性，斑点噪声的小波系数被认为服从瑞利分布

$p_N\left(n\right)=\frac{n}{{\mathrm{\&sigma;}}_n^2}\exp (-\frac{n^2}{2{\mathrm{\&sigma;}}_n^2})---\left(7\right)$

式中σ_n为小波域内噪声的标准差；

步骤3)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理；

在小波去噪方法中，阈值函数的选择会直接影响到最终的图像去噪结果；当阈值选择较小时，一部分大于该阈值的噪声系数会被当作有用信号保留下来，这就导致去噪后的图像依然存在大量噪声；当阈值选择较大时，会将很多系数很小的有用信息当作噪声而置零，这将使得去噪后的图像变得很平滑，损失很多细节信息；因此选择恰当的小波阈值函数非常重要；

Donoho等人提出了一种典型的阈值选取方法，并且从理论上证明了该阈值与噪声的标准差成正比，改阈值函数又称为统一阈值函数，其公式如下

$T={\mathrm{\&sigma;}}_n\sqrt{2\log M}---\left(8\right)$

其中，M即是对应小波域内小波系数的总体个数，σ_n是噪声的标准差；在这种阈值函数中，阈值T受小波系数的个数影响较大，即当M过大时，较大的阈值可能会平滑掉那些系数较小的有用信息；

在式(8)的基础之上，提出了一种更加适合超声图像的阈值函数，其公式如下

$T={\mathrm{\&alpha;}}_j{\mathrm{\&sigma;}}_n\sqrt{2\log M}---\left(9\right)$

其中，σ_n是噪声的标准差，a_j代表j层的自适应参数；这是种常见的阈值改进的方法，a_j的选取是根据实验决定的，如果选择恰当将会得到更好的效果，试验中选择 $a_j=\frac{1}{ln(1+j)};$

在小波去噪方法中，首先选定一个给定阈值，然后按照一定的规则对小波系数进行收缩，便完成了对小波系数的去噪；即给定一个阈值，所有绝对值小于这个阈值的系数被当作噪声，然后对其作置零处理；对绝对值大于阈值的小波系数用一定的方法进行缩减，然后得到缩减后的新值；

经典的小波收缩方法有软阈值法和硬阈值法，但是在软阈值法中，较大的小波系数总是被阈值缩减，因此收缩后的信号的数学期望与收缩之前不同，所以处理后的图像相对平滑一些；硬阈值法的缺点是在零值域附近的小波系数被突然置零，导致了小波数据的不连续性，并且这使得信号的方差更大了，这些变换对于图像中的细节影响较大；但是在实际应用中，特别是噪声水平很高时，硬阈值法处理后的图像在不连续点周围会产生震荡，影响图像的去噪效果；

由于经典的阈值收缩方法不能满足对医学超声图像去噪的要求，所以对收缩方法做了改进；

无噪信号的小波系数服从广义拉普拉斯分布，小波域内的斑点噪声部分服从瑞利分布；为了简化计算，选择β＝1，则式(6)变为拉普拉斯分布

$p_G\left(g\right)=\frac{1}{\sqrt{2}{\mathrm{\&sigma;}}_g}\exp (-\frac{\sqrt{2}\left|g\right|}{{\mathrm{\&sigma;}}_g})---\left(10\right)$

为了得到小波域内的信号估计值，使用贝叶斯最大后验估计的方法；在后验概率的计算过程中，使用贝叶斯公式如下

$\begin{array}{c}{p}_{G|F}\left(g|f\right)=\frac{1}{p_F\left(f\right)}{p}_{F|G}\left(f|g\right)\&CenterDot;p_G\left(g\right)=\\ \frac{1}{p_F\left(f\right)}{p}_N\left(f-g\right)\&CenterDot;p_G\left(g\right)\end{array}---\left(11\right)$

将式(7)、式(10)带入上式(11)，得到

$\begin{array}{c}{p}_{G|F}\left(g|f\right)=\frac{1}{p_F\left(f\right)}\&CenterDot;\frac{f-g}{\sqrt{2}{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_g}\&times;\\ \exp \left(-\frac{2\sqrt{2}{\mathrm{\&sigma;}}_N^2\left|g\right|+{\mathrm{\&sigma;}}_g{\left(f-g\right)}^2}{2{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_g}\right)\end{array}---\left(12\right)$

为了得到最大后验概率，将ln(p_G|F(g|f))对g求一次导数的方程置零，最后得到

$\hat{g}=sign\left(f\right)\&CenterDot;max\left(\right|f|-\frac{{\mathrm{\&sigma;}}_n^2+\sqrt{{\mathrm{\&sigma;}}_n^4+2{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_g^2}}{\sqrt{2}{\mathrm{\&sigma;}}_g},0)---\left(13\right)$

为g的估计，并且假设包络信号f和经过对数变换后的无噪信号g同号；这样就得到新的收缩方法

$\hat{g}=\left\{\begin{array}{cc}0& f\&le;T_j\\ sign\left(f\right)\&CenterDot;max\left(\right|f|-\frac{{\mathrm{\&sigma;}}_n^2+\sqrt{{\mathrm{\&sigma;}}_n^4+2{\mathrm{\&sigma;}}_n^2{\mathrm{\&sigma;}}_g^2}}{\sqrt{2}{\mathrm{\&sigma;}}_g},0)& f>T_j\end{array}\right.---\left(14\right)$

斑点噪声的标准差由小波域内最高频的小波系数(即式5中最高频子带)绝对值的中值计算得到

${\widehat{\mathrm{\&sigma;}}}_N=\frac{median\left(\right|F_{l,k}^{{\mathrm{HH}}^1}\left|\right)}{0.6745}---\left(15\right)$

由式(5)可以得到j层小波系数的方差关系式

${\mathrm{\&sigma;}}_{F,j}^2={\mathrm{\&sigma;}}_{G,j}^2+{\mathrm{\&sigma;}}_N^2---\left(16\right)$

由于和都被建模为零均值的式(6)和式(7)模型，因此小波系数在j层的标准差σ_F，j可由小波系数的标准差计算得到

${\widehat{\mathrm{\&sigma;}}}_{F,j}=\sqrt{\frac{1}{N}\underset{l,k=1}{\overset{m,n}{\&Sigma;}}{\left(F_{l,k}^j\right)}^2}---\left(17\right)$

这里的N＝m×n即是对应小波域内小波系数的总个数，m和n分别是对应小波域内系数的行数和列数；

由式(15)、式(16)和式(17)可以得到

${\widehat{\mathrm{\&sigma;}}}_{G,j}^2=max({\widehat{\mathrm{\&sigma;}}}_{F,j}^2-{\widehat{\mathrm{\&sigma;}}}_N^2,0)---\left(18\right)$

小波收缩函数在曲线图像上表现的更加平滑，尤其当小波系数大于小波阈值的区间范围内；

但是就传统的噪声方差估计方法(式15)，实验表明当图像噪声比较弱的时候，在小波域内最高频分区内，边缘细节的比重相对于噪声的分量要大的多，所以按照传统的噪声方差估计方法的到的σ_n与实际的噪声方差有很大的误差，这直接影响到阈值选取的和收缩算法执行的正确性；另外这种方法是建立在噪声绝对服从加性高斯分布假设的基础上，因而对高斯白噪声会有很好的去噪效果，但是对其他噪声(如脉冲噪声)会产生很大的误差，前面的分析中，医学超声图像中的噪声是服从瑞利分布的，并不是严格意义上的高斯分布；以两个Laplacian算子之差作为滤波器对图像进行滤波通过滤波结果的方差估计噪声方差；

步骤4)利用引导滤波器对最后一层的低频部分(LL^J)中的小波系数做滤波处理；

为了更有效地滤除低频域内的斑点噪声，选择引导滤波器对低频域内的小波系数作滤波处理；

引导图像滤波由局部线性模型发展而来，该法基本原理如下式所示

$q_i=\underset{j}{\&Sigma;}{W}_{ij}\left(I\right)p_j---\left(19\right)$

式(8)中，I为引导图像，p为输入图像，q为输出图像，W_ij为关于引导图像I的函数，i和j为像素点的位置，I由具体问题确定，可以令I＝p；

假设在窗口w_k内，中心点为后k，q为输出图像，如式(20)所示

$q_i=a_k{I}_i+b_k,\&ForAll;i\&Element;w_k---\left(20\right)$

在图像滤波中，希望能在达到滤波效果的前提下最小化输入图像和输出图像的差异，减小原始图像细节的损失，故通过最小化p和q的差来确定系数a_k和b_k，即使式(21)最小

$E(a_k,b_k)=\underset{i\&Element;w_k}{\&Sigma;}\&lsqb;{(a_k{I}_i+b_k-p_i)}^2+{\mathrm{\&epsiv;a}}_k^2\&rsqb;---\left(21\right)$

式(21)中，ε是正则化参数，目的是为了防止a_k过大；求解式(21)，得

$a_k=\frac{\frac{1}{\left|w\right|}\underset{i\&Element;w_k}{\&Sigma;}{I}_i{p}_i-{\mathrm{\&mu;}}_k{\stackrel{\&OverBar;}{p}}_k}{{\mathrm{\&sigma;}}_k^2+\mathrm{\&epsiv;}}$

$b_k={\stackrel{\&OverBar;}{p}}_k-a_k{\mathrm{\&mu;}}_k---\left(22\right)$

${\stackrel{\&OverBar;}{p}}_k=\frac{1}{\left|w\right|}\underset{i\&Element;w_k}{\&Sigma;}{p}_i$

式中，μ_k和分别为I在w_k中的均值和方差；|w|为w_k中的像素个数，是输入图像p在w_k中的均值；求出该线性模型后，带入整幅图像，由于每一个像素点会有多个包含该像素点的窗口w_k，所以当在不同窗口w_k计算时，q_i的值会不同；故需要对其进行平均处理

$q_i=\frac{1}{\left|w\right|}\underset{k,i\&Element;w_k}{\&Sigma;}(a_k{I}_i+b_k)={\stackrel{\&OverBar;}{a}}_i{I}_i+\stackrel{\&OverBar;}{b_i}---\left(23\right)$

式中，系数a的平均值系数b的平均值

综上所述，核函数W_ij可以定义如下

$W_{ij}=\frac{1}{|w{|}^2}\underset{k,(i,j)\&Element;w_k}{\&Sigma;}(1+\frac{(I_i-{\mathrm{\&mu;}}_k)(I_j-{\mathrm{\&mu;}}_k)}{{\mathrm{\&sigma;}}_k^2+\mathrm{\&epsiv;}})---\left(24\right)$

由以上原理可知，引导滤波器去噪的过程如下：

(51)输入图像p；

(52)输入滤波窗口w_k的尺寸和正则化参数ε；

(53)计算I,p和I*p的均值；

(54)计算(I,p)的协方差；

(55)计算(I*I)的均值并计算I的方差；

(56)计算系数a,b；

(57)分别计算a和b的均值；

(58)得到输出图像q；

利用简单的高斯滤波去对图像进行预处理，将得到的图像作为引导图像，再利用上述引导滤波器去噪的过程，去噪效果大大改善，解决了以图像自身作为引导图像时去噪效果不佳的问题；

步骤5)作小波逆变换处理，得到去噪后的医学超声图像。