CN106901776A - 基于可变滤波器长度的超声弹性成像方法 - Google Patents

Abstract

本发明提供了一种基于可变滤波器长度的超声弹性成像方法，首先根据变形前后的两帧RF信号，利用数值优化算法，计算出位移；然后通过传统固定滤波器长度的方法计算出初应变；再对该应变场进一步差分运算求出应变梯度，并对此线性变换得到可变的SG差分滤波器长度分布；最后，结合前面计算出的位移值与滤波器长度分布计算得到最终应变场。本发明获取的应变图中的运动边界更清晰，计算域内CNR值更高，具有更优的目标检测能力。

Description

基于可变滤波器长度的超声弹性成像方法

技术领域

本发明涉及超声弹性成像技术领域，具体涉及一种基于可变滤波器长度的超声弹性成像方法。

背景技术

生物组织的弹性特性，与组织的病理学状态紧密相关，具有重要的临床应用。超声弹性成像能无创定性定量生物组织的弹性，弥补了传统医学成像模态的不足，具有广阔的应用前景。应变计算是超声弹性成像技术的关键步骤，计算准确度将直接影响成像质量。从数学理论的角度，位移到应变的计算相当于一个数值微分的过程，即应变是位移的梯度。数值微分操作被当作一个数字差分滤波器。然而，数值微分操作是一个不稳定和危险的操作，因为它对信号的噪声或误差非常敏感，具有明显的放大作用，即很小的噪声或误差就会引起数值微分结果的很大误差。

Savitzky-Golay差分滤波器(SG差分滤波器)的基本思想是移动、滑动、局部、分段、逐点多项式拟合，具有计算过程清晰简洁，滤波器系数容易获取，滤波器长度可以设置为任意值等诸多优点，能很好的应用于弹性成像中的应变计算。其中，SG差分滤波器长度是影响应变计算结果的一个关键参数。如果滤波器长度设置过短，则会增加应变计算错误率。反之，则会降低应变弹性图的分辨率，即会降低目标检测能力。因此，设置一个最优的滤波器长度是在超声弹性成像中准确计算应变的一个关键因素。传统利用SG差分滤波器计算应变，对滤波器长度都是设置为固定值，这在力学线性环境下，能取得较为满意的结果。但实际上生物组织是力学非线性，特别是对于含有异物(如肿瘤)的生物组织，周围组织与异物之间弹性差值更明显，射频数据采集环境复杂，获取的被测区域射频数据在各处的噪声含量都不一样。此时，采用固定的滤波器长度，成像质量会不佳。

发明内容

本申请通过提供一种基于可变滤波器长度的超声弹性成像方法，以解决采用固定滤波器长度而造成的成像质量不佳的技术问题。

为解决上述技术问题，本申请采用以下技术方案予以实现：

一种基于可变滤波器长度的超声弹性成像方法，包括如下步骤：

S1：根据变形前后两帧超声射频信号，利用数值优化算法计算出位移场；

S2：利用滤波器长度为固定值的SG差分滤波器对位移场进行差分运算，得出初应变场；

S3：对初应变场进行差分运算求出含有异物边界信息的应变梯度场；

S4：通过映射函数对应变梯度场的绝对值进行线性变换，得出与异物边界信息相关的滤波器长度分布场；

S5：按照滤波器长度分布场，对步骤S1中求出的位移场进行SG差分滤波器差分运算，得出最终的应变场。

进一步地，步骤S1中的数值优化算法具体为：

假设图像窗口在运动前后满足仿射变换关系，运动前窗口中的任意一点P(x,y)对应于运动后窗口中点Q(x′,y′)，两者的坐标满足：

式中，向量为窗口的移动和变形，u，v为刚体位移，为窗口的应变，即窗口的一阶变形，将目标函数Ob(V_rs)最小化即可同时得到向量V_rs中的6个参数值，其中，目标函数

