“HIFU Beam：”用于预测分层介质中聚焦换能器产生的轴对称非线性声场的模拟器

A. 单元件聚焦换能器，在水中传播

在第一个示例中，声场由球形碗形状的单元件源生成。仿真代表一个直径为 100 毫米、焦距为 90 毫米、工作频率为 1 MHz 且以水为传播介质的换能器。水参数从默认 GUI 设置中选择：c 0 = 1500 m/s，ρ 0 = 1000 kg/m 3，β = 3.5，δ = 4.33 mm 2 /s。我们之前的论文也考虑了这种典型的实验室超声源[ 60 , 62]。这里，给出了在 KZK 模式下对等效源模型进行的模拟结果，以证明与激波形成聚焦条件相关的临界非线性波效应，并比较了使用求解器的 KZK 和 WAPE 模式进行建模的结果。

当选择KZK模式时，用户在界面中指示源的参数，然后在软件内部重新计算初始平面z = 0中定义的等效源的相应参数[ 59 ]。对于所考虑的单元件换能器，等效源的参数为：直径114.06 mm，焦距98.27 mm，源压力幅度的比例因子p 0 e q / p 0 = 0.92。

图2(a)所示是针对不同源压力幅度p 0使用KZK模式在焦点处获得的归一化压力波形。在非常低的源压力幅度（p 0 = 0.015 MPa）下，观察到线性传播，并且波谱中仅存在基频，如图2（b）所示。随着源处的压力p 0增加，非线性效应变得明显并且焦点处的压力波形开始以不对称的方式扭曲。按照总波强度的 10% 分布在基频谐波上的标准，准线性波形失真的阈值对应于p 0 = 0.15 MPa（图2（a））[ 74 ]。对于该波形，二次谐波的压力振幅达到基频振幅的约1/3（图2（b））。

将p 0增加到超出准线性条件阈值p 0 = 0.15 MPa 后，在一定的源压力水平下会形成激波前沿。首先，冲击出现在波形的正峰值附近。p 0的进一步增加导致激波的底部边缘向零压力水平移动。当归一化为源压力p 0 的冲击振幅达到最大值时，非线性失真水平可以描述为对应于已发展的冲击的形成[ 60]。在这个水平上，震源波形中激波前沿的底部边缘位于零压力水平，因此激波幅度等于波形的峰值正压力[ 60 ]。在所考虑的单元件传感器示例中，产生的冲击在p 0 = 0.45 MPa 时形成（图 2(a)）。

随着p 0进一步增加，超过震源波形中已形成激波形成的水平（p 0 = 0.45 MPa），激波的底部边缘继续向下移动到零压力水平以下。然而，由于激波处强烈的能量吸收，激波前沿振幅的增长速度减慢（图2（a），0.7 MPa和1.45 MPa的波形），激波开始在焦前形成。对于这些饱和冲击，波形变得更接近非对称锯齿形状（图2（b））。请注意，“HIFU 光束”默认界面不仅允许显示焦点处的压力波形，还允许沿轴向坐标z （ r处的轴上）显示压力波形。= 0 mm) 并沿着焦平面 ( z = F ) 处的横向坐标r。

由非线性效应引起的主焦瓣形状的特征变化如图3所示，其中归一化峰值的轴向（图3（a，c））和横向（图3（b，d） ）分布正p +和峰值负p −压力（图 3 (a,b)）以及归一化强度（图 3 (c,d) ）针对增加的源压力幅度p 0给出。对于线性聚焦条件（p 0 = 0.015 MPa），峰值正压和峰值负压的分布相对于零对称（图3（a，b）），焦压增益等于 63.4。非线性和衍射的综合效应导致较高输出水平的分布显着不对称。在准线性聚焦情况下（p 0 = 0.15 MPa），峰值正压= 的聚焦增益为92.6，峰值负压的聚焦增益几乎低两倍，为49.6。在焦点处形成激波的情况下（p 0 = 0.45 MPa），峰值正压的聚焦增益为218.7，比线性情况高3.5倍；对于峰值负压，聚焦增益为 35.8，比线性情况低 1.7 倍。p +和p −的分布彼此之间变得强烈不对称，峰值压力水平相差 6 倍。在更高的源压力下，在饱和聚焦条件下，当p 0 = 0.7 MPa 和 1.45 MPa时，峰值正压的聚焦增益分别为 165.5 和 91.7 ；峰值负压对应的焦点增益为29.7和20.8。

