基于FDTD模拟仿真聚焦超声温度场的方法、模拟器、系统

摘要

本发明提供一种基于FDTD模拟仿真聚焦超声温度场的方法、模拟器、系统，所述的方法包括如下步骤：S1：通过获取聚焦超声温度场的数值模拟方程；其中，所述的聚焦超声温度场的数值模拟方程采用Pennes生物热传导方程表示；S2：采用空间二阶精度、时间一阶精度的时域有限差分法求解S1中的数值模拟方程，实现模拟聚焦超声引发的温升效应。本发明利用空间二阶精度、时间一阶精度的时域有限差分算法求解Pennes生物热传导方程，从而建立了一种基于FDTD模拟仿真聚焦超声温度场的方法。聚焦超声温度场的数值模拟精度确实对应空间二阶精度、时间一阶精度的时域有限差分法，确保了模拟仿真的有效性和准确性。

权利要求

1.一种基于FDTD模拟仿真聚焦超声温度场的方法，其特征在于：所述的方法包括如下步骤：S1：通过获取聚焦超声温度场的数值模拟方程；其中，所述的聚焦超声温度场的数值模拟方程采用Pennes生物热传导方程表示；

S2：采用空间二阶精度、时间一阶精度的时域有限差分法求解S1中的数值模拟方程，实现模拟聚焦超声引发的温升效应；

所述的Pennes生物热传导方程：

其中，T表示生物组织温度；κt表示导热系数；Ct表示生物组织比热；Cb表示血管比热；Wb表示生物组织内毛细血管的血液灌注率；Ta表示初始温度；Qv表示单位体积单位时间内组织吸收的热量；Qm表示生物代谢生热率；

取 表示温度的升高量，忽略Qm，一般起始温度是常数温度，将Pennes生物热传导方程表示为：对于聚焦超声温度场模拟，使用单频稳态聚焦超声场，Qv表示如下：

2 2

Qv＝αIav＝αρ0c0＝1/2×αρ0c0V0

其中，α表示声吸收系数；Iav表示平均声强；c0表示波速；V表示质点振动速度；<*>表示对里面的变量求周期平均值；V0表示质点振动速度幅值。

2.根据权利要求1所述的基于FDTD模拟仿真聚焦超声温度场的方法，其特征在于：所述的空间二阶精度的时域有限差分法具体采用空间网格点对模拟仿真区域进行划分，相邻的空间网格点的间距为0.1λ，其中，λ表示超声波的波长。

3.根据权利要求2所述的基于FDTD模拟仿真聚焦超声温度场的方法，其特征在于：所述的空间二阶精度的时域有限差分法，具体求解如下：对于一元函数f(x)，其中x是自变量，其二阶导数的二阶精度中心差分格式为：(2)

其中，f |i中的上标(2)表示一元函数f(x)的二阶导数、f|i中的下标i表示第i个网格点对应的自变量的值(x0+iΔx)，x0是自变量初始值，Δx是相邻网格点之间自变量的差值；

不同时刻、不同空间网格点位置的温升分别定义为 上标n表示时间坐标(t0+nΔt)，其中t0表示起始时间，Δt表示数值模拟的时间步长；下标i、j、k表示笛卡尔空间坐标(x0+iΔx,y0+jΔy,z0+kΔz)，其中(x0,y0,z0)表示起始位置，Δx、Δy、Δz分别表示x、y、z方向的空间步长。

4.根据权利要求3所述的基于FDTD模拟仿真聚焦超声温度场的方法，其特征在于：所述的时间一阶精度的时域有限差分法，具体求解如下：时间离散格式为：

方程(4)中二阶空间导数 的差分格式采用方程(3)，从而有：

5.一种基于权利要求1～4任一项所述的方法的模拟器，其特征在于：所述的模拟器包括获取模块、空间二阶精度的时域有限差分法求解模块、时间一阶精度的时域有限差分法求解模块；

所述的获取模块，用于获取Pennes生物热传导方程；

