k-wave伪谱法在生物医学超声仿真中的研究进展

2021-10-14陈思瑶杜宇翎李发琪李成海

科技创新与应用 2021年28期

陈思瑶，杜宇翎，周野，李发琪，李成海*

（1.重庆医科大学生物医学工程学院超声医学工程国家重点实验室，重庆 400016；2.重庆市生物医学工程重点实验室，重庆 400016）

超声在生物组织中的非线性传播模拟对于许多应用具有重要作用，包括超声波换能器的设计，聚焦超声治疗计划的制定，诊断超声成像等[1-2]。然而，模拟非线性波的传播是一项计算困难的任务，因为介质离散必须足够精细，以捕捉在超声波传播时产生的谐波。对于强激波，可能会产生几百次谐波，因此需要非常密集的计算网格[3]。目前，多数介质中的非线性超声模拟主要基于Khokhlov-Zabolotskaya-Kuznetsov（KZK）方程[4]或Burgers方程[5]。虽然这些模型在许多实际情况下是准确的[6]，但它们仅限于模拟定向声束，并且通常只考虑均匀介质中的单向传播[7]。有限差分或有限元方法常被用来求解控制方程，但由于每个最小波长需要大量的网格点来避免数值色散[8]，增加了计算的复杂程度。为了减少这种计算负担，基于傅里叶伪谱方法的非线性超声仿真模型被提出[9-10]。与有限差分方法相比，这种方法可以减少所需网格点的总数，最多可减小两个数量级[11-12]。基于k-space伪谱方法的全波非线性声学模拟作为开源k-wave声学工具箱的一部分[11，13]。该方法可用于描述声波在非均匀介质中的声传播，考虑了幂律吸收，对波的方向性没有限制。本文将从k-wave声学工具箱的发展以及k-wave在超声仿真应用中的研究进行综述，为超声在生物组织中非线性传播的模拟研究提供参考。

1 k-wave声学工具箱的发展

k-wave是基于MATLAB和C++编写的开源声学工具箱，由Bradley Treeby、Ben Cox和Jiri Jaros开发。该工具箱使用Westervelt方程的广义模型来模拟超声波在软组织中的传播，考虑了声传播的非线性、组织的不均匀性和幂律吸收。可以实现复杂多层生物组织介质中的时域声学和超声模拟。k-wave工具箱最初在2009年由伦敦大学光声成像小组开发，主要解决无损介质[13]中光声波场的模拟和重建问题。在后续的版本中扩展了包括时变压力和速度源、声吸收、非线性、弹性材料和超声波换能器模型的功能。

求解非线性Westervelt方程常见的方法为有限差分法、有限元法和谱法[14]。有限差分法和有限元法被称为局部法，因为所关注的波动传播方程仅根据附近点的条件在每一点求解。相比而言，谱方法，如k-space方法[15]和伪谱方法[16-17]，被认为是全局的，因为在每一点上都是利用整个波场的信息来求解波的传播方程。k-wave工具箱的控制方程采用k-space伪谱方法求解，其中傅里叶配置谱方法用于计算空间梯度，k-space校正有限差分格式用于在时间上向前积分，具有快速、简便的特点。由于其全局性质，谱方法可以比局部方法更精确，例如，应用于周期问题的伪谱方法已被证明与无限阶有限差分方法等价[18]。非线性k-space仿真的优势已经由临床超声换能器在均匀和非均匀介质中的波束方向图的三维模拟进行验证。通过将3D中的完整空间域离散化，并以时间步进的方式计算域内各处的压力和粒子速度场的解。与以前基于KZK方程的超声模拟相比，此方法对声波的方向性或空间变化没有任何限制。与有限差分法和有限元法相比，这种方法在相同精度下允许更粗的网格间距和更大的时间步长[11-12]。

k-space计算声场的物理模型是基于运动方程、连续性方程以及物态方程构建的Westervelt方程，运动方程、连续性方程以及物态方程的描述如下：

u是声粒子速度，p是声压，ρ是声学介质密度，ρ0是环境（或平衡）介质密度，C0是等熵的声波的速度，d为声粒子位移。这里的B/A是表征有限振幅效应对声速的相对贡献的非线性参数。与线性情况相比，质量守恒方程包括一个附加项，它解释了一个对流非线性，其中粒子速度对波速有贡献。

压力-密度关系中的L算子是一个线性积分-微分算子，它解释了声波的吸收和扩散，遵循频率幂律：

这里τ、η分别是吸收和色散比例系数：

其中α0是幂律吸收因子，单位为Np（rad/s）-ym-1，y是幂律指数。该方程考虑了二阶声学非线性、幂律吸声和材料特性（声速、密度、非线性和吸收系数）的非均匀分布。与传统的有限差分方法相比，这提高了梯度计算的精度，从而降低了对精细计算网格的要求。

