基于贝叶斯模型的超声图像斑点噪声滤波方法

摘要

本发明公开了一种基于贝叶斯模型的超声图像斑点噪声滤波方法，具体包括：S1、建立基于贝叶斯模型的非局部滤波模型；S2、采用Gamma分布拟合经对数压缩后超声图像中的斑点噪声，进而得出基于Pearson统计距离的权重函数；S3、求解非局部滤波模型中的变量；S4、采用求解后的非局部滤波模型对超声图像中斑点噪声进行滤波处理。本发明能够有效的滤除超声图像中的斑点噪声，提高超声图像的成像清晰度。

说明书

技术领域

本发明涉及超声图像滤波技术领域，尤其涉及一种基于贝叶斯模型的超声 图像斑点噪声滤波方法。

背景技术

在超声信号的采集过程中，由于斑点噪声的存在，使得重建的超声图像质 量明显比CT、MRI等成像模态的差。为了提高超声的图像质量，通常需要对斑 点噪声进行滤波去噪。但由于斑点噪声对人体组织的依赖性，导致对斑点噪声 的建模和滤除非常困难。对超声图像的斑点噪声滤波问题，虽然已有大量的研 究报道，但经典的滤波方法在抑制噪声的同时会丢失图像中的边缘等细节信息。 近些年来，基于偏微分方程的各向异性扩散滤波模型开始应用于图像的平滑去 噪。各向异性扩散模型的优点是在对图像进行平滑去噪的同时会很好的保留图 像中的边缘等细节信息。

最简单的扩散方程是把图像平滑看作是各向同性的热传导方程的解。为了 获得对加性噪声模型的平滑滤波结果，一个直观的想法就是最小化由估计的 图像变化值。1963年Tikhonov提出了对应的变分问题，即在物理学领域著名的 热扩散方程。热扩散方程的解等价于观察图像与高斯核函数的卷积，因为卷积 算子是线性的，所以热扩散方程的解也是线性的。由于一幅图像与高斯核地卷 积等价于它的Fourier变换与另一个高斯核的乘积，所以各向同性的扩散过程本 质上是一个低通滤波的过程，它将抑制图像中的高频信号。不幸的是，图像的 许多特征，如边缘、纹理等，通常表现为高频信号。因此边缘、纹理等有用信 息也将与噪声一起被各向同性的滤波器滤除掉。

为了克服上述线性滤波方法导致的各向同性平滑，Perona和Malik对热扩散 方程进行了扩展，提出了著名的P-M非线性滤波方法。P-M方程通过把热扩散 方程改写成散度的形式，并在散度算子内添加递减函数作为图像特征检测算子， 以允许对图象滤波过程进行更加精确的控制。但P-M方程在一些图像点上产生 逆扩散。而我们知道，逆扩散虽然会增强图像边缘特征，但噪声同样也有可能 得到增强，因此是一个不稳定的扩散过程。但经典的扩散方程只是基于图像局 部边缘信息的，没有从全局图像去考虑。

Gilboa和Osher引入一种非局部平滑各向异性扩散方程，该方法基本思想是 利用图像的冗余信息以及自相似性对图像进行去噪。非局部滤波方法将在全局 图像中搜索与当前像素块相似的像素块，最后将所有这些相似像素块的加权平 均值作为噪声点处的像素值。但非局部扩散方程只适用于高斯白噪声模型的平 滑去噪，但并不适用于服从Gamma统计分布的超声图像斑点噪声滤除。

因此，针对上述技术问题，有必要提供一种基于贝叶斯模型的超声图像斑 点噪声滤波方法。

发明内容

有鉴于此，本发明的目的在于提供一种基于贝叶斯模型的超声图像斑点噪 声滤波方法。

为了达到上述目的，本发明实施例提供的技术方案如下：

一种基于贝叶斯模型的超声图像斑点噪声滤波方法，所述方法包括：

S1、建立基于贝叶斯模型的非局部滤波模型

$J_{\mathrm{Bayesian}-\mathrm{NL}-\mathrm{TV}}\left(u\right)={\int }_{\Omega }\left|\right|{▿}_{\mathrm{NL}}u\left|\right|+\frac{\lambda }{2}{(u-u_0)}^2\mathrm{dx}$

$={\int }_{\Omega }\sqrt{{\int }_{\Omega }{(u\left(y\right)-u\left(x\right))}^{2}w(x,y)\mathrm{dy}}+\frac{\lambda }{2}{(u-u_0)}^2\mathrm{dx},$

