CN102920485A - 一种超声弹性成像中生物组织二维位移场的估计方法 - Google Patents

Abstract

本发明公开了一种超声弹性成像中生物组织二维位移场的估计方法，包括：(1)采集超声图像并估算纵向位移分量；(2)通过生物组织力学模型建立系统离散方程；(3)通过H∞滤波估算出生物组织的二维的位移场。本发明提出二维位移场估计方法是将现有所得的纵向位移测量值在生物力学模型约束下经H∞滤波算法得到位移场的最优估计，在有效滤除系统内外部噪声的同时，恢复出高质量的侧向位移分量，用于应变计算和弹性参数的重建；通过与现有方法进行比较，本发明方法在超声弹性成像侧向位移估计方面将会起到非常有效的作用。

Description

一种超声弹性成像中生物组织二维位移场的估计方法

技术领域

本发明属于超声成像技术领域，具体涉及一种超声弹性成像中生物组织二维位移场的估计方法。

背景技术

超声弹性成像技术是近年兴起的用以反应生物组织内部弹性特征的新型超声成像技术，它利用超声波RF信号来追踪生物组织在受到外力挤压前后的形变和运动信息，通过一定的算法得到生物组织的应变、弹性模量等力学特征。由于生物组织病变与其力学特征密切相关，如常见的乳腺癌病变区域通常要比正常组织硬，因此超声弹性成像可以为临床医学诊断提供重要的辅助信息，具有非常广阔的应用前景。

超声准静态弹性成像的概念最早由Ophir等人于1991年首次提出，其目的是借助超声成像仪对生物组织的硬度信息、弹性参数进行成像显示。它的主要步骤包括：(1)用超声探头记录待测生物组织在某一时刻的一帧超声信号；(2)利用超声探头挤压待测组织，同时记录组织受挤压后的超声信号；(3)比对所记录的两帧超声信号，采用位移估计算法(如经典的互相关算法)计算出待测组织在受挤压后的位移场；(4)通过对组织位移场的分析，计算出最终的弹性图(elastogram，可以是应变图也可以是弹性模量图)。

超声弹性成像的一般方法是利用超声探头挤压生物组织，同时记录返回的超声RF信号进行成像。由于超声信号在声波传播的方向(纵向)分辨率要明显高于垂直其传播方向(侧向)，因此传统超声弹性成像中，都仅仅计算位移场中的纵向分量，而侧向位移分量往往由于分辨率低、噪声大而被忽略。但是生物组织在受到探头挤压的运动并不是简单的一维运动，在超声成像平面中，它还伴随着不可忽略的侧向运动，因此如何准确地计算出侧向位移分量也开始收到研究人员的关注。为了使一维位移分量扩展到二维平面位移场，通常将用于一维纵向位移估计的方法(主要是基于信号窗口追踪的方法)应用到二维追踪，这些方法最早应用于血流成像技术。但是在超声弹性成像中，由于超声信号在侧向的空间分辨率明显低于纵向空间分辨率，使得在侧向方向需要高强度的数值插值。因此到目前为止，二维位移场的估计仍然受到超声探头硬件条件的制约，以至于侧向位移分量的准确度和分辨率并不能得到非常满意的结果。

发明内容

针对现有技术所存在的上述技术缺陷，本发明提供了一种超声弹性成像中生物组织二维位移场的估计方法，能够在有效滤除系统内外部噪声的同时，恢复出高质量的侧向位移分量，得到具有较好分辨率和信噪比的二维位移场。

一种超声弹性成像中生物组织二维位移场的估计方法，包括如下步骤：

(1)采用超声探头通过挤压对生物组织进行探测，连续采集得到多帧超声图像，并估算出每帧超声图像对应的一维纵向位移分量；

(2)根据生物组织的材质特性，建立生物组织力学模型的动态方程；对所述的动态方程进行状态空间转换以及离散化，得到生物组织的系统离散方程；

(3)引入与所述的系统离散方程对应的测量方程，根据系统离散方程、测量方程以及各帧超声图像的一维纵向位移分量，通过H滤波迭代算法估算出生物组织的二维位移场。

