CN111388008B - 一种超声热成像技术的自适应运动补偿方法及系统 - Google Patents

Abstract

本发明提供一种超声热成像技术的自适应运动补偿方法及系统。其中，涉及的一种超声热成像技术的自适应运动补偿方法，包括步骤：S1.对超声热成像的目标区域进行分块，并选取运动补偿的局域参考点；S2.计算运动补偿的滤波系数；S3.根据所述选取的参考点以及计算得到的运动补偿的滤波系数，对温度分布进行运动补偿。本发明与已公开的超声热成像技术的运动补偿方案相比，充分考虑了人体自由呼吸和血液脉动等因素对温度估计的影响，通过对热消融目标区域的分块化处理和设置运动补偿滤波系数具有时变特性，提高了热消融过程中温度估计的准确性。此外，通过采用最速梯度下降算法进一步优化了算法的运行效率。

Description

一种超声热成像技术的自适应运动补偿方法及系统

技术领域

本发明涉及超声监控成像的技术领域，尤其涉及一种超声热成像技术的自适应运动补偿方法及系统。

背景技术

超声热成像技术是热消融过程中无创监测和温度控制的一种有效手段。因热膨胀引起的局部温度变化会导致脉冲回波产生移位。具体原因在于温度变化会引起声速变化，而局部声速的变化会最终引起回波移位。

基于热膨胀的交换效应和物理位移之间的关系，Bayat等人(Bayat M,Ballard JR,Ebbini E S.In Vivo Ultrasound Thermography in Presence of TemperatureHeterogeneity and Natural Motions[J].IEEE Transactions on BiomedicalEngineering,2015,62(2):450-457)对温度场进行了重建，并且表明温度重建滤波器具有微分-积分形式。积分器为低通滤波器是由局部组织特性决定的，可以减少温度波动和热效应伪像。最终推导得到的回波位移与温度之间的关系为：

其中，

表示回波位移的变化量对时间求一阶导数；/>

表示温度的变化量对时间求一阶导数。t表示当前纵向深度z处对应的超声传输延时；c表示声速，α和β分别表示由介质决定的热传导系数(Claudio Simon,Philip VanBaren,Emad S.Ebbini.Two-dimensional Temperature Estimation Using Diagnostic Ultrasound[J].IEEETransac-tions on Ultrasonics,Ferroelectrics,and Frequency Control,1998,45(4):1088-1099)。

此外，在估计体内温度时，自由呼吸、血液脉动等自然运动引起的变形也会引入非热位移。利用运动应变的全局有效性和自适应滤波技术，Bayat等人提出了一种可以显著减少运动伪像的方法。具体为使用热消融区域外的应变信息去校正热消融区域内的估计误差。在热消融进行前，假设自然运动引起的误差可以表示为：

其中，(m,n)表示纵向横向交织的空间网格点；k表示离散时间序列的索引；a_i(m,n)为(m,n)对应的自适应滤波器系数；N表示构建自适应滤波器的参考点数；θ_a(m_i，n_i，k)表示热消融区域外参考点对应的温度。使用归一化的最小均方法(Normalized Least MeanSquares,NLMS)对a(m,n)＝[a1(m,n),...,aN(m,n)]进行迭代求解：

其中，t_foci(k)＝[θ_a(m₁，n₁，k)，θ_a(m₂，n₂，k)，...θ_a(m_N，n_N，k)]^T；a^(k)(m,n)表示上一次的迭代结果；a^(k+1)(m,n)表示当前次的迭代结果。选取梯度下降算法训练a(m,n)，当温度值下降到某个阈值之下时迭代终止得到最优解a*(m,n)。因而运动补偿后的温度为：

然而，上述运动补偿方案存在以下问题：首先，不同空间位置处受到的自由呼吸和血液脉动程度不同，在运动补偿参考点也应该在一定局域范围内进行选取，即应变分布的不均匀性。其次，自由呼吸和血液脉动具有周期性，而上述过程的自适应滤波系数并未考虑时变特性。最后，梯度下降法求解滤波系数的步长是固定的，求解效率不高。

发明内容

本发明的目的是针对现有技术的缺陷，提供了一种超声热成像技术的自适应运动补偿方法及系统，减小自由呼吸和血液脉动等非热因素引起的回波位移，从而提高温度估计的准确性。