其中，第一项λ为平滑项，第二项(u-u₀)²为数据保真项，λ是一个正常数， 用于控制平滑项与数据保真项之间的平衡；数据保真项用于使求解的值u不会 偏移原始观测值u₀太远；

S2、采用Gamma分布拟合经对数压缩后超声图像中的斑点噪声，进而得出 基于Pearson统计距离的权重函数；

S3、求解非局部滤波模型中的变量u和d，其中

S4、采用求解后的非局部滤波模型对超声图像中斑点噪声进行滤波处理。

作为本发明的进一步改进，所述步骤S2包括：

采用Gamma分布拟合经对数压缩后超声图像中的斑点噪声：

u(x)＝v(x)+v^γ(x)η(x)，

其中，γ是一个依赖于超声设备并与后续成像过程相关的参数。

作为本发明的进一步改进，所述步骤S2中γ为0.5。

作为本发明的进一步改进，所述步骤S2中：

两两像素块之间相似度的Pearson统计距离为：

$d_p(x,y)={\int }_{\Omega }\frac{{(u(x+·)-u(y+·))}^2}{u(y+·)}\mathrm{dy},$

基于Pearson统计距离的权重函数为：

$w(x,y)=\frac{1}{Z}\exp (-\frac{1}{h^2}{\int }_{\Omega }\frac{{(u(x+·)-u(y+·))}^2}{u(y+·)}\mathrm{dy}),$

其中，u(x+·)和u(y+·)为图像中的两个像素块。

作为本发明的进一步改进，所述步骤S3中非局部滤波模型的求解采用 Split-Bregman分解法、最陡梯度下降法、或对偶投影法。

作为本发明的进一步改进，所述步骤S3中的非局部滤波模型为：

$\left(\begin{array}{c}\underset{u,d}{\arg \min }{\int }_{\Omega }\left|d\right|\mathrm{dx}+\frac{\lambda }{2}{\int }_{\Omega }{(u-u^0)}^2\mathrm{dx}\\ \mathrm{subjecttod}={▿}_{\mathrm{NL}}u\end{array}\right).$

作为本发明的进一步改进，所述步骤S3中非局部滤波模型的求解通过求解 如下以离散形式给出的迭代方程获得：

$\underset{u,d}{\arg \min }\underset{i,j}{\Sigma }\left|d_{\mathrm{ij}}\right|+\frac{\lambda }{2}\underset{i,j}{\Sigma }{(u_{\mathrm{ij}}-u_{0,\mathrm{ij}})}^2+\frac{\beta }{2}\underset{i,j}{\Sigma }{|d_{\mathrm{ij}}-{▿}_{\mathrm{NL}}{u}_{\mathrm{ij}}-b_{\mathrm{ij}}|}^2,$

其中，b是与Bregman迭代算法相关的一个变量，也称作Bregman变量，β 是一个大于零的常数。

作为本发明的进一步改进，所述步骤S3中非局部滤波模型中变量u的求解 具体包括：

固定变量d，则非局部滤波模型简化为求关于变量u的最小化子问题：

$\underset{u,d}{\arg \min }\frac{\lambda }{2}\underset{i,j}{\Sigma }{(u_{\mathrm{ij}}-u_{0,\mathrm{ij}})}^2+\frac{\beta }{2}\underset{i,j}{\Sigma }{|d_{\mathrm{ij}}-{▿}_{\mathrm{NL}}{u}_{\mathrm{ij}}-b_{\mathrm{ij}}|}^2,$

通过Eular-Lagrange方程，非局部滤波模型的最小解满足：

$\lambda (u-u_0)-\beta {\mathrm{div}}_{\mathrm{NL}}(d^k-{▿}_{\mathrm{NL}}u-b^k)=0,$

使用高效的Gauss-Seidel迭代法进行求解，即对n≥0有：

$\left(\begin{array}{c}{u}_{0,i}^{k+1}=u_i^{k+1}\\ u_i^{k+1,n+1}=\frac{1}{\lambda +\beta {\Sigma }_j{w}_{\mathrm{ij}}}\left(\begin{array}{c}\beta {\Sigma }_j{w}_{\mathrm{ij}}{u}_j^{k+1,n}+{\mathrm{\lambda u}}_{0,i}-\\ \beta {\Sigma }_j\sqrt{w_{\mathrm{ij}}}(d_{\mathrm{ij}}^{k+1,n}-d_{\mathrm{ji}}^{k+1,n}-b_{\mathrm{ij}}^{k+1,n}-b_{\mathrm{ji}}^{k+1,n})\end{array},\right)\end{array}\right),$