所述的步骤(1)中，采用时域位移估计法或相位域位移估计法估算出每帧超声图像对应的一维纵向位移分量。时域位移估计法是通过计算两个信号的时延值来计算位移的，而相位域算法则是基于相位差的来估计位移的。

所述的生物组织力学模型的动态方程如下：

MZ+CV+KU＝R

其中：M为生物组织的质量矩阵(其与生物组织的密度有关)，C为生物组织的阻尼矩阵且C＝αM+βK，α和β均为权重系数(均为给定的实际经验值且取值范围为0.01～0.05)，K为生物组织的刚度矩阵(其与生物组织的材质有关)，R为外力载荷向量(其与超声探头的挤压用力有关)，U为生物组织的位移矩阵，V为生物组织的速度矩阵，Z为生物组织的加速度矩阵。

所述的步骤(2)中，对动态方程进行状态空间转换以及离散化的过程为：首先，对动态方程进行状态空间转换，得到生物组织的时域状态方程如下：

$\stackrel{\·}{x}\left(t\right)=a\left(t\right)x\left(t\right)+\mathrm{bw}\left(t\right)$

其中： $\stackrel{\·}{x}\left(t\right)=\left[\begin{array}{c}V\left(t\right)\\ Z\left(t\right)\end{array}\right],$ $x\left(t\right)=\left[\begin{array}{c}U\left(t\right)\\ V\left(t\right)\end{array}\right],$ $a\left(t\right)=\left(\begin{array}{cc}0& I\\ -M^{-1}K\left(t\right)& -M^{-1}C\left(t\right)\end{array}\right),$ $b=\left(\begin{array}{cc}0& 0\\ 0& M^{-1}\end{array}\right),$ $w\left(t\right)=\left[\begin{array}{c}0\\ R\left(t\right)\end{array}\right],$ C(t)为t时刻生物组织的阻尼矩阵，K(t)为t时刻生物组织的刚度矩阵，w(t)为t时刻生物组织的力学向量，R(t)为t时刻生物组织的外力载荷向量，U(t)为t时刻生物组织的位移矩阵，V(t)为t时刻生物组织的速度矩阵，Z(t)为t时刻生物组织的加速度矩阵，I为单位矩阵；

然后，对生物组织的时域状态方程进行离散化，得到生物组织的系统离散方程如下：

x_i+1＝A_ix_i+B_iw_i+n_i

其中：x_i和x_i+1分别为第i帧和第i+1帧超声图像对应的位移矩阵，A_i和B_i均为第i帧超声图像采样时刻的参数矩阵，且A_i＝e^a(t)T，B_i＝a(t)^-1(A_i-I)b，T为采样间隔，n_i和w_i分别为第i帧超声图像采样时刻的过程噪声矩阵和力学向量。

所述的测量方程为y_i＝Dx_i+v_i；其中，y_i为第i帧超声图像对应的一维纵向位移分量，D为测量矩阵(其与超声探头设备参数有关)，v_i为第i帧超声图像采样时刻的测量噪声矩阵。

所述的步骤(3)中，所述的H∞滤波迭代算法基于以下公式：

x_i+1＝A_ix_i+B_iw_i+H_i(y_i-Dx_i)

$H_i=A_i{P}_i{S}_i{D}^T{v}_i^{-1}$

$S_i={(I-{\mathrm{\γ}}^{-2}{Q}_i{P}_i+D^T{v}_i^{-1}{\mathrm{DP}}_i)}^{-1}$

$P_{i+1}=A_i{P}_i{S}_i{A}_i^T+n_i$

其中，Q_i和P_i分别为第i帧超声图像采样时刻估计误差和初始误差的权重矩阵，γ为约束衰减系数(为给定的实际经验值)。Q_i为给定的实际经验值，P₀和x₀均为随机初始化确定。

所述的步骤(3)中，所述的H∞滤波迭代算法的收敛条件如下：

$\left|\right|\frac{x_{i+1}-x_i}{x_i}\left|\right|<\mathrm{\&rho;}$

其中：ρ为收敛阈值(为给定的实际经验值)。

本发明的有益技术效果：

(1)本发明由于采用生物力学模型作为约束条件，使得估计结果是具有生理意义；

(2)本发明由于采用H∞无穷滤波算法，不在乎具体的噪声分布，使得求解具备鲁棒性，更符合实际应用情况；