随着波谱的非线性变化和冲击形成时能量沉积的增加，强度的空间分布随着源振幅的增加而变化。开始时，由于光束中产生的高次谐波的聚焦更加有效，强度聚焦效率随着源压力p 0 的增加而增加 [ 58 , 75 ]。形成发达激波后，由于聚焦前发生的激波处的波能被吸收，聚焦效率下降（图3（c，d））。当p 0从0.45 MPa（发展冲击）增加到1.45 MPa（饱和冲击）时，强度的聚焦增益下降了近4倍（图3（c，d）））。

与聚焦增益值类似，峰值正压的聚焦区域的轴向和横向尺寸（-3 dB 水平）非单调变化。最初，随着源压力的增加，焦点区域的轴向尺寸减小，当发生激波形成时达到最小值；然后，随着饱和激波的形成，这些轴向尺寸开始缓慢增长。图4提供了峰值正压和峰值负压的焦点区域的形状演变的直观表示。-3 dB 级峰值正压焦点区域的尺寸从 7.8 mm×1.32 mm 沿z和r开始分别是线性情况下的坐标。在准线性情况下这些尺寸缩小到6.0毫米×0.95毫米（图4（a）），在发展的激波形成情况下达到最小值5.1毫米×0.52毫米（图4（b）），然后在p 0 = 1.45 MPa的饱和状态下生长到11.7mm×0.72 mm （图4（c））。在最后一种情况下，峰值正压绝对水平达到约 150 MPa。请注意，横向r上焦点区域的宽度如此之小（ p 0为 0.52 mm）= 0.45 MPa）与实验中水听器尖端的尺寸相当，因此会受到平均效应的影响，并且对尖端的精确定位非常敏感。这是在非常高的源输出水平下非线性场的水听器测量中峰值正压有时被低估的主要原因[ 6 ]。这种平均效应强调了建模工具在表征具有由于非线性传播效应而产生的精细空间和时间结构的声场时的好处和重要性。

最大峰值正压p +的位置随着源压力p 0 的增加而变化。在线性情况下， p +的最大值位于z = 89.8 mm 处，比几何焦点z = F = 90 mm稍微靠近光源。然后， p +的最大值远离光源移动到z = 90.1 mm，甚至超出几何焦点。此外，在p 0 = 0.45 MPa 时，它在z = 89.5 mm 处返回到更靠近源的位置，并且对于饱和冲击，由于非线性散焦效应 ( zp 0 = 0.7 MPa 和 1.45 MPa 时，z = 90.5 mm 和z = 93.3 mm 。焦点移动现象对于解释点水听器测量非常重要，它表明在哪个源输出处确定了p +的最大值[ 5 , 6 ]。

与峰值正压p +不

同，峰值负压p -的焦点区域的尺寸随着p 0的增加而单调变化：其在z和r两个方向上的尺寸略微增加并且焦点最大值向换能器移动。在准线性情况下，聚焦区域沿着z和r坐标分别从线性光束的7.8mm×1.32mm变化到9.2mm×1.58mm （图4（d））。在激波形成的情况下，相应的尺寸增加到10.7 mm×1.86 mm（图4（e）），在饱和状态下变成12.5 mm×2.1 mm（图4(f) )。

在与焦点处的饱和效应相对应的源压力水平p 0处，使用KZK模拟模型时观察到p +的轴向分布上的特定双尖峰结构。事实上，这种结构在图 3(a)中在p 0 = 0.7 MPa 和p 0 = 1.45 MPa处的p +曲线上清晰可见为焦前尖峰。为了确定形成这种尖峰的物理原因，分析其附近的压力波形是有启发性的（图5）。