通过至少两次Gauss-Seidel迭代收敛到方程的最小解处。

作为本发明的进一步改进，所述步骤S3中非局部滤波模型中变量d的求解 具体包括：

固定变量u，则非局部滤波模型简化为求关于变量d的最小化子问题：

$\underset{u,d}{\arg \min }\underset{i,j}{\Sigma }\left|d_{\mathrm{ij}}\right|+\frac{\beta }{2}\underset{i,j}{\Sigma }{|d_{\mathrm{ij}}-{\mathrm{\Delta u}}_{\mathrm{ij}}{-b}_{\mathrm{ij}}|}^2,$

采用软阈值法，则非局部滤波模型的Eular-Lagrange方程为：

$d_{\mathrm{ij}}^{k+1}=\frac{{▿}_{\mathrm{NL}}{u}_{\mathrm{ij}}^{k+1}+b_{\mathrm{ij}}^k}{|{▿}_{\mathrm{NL}}{u}_{\mathrm{ij}}^{k+1}+b_{\mathrm{ij}}^k|}\max \left(\right|{▿}_{\mathrm{NL}}{u}_{\mathrm{ij}}^{k+1}+b_{\mathrm{ij}}^k|-\beta ,0)$

$=\frac{\sqrt{w_{\mathrm{ij}}}(u_j^{k+1}-u_i^{k+1})+b_{\mathrm{ij}}^k}{\sqrt{{\Sigma }_j{w}_{\mathrm{ij}}{(u_j^{k+1}-u_i^{k+1})}^2+{b_{\mathrm{ij}}^k}^2}}\max (\sqrt{{\Sigma }_j{w}_{\mathrm{ij}}{(u_j^{k+1}-u_i^{k+1})}^2+{b_{\mathrm{ij}}^k}^2}-\beta ,0),$

根据上述公式求解变量d；

最后，变量b可根据如下公式进行更新：

$b_{\mathrm{ij}}^{k+1}=b_{\mathrm{ij}}^k+\sqrt{w_{\mathrm{ij}}}(u_j^{k+1}-u_i^{k+1})-d_{\mathrm{ij}}^{k+1}.$

本发明具有以下有益效果：

本发明根据经对数变换后的超声图像中斑点噪声数理统计特性，获得了基 于Pearson统计距离的权重计算函数，是一种非参数统计的方法，避免了传统扩 散方程中对参数的复杂设计，对斑点噪声取得了较好的滤波效果。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施 例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述 中的附图仅仅是本发明中记载的一些实施例，对于本领域普通技术人员来讲， 在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明基于贝叶斯模型的超声图像斑点噪声滤波方法的流程图。

具体实施方式

为了使本技术领域的人员更好地理解本发明中的技术方案，下面将结合本 发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述， 显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基 于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获 得的所有其他实施例，都应当属于本发明保护的范围。

参图1所示，本发明公开了一种基于贝叶斯模型的超声图像斑点噪声滤波 方法，具体包括：

S1、建立基于贝叶斯模型的非局部滤波模型；

S2、采用Gamma分布拟合经对数压缩后超声图像中的斑点噪声，进而得出 基于Pearson统计距离的权重函数；

S3、求解非局部滤波模型中的变量；

S4、采用求解后的非局部滤波模型对超声图像中斑点噪声进行滤波处理。

以下结合具体实施方式对本发明作进一步说明。

Bresson的研究表明，在现有的非局部扩散滤波模型中非局部全变分扩散模 型具有更好的平滑滤波效果。根据贝叶斯定理，给非局部扩散项添加数据保真 项(fidelityterm)的限制：

$J_{\mathrm{Bayesian}-\mathrm{NL}-\mathrm{TV}}\left(u\right)={\int }_{\Omega }\left|\right|{▿}_{\mathrm{NL}}u\left|\right|+\frac{\lambda }{2}{(u-u_0)}^2\mathrm{dx}$

$={\int }_{\Omega }\sqrt{{\int }_{\Omega }{(u\left(y\right)-u\left(x\right))}^{2}w(x,y)\mathrm{dy}}+\frac{\lambda }{2}{(u-u_0)}^2\mathrm{dx}.---\left(1\right)$