进一步地，步骤S2中采用1阶或2阶多项式拟合的SG差分滤波器进行应变计算，假设SG差分滤波器的长度为2M+1，则1阶或2阶多项式拟合的SG差分滤波器的输入输出关系为：

进一步地，步骤S4中的映射函数为：FL_var＝round{ζ[1-(1-σ)Φ]}，FL_var为滤波器长度，round()为四舍五入函数，式中，ζ为一常数，σ为调整因子，Φ为应变梯度归一化函数，式中，MB为应变梯度场，MB_max为应变梯度场的最大值，MB_min为应变梯度场的最小值。

可选地，调整因子σ取值为0.05～0.95。

与现有技术相比，本申请提供的技术方案，具有的技术效果或优点是：本发明获取的应变图中的运动边界更清晰，计算域内CNR值更高，具有更优的目标检测能力。

附图说明

图1(a)为弹性图轴线上的位移曲线图；

图1(b)为不同滤波器长度下的SG差分滤波器对应图1(a)的应变曲线图；

图2(a)为物理模型中A-line的位置示意图；

图2(b)为对应A-line的位移曲线图；

图2(c)为对应A-line的理想应变曲线图；

图2(d)为对应A-line的理想应变梯度曲线图；

图3为本发明的流程图；

图4(a)为单异物体模的侧向轴线上的应变曲线图；

图4(b)为单异物体模的轴向轴线上的应变曲线图；

图4(c)为单异物体模的CNR图；

图5(a)为离体组织的侧向轴线上的应变曲线图；

图5(b)为离体组织的轴向轴线上的应变曲线图；

图5(c)为离体组织的CNR图。

具体实施方式

本申请实施例通过提供一种基于可变滤波器长度的超声弹性成像方法，以解决采用固定滤波器长度而造成的成像质量不佳的技术问题。

为了更好的理解上述技术方案，下面将结合说明书附图以及具体的实施方式，对上述技术方案进行详细的说明。

实施例

噪声或误差一般表现为高频成分，理论分析与实验研究表明，采用低通数字差分器SG差分滤波器进行应变计算，能够取得较好的效果。

SG差分滤波器的属性与多项式的阶次以及滤波器长度相关。滤波器长度越长，多项式的阶次越低，滤波器的噪声放大因子(滤波器系数的平方和)越小。

本实施例采用1或2阶多项式拟合的SG差分滤波器进行应变计算研究。假设滤波器的长度为2M+1，1或2阶多项式拟合的SG差分滤波器的输入输出关系可以表示为：

采用滤波器长度较长的SG差分滤波器能有效的去除噪声，但同时降低了目标检测能力。假设一弹性图轴线上的位移曲线如图1(a)所示，利用不同滤波器长度下的SG差分滤波器对该曲线进行应变计算，结果如图1(b)所示。从图1(b)可看出，随着滤波器长度的增加，应变曲线的运动边界(即背景组织与异物之间的边界)就越平滑，即目标检测能力就越低。

从以上分析可知，SG差分滤波器长度是影响弹性成像质量的一个关键因素。不同于以往通过设置固定滤波器长度的方法计算应变，本发明通过获取的应变梯度作为反馈参数，并利用一个映射函数得到一个可变的滤波器长度分布场，再据此进行应变计算。所谓应变梯度，即是对位移的二次差分结果。定义超声数据图像轴向上的一条轴线为“A-line”，则A-line对应的位移、应变、应变梯度的演化过程如图2所示。位移、应变、应变梯度都包含有异物的位置信息，但应变梯度更能量化和定位异物的边界。

基于上述理论，本发明提出了一种基于可变滤波器长度的超声弹性成像方法，如图3所示，包括如下步骤：

S1：根据变形前后两帧超声射频信号，利用数值优化算法计算出位移场；

数值优化算法具体为：

假设图像窗口在运动前后满足仿射变换关系，运动前窗口中的任意一点P(x,y)对应于运动后窗口中点Q(x′,y′)，两者的坐标满足：

