基于交错网格时域有限差分法的经颅聚焦超声模拟器实现

王祥达 滕世国 李建霖 吕根品

（1.乳源瑶族自治县东阳光实业发展有限公司 广东省韶关市 512721）（2.乳源东阳光机械有限公司 广东省韶关市 512721 3.韶关东阳光自动化设备有限公司 广东省韶关市 512721）

经颅聚焦超声是一项蓬勃发展的脑部疾病无创治疗手段。其在颅内肿瘤消融、神经调控、血脑屏障打开、颅内靶向给药、颅内血栓溶解等方面具有应用前景。当前，经颅聚焦超声的主要障碍是颅骨与其周围介质较大的声阻抗差异、颅骨较强的非均匀性及声吸收特性所带来的颅骨及其周围组织过热和焦点偏移或散焦。目前主要通过数值模拟来预测颅骨给经颅聚焦超声疗法带来的困难，并给出相应的解决办法。传统的经颅聚焦超声模拟采用的控制方程主要是HIFU 中常用的粘性流体Westervelt 方程。该方程只考虑纵波传播，并未考虑颅骨中存在的横波。近年来，经颅聚焦超声模拟对颅骨横波的考虑逐渐增多，包括利用交错网格时域有限差分法数值求解Biot 方程、利用时域伪谱法数值求解Kelvin-Voigt 方程所建立的考虑横波的经颅聚焦超声模拟器。本文则利用交错网格时域有限差分法数值求解Kelvin-Voigt 方程，实现一套经颅聚焦超声场模拟器。

1 基本方程

1.1 各向同性非均匀粘弹性介质中超声波的传播

对各项同性非均匀粘弹性介质，应力和应变之间的能量守恒关系为：

T、ε、V 是应力、应变、振速。λ、μ 是拉密常数，χ、η 是纵、横波粘性系数，满足：

ρ0是密度，cp、cs是纵、横波波速，αp、αs是纵、横波吸收系数，ω=2πf0、f0是角频率、频率。

为刻画弹性波传播，还必须结合应力-应变之间的动量守恒关系：

方程（1）和（3）构成考虑颅骨横波的经颅聚焦超声场所用的Kelvin-Voigt 控制方程。

1.2 完美匹配层（PML）吸收边界

本文使用不分裂卷积完全匹配层(CPML)作为模拟吸收边界。对空间微分算子操作如下：

参数κx≥1、αp≥，dx是衰减因子，i 是虚数单位。满足迭代式：

2 数值方法

2.1 三维颅脑模型的建立

三维颅骨模型的建立是基于颅骨CT 扫描文件，如图1所示。

图1：基于颅脑CT 扫描的三维颅脑模型建立

颅骨声参数具有较强的非均匀性，其与CT 值H 之间具有如下关系：

φ是颅骨孔隙率，ρskull,min、cp,skull,min、αp,skull,min是颅骨中最小的密度、纵波声速、纵波吸收系数；ρskull,max、cp,skull,max、αp,skull,max是颅骨中最大的密度、纵波声速、纵波吸收系数。对于颅骨内横波声速和吸收系数，基于经验计算，其值可以设置为：cs=4/7cp，αs=90/85αp。

2.2 Kelvin-Voigt方程的交错网格时域有限差分格式

本文使用空间四阶精度、时间二阶精度的交错网格时域有限差分法数值求解Kelvin-Voigt 方程。不同时刻应力和振速分别定义为Tij|n、Vi|n+1/2，下标i、j 代表笛卡尔坐标，上标n 表示时间坐标(t0+nΔt)，Δt 为时间间隔。方程（1）、（3）结合CPML 的时间离散格式为：

Tφij、Vφij是CPML 吸收边界引入的吸收应力、振速的辅助变量。

3 模拟器的验证

3.1 周期性波的传播

利用周期性波的传播上述模拟器。忽略，方程（1）、（3）退化为无黏流体声波方程：