其中，λ是一个正的常数，用于控制平滑项与数据保真项之间的平衡；第二 项(u-u₀)²为数据保真项，它的作用是使求解的值u不会偏移原始观测值u₀太远， 否则也容易导致对图像的过度平滑滤波。从数理统计的角度看，公式(1)中的 平滑项对应贝叶斯定理中有关噪声模型的先验项(prior term)，而数据保真项则 对应其中的似然项(likelihood term)，因此该扩散方程的求解过程也就是使后验 概率(posterior)最大化的过程。

对超声图像的研究表明，Gamma分布能更好的对经对数压缩后超声图像中 的斑点噪声进行近似：

u(x)＝v(x)+v^γ(x)η(x)。   (2)

其中，γ是一个依赖于超声设备并与后续成像过程相关的参数。与Rayleigh 分布相比，Gamma分布更加灵活并能可靠地捕捉到超声图像的统计特征；与高 斯白噪声模型相比，伽马噪声模型是与图像相关的。

通过计算经对数压缩变换后图像中均值和标准差的比值，Loupas等人的实 验结果表明当参数γ取0.5时，伽马分布对超声图像中斑点噪声的拟合要优于瑞 利分布模型和乘性噪声模型。因此，本实施方式中将使用参数γ＝0.5的伽马噪声 模型对超声图像中的斑点噪声进行建模和处理。

针对超声图像中斑点噪声的Gamma分布统计特征，本发明将引入基于 Pearson统计距离的权重函数，以取代经典的基于高斯加权的权重函数。已知超 声图像中的斑点噪声服从Gamma分布，那么对采样图像中每一个像素的观察值 有：u(x)|v(x)～N(v(x),v(x)σ²)，其条件概率密度函数可表述为：

$p\left(u\left(x\right)\right|v\left(x\right))\propto \exp (-\frac{{(u\left(x\right)-v\left(x\right))}^2}{2v\left(x\right)})..---\left(3\right)$

假设像素间的噪声分布是相互独立的，根据概率链式法则，图像中像素块 u(x+·)和u(y+·)间的条件概率密度函数为：

$p\left(u(x+·)\right|u(y+·))=\underset{i=1}{\overset{N}{\Pi }}p\left(u^{\left(i\right)}\left(x\right)\right|u^{\left(i\right)}\left(y\right))$

$\propto \exp (-{\int }_{\Omega }\frac{{(u(x+·)-u(y+·))}^2}{u(y+·)}\mathrm{dy})..---\left(4\right)$

因此，从公式(4)可得如下用于度量两两像素块之间相似度的Pearson统 计距离：

$d_p(x,y)={\int }_{\Omega }\frac{{(u(x+·)-u(y+·))}^2}{u(y+·)}\mathrm{dy}.---\left(5\right)$

用Pearson统计距离代替公式(1)中的欧拉距离，可得如下基于噪声统计 特性的权重距离函数：