式中，向量为窗口的移动和变形，u，v为刚体位移，为窗口的应变，即窗口的一阶变形，将目标函数Ob(V_rs)最小化即可同时得到向量V_rs中的6个参数值，其中，目标函数

S2：利用滤波器长度为固定值的SG差分滤波器对位移场进行差分运算，得出初应变场；

即采用1阶或2阶多项式拟合的SG差分滤波器进行应变计算，假设SG差分滤波器的长度为2M+1，则1阶或2阶多项式拟合的SG差分滤波器的输入输出关系为：

S3：对初应变场进行差分运算求出含有异物边界信息的应变梯度场；

S4：通过映射函数对应变梯度场的绝对值进行线性变换，得出与异物边界信息相关的滤波器长度分布场；

映射函数为：FL_var＝round{ζ[1-(1-σ)Φ]}，FL_var为滤波器长度，round()为四舍五入函数，式中，ζ为一常数，σ为调整因子，Φ为应变梯度归一化函数，式中，MB为应变梯度场，MB_max为应变梯度场的最大值，MB_min为应变梯度场的最小值调整因子σ取值为0.05～0.95。

S5：按照滤波器长度分布场，对步骤S1中求出的位移场进行SG差分滤波器差分运算，得出最终的应变场。

此时的滤波器长度为可变值，即计算矩阵中，每个计算位置所对应的滤波器长度值为S4中求出的滤波器长度分布场中相对应位置的滤波器长度值。

由于对比度噪声比(contrast-to-noise ratio,CNR)结合了弹性图中的噪声和对比状态，本实施例选用该参数来评估目标检测能力。在背景组织包含有异物的模型的CNR可表示为：式中，μ_sb为背景组织的弹性值的均值，μ_st为异物的弹性值的均值，为背景组织的弹性值的标准差，为异物的弹性值的标准差。

为了进一步验证本发明的显著效果，接下来将通过体模实验和离体组织实验进行验证。

体模实验是在一个仿生物组织弹性体模上进行(model 049A,CIRS Inc.,Norfolk,VA,USA)。体模中背景组织的弹性模量为25±6kPa(均值±标准差)，内嵌异物的弹性模量为80±12kPa(均值±标准差)，内嵌异物的直径为6.5mm。

采用一个Philips iU22超声系统(Philips Medical Systems,Bothell,WA,USA)采集超声射频数据，采样频率为32MHz，轴向分辨率为0.0241mm，侧向分辨率为0.1187mm。获取持续压缩体模过程中的5帧射频数据，并选取其中两帧数据用于应变计算。

本实施例方法与传统固定滤波器长度方法在轴线上的应变曲线如图4所示。由图4(a)侧向轴线上的应变曲线和图4(b)轴向轴线上的应变曲线可看出，相比于固定滤波器长度，本实施例方法获取的应变曲线更平滑，边界处更锐利。图4(c)为分别基于固定滤波器长度方法和本实施例方法的CNR计算结果，基于可变滤波器长度的SG差分滤波器方法的CNR高于固定滤波器长度方法的计算值，图中固定滤波器长度为2M+1。

离体组织实验在离体的猪肝上进行，利用无水乙醇可以使蛋白质变性的原理来使肝脏组织产生损伤。本次实验中注射的无水乙醇(Bei Hua Fine Chemicals Co,Ltd,Beijing,China)浓度为99.9％，注射剂量约2mL，并认为注射约5分钟后，注射酒精部位肝脏组织蛋白质变性成功，形成损伤。

一个带有VF10-5探头的Sonoline Antares系统(Siemens Medical SolutionsUSA,Inc,Mountain View,CA)被用于本次实验。通过该系统的超声开放接口(ultrasoundresearch interface,URI)采集原始射频数据。在挤压的过程中同时以40MHz的采样率采集9帧连续的射频数据。本次试验中选取其中的两帧数据用于应变计算。