为了实现以上目的，本发明采用以下技术方案：

一种超声热成像技术的自适应运动补偿方法，包括步骤：

S1.对超声热成像的目标区域进行分块，并选取运动补偿的局域参考点；

S2.计算运动补偿的滤波系数；

S3.根据所述选取的参考点以及计算得到的运动补偿的滤波系数，对温度分布进行运动补偿。

进一步的，所述步骤S1具体为：

S11.对超声热成像的目标区域进行纵向分块；

S12.对进行纵向分块的目标区域中的小块进行边界扩展，得到有效邻域；

S13.在所述得到的有效区域内选取数个点作为当前块的运动补偿参考点；

S14.判断所有小块是否全部处理完成，若否，则继续执行步骤S12-S13；若是，则停止处理。

进一步的，所述步骤S12中得到的有效区域为属于扩展区域但不属于目标区域的数据集合。

进一步的，所述步骤S2中计算运动补偿的滤波系数是基于最速梯度下降法进行计算的。

进一步的，所述步骤S2具体包括：

S21.设定迭代精度eps并令离散时间序列的索引k＝0；

S22.设滤波系数的初始值为a^(t)(m,n,k)，并令迭代次数t＝0；

S23.计算当前误差e(m,n,k)相对于a^(t)(m,n,k)的一阶导数，所述一阶导数记为

判断/>

的绝对值是否小于等于eps，若是，则迭代过程终止，并执行步骤S26；否则，则令/>

S24.根据线性搜索计算步长u^(t)＝argmin[e(a^(t)(m,n,k)+u^(t)d^(t))]；

S25.令a^(t+1)(m,n,k)＝a^(t)(m,n,k)+u^(t)d^(t)，t＝t+1，则继续执行步骤S23-S25；

S26.令k＝k+1，则继续执行步骤S22-S26；

S27.根据各个时刻的滤波系数解a^*(m,n,k)确定整个热消融过程的时变滤波系数；其中：

进一步的，所述步骤S3具体为：

S31.根据回波位移和温度估计之间的对应关系，计算当前时刻数据点的温度；

S32.基于滤波系数解计算运动补偿后的时空温度分布；

S33.根据计算得到的运动补偿后的时空温度分布确定运动补偿后整个热消融过程的温度分布。

进一步的，所述步骤S32中基于滤波系数解计算运动补偿后的时空温度分布，表示为：

其中，θ(m,n,k)表示当前时刻k数据点(m,n)的温度。

相应的，还提供一种超声热成像技术的自适应运动补偿系统，包括：

选取模块，用于对超声热成像的目标区域进行分块，并选取运动补偿的局域参考点；

第一计算模块，用于计算运动补偿的滤波系数；

运动补偿模块，用于根据所述选取的参考点以及计算得到的运动补偿的滤波系数，对温度分布进行运动补偿。

进一步的，所述第一计算模块中计算运动补偿的滤波系数是基于最速梯度下降法进行计算的。

进一步的，所述运动补偿模块包括：

第二计算模块，用于根据回波位移和温度估计之间的对应关系，计算当前时刻数据点的温度；

第三计算模块，用于基于滤波系数解计算运动补偿后的时空温度分布；

确定模块，用于根据计算得到的运动补偿后的时空温度分布确定运动补偿后整个热消融过程的温度分布。

与现有技术相比，本发明具有以下优点：

1.考虑到自由呼吸和血液脉动等因素引起的应变分布的不均匀性，对热消融目标区域进行分块处理，不同区域选取不同的运动补偿参考点。

2.考虑到自由呼吸和血液脉动具有周期性，因而设置自适应运动补偿滤波系数随时间变化，即滤波系数满足时变特性。

3.采用最速梯度下降法求解运动补偿滤波系数，迭代步长自适应变化，提高了算法的运行效率。

附图说明

图1是实施例一提供的一种超声热成像技术的自适应运动补偿方法流程图；

图2是实施例二提供的一种超声热成像技术的自适应运动补偿系统结构图。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

本发明的目的是针对现有技术的缺陷，提供一种超声热成像技术的自适应运动补偿方法及系统。

实施例一

本实施例提供一种超声热成像技术的自适应运动补偿方法，如图1所示，包括步骤：

S1.对超声热成像的目标区域进行分块，并选取运动补偿的局域参考点；

S2.计算运动补偿的滤波系数；