$w(x,y)=\frac{1}{Z}\exp (-\frac{1}{h^2}{\int }_{\Omega }\frac{{(u(x+·)-u(y+·))}^2}{u(y+·)}\mathrm{dy}.---\left(6\right)$

相比于高斯加权函数，公式(6)的最大优点是该加权函数是完全基于当前 图像的观测值，是一种非参数估计的方法，更适用于超声图像的斑点去除。

为了求解非局部滤波方程(1)，本发明采用Split-Bregman分解的数值求解 方法。通过引入一个辅助变量d:Ω×Ω→R，并对该辅助变量加入限制那么扩散方程(1)可重新表述为：

$\left(\begin{array}{c}\underset{u,d}{\arg \min }{\int }_{\Omega }\left|d\right|\mathrm{dx}+\frac{\lambda }{2}{\int }_{\Omega }{(u-u^0)}^2\mathrm{dx}\\ \mathrm{subjecttod}={▿}_{\mathrm{NL}}u\end{array}\right).---\left(7\right)$

根据Bregman迭代方法，公式(7)的解可通过求解如下以离散形式给出的 迭代方程获得：

$\underset{u,d}{\arg \min }\underset{i,j}{\Sigma }\left|d_{\mathrm{ij}}\right|+\frac{\lambda }{2}\underset{i,j}{\Sigma }{(u_{\mathrm{ij}}-u_{0,\mathrm{ij}})}^2+\frac{\beta }{2}\underset{i,j}{\Sigma }{|d_{\mathrm{ij}}-{▿}_{\mathrm{NL}}{u}_{\mathrm{ij}}-b_{\mathrm{ij}}|}^2.---\left(8\right)$

其中b是与Bregman迭代算法相关的一个变量，也称作Bregman变量，β是 一个大于零的常数。方程(8)是一个典型的多元变量函数，通常采用分解法将 该方程分解为基于单个变量的子问题。

*关于变量u的子问题

固定变量d，则方程(8)简化为求关于变量u的最小化子问题：

$\underset{u,d}{\arg \min }\frac{\lambda }{2}\underset{i,j}{\Sigma }{(u_{\mathrm{ij}}-u_{0,\mathrm{ij}})}^2+\frac{\beta }{2}\underset{i,j}{\Sigma }{|d_{\mathrm{ij}}-{▿}_{\mathrm{NL}}{u}_{\mathrm{ij}}-b_{\mathrm{ij}}|}^2.---\left(9\right)$

通过Eular-Lagrange方程，方程(9)的最小解满足：

$\lambda (u-u_0)-\beta {\mathrm{div}}_{\mathrm{NL}}(d^k-{▿}_{\mathrm{NL}}u-b^k)=0.---\left(10\right)$

使用高效的Gauss-Seidel迭代法对式(10)进行求解，即对n≥0有：

$\left(\begin{array}{c}{u}_{0,i}^{k+1}=u_i^{k+1}\\ u_i^{k+1,n+1}=\frac{1}{\lambda +\beta {\Sigma }_j{w}_{\mathrm{ij}}}\left(\begin{array}{c}\beta {\Sigma }_j{w}_{\mathrm{ij}}{u}_j^{k+1,n}+{\mathrm{\lambda u}}_{0,i}-\\ \beta {\Sigma }_j\sqrt{w_{\mathrm{ij}}}(d_{\mathrm{ij}}^{k+1,n}-d_{\mathrm{ji}}^{k+1,n}-b_{\mathrm{ij}}^{k+1,n}-b_{\mathrm{ji}}^{k+1,n})\end{array},\right)\end{array}\right).---\left(11\right)$

实验结果表明，两次的Gauss-Seidel迭代就足够收敛到方程的最小解处。

*关于变量d的子问题

固定变量u，则方程(8)简化为求关于变量d的最小化子问题：

$\underset{u,d}{\arg \min }\underset{i,j}{\Sigma }\left|d_{\mathrm{ij}}\right|+\frac{\beta }{2}\underset{i,j}{\Sigma }{|d_{\mathrm{ij}}-{\mathrm{\Delta u}}_{\mathrm{ij}}{-b}_{\mathrm{ij}}|}^2.---\left(12\right)$

采用软阈值法，则方程(12)的Eular-Lagrange方程为：

$d_{\mathrm{ij}}^{k+1}=\frac{{▿}_{\mathrm{NL}}{u}_{\mathrm{ij}}^{k+1}+b_{\mathrm{ij}}^k}{|{▿}_{\mathrm{NL}}{u}_{\mathrm{ij}}^{k+1}+b_{\mathrm{ij}}^k|}\max \left(\right|{▿}_{\mathrm{NL}}{u}_{\mathrm{ij}}^{k+1}+b_{\mathrm{ij}}^k|-\beta ,0)$

$=\frac{\sqrt{w_{\mathrm{ij}}}(u_j^{k+1}-u_i^{k+1})+b_{\mathrm{ij}}^k}{\sqrt{{\Sigma }_j{w}_{\mathrm{ij}}{(u_j^{k+1}-u_i^{k+1})}^2+{b_{\mathrm{ij}}^k}^2}}\max (\sqrt{{\Sigma }_j{w}_{\mathrm{ij}}{(u_j^{k+1}-u_i^{k+1})}^2+{b_{\mathrm{ij}}^k}^2}-\beta ,0).---\left(13\right)$

最后，Bregman变量b可根据如下公式进行更新：

$b_{\mathrm{ij}}^{k+1}=b_{\mathrm{ij}}^k+\sqrt{w_{\mathrm{ij}}}(u_j^{k+1}-u_i^{k+1})-d_{\mathrm{ij}}^{k+1}.---\left(14\right)$

本发明对模拟数据和临床超声真实数据进行了滤波实验，结果表明本发明 能够有效的滤除超声图像中的斑点噪声，提高超声图像的成像清晰度。

综上所述，本发明根据经对数变换后的超声图像中斑点噪声数理统计特性， 获得了基于Pearson统计距离的权重计算函数，是一种非参数统计的方法，避免 了传统扩散方程中对参数的复杂设计，对斑点噪声取得了较好的滤波效果。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节， 而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现 本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非 限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落 在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权 利要求中的任何附图标记视为限制所涉及的权利要求。

此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施 方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见， 本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经 适当组合，形成本领域技术人员可以理解的其他实施方式。