所述的空间二阶精度的时域有限差分法求解模块、时间一阶精度的时域有限差分法求解模块分别对Pennes生物热传导方程进行求解Pennes生物热传导方程，实现模拟聚焦超声引发的温升效应。

6.一种计算机系统，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，其特征在于：所述的处理器执行所述的计算机程序时，实现如权利要求1～4任一项所述的方法的步骤。

说明书

基于FDTD模拟仿真聚焦超声温度场的方法、模拟器、系统

技术领域

[0001] 本发明涉及聚焦超声温度场仿真技术领域，更具体地，涉及一种基于FDTD模拟仿真聚焦超声温度场的方法、模拟器、系统。

背景技术

[0002] 高强度聚焦超声波(HIFU)肿瘤消融技术是近年比较开始流行的新型肿瘤治疗方法。其最大特点是利用体外超声波换能器发射的高强度超声波无创地穿过皮肤并在体内的肿瘤形成超声波焦点。焦点处的温度在瞬间升高，从而令周围的肿瘤组织发生凝固坏死，最后形成伤疤或被新陈代谢吸收。如中国专利公开号：CN101108268A，公开日：2008‑01‑23，公开了一种生物医学工程领域的相控阵聚焦超声的多模式热场形成方法，包括如下步骤：(1)采用多元阵列球面开孔不等间距环形阵列治疗头；(2)采用每种声场模式与一组驱动信号的相位与幅值相对应；(3)利用试探法确定目标靶组织中可用的高热模式声强；(4)温热模式热场声强的调控，以保证加热区域冷点温度上升速率保持为预定值，最终达到形成任意预期均匀温热模式热场的目的。

[0003] 但目前，现有技术均存在焦点的余热或者焦点聚焦位置不准确等等原因，周围的健康组织也会很容易被杀死，造成不必要的损伤。因此，通过仿真对聚焦超声导致的温度场进行研究，对于聚焦超声设备的设计、聚焦超声疗法的安全性和有效性的保障、聚焦超声疗法的合适策略的制定等方面具有重要的意义。

发明内容

[0004] 本发明为克服现有技术模拟聚焦超声导致的温度场不够准确的问题，提供了一种基于FDTD模拟仿真聚焦超声温度场的方法、模拟器、系统，其能有效性和准确性模拟仿真聚焦超声温度场。

[0005] 为解决上述技术问题，本发明的技术方案如下：一种基于FDTD模拟仿真聚焦超声温度场的方法，所述的方法包括如下步骤：

[0006] S1：通过获取聚焦超声温度场的数值模拟方程；其中，所述的聚焦超声温度场的数值模拟方程采用Pennes生物热传导方程表示；

[0007] S2：采用空间二阶精度、时间一阶精度的时域有限差分法求解S1中的数值模拟方程，实现模拟聚焦超声引发的温升效应。

[0008] 优选地，所述的Pennes生物热传导方程：

[0009]

[0010] 其中，T表示生物组织温度；κt表示导热系数；Ct表示生物组织比热；Cb表示血管比热；Wb表示生物组织内毛细血管的血液灌注率；Ta表示初始温度；Qv表示单位体积单位时间内组织吸收的热量；Qm表示生物代谢生热率。

[0011] 进一步地，取 表示温度的升高量，忽略Qm，一般起始温度是常数温度，将Pennes生物热传导方程表示为：

[0012]

[0013] 对于聚焦超声温度场模拟，Qv即聚焦超声给生物组织带来的外部热源输入；对于单频稳态聚焦超声场，超声热源Qv表示如下：

[0014] Qv＝αIav＝αρ0c0＝1/2×αρ0c0V02

[0015] 其中，α表示声吸收系数；Iav表示平均声强；c0表示波速；V表示质点振动速度；<*>表示对里面的变量求周期平均值；V0表示质点振动速度幅值。

[0016] 再进一步地，所述的空间二阶精度的时域有限差分法具体采用空间网格点对模拟仿真区域进行划分，相邻的空间网格点的间距为0.1λ，其中，λ表示超声波的波长。

[0017] 再进一步地，所述的空间二阶精度的时域有限差分法，具体求解如下：