S3.根据所述选取的参考点以及计算得到的运动补偿的滤波系数，对温度分布进行运动补偿。

在步骤S1中，对超声热成像的目标区域进行分块，并选取运动补偿的局域参考点。具体包括：

S11.对超声热成像的目标区域进行纵向分块；

具体为，对超声热成像的目标区域进行纵向分块，纵向分块数AxialNum＝round(M/m)。M表示目标区域的纵向长度；m表示单个小块的纵向长度；round表示对运算结果四舍五入。考虑到应变分布的纵向特征比横向特征显著，因而上述方案仅对目标区域进行纵向分块。

在本实施例中，m值取决于运动补偿的精度和算法的运行效率，设置m为0.25cm。需要说明的是，本实施例对m值不做限定，可根据实际情况设置。

S12.对进行纵向分块的目标区域中的小块进行边界扩展，得到有效邻域；

具体为，对第一个小块进行边界扩展，上下各扩展y个数据，左右各扩展x个数据，取属于扩展区域但不属于目标区域的数据集合为有效邻域。

在本实施例中，x和y的取值在一定范围内，值太大则有效邻域不能准确反映当前块的非热位移，值太小则运动补偿结果的信噪比下降。x和y均设置为5个像素。需要说明的是，本实施例对x和y值不做限定，可根据实际情况设置。

S13.在所述得到的有效区域内选取数个点作为当前块的运动补偿参考点；

具体为，在有效邻域中选取N个点作为当前块的运动补偿参考点。参考点可以随机选取也可以按照指定方式选取。例如：N设置为25。

S14.判断所有小块是否全部处理完成，若否，则继续执行步骤S12-S13；若是，则停止处理。

具体为，对下一个小块重复步骤S12-S13，直到所有小块处理完成。

在步骤S2中，计算运动补偿的滤波系数，其中，计算运动补偿的滤波系数是基于最速梯度下降法进行计算的。具体包括：

S21.设定迭代精度eps并令离散时间序列的索引k＝0；

具体为，设定迭代精度eps并令离散时间序列的索引k＝0，即开始计算第一帧的滤波系数。

在本实施例中，设eps＝1e-3，k＝0。需要说明的是，本实施例对eps和k值不做限定，可根据实际情况设置。

S22.设滤波系数的初始值为a^(t)(m,n,k)，并令迭代次数t＝0；

具体为，给定滤波系数的初始值a^(t)(m,n,k)，并令迭代次数t＝0。滤波系数a中包含离散时间序列索引k表示滤波系数考虑了时变特性。

S23.计算当前误差e(m,n,k)相对于a^(t)(m,n,k)的一阶导数，所述一阶导数记为

判断/>

的绝对值是否小于等于eps(1e-3)，若是，则迭代过程终止，并执行步骤S26；否则，则令/>

S24.根据线性搜索计算步长u^(t)＝argmin[e(a^(t)(m,n,k)+u^(t)d^(t))]；

S25.令a^(t+1)(m,n,k)＝a^(t)(m,n,k)+u^(t)d^(t)，t＝t+1，则继续执行步骤S23-S25；

S26.令k＝k+1，则继续执行步骤S22-S26，即开始计算下一帧的滤波系数；

S27.根据各个时刻的滤波系数解a^*(m,n,k)确定整个热消融过程的时变滤波系数；其中：

在步骤S3中，根据所述选取的参考点以及计算得到的运动补偿的滤波系数，对温度分布进行运动补偿。具体包括：

S31.根据回波位移和温度估计之间的对应关系，计算当前时刻数据点的温度；

具体为，根据回波位移和温度估计之间的对应关系，计算当前时刻k数据点(m,n)的温度θ(m,n,k)。

S32.基于滤波系数解计算运动补偿后的时空温度分布，表示为：

其中，θ(m,n,k)表示当前时刻k数据点(m,n)的温度。

S33.根据计算得到的运动补偿后的时空温度分布θ_c(m,n,k)确定运动补偿后整个热消融过程的温度分布。

本实施例与已公开的超声热成像技术的运动补偿方案相比，充分考虑了人体自由呼吸和血液脉动等因素对温度估计的影响，通过对热消融目标区域的分块化处理和设置运动补偿滤波系数具有时变特性，提高了热消融过程中温度估计的准确性。此外，通过采用最速梯度下降算法进一步优化了算法的运行效率。