这里考虑的单元件换能器是产生直波和边波的活塞源。在光束的焦点波瓣中，就在焦点附近形成尖峰之前，压力波形包含两个前沿（图 5，z = 86.2 和 86.7 mm 处的波形）。波形中的左（下）冲击对应于直达波，而右（上）冲击来自活塞源的边波[ 75 ]。由于激波前沿的能量吸收与激波振幅的立方成正比𝐴3sh，具有两次冲击的波的能量吸收与两次冲击幅度的立方之和成正比(𝐴3sh1+𝐴3sh2)。传播到靠近焦点时，直达波和边波的激波前沿在一定距离处合并（图 5 ， z = 87.2 mm处的波形）。尖峰的顶部正好位于这个距离处。这两种冲击合并后，能量吸收与(𝐴sh1+𝐴sh2)3并且由于不平等而急剧增加(𝐴sh1+𝐴sh2)3>(𝐴3sh1+𝐴3sh2)导致合并冲击的振幅和波形的峰值压力下降（图5 ， z = 87.4 mm处的波形）。由于这种冲击-冲击相互作用发生在焦前，由于焦点附近的衍射放大， p +的急剧下降又恢复到增加。之前已经观察到类似的冲击-冲击相互作用[ 76-78 ]。

在气动声学中，当马赫干形成时，激波前沿的类似相互作用会发生在反射边界附近[ 79 ]。在实验中使用纹影技术来可视化马赫干是很常见的。在这里，提出了梁纹影图像的数值模拟，以显示冲击锋的结构和几何形状。图 6(a)所示为距声束轴不同横向距离r处的压力波形，距声源的距离为z = 86.7 mm。在这个距离处，边缘波和中心波还没有融合。这些压力波形的时间导数的二维分布（图6（b））是纹影图像的类比。该图像上的深色条纹对应于中心波中激波前沿的位置，而边缘波的前沿由较亮的灰色线表示。可以清楚地看到，中心波和边缘波的激波的合并发生在这些锋面之间的中心部分形成的三角形内部（图6（b））。合并后，马赫干垂直于梁轴（此处未显示）。

一个波周期内两个激波形成和合并的物理机制是衍射，特别是边缘波的到来。激波与激波相互作用的现象非常难以研究，因为应该分析激波锋面的精细结构以进行正确的描述和解释。由于 KZK 方程在抛物线近似中考虑了衍射效应，因此这里的一个逻辑问题是，这种由冲击-冲击相互作用引起的空间结构是否可以在抛物线近似中正确描述。

为了测试这一点，在“HIFU 光束”的 WAPE 模式下进行了额外的模拟。对于KZK 和WAPE 模式，对于选择作为界面中提供的推荐网格参数的相同p 0值，数值网格是相同的。所有计算的径向网格步长均为 0.0125 mm，最大谐波数为 1000，线性传播的轴向网格步长为 0.08 mm，准线性和发达激波形成情况为 0.04 mm，p 0 = 0.7 MPa 时为 0.025 mm，p 0 = 0.7 MPa 时为 0.0125毫米，p 0 = 1.45 MPa。模拟示例的运行时间范围从 2 分钟（线性情况）到 1 小时 17 分钟（p 0处的饱和冲击）= 1.45 MPa) 在 KZK 模式下，在 WAPE 模式下对应情况为 6 分钟到 2 小时 10 分钟。模拟是在配备 AMD Ryzen 3800X 8 物理核心处理器的台式电脑上进行的。

总体而言，这两个模型获得的结果在几个百分点内非常一致，表明 KZK 方程与等效源模型的组合提供了高精度声学非线性场的表征。例如，几何焦点处的峰值正压和峰值负压之间的差异不超过3%。在尖峰附近观察到最大的差异，p 0 = 0.7 MPa，其中冲击-冲击峰值形成的影响导致 KZK 模型的峰值正压比 WAPE 模型预测的高约 8%（图 1）。7）。我们的建议是从峰值正压最大值的评估中排除该尖峰，或者如果不偏好模拟 KZK 方程，则使用“HIFU 光束”的 WAPE 模式。