[0018] 对于一元函数f(x)，其中x是自变量，其二阶导数的二阶精度中心差分格式为：

[0019]

[0020] 其中，f(2)|i中的上标(2)表示一元函数f(x)的二阶导数、f|i中的下标i表示第i个网格点对应的自变量的值(x0+iΔx)，x0是自变量初始值，Δx是相邻网格点之间自变量的差值；

[0021] 不同时刻、不同空间网格点位置的温升分别定义为 上标n表示时间坐标(t0+nΔt)，其中t0表示起始时间，Δt表示数值模拟的时间步长；下标i、j、k表示笛卡尔空间坐标(x0+iΔx,y0+jΔy,z0+kΔz)，其中(x0,y0,z0)表示起始位置，Δx、Δy、Δz分别表示x、y、z方向的空间步长。

[0022] 再进一步地，所述的时间一阶精度的时域有限差分法，具体求解如下：

[0023] 方程(2)的显式的时间离散格式为：

[0024]

[0025] 方程(4)中二阶空间导数 的差分格式采用方程(3)，从而有：

[0026]

[0027] 一种模拟器，所述的模拟器包括获取模块、空间二阶精度的时域有限差分法求解模块、时间一阶精度的时域有限差分法求解模块；

[0028] 所述的获取模块，用于获取Pennes生物热传导方程；

[0029] 所述的空间二阶精度的时域有限差分法求解模块、时间一阶精度的时域有限差分法求解模块分别对Pennes生物热传导方程进行求解Pennes生物热传导方程，实现模拟聚焦超声引发的温升效应。

[0030] 优选地，所述的空间二阶精度的时域有限差分法求解模块实现的求解如下：

[0031] 对于一元函数f(x)，其中x是自变量，其二阶导数的二阶精度中心差分格式为：

[0032]

[0033] 其中，f(2)|i中的上标(2)表示一元函数f(x)的二阶导数、f|i中的下标i表示第i个网格点对应的自变量的值(x0+iΔx)，x0是自变量初始值，Δx是相邻网格点之间自变量的差值；

[0034] 不同时刻、不同空间网格点位置的温升分别定义为 上标n表示时间坐标(t0+nΔt)，其中t0表示起始时间，Δt表示数值模拟的时间步长；下标i、j、k表示笛卡尔空间坐标(x0+iΔx,y0+jΔy,z0+kΔz)，其中(x0,y0,z0)表示起始位置，Δx、Δy、Δz分别表示x、y、z方向的空间步长。

[0035] 进一步地，时间一阶精度的时域有限差分法求解模块实现的求解如下：

[0036] 方程(2)的显式的时间离散格式为：

[0037]

[0038] 方程(4)中二阶空间导数 的差分格式采用方程(3)，从而有：

[0039]

[0040] 一种计算机系统，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，所述的处理器执行所述的计算机程序时，实现如权利要求1～6任一项所述的方法的步骤。

[0041] 与现有技术相比，本发明技术方案的有益效果是：

[0042] 本发明利用空间二阶精度、时间一阶精度的时域有限差分算法求解Pennes生物热传导方程，从而建立了一种基于FDTD模拟仿真聚焦超声温度场的方法。将聚焦超声温度场的数值模拟精度确实对应空间二阶精度、时间一阶精度的时域有限差分法，提高了模拟仿真的有效性和准确性。

附图说明

[0043] 图1是本实施例所述的方法的步骤流程图。

[0044] 图2是本实施例归一化温升的误差的2‑范数和∞‑范数随着空间步长的变化。

[0045] 图3是本实施例归一化温升的误差的2‑范数和∞‑范数随着时间步长的变化。

[0046] 图4是本实施例区域中心温升随着时间的变化。

具体实施方式

[0047] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，仅用于示例性说明，不能理解为对本专利的限制。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0048] 下面结合附图和实施例对本发明的技术方案做进一步的说明。

[0049] 实施例1

[0050] 如图1所示，一种基于FDTD模拟仿真聚焦超声温度场的方法，所述的方法包括如下步骤：