(3)本发明由于采用了状态空间解法，其突破了传统的方法在准静态超声弹性成像技术中的能够得到具有较高分辨率和准确度的二维位移场。

附图说明

图1为本发明估计方法的步骤流程示意图。

图2为二维矩阵材料模型的示意图。

图3(a)为纵向位移分量的参考真值图。

图3(b)为侧向位移分量的参考真值图。

图4(a)为采用最小二乘位移估计法得到的侧向位移分量的参考真值图。

图4(b)为采用本发明方法得到的侧向位移分量的参考真值图。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

如图1所示，一种超声弹性成像中生物组织二维位移场的估计方法，包括如下步骤：

(1)采集超声图像并估算纵向位移分量。

采用超声探头(一般采用线性探头)缓慢挤压生物组织表面，对其进行探测，采集到一系列连续多帧超声图像以记录生物组织的内部运动过程；本实施方式中，超声探头的中心频率为5MHz，采样频率为20MHz。

本实施方式采用时域位移估计法对超声图像纵向位移分量进行估计：将每一帧超声图像按照A-line线在横向方向分成N(A-line线数)条一维信号。对于每条A-line线，我们分成一定数量的相互之间有重叠的信号段，称之为信号窗口。由于组织在受挤压后会产生一定的位移，因此挤压前后两帧超声信号的两个对应信号窗口会在时域上产生时延τ(t)。假定超声在生物组织中的传播速度为c，那么生物组织任一网格节点(信号窗口所在的物理坐标)纵向位移u(t)可以通过公式u(t)＝cτ(t)/2得到。

故结合上述位移估计算法，从所采集的超声图像序列中可估算出生物组织的内部纵向位移场，即每帧超声图像对应的一维纵向位移分量。

(2)通过生物组织力学模型建立系统离散方程。

根据的生物组织特性，设定材料参数，建立生物组织的力学模型如下，模拟生物组织内部运动状态；

σ_ij＝2με_ij+λδ_ijε_nn

其中：μ＝E/2(1+υ)，λ为Lame常数且λ＝Eυ/(1+υ)(1-2υ)，E是材料的杨氏模量，υ是泊松比；ε_ij和ε_nn都是无穷小应变张量的分量，σ_ij是应力张量分量，δ_ij是Kronecher’s delta函数。

根据生物组织力学模型，通过有限元方法建立生物组织力学模型的动态方程：MZ+CV+KU＝R

其中：M为生物组织的质量矩阵其与生物组织的密度有关，C为生物组织的阻尼矩阵且C＝αM+βK，α和β均为权重系数，本实施方式中α和β均为0.01，K为生物组织的刚度矩阵其与生物组织的材质有关，R为外力载荷向量其与超声探头的挤压用力有关，U为生物组织的位移矩阵，V为生物组织的速度矩阵，Z为生物组织的加速度矩阵。

对生物组织力学模型的动态方程进行状态空间转换，得到生物组织的时域状态方程如下：

$\stackrel{\&CenterDot;}{x}\left(t\right)=a\left(t\right)x\left(t\right)+\mathrm{bw}\left(t\right)$

其中： $\stackrel{\&CenterDot;}{x}\left(t\right)=\left[\begin{array}{c}V\left(t\right)\\ Z\left(t\right)\end{array}\right],$ $x\left(t\right)=\left[\begin{array}{c}U\left(t\right)\\ V\left(t\right)\end{array}\right],$ $a\left(t\right)=\left(\begin{array}{cc}0& I\\ -M^{-1}K\left(t\right)& -M^{-1}C\left(t\right)\end{array}\right),$ $b=\left(\begin{array}{cc}0& 0\\ 0& M^{-1}\end{array}\right),$ $w\left(t\right)=\left[\begin{array}{c}0\\ R\left(t\right)\end{array}\right],$ C(t)为t时刻生物组织的阻尼矩阵，K(t)为t时刻生物组织的刚度矩阵，w(t)为t时刻生物组织的力学向量，R(t)为t时刻生物组织的外力载荷向量，U(t)为t时刻生物组织的位移矩阵，V(t)为t时刻生物组织的速度矩阵，Z(t)为t时刻生物组织的加速度矩阵，I为单位矩阵。