B. 环形阵列，水中传播

典型 HIFU 源的下一个例子是环形阵列，用于开发小鼠肿瘤组织消融的热和机械方法 [ 80 , 81 ]。源是一个 16 元件环形阵列（3 MHz 频率、48 毫米直径、35 毫米曲率半径以及元件之间的 0.15 毫米空间间隙）。在水中声场的计算中，使用了 WAPE 模式，并说明了与电子聚焦控制相关的特定非线性效应。

对线性、准线性和发展的激波形成情况进行了模拟。在第一个系列的模拟中，没有对阵列元素进行额外的定相。当聚焦于阵列的几何焦点时，准线性情况对应于p 0 = 0.14 MPa，并且在p 0 = 0.425 MPa 的焦点处形成激波。模拟对应于饱和度的情况，p 0 = 0.75 MPa。在第二系列模拟中，通过在阵列元素上添加相位将焦点转向距离z = 30 mm。在这种情况下，转向焦点处的波形在p 0处呈准线性= 0.17 MPa，并且在p 0 = 0.525 MPa时，在转向焦点处形成了发达的激波。还模拟了p 0 = 0.75 MPa的饱和情况，以显示在有和没有焦点控制的情况下相同源输出产生的差异。两个系列的模拟均选择相同的网格步长：径向r为 0.0025 毫米，轴向z为 0.025 毫米，最大谐波数为 1000。运行时间为 18 分钟（线性情况）至 2 小时（形成冲击波的情况）在上面提到的台式电脑上。

图8所示的压力波形对应于聚焦在几何焦点（图8（a））并转向5mm时焦点处的三个程度的非线性波形畸变（准线性、形成冲击和饱和）更接近源（图8（b））。源的小尺寸及其相对较高的工作频率（3 MHz）允许在焦点处产生非常高振幅（约 200 MPa）的冲击锋。将焦点转向震源会增加聚焦角，这是控制非线性效应的震源最重要的参数、激波形成条件的特征震源输出以及焦点处形成的震波的特征振幅[ 60 , 82]]。当焦点转向z = 30 mm 时，光束比没有转向时具有更高的聚焦角，在更高的压力水平下形成冲击波（p 0 = 0.525 MPa 与 0.425 MPa），并且具有更大的振幅（图 8）。

焦瓣的大小也随着聚焦角度的变化而变化。弱聚焦光束具有较长的焦瓣，如峰值压力分布中清楚可见：转向情况下的焦瓣在z和r方向上都比无转向情况下窄（图 9））。因此，对于所考虑的 16 元件阵列，将焦点转向靠近源 5 mm 会导致以下相对特征：主焦瓣尺寸较小；形成更大幅度、更高震源压力的发达激波锋；在聚焦饱和状态下，受控焦点处的压力水平更高。在设计环形阵列、表征其非线性场以及在开发冲击波暴露协议时测试非线性转向能力时，应考虑这些效应，并可以使用“HIFU 光束”进行分析。

C. 单元件球形 HIFU 源，多层介质