[0051] S1：通过获取聚焦超声温度场的数值模拟方程；其中，所述的聚焦超声温度场的数值模拟方程采用Pennes生物热传导方程表示；

[0052] S2：采用空间二阶精度、时间一阶精度的时域有限差分法求解S1中的数值模拟方程，实现模拟聚焦超声引发的温升效应。

[0053] 在一个具体的实施例中，所述的Pennes生物热传导方程：

[0054]

[0055] 其中，T表示生物组织温度；κt表示导热系数；Ct表示生物组织比热；Cb表示血管比热；Wb表示生物组织内毛细血管的血液灌注率；Ta表示初始温度；Qv表示单位体积单位时间内组织吸收的热量；Qm表示生物代谢生热率，计算时Qm一般可忽略。

[0056] 取 表示温度的升高量，忽略Qm，一般起始温度是常数温度，将Pennes生物热传导方程表示为：

[0057]

[0058] 对于聚焦超声温度场模拟，Qv即聚焦超声给生物组织带来的外部热源输入；对于单频稳态聚焦超声场，超声热源Qv表示如下：

[0059] Qv＝αIav＝αρ0c0＝1/2×αρ0c0V02

[0060] 其中，α表示声吸收系数；Iav表示平均声强；c0表示波速；V表示质点振动速度；<*>表示对里面的变量求周期平均值；V0表示质点振动速度幅值。

[0061] 在一个具体的实施例中，所述的空间二阶精度的时域有限差分法具体采用空间网格点对模拟仿真区域进行划分，相邻的空间网格点的间距为0.1λ，其中，λ表示超声波的波长。

[0062] 再进一步地，所述的空间二阶精度的时域有限差分法，具体求解如下：

[0063] 对于一元函数f(x)，其二阶导数的二阶精度中心差分格式为：

[0064]

[0065] 其中，f(2)|i中的上标(2)表示一元函数f(x)的二阶导数、f|i中的下标i表示第i个网格点对应的自变量的值(x0+iΔx)，x0是自变量初始值，Δx是相邻网格点之间自变量的差值；

[0066] 不同时刻、不同空间网格点位置的温升分别定义为 上标n表示时间坐标(t0+nΔt)，其中t0表示起始时间，Δt表示数值模拟的时间步长；下标i、j、k表示笛卡尔空间坐标(x0+iΔx,y0+jΔy,z0+kΔz)，其中(x0,y0,z0)表示起始位置，Δx、Δy、Δz分别表示x、y、z方向的空间步长。

[0067] 在一个具体的实施例中，所述的时间一阶精度的时域有限差分法，具体求解如下：

[0068] 方程(2)的显式的时间离散格式为：

[0069]

[0070] 方程(4)中二阶空间导数 的差分格式采用方程(3)，从而有：

[0071]

[0072] 为了证实本实施例所述的方法的可靠行，进行如下验证：

[0073] 1.有限区域的温升初值问题

[0074] 不考虑生物热传导方程(2)右边的热源项Qv和与血管有关的项 则方程(2)具有如下形式：

[0075]

[0076] 对于边长分别为a、b、c的长方体有限区域，假设该区域表面温升保持不变：

[0077]

[0078] 初始温升分布为：

[0079]

[0080] 此后温升随时空的变化为：[9]

[0081]

[0082] 数值模拟设置如下。长方体边长相等，a＝b＝c＝10cm，密度ρ0＝1000kg/m3，导热系数κt＝0.6W/(m·K)，比热Ct＝4180J/(kg·K)，

[0083] 2.利用2‑范数和∞‑范数表征数值计算误差，具体如下：

[0084]

[0085] 其中， 表示温升数值解， 表示方程(9)的温升解析解。

[0086] 设置时间步长满足Δt＝0.0001s并保持不变，分析空间步长对温升的数值计算结果的影响。温升误差的2‑范数和∞‑范数随着空间步长的变化如图1所示，温升的数值计算误差随着空间步长的增大而增大。根据范数表征的误差随着空间步长分布的拟合曲线也在图1中画出，拟合曲线的结果表明，温升的数值计算结果在空间上确实是二阶精度的。

