CN111815527B - 基于Weibull分布的混合高阶变分超声图像去噪方法 - Google Patents

Abstract

基于Weibull分布的混合高阶变分超声图像去噪方法，对超声图像进行噪声分布拟合，采用Weibull概率密度函数，对选定的图像局部均匀区域进行灰度直方图拟合，并对拟合分布参数采用最大似然法进行估计；根据拟合的噪声分布，推导对数真实无噪图像的似然函数，将似然函数进行负的对数变换，从而得到最小化能量函数的数据保真项；将对数真实无噪图像一阶正则化表达式和二阶正则化表达式进行加权融合，构建最小化能量函数的混合高阶正则化项；将数据保真项和混合高阶正则化项形成最小化能量函数模型，采用Split‑Bregman迭代方法进行快速求解，迭代收敛后，利用指数变换得到去噪超声图像。本发明方法能有效降低散斑噪声的同时更好地保留图像的细节和边缘信息，从而增强超声图像的视觉判读性。

Description

基于Weibull分布的混合高阶变分超声图像去噪方法

技术领域

本发明涉及图像处理技术领域，具体涉及一种基于Weibull分布的混合高阶变分超声图像去噪方法。

背景技术

目前，超声影像已成为临床多种疾病诊断的首选方法，并已发展成为一种能提供多种参数的系列诊断技术。超声成像方法常用来判断脏器的位置、大小、形态，确定病灶的范围和物理性质，提供一些腺体组织的解剖图，鉴别胎儿的正常与异常，在妇产科及心血管系统、消化系统、泌尿系统的应用十分广泛。超声成像采用回声原理，由仪器探头向人体发射一束超声进入体内，并进行线形、扇形或其它形式的扫描，遇到不同声阻抗的两种组织的交界面，即有超声反射回来，由探头接收后，经过信号放大或信息处理，形成人体的断层图像。除了反射回波外，人体组织的微小结构作为散射体使声波产生散射后形成的散射回波也会被接收，但这部分回波相互干涉通常产生所谓的散斑噪声。由于该类噪声的存在，使得图像分辨率和对比度严重下降，影响计算机辅助的临床诊断，并为后续的组织或病灶体分割、识别、配准或三维重建等造成极大困难。因此超声图像散斑噪声去除一直是医学成像领域研究的热点问题。

多年来，国内外学者致力于超声图像的去噪研究，提出了多种不同类型的去噪方法，包括局部空域滤波、各项异性扩散滤波、非局部均值滤波、总变分方法、多尺度方法以及同态滤波方法等。这些方法，有的去噪效果不够理想，有的去噪算法复杂度较高。因此，寻求一种算法简单有效且去噪效果较理想的方案，是目前超声图像实时处理的必要条件，也是去噪算法移植到实际医疗设备的先决条件。

发明内容

本发明提出了一种基于Weibull分布的混合高阶变分超声图像去噪方法，通过拟合超声图像均匀区域得到整体噪声分布参数，以此构建最小化能量函数的数据保真项；并设计自适应参数以及混合高阶变分得到正则化项。两者共同作用，能有效降低散斑噪声的同时更好地保留图像的细节和边缘信息，从而增强超声图像的视觉判读性。

本发明采取的技术方案为：

基于Weibull分布的混合高阶变分超声图像去噪方法，包括以下步骤：

步骤1：对超声图像进行噪声分布拟合，采用Weibull概率密度函数，对选定的图像局部均匀区域进行灰度直方图拟合，并对拟合分布参数采用最大似然法进行估计；

步骤2：根据拟合的噪声分布，推导对数真实无噪图像的似然函数，将似然函数进行负的对数变换，从而得到最小化能量函数的数据保真项；

步骤3：将对数真实无噪图像一阶正则化表达式和二阶正则化表达式进行加权融合，构建最小化能量函数的混合高阶正则化项，且权值能根据图像边缘特征自适应调整；

步骤4：将步骤2的数据保真项和步骤3的混合高阶正则化项形成最小化能量函数模型，采用Split-Bregman迭代方法进行快速求解，迭代收敛后，得到对数去噪超声图像，最后通过指数变换得到去噪超声图像。

本发明一种基于Weibull分布的混合高阶变分超声图像去噪方法，技术效果如下：

1)、本发明利用图像的噪声分布参数构建数据保真项，该项可根据超声图像的噪声强度自适应改变；采用混合高阶变分项并设计自适应参数，在有效去除噪声时能更精确地描述边缘信息以及微小细节信息，同时避免一阶变分计算时产生的阶梯块效应，且引入更少的参数以减少算法复杂度。综合两项的优势使去噪算法能够在均匀区域有效降低散斑噪声，在细节区域有效保留或增强对比度。

2)、本发明通过图像自身灰度统计、自适应参数设计、高阶变分混合的形式对图像进行滤波，引入更少的参数降低算法难度。且有效降低噪声以及一阶变分框架阶梯块效应，更好地保留图像细节和边缘信息，从而增强超声图像的视觉判读性。