最后一个示例是单元件球形聚焦传感器，其参数与临床 Sonalleve V1 MRgHIFU 系统（加拿大安大略省米西索加市 Profound Medical Corp.）中使用的 256 元件阵列的参数相当[ 5 ]。具体来说，我们考虑工作频率f 0 = 1.2 MHz、焦距F = 120 mm、直径D = 120 mm的单元件源。考虑在由水 (0 ≤ z ≤ 80 mm)、肌肉 (80 mm ≤ z ≤ 105 mm) 和肾 ( z ≥ 105 mm)组成的分层介质中进行繁殖。组织参数选自默认的“HIFU 光束”组织轮廓 [ 63 ]：c0 = 1585 m/s，ρ 0 = 1060 kg/m 3，β = 4.8，δ = 4.33 mm 2 /s，α 0 = 0.12 Np/cm/MHz，对于肌肉，η = 1.1， c 0 = 1570 m/ s，ρ 0 = 1050 kg/m 3，β = 4.7，δ = 4.33 mm 2 /s，α 0 = 0.1 Np/cm/MHz，η = 1.1（对于肾脏）。在 WAPE 模式下对线性、准线性 ( p 0 = 0.15 MPa) 进行了模拟，并开发了激波形成 ( p0 = 0.375 MPa）情况。网格步长在径向r上选择为 0.0025 mm，在轴向z上选择为 0.025 mm ，最大谐波数为 1000。在上述台式电脑上，运行时间为 16 分钟（线性情况）到 1 小时（开发的冲击形成情况）。

在存在组织层的情况下，峰值正压p +的最大值和峰值负压 | 的绝对值。p- |由于折射效应，比在水中更接近源。对于焦距F = 120 mm 的传感器，水中的最大线性压力幅度位于z = 119.8 mm（此处未显示）。在水-肌肉-肾脏介质中，峰值移动 2.6 mm，更靠近z = 117.2 mm处的源（图 10(a)）。水中的线性焦点增益为80.8，但在吸收性层状介质中降低至约50（图10（a））。

与前面仅在水中传播的示例类似，在组织层中聚焦的非线性光束表现出p +聚焦增益的增加，直到发生激波形成。在准线性情况下，p +的焦点增益达到70（图10）。当在焦点处形成激波时，p +的焦点增益为 165。然而，由于组织中非线性参数的较高值导致激波在较低压力水平下形成，所以产生的激波的振幅仅为 60 MPa 左右。

随着p 0的增加超过了激波形成的水平，激波不仅出现在焦点处，而且出现在邻近区域。此外，聚焦增益开始减小，聚焦区域的尺寸增大（此处未显示）。峰值负压p −的焦点波瓣的行为方式也与之前的情况相同：随着p 0 的增加，它沿轴向z方向增长，并且其最大值向源移动（图 10））。因此，组织层的存在不会定性地改变靠近焦点的非线性效应的表现。然而，实现可比的冲击焦点波形需要更高的源输出水平，以补偿组织中的吸收效应以及与组织的更高非线性相关的改变的冲击形成条件。当重新计算在水中获得的声场参数以在组织中传播时，这些效应被用于非线性降额方法 [ 57 , 83 ]。

除了表征非线性场的声学特性之外，“HIFU Beam”软件还可以评估超声波聚焦到组织中的热效应。用户可用的选项是总声束功率W随传播距离的变化以及热沉积速率的一维和二维分布的可视化。

在线性情况下，组织中波束的总声功率呈指数下降（图11（a））。激波前沿的形成导致额外的热沉积，与激波振幅的立方成正比𝐴3sh。如图11(a)清楚地所示，在激波形成的情况下，总束功率的下降在“肌肉-肾”边界之后靠近焦点处显着偏离指数定律。在焦点区域，由激波形成引起的功率损失急剧增加（图11（b））。尽管线性光束聚焦和焦点中产生激波的非线性状态之间，焦点中光束功率的归一化总损耗仅相差约两倍（图11（b）），但这些损耗分布在不同的区域。体积（图12）：当激波形成时，热源焦点区域的轴向和径向尺寸缩小3倍（图12））。这导致焦点处局部热沉积的显着差异。在线性情况下，计算的功率与发达激波形成时的功率相同（568W，p 0 = 0.375 MPa），最大热沉积速率为 2.6 W/mm 3，而对于非线性传播，则高出约 20 倍（50.6 W/毫米3）。

因此，冲击波暴露方案可以在非常局部的体积内提供非常快速的组织加热。这种快速加热已用于沸腾组织解剖方法以及现有临床 HIFU 系统中冲击波快速体积热组织消融方法的开发 [ 61 , 83 – 86 ]。