对生物组织的时域状态方程进行离散化，得到生物组织的系统离散方程如下：

x_i+1＝A_ix_i+B_iw_i+n_i

其中：x_i和x_i+1分别为第i帧和第i+1帧超声图像对应的位移矩阵，A_i和B_i均为第i帧超声图像采样时刻的参数矩阵，且A_i＝e^a(t)T，B_i＝a(t)^-1(A_i-I)b，T为采样间隔，n_i和w_i分别为第i帧超声图像采样时刻的过程噪声矩阵和力学向量。

(3)通过H∞滤波估算出生物组织的二维的位移场。

引入与系统离散方程对应的测量方程：y_i＝Dx_i+v_i；其中，y_i为第i帧超声图像对应的一维纵向位移分量，D为测量矩阵其与超声探头设备参数有关，v_i为第i帧超声图像采样时刻的测量噪声矩阵。

根据生物组织的系统离散方程和测量方程 $\left\{\begin{array}{c}{x}_{i+1}=A_i{x}_i+B_i{w}_i+n_i\\ y_i={\mathrm{Dx}}_i+v_i\end{array}\right.$ 以及各帧超声图像对应的一维纵向位移分量，通过以下H∞滤波迭代算法估算出生物组织的二维位移场。

x_i+1＝A_ix_i+B_iw_i+H_i(y_i-Dx_i)

$H_i=A_i{P}_i{S}_i{D}^T{v}_i^{-1}$

$S_i={(I-{\mathrm{\&gamma;}}^{-2}{Q}_i{P}_i+D^T{v}_i^{-1}{\mathrm{DP}}_i)}^{-1}$

$P_{i+1}=A_i{P}_i{S}_i{A}_i^T+n_i$

其中，Q_i和P_i分别为第i帧超声图像采样时刻估计误差和初始误差的权重矩阵，γ为约束衰减系数，本实施方式中γ＝1。Q_i为根据实际经验给定且本实施方式中每次迭代均采用相同的权重矩阵Q_i，P₀和x₀均为随机初始化确定。

当满足以下迭代收敛条件，终止迭代算法，输出x_i+1作为生物组织的二维位移场；

$\left|\right|\frac{x_{i+1}-x_i}{x_i}\left|\right|<\mathrm{\&rho;}$

其中：ρ为收敛阈值，本实施方式中ρ＝10^-5。

以下通过计算机模拟实验将本实施方式与基于最小二乘的位移估计方法(其利用组织不可压缩性特征，采用最小二乘进行侧向位移估计的方法)进行了比较。通过有限元分析软件ABAQUS建立一个20mm×20mm的二维矩阵材料模型，如图2所示，中心处包含一个直径为7.5mm的圆形硬块，中心圆硬度为80kPa，背景材料硬度为80kPa。设定合理的外力载荷等参数，导出模型的位移场(包括纵向位移分量和侧向位移分量)作为实验的参照真值，如图3所示了纵向位移分量和侧向位移分量的参照真值。实验中，通过在纵向位移场参照真值上加入不同程度不同类型的噪声，作为位移测量值使用。

应用本实施方式估计模型收超声探头挤压下的位移场，然后将位移场减去纵向位移分量得到模型的侧向位移分量，如图4(b)所示；而采用最小二乘估计方法计算得到的侧向位移分量，如图4(a)所示；可以看到应用本实施方式技术所得到的实验结果相对更加理想。

Claims (7)

1.一种超声弹性成像中生物组织二维位移场的估计方法，包括如下步骤：

(1)采用超声探头通过挤压对生物组织进行探测，连续采集得到多帧超声图像，并估算出每帧超声图像对应的一维纵向位移分量；

(2)根据生物组织的材质特性，建立生物组织力学模型的动态方程；对所述的动态方程进行状态空间转换以及离散化，得到生物组织的系统离散方程；

(3)引入与所述的系统离散方程对应的测量方程，根据系统离散方程、测量方程以及各帧超声图像的一维纵向位移分量，通过H∞滤波迭代算法估算出生物组织的二维位移场。

2.根据权利要求1所述的超声弹性成像中生物组织二维位移场的估计方法，其特征在于：所述的步骤(1)中，采用时域位移估计法或相位域位移估计法估算出每帧超声图像对应的一维纵向位移分量。

3.根据权利要求1所述的超声弹性成像中生物组织二维位移场的估计方法，其特征在于：所述的生物组织力学模型的动态方程如下：

MZ+CV+KU＝R

其中：M为生物组织的质量矩阵，C为生物组织的阻尼矩阵且C＝αM+βK，α和β均为权重系数，K为生物组织的刚度矩阵，R为外力载荷向量，U为生物组织的位移矩阵，V为生物组织的速度矩阵，Z为生物组织的加速度矩阵。

4.根据权利要求1所述的超声弹性成像中生物组织二维位移场的估计方法，其特征在于：所述的步骤(2)中，对动态方程进行状态空间转换以及离散化的过程为：首先，对动态方程进行状态空间转换，得到生物组织的时域状态方程如下：

$\stackrel{\&CenterDot;}{x}\left(t\right)=a\left(t\right)x\left(t\right)+\mathrm{bw}\left(t\right)$

其中： $\stackrel{\&CenterDot;}{x}\left(t\right)=\left[\begin{array}{c}V\left(t\right)\\ Z\left(t\right)\end{array}\right],$ $x\left(t\right)=\left[\begin{array}{c}U\left(t\right)\\ V\left(t\right)\end{array}\right],$ $a\left(t\right)=\left(\begin{array}{cc}0& I\\ -M^{-1}K\left(t\right)& -M^{-1}C\left(t\right)\end{array}\right),$ $b=\left(\begin{array}{cc}0& 0\\ 0& M^{-1}\end{array}\right),$ $w\left(t\right)=\left[\begin{array}{c}0\\ R\left(t\right)\end{array}\right],$ C(t)为t时刻生物组织的阻尼矩阵，K(t)为t时刻生物组织的刚度矩阵，w(t)为t时刻生物组织的力学向量，R(t)为t时刻生物组织的外力载荷向量，U(t)为t时刻生物组织的位移矩阵，V(t)为t时刻生物组织的速度矩阵，Z(t)为t时刻生物组织的加速度矩阵，I为单位矩阵；

然后，对生物组织的时域状态方程进行离散化，得到生物组织的系统离散方程如下：

x_i+1＝A_ix_i+B_iw_i+n_i

其中：x_i和x_i+1分别为第i帧和第i+1帧超声图像对应的位移矩阵，A_i和B_i均为第i帧超声图像采样时刻的参数矩阵，且A_i＝e^a(t)T，B_i＝a(t)^-1(A_i-I)b，T为采样间隔，n_i和w_i分别为第i帧超声图像采样时刻的过程噪声矩阵和力学向量。

5.根据权利要求1所述的超声弹性成像中生物组织二维位移场的估计方法，其特征在于：所述的测量方程为y_i＝Dx_i+v_i；其中，y_i为第i帧超声图像对应的一维纵向位移分量，D为测量矩阵，v_i为第i帧超声图像采样时刻的测量噪声矩阵。

6.根据权利要求1所述的超声弹性成像中生物组织二维位移场的估计方法，其特征在于：所述的步骤(3)中，所述的H∞滤波迭代算法基于以下公式：

x_i+1＝A_ix_i+B_iw_i+H_i(y_i-Dx_i)

$H_i=A_i{P}_i{S}_i{D}^T{v}_i^{-1}$

$S_i={(I-{\mathrm{\&gamma;}}^{-2}{Q}_i{P}_i+D^T{v}_i^{-1}{\mathrm{DP}}_i)}^{-1}$

$P_{i+1}=A_i{P}_i{S}_i{A}_i^T+n_i$

其中：Q_i和P_i分别为第i帧超声图像采样时刻估计误差和初始误差的权重矩阵，γ为约束衰减系数。

7.根据权利要求6所述的超声弹性成像中生物组织二维位移场的估计方法，其特征在于：所述的H∞滤波迭代算法的收敛条件如下：

$\left|\right|\frac{x_{i+1}-x_i}{x_i}\left|\right|<\mathrm{\&rho;}$

其中：ρ为收敛阈值。

Images