本实施例方法与传统固定滤波器长度方法在含乙醇损伤的离体猪肝上的轴线上的应变曲线以及计算域内的CNR计算结果如图5所示。由图5(a)侧向轴线上的应变曲线和图5(b)轴向轴线上的应变曲线可看出，相比于固定滤波器长度，本实施例方法获取的应变曲线更平滑，边界处更锐利。图5(c)为分别基于固定滤波器长度方法和本实施例方法的CNR计算结果，本实施例方法的CNR高于固定滤波器长度方法的计算值。

本申请的上述实施例中，通过提供一种基于可变滤波器长度的超声弹性成像方法，首先根据变形前后的两帧RF信号，利用数值优化算法，计算出位移；然后通过传统固定滤波器长度的方法计算出初应变；再对该应变场进一步差分运算求出应变梯度，并对此线性变换得到可变的SG差分滤波器长度分布；最后，结合前面计算出的位移值与滤波器长度分布计算得到最终应变场。本发明获取的应变图中的运动边界更清晰，计算域内CNR值更高，具有更优的目标检测能力。

应当指出的是，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的普通技术人员在本发明的实质范围内所做出的变化、改性、添加或替换，也应属于本发明的保护范围。

Claims (5)

1.一种基于可变滤波器长度的超声弹性成像方法，其特征在于，包括如下步骤：

S1：根据变形前后两帧超声射频信号，利用数值优化算法计算出位移场；

S2：利用滤波器长度为固定值的SG差分滤波器对位移场进行差分运算，得出初应变场；

S3：对初应变场进行差分运算求出含有异物边界信息的应变梯度场；

S4：通过映射函数对应变梯度场的绝对值进行线性变换，得出与异物边界信息相关的滤波器长度分布场；

S5：按照滤波器长度分布场，对步骤S1中求出的位移场进行SG差分滤波器差分运算，得出最终的应变场。

2.根据权利要求1所述的基于可变滤波器长度的超声弹性成像方法，其特征在于，步骤S1中的数值优化算法具体为：

假设图像窗口在运动前后满足仿射变换关系，运动前窗口中的任意一点P(x,y)对应于运动后窗口中点Q(x′,y′)，两者的坐标满足：

$\left\{\begin{array}{c}{x}^{\′}=x+\mathrm{\ξ}\left(x,y,V_{rs}\right)=x+u+\frac{\∂u}{\∂x}x+\frac{\∂u}{\∂y}y\\ y^{\′}=y+\mathrm{\η}\left(x,y,V_{rs}\right)=y+v+\frac{\∂v}{\∂x}x+\frac{\∂v}{\∂y}y\end{array}\right.$

式中，向量为窗口的移动和变形，u，v为刚体位移，为窗口的应变，即窗口的一阶变形，将目标函数Ob(V_rs)最小化即可同时得到向量V_rs中的6个参数值，其中，目标函数

$Ob\left(V_{rs}\right)=\underset{y=1}{\overset{L}{\Σ}}\underset{x=1}{\overset{N}{\Σ}}\[P(x,y)-Q(x+\mathrm{\ξ}(x,y,V_{rs}),y+\mathrm{\η}(x,y,V_{rs}\left)\right)\].$

3.根据权利要求1所述的基于可变滤波器长度的超声弹性成像方法，其特征在于，步骤S2中采用1阶或2阶多项式拟合的SG差分滤波器进行应变计算，假设SG差分滤波器的长度为2M+1，则1阶或2阶多项式拟合的SG差分滤波器的输入输出关系为：

4.根据权利要求1所述的基于可变滤波器长度的超声弹性成像方法，其特征在于，步骤S4中的映射函数为：FL_var＝round{ζ[1-(1-σ)Φ]}，FL_var为滤波器长度，round()为四舍五入函数，式中，ζ为一常数，σ为调整因子，Φ为应变梯度归一化函数，式中，MB为应变梯度场，MB_max为应变梯度场的最大值，MB_min为应变梯度场的最小值。

5.根据权利要求4所述的基于可变滤波器长度的超声弹性成像方法，其特征在于，调整因子σ取值为0.05～0.95。