[0087] 同理，设置空间步长为Δx＝0.001m并保持不变，分析时间步长对温升的数值计算结果的影响。归一化温升的误差的2‑范数和∞‑范数随着时间步长的变化如图2所示，温升的数值计算误差同样也是随着时间步长的增大而增大，而且，拟合曲线的结果也证实了温升的数值计算结果在时间上确实是一阶精度的。

[0088] 设置空间步长和时间步长分别为Δx＝0.005m和Δt＝0.1s，总的计算时间为2500s。图3给出了数值模拟计算得到的区域中心温升随时间变化的曲线和利用方程(9)从理论上计算得到的区域中心温升随时间变化的曲线，从图中可以很清晰的看到两条曲线的吻合度非常高，数值模拟结果和解析结果具有很好的一致性。

[0089] 实施例2

[0090] 一种模拟器，所述的模拟器包括获取模块、空间二阶精度的时域有限差分法求解模块、时间一阶精度的时域有限差分法求解模块；

[0091] 所述的获取模块，用于获取Pennes生物热传导方程；

[0092] 所述的Pennes生物热传导方程：

[0093]

[0094] 其中，T表示生物组织温度；κt表示导热系数；Ct表示生物组织比热；Cb表示血管比热；Wb表示生物组织内毛细血管的血液灌注率；Ta表示初始温度；Qv表示单位体积单位时间内组织吸收的热量；Qm表示生物代谢生热率，计算时Qm一般可忽略。

[0095] 取 表示温度的升高量，忽略Qm，一般起始温度是常数温度，将Pennes生物热传导方程表示为：

[0096]

[0097] 对于聚焦超声温度场模拟，Qv即聚焦超声给生物组织带来的外部热源输入；对于单频稳态聚焦超声场，超声热源Qv表示如下：

[0098] Qv＝αIav＝αρ0c0＝1/2×αρ0c0V02

[0099] 其中，α表示声吸收系数；Iav表示平均声强；c0表示波速；V表示质点振动速度；<*>表示对里面的变量求周期平均值；V0表示质点振动速度幅值。

[0100] 所述的空间二阶精度的时域有限差分法求解模块、时间一阶精度的时域有限差分法求解模块分别对Pennes生物热传导方程进行求解耦合，实现模拟聚焦超声引发的温升效应。

[0101] 其中，所述的空间二阶精度的时域有限差分法求解模块实现的求解如下：

[0102] 对于一元函数f(x)，其中x是自变量，其二阶导数的二阶精度中心差分格式为：

[0103]

[0104] 其中，f(2)|i中的上标(2)表示一元函数f(x)的二阶导数、f|i中的下标i表示第i个网格点对应的自变量的值(x0+iΔx)，x0是自变量初始值，Δx是相邻网格点之间自变量的差值；

[0105] 不同时刻、不同空间网格点位置的温升分别定义为 上标n表示时间坐标(t0+nΔt)，其中t0表示起始时间，Δt表示数值模拟的时间步长；下标i、j、k表示笛卡尔空间坐标(x0+iΔx,y0+jΔy,z0+kΔz)，其中(x0,y0,z0)表示起始位置，Δx、Δy、Δz分别表示x、y、z方向的空间步长。

[0106] 所述的时间一阶精度的时域有限差分法求解模块实现的求解如下：

[0107] 方程(2)的显式的时间离散格式为：

[0108]

[0109] 方程(4)中二阶空间导数 的差分格式采用方程(3)，从而有：

[0110]

[0111] 实施例3

[0112] 一种计算机系统，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，所述的处理器执行所述的计算机程序时，实现的方法的步骤如下：

[0113] S1：通过获取聚焦超声温度场的数值模拟方程；其中，所述的聚焦超声温度场的数值模拟方程采用Pennes生物热传导方程表示；

[0114] S2：采用空间二阶精度、时间一阶精度的时域有限差分法求解S1中的数值模拟方程，实现模拟聚焦超声引发的温升效应。

[0115] 显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。