3)、本发明方法采用Split-Bregman迭代求解算法，收敛速度快，适宜进行硬件固化。

附图说明

图1为本发明设计实现流程图。

图2为本发明实施例中的待去噪超声图。

图3为本发明实施例中的去噪结果图。

图4为本发明方法对比三幅不同噪声强度测试图像的去噪结果的峰值信噪比(Peak Signal to Noise Ratio,PSNR)随数据保真项权系数α的变化曲线图。

图5为本发明方法对比三幅不同噪声强度测试图像去噪结果的峰值信噪比PSNR随调节系数θ₁的变化曲线图。

图6为本发明方法对比三幅不同噪声强度测试图像去噪结果的峰值信噪比PSNR随调节系数θ₂的变化曲线图。

图7为去噪结果对比的局部区域示意图。

图8(1)为原始超声图像的局部区域；

图8(2)为Exp超声图像去噪方法的处理结果图；

图8(3)为MIDAL超声图像去噪方法的处理结果图；

图8(4)为NTV超声图像去噪方法的处理结果图；

图8(5)为Huang-Yang超声图像去噪方法的处理结果图；

图8(6)为TM超声图像去噪方法的处理结果图；

图8(7)为本发明WHHVM超声图像去噪方法的处理结果图。

具体实施方式

基于Weibull分布的混合高阶变分超声图像去噪方法，包括以下步骤：

步骤1：对超声图像进行噪声分布拟合，采用Weibull概率密度函数，对选定的图像局部均匀区域进行灰度直方图拟合，并对拟合分布参数采用最大似然法进行估计。

所述步骤1中，真实无噪图像v，受到乘性噪声n的干扰，生成观察到的超声图像f，且满足超声图像成像模型f＝vn，其中，乘性噪声n在图像的局部灰度均匀区域近似为Weibull分布，即噪声的概率密度函数满足：

其中，exp()为指数函数；λ和k分别为Weibull分布的尺度参数和形状参数，选取超声图像的局部灰度均匀区域，获得直方图拟合的样本数据，利用最大似然法对尺度参数λ和形状参数k进行估计。

步骤2：根据拟合的噪声分布，推导对数真实无噪图像的似然函数，将似然函数进行负的对数变换，从而得到最小化能量函数的数据保真项。

所述步骤2中，对超声图像成像模型f＝vn，进行对数变换可得logf＝logv+logn。令g＝logf，u＝logv和w＝logn，则有g＝u+w。根据乘性噪声n的概率密度函数p_N(n)推导对数变换后的噪声w满足：

进一步，求得对数变换后的对数真实无噪图像u的似然函数为：

将似然函数进行负的对数变换，且忽略常数部分，可得：

-logp_G|U(g|u)＝-k(g-u)+λ^-ke^k(g-u)；

假定图像区域Ω内的所有像素g(x,y)∈Ω满足独立同分布，则能构建最小化能量函数的数据保真项如下：

步骤3：将对数真实无噪图像一阶正则化表达式和二阶正则化表达式进行加权融合，构建最小化能量函数的混合高阶正则化项，且权值能根据图像边缘特征自适应调整。

所述步骤3中，将对数真实无噪图像u的一阶正则化表达式和二阶正则化表达式/>进行加权组合，得到混合高阶正则化项为：

其中，分别代表x和y方向的偏导数，·代表向量或矩阵的Euclid范数，加权系数γ₁、γ₂根据对数真实无噪图像u的边缘特征自适应变化，具体表达式为：

其中，η为控制函数平滑程度的阈值，取值范围为0.1～0.9，本发明实施例中取值为0.65。步骤4：将步骤2的数据保真项和步骤3的混合高阶正则化项形成最小化能量函数模型，采用Split-Bregman迭代方法进行快速求解，迭代收敛后，得到对数去噪超声图像，最后通过指数变换得到去噪超声图像。

所述步骤4中，将步骤2的数据保真项和步骤3的混合高阶正则化项组合，引入数据保真项权系数α，得到最小化能量函数模型为：

采用Split-Bregman迭代方法，对最小化能量函数模型进行求解，具体过程如下：

a)：引入辅助分离变量以及与其分别对应解析空间下的迭代参量：

则最小化能量函数模型E(u)可变化为：

其中，θ₁、θ₂是算法求解时引入的调节系数。

b)：为求u的最优解，首先将函数对u求偏导，得出欧拉-拉格朗日方程：

令偏导为0，整理后有：

其中，div(·)、div²(·)分别为一、二阶散度算子。

根据图像偏微分以及网格离散差分运算，有：

其中，上标“+”代表前向差分格式，上标“-”代表后向差分格式。

将上述离散化表达式，代入偏微分方程中，利用Gauss-Seidel迭代算法，可以得出：

其中，等的上标表示迭代次数m，下标表示图像u在像素点(i,j)及其周围的灰度值。本发明实施例中数据保真项权系数取值为α＝0.2，调节系数分别取值为θ₁＝200，θ₂＝30。