k-wave声学工具箱包括了以下几点优势：提供了较为详细的用户手册，包括对控制方程和数值方法的一般介绍，还提供了软件架构的基本概述和许多易于遵循的教程示例，以说明工具箱的功能；与基于时域有限差分（FDTD）格式的模型相比，k-wave数值模型的主要优点是精确模拟所需的空间和时间网格点较少，这意味着这些模型运行得更快，使用的内存更少；能够模拟压力和速度源，包括光声源，以及诊断和治疗超声换能器；可以使用定向元素指定任意检测表面的功能，并可以选择记录声压，粒子速度和声强；优化的C++程序代码，最大限度地提高了大型模拟的计算性能；选择使用前向模型作为一个灵活的时间反转图像重建算法的光声层析与任意测量表面；能快速进行光声图像重建算法，用于记录在线性（2D）或平面（3D）测量面上的数据。

2 k-wave在超声仿真应用中的研究

k-wave作为一种开源的声学工具箱，是实现生物组织的声传播仿真的重要工具，目前国内外已报道基于kwave开展了相关的研究，应用领域主要包括声波在非均匀生物组织介质中传播的数值模拟、光声图像重建以及无损探伤等。

声波在非均匀生物组织介质中传播的数值模拟对开展生物医学超声研究具有重要意义。例如，在超声治疗应用中，数值模拟可用于研究磁共振引导的聚焦超声手术中的相位像差[19-20]，并改善治疗结果。对于诊断超声，数值建模已被用作图像重建以及理解超声成像中图像退化来源的重要工具[21-23]。Suomi等[24]人使用的开源k-wave声学模拟工具箱，基于三个不同患者分割的三维CT数据集预测高强度聚焦超声（High Intensity Focussed Ultrasound，HIFU）治疗剂量，模拟了三位患者肾脏前方组织层的衰减、反射和折射的综合效应如何影响超声场的强度和形状。模拟结果表明，超声场的强度平均下降了11.1 dB，并且发现强度损失可以在衰减和折射之间大致平均分配，为临床HIFU治疗肾癌提供了一定的参考。Kittiphot等人[25]基于k-wave工具包模拟了高强度聚焦超声热疗在肿瘤组织中疗效，模拟结果可见，聚焦区域周围声压的差异随着聚焦深度的增加而增大。进一步使用数据可视化装置，获得相应的温度分布2D图像，以显示乳腺癌的消融疗效，观察到除了焦点区域的损伤，周围的健康组织未受影响。此外，姜翔飞等人[26]利用此工具包，对基于时间反演的经颅二维超声进行仿真，研究了不同换能器阵列数和不同的超声发射频率对经颅超声聚焦效果的影响。张值豪[27]也利用k-wave工具箱对超声断层成像透射过程进行仿真，在环形探头上实现了对平行束以及扇形束重建效果的仿真。

为了精确求解光声反演问题，需要建立光致声波在生物软组织中产生和传播的数值模型。有大量的文献描述了超声在生物软组织中的传播，部分理论可直接适用于光声学，同时开发了相应的方法来数值求解偏微分方程。其中基于k-space伪谱法的数值仿真具有较大优势[28]。有限元（FE）和有限差分（FD）建模技术可以用于建模光致声波的传播，但这两种方法精确的仿真需要网格单元的尺寸约为十分之一波长。对于较高频率条件下的数值仿真，这一要求加大了对计算硬件和算法的要求。COX等人[29]利用kwave进行建模研究了二维传感器元件记录表面平均压力值的理论方向性，然后针对时间反转图像重建，研究了传感器方向对光声层析成像的影响。研究发现即使在使用完整连续的测量曲面时，传感器指向性的加入也会带来图像重建的假象，研究明晰了传感器元件阵列的使用对重建光声图像的结构保真度和定量精度的影响。

除此之外，Prieur等人[30]利用k-wave工具箱模拟横波弹性成像，将均匀和各向同性介质中的模拟结果与解析解和有限元建模仿真的结果进行比较。从k-wave获得的剪切位移与解析解非常相似，在短传播时间和大传播时间下，聚焦区域周围的均方根误差分别低于9%和21%，说明k-wave是模拟声辐射力和剪切波传播的一种有效工具。另外，Gu等人[31]模拟声波在强非均匀介质中传播时，使用k-wave的结果用作比较和验证的基准，系统地评估了改进混合域方法。吴其洲等人[32]在圆柱体缺陷构件的超声探伤重构研究中，利用k-wave工具箱建立仿真平台，在圆柱体周围设置360个传感器，同时接受散射信号，实现缺陷信号的采集，为缺陷重构提供数据支撑。孙明健等人[33]在基于光声信号的高铁钢轨表面缺陷检测方法中，使用有限元及k-wave方法对钢轨表面的光声图像进行了重建，结果表明，采用这种方法能较好的检测到表面微裂纹，此结果在钢轨探伤领域有较大的可行性及发展潜力。张思思等人[34]基于k-wave建立了表面缺陷的超声检测模型，结果为表面缺陷检测、数值模拟、超声回波信号处理以及特征提取等技术提供了理论依据。综上可见，基于k-wave进行的声传播模拟在较多的领域都具有重要的应用价值。