与现有技术相比，本实施例具有以下优点：

1.考虑到自由呼吸和血液脉动等因素引起的应变分布的不均匀性，对热消融目标区域进行分块处理，不同区域选取不同的运动补偿参考点。

2.考虑到自由呼吸和血液脉动具有周期性，因而设置自适应运动补偿滤波系数随时间变化，即滤波系数满足时变特性。

3.采用最速梯度下降法求解运动补偿滤波系数，迭代步长自适应变化，提高了算法的运行效率。

实施例二

本实施例提供一种超声热成像技术的自适应运动补偿系统，如图2所示，包括：

选取模块11，用于对超声热成像的目标区域进行分块，并选取运动补偿的局域参考点；

第一计算模块12，用于计算运动补偿的滤波系数；

运动补偿模块13，用于根据所述选取的参考点以及计算得到的运动补偿的滤波系数，对温度分布进行运动补偿。

进一步的，所述第一计算模块12中计算运动补偿的滤波系数是基于最速梯度下降法进行计算的。

进一步的，所述运动补偿模块13包括：

第二计算模块，用于根据回波位移和温度估计之间的对应关系，计算当前时刻数据点的温度；

第三计算模块，用于基于滤波系数解计算运动补偿后的时空温度分布；

确定模块，用于根据计算得到的运动补偿后的时空温度分布确定运动补偿后整个热消融过程的温度分布。

需要说明的是，本实施例提供的一种超声热成像技术的自适应运动补偿系统与实施例一类似，在此不多做赘述。

与现有技术相比，本实施例具有以下优点：

1.考虑到自由呼吸和血液脉动等因素引起的应变分布的不均匀性，对热消融目标区域进行分块处理，不同区域选取不同的运动补偿参考点。

2.考虑到自由呼吸和血液脉动具有周期性，因而设置自适应运动补偿滤波系数随时间变化，即滤波系数满足时变特性。

3.采用最速梯度下降法求解运动补偿滤波系数，迭代步长自适应变化，提高了算法的运行效率。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (1)

1.一种超声热成像技术的自适应运动补偿系统，其特征在于，包括：

选取模块，用于对超声热成像的目标区域进行分块，并选取运动补偿的局域参考点；具体过程包括：S11.对超声热成像的目标区域进行纵向分块；S12.对进行纵向分块的目标区域中的小块进行边界扩展，得到有效邻域；S13.在所述得到的有效区域内选取数个点作为当前块的运动补偿参考点；S14.判断所有小块是否全部处理完成，若否，则继续执行步骤S12-S13；若是，则停止处理；所述步骤S12中得到的有效区域为属于扩展区域但不属于目标区域的数据集合；

第一计算模块，用于计算运动补偿的滤波系数；计算运动补偿的滤波系数是基于最速梯度下降法进行计算，具体包括以下过程：S21.设定迭代精度eps并令离散时间序列的索引k＝0；

S22.设滤波系数的初始值为a^(t)(m,n,k)，并令迭代次数t＝0；

S23.计算当前误差e(m,n,k)相对于a^(t)(m,n,k)的一阶导数，所述一阶导数记为

判断/>

的绝对值是否小于等于eps，若是，则迭代过程终止，并执行步骤S26；否则，则令/>

S24.根据线性搜索计算步长u^(t)＝argmin[e(a^(t)(m,n,k)+u^(t)d^(t))]；

S25.令a^(t+1)(m,n,k)＝a^(t)(m,n,k)+u^(t)d^(t)，t＝t+1，则继续执行步骤S23-S25；

S26.令k＝k+1，则继续执行步骤S22-S26；

S27.根据各个时刻的滤波系数解a^*(m,n,k)确定整个热消融过程的时变滤波系数；其中：

所述第一计算模块中计算运动补偿的滤波系数是基于最速梯度下降法进行计算的；

所述运动补偿模块包括：

第二计算模块，用于根据回波位移和温度估计之间的对应关系，计算当前时刻数据点的温度；

第三计算模块，用于基于滤波系数解计算运动补偿后的时空温度分布，表示为：

其中，θ(m,n,k)表示当前时刻k数据点(m,n)的温度；

确定模块，用于根据计算得到的运动补偿后的时空温度分布确定运动补偿后整个热消融过程的温度分布。

Images