c)：根据b)中得到的结果，更新步骤3中的自适应加权系数γ₁、γ₂。

d)：求解最小化能量函数模型中分离变量q,s，参考步骤b)，分别将函数对q,s求偏导并置零：

采用阈值收缩算法求得：

e)：得到q^m+1,s^m+1后，即可对迭代参量进行更新：

f)：得到第m+1次迭代图像u^m+1后，根据前后两次迭代的滤波图像均方差小于给定的阈值，即||u^m+1-u^m||≤ξ来判定迭代收敛，阈值ξ典型取值为0.01；迭代收敛的图像u^m+1为对数变换后的对数去噪超声图像，最后利用指数变换得到去噪超声图像v，即v＝exp(u^m+1)。

图2为本发明实施例所用的肌肉组织超声图像，可以看到由于散斑噪声的影响，肌肉组织的细节比较模糊。图3为本发明实施例中的肌肉组织超声图像去噪结果，可以看到去噪后的图像有效减弱了散斑噪声，同时保留了肌肉组织结构以及纹理细节，视觉上获得有效增强。

由图4～图6可以看出本发明方法在三种不同强度噪声情况下，模型中的数据保真项权系数α对模型影响较大，噪声越强α取值应该适当增大以增大数据保真项权重；调节系数θ₁、θ₂对模型影响较小。

图8(2)～图8(7)为不同超声图像去噪方法的处理结果对比图：包括Exp(Exponential Model,指数模型)、MIDAL(Multiplicative Image Denoising byAugmented Lagrangian，增广拉格朗日乘性图像去噪)、NTV(Nakagami Total Variation，Nakagami变分去噪)、Huang-Yang(Huang-Yang变分模型)、TM(Taylor’s expansionapproximation Model,泰勒展开近似模型)以及本发明的WHHVM(Weibull-based HybridHigh-order Variational Model，)基于Weibull分布的混合高阶变分模型方法。由图8(2)～图8(7)可以看出：MIDAL、Huang-Yang和TM的滤波结果仍有少量的散斑噪声残留，Exp和NTV处理出现了伴生的伪影现象，而本发明WHHVM方法去噪同时保留了完整的细节，具有最好的视觉效果。

Claims (1)

1.基于Weibull分布的混合高阶变分超声图像去噪方法，其特征在于包括以下步骤：

步骤1：对超声图像进行噪声分布拟合，采用Weibull概率密度函数，对选定的图像局部均匀区域进行灰度直方图拟合，并对拟合分布参数采用最大似然法进行估计；

步骤2：根据拟合的噪声分布，推导对数真实无噪图像的似然函数，将似然函数进行负的对数变换，从而得到最小化能量函数的数据保真项；

步骤3：将对数真实无噪图像一阶正则化表达式和二阶正则化表达式进行加权融合，构建最小化能量函数的混合高阶正则化项，且权值能根据图像边缘特征自适应调整；

步骤4：将步骤2的数据保真项和步骤3的混合高阶正则化项形成最小化能量函数模型，采用Split-Bregman迭代方法进行快速求解，迭代收敛后，得到对数去噪超声图像，最后通过指数变换得到去噪超声图像；

所述步骤1中，真实无噪图像v，受到乘性噪声n的干扰，生成观察到的超声图像f，且满足超声图像成像模型f＝vn，其中，乘性噪声n在图像的局部灰度均匀区域近似为Weibull分布，即噪声的概率密度函数满足：

其中，exp()为指数函数；λ和k分别为Weibull分布的尺度参数和形状参数，选取超声图像的局部灰度均匀区域，获得直方图拟合的样本数据，利用最大似然法对尺度参数λ和形状参数k进行估计；

所述步骤2中，对超声图像成像模型f＝vn，进行对数变换可得logf＝logv+logn；令g＝logf，u＝logv和w＝logn，则有g＝u+w；根据乘性噪声n的概率密度函数p_N(n)推导对数变换后的噪声w满足：

进一步，求得对数变换后的对数真实无噪图像u的似然函数为：

将似然函数进行负的对数变换，且忽略常数部分，可得：

-logp_G|U(g|u)＝-k(g-u)+λ^-ke^k(g-u)；

假定图像区域Ω内的所有像素g(x,y)∈Ω满足独立同分布，则能构建最小化能量函数的数据保真项如下：

所述步骤3中，将对数真实无噪图像u的一阶正则化表达式和二阶正则化表达式/>进行加权组合，得到混合高阶正则化项为：

其中， 分别代表x和y方向的偏导数，||·||代表向量或矩阵的Euclid范数，加权系数γ₁、γ₂根据对数真实无噪图像u的边缘特征自适应变化，具体表达式为：

γ₂＝min(γ₁)+max(γ₁)-γ₁；

其中，η为控制函数平滑程度的阈值，目标轮廓越明显取值越大；

所述步骤4中，将步骤2的数据保真项和步骤3的混合高阶正则化项组合，引入数据保真项权系数α，得到最小化能量函数模型为：