p=-Txx=-Tyy=-Tzz是声压，Txy=Tyz=Tzz=0。

在均匀流体中，方程（8）一组周期性波的解析解为：

i、j、k 是x、y、z 方向单位向量，k 是波数，p0和V0满足：p0+ρ0cpV0=0。

方程（9）显示声压和振速具有周期性空间分布，可将数值域设置为棱长为一个波长的正方体，同时用周期性边界条件取代CPML 吸收边界。f0=1MHz，ρ0=1000kg/m3，cp=1500m/s，V0=1m/s，网格点间距Δx=Δy=Δz。在超过半个周期的计算之后，检验归一化声压误差(pn-pa)/ ρ0cp，pn是声压数值解，pa是方程（9）中的声压解析解。以整个数值域中所有时刻[1,nt]最大2-范数和∞-范数表征该误差，其定义分别为：

设置时间步长满足Δt/Tp=0.0001 并保持不变，Tp=1/f0是周期。归一化声压误差的2-范数和∞-范数随空间步长的变化及其拟合曲线如图2所示。声压的数值计算误差随着空间步长的增大而增大，声压的数值计算结果在空间上确实是四阶精度的。

图2：归一化声压的误差的2-范数和∞-范数随着空间步长的变化

同理，设置空间步长满足Δx/λ=0.01并保持不变，λ= cp/f0是波长。归一化声压误差的2-范数和∞-范数随时间步长的变化及其拟合曲线如图3所示。声压的数值计算误差同样也是随着时间步长的增大而增大，声压的数值计算结果在时间上也确实是二阶精度的。

图3：归一化声压的误差的2-范数和∞-范数随着时间步长的变化

3.2 圆形活塞轴线声压分布

无限大平面障板上半径为a 的圆形活塞声源，圆心位于坐标原点，活塞表面以Vz|z=0,ρ≤a=Vacos(ωt)均匀振动并向+ 轴辐射声场。(ρ, φ,z)为柱坐标系坐标，Va是活塞表面振速幅值。则该圆形活塞声源在轴线上(x=y=0，z ≥0)的归一化（归一化因子为2ρ0cpVa）声压分布为：

圆形活塞声源辐射声场的数值计算域如图4所示。红线包围区域为整个计算域，绿线和红线之间区域为PML 域。数值模拟设置如下：a=2λ，2 个数值域Lx*Ly*Lz=12λ*12λ*27λ & 12λ*12λ*31λ}，空间步长Δx=Δy=Δz=λ/10，PML厚度3λ，f0=1MHz，ρ0=1000kg/m3，cp=1500m/s，Va=1m/s。为使计算结果收敛，时间步长设置满足cmaxΔt/Δx=0.1，cmax是最大声速。

图4：圆形活塞声源辐射超声场的数值计算域示意图

Lz=27λ 和Lz=31λ 的计算域数值模拟得到的圆形活塞声源轴线上的归一化声压幅值分布与方程（11）的理论解析分布结果对比如图5所示。数值模拟结果和解析结果具有很好的一致性，PML 吸收边界对边界内声场的吸收效果非常好，几乎没有数值反射波。

图5：圆形活塞声源的轴线上的归一化声压幅值

4 结语

本文基于大脑的CT 扫描文件重建三维颅脑模型，分别利用空间四阶精度、时间二阶精度的交错网格时域有限差分算法求解Kelvin-Voigt 方程，从而建立了经颅聚焦超声场模拟器。利用周期性声波的传播问题和圆形活塞辐射声场问题对经颅聚焦超声场模拟器进行了核实。结果表明，模拟器确实对应空间四阶精度、时间二阶精度的时域有限差分法；圆形活塞声源轴线上声压分布的模拟值与解析解的一致性很好；CPML 吸收边界很好的消除了声波在计算区域边界的反射。从而证实了经颅聚焦超声场准确性和有效性。