基于蒙特卡洛法的射频消融快速仿真方法及系统与流程

1.本发明涉及医疗领域，具体地，涉及一种基于蒙特卡洛法的射频消融快速仿真方法及系统。


背景技术：

2.在射频消融手术中，常常使用仿真方法对消融区域及温度场分布进行估计，需要对电热耦合场进行迭代计算，目前最常用的有限元算法具有耗时高、内存占用高等缺点。
3.专利文献为cn111027225a的发明专利公开了一种基于数据拟合的消融仿真方法，该方法包括：计算同一消融功率、消融时间下的消融边界的短径、长径数据的实测平均值，并对靶区进行热场仿真，以短径数据为基准，将温度阈值所在等高线对应的长径为长径仿真值；在低消融功率下，建立温度阈值与消融功率关系模型；在高消融功率下，建立温度阈值与消融时间关系模型；建立修正系数与消融功率关系模型；输入消融功率、消融时间，根据消融功率高低选择温度阈值与消融功率关系模型或温度阈值与消融时间关系模型，得到温度阈值，根据温度阈值得到对应的消融仿真图像，利用修正系数与消融功率关系模型，对消融仿真图像的长径进行修正。但是上述方案基于经验数据统计，缺乏对个体化差异的适应；温度阈值难以确定，无法与真实消融过程建立直接联系。
4.唐茹冰发表的《射频加热治疗肿瘤体模实验的数值模拟》一文中提出的肿瘤体模的有限元稳态热仿真方法，虽然能够计算稳态导热，但是对于非稳态以及非射频加热的场景缺乏适应性。mohammad h.n.naraghi发表的《a boundary-dispatch monte carlo(exodus)method for analysis of conductive heat transfer problems》一文中提出的边界派发式的蒙特卡洛热仿真方法，对于边界不规则的情况适应程度不佳，无法处理内部热源点。


技术实现要素：

5.针对现有技术中的缺陷，本发明的目的是提供一种基于蒙特卡洛法的射频消融快速仿真方法及系统。
6.根据本发明提供的一种基于蒙特卡洛法的射频消融快速仿真方法，包括如下步骤：
7.基础仿真模型构建步骤：输入ct/mr形式的影像数据，根据影像数据生成体素网格，构建基础仿真模型；
8.参数估计步骤：输入基础仿真模型，基于影像数据的分割结果，细化模型边界，对边界体素进行细分，对于基础仿真模型计算域中的每一处，根据影像数据的数值，估算差异化参数；
9.计算域约化步骤：基于多维度条件对计算域进行聚类约化；
10.电/热物理场耦合计算步骤：使用蒙特卡洛法对电物理场耦合计算进行简化，使用有限差分法对热物理场进行计算，迭代进行快速的热/电物理场耦合计算；
11.细胞凋亡计算步骤：计算细胞凋亡结果，判断是否满足停止条件，若判断结果为是，则输出消融区域；若判断结果为否，则跳转至计算域约化步骤。
12.优选地，每个体素单元之间使用pennes生物导热方程的差分形式确定热流传导方向和关系。
13.优选地，所述计算域约化步骤包括：
14.基于影像数据的分割结果确定组织类型，根据组织边界调整计算域，将目标组织边界外的部位排除出计算域；
15.根据消融探针位置，将超出距离阈值以及不在手术范围内的区域排除出计算域；
16.依据仿真中间结果，对相变、凋亡等动态过程造成的参数变化进行响应；
17.对已经变性、烧焦的区域调整模型参数，根据对称性条件，对模型计算域中电/热属性以及距离构成的计算条件相似的位置进行聚类约化。
18.优选地，所述差异化参数包括差异化的导热率、电阻率、血流灌注率。
19.优选地，使用蒙特卡洛法对电物理场耦合计算包括：
20.步骤s1：确定射频发射位置为起始位置，根据预设参数生成数量为n的模拟粒子，有比例为tb/th的粒子为高能粒子，其余为低能粒子，其中tb代表加热的目标温度，th代表组织细胞的最高温度；
21.步骤s2：计算粒子的随机游走，每运行一个时间步长，模拟粒子根据模型的电传导可行方向游走m步，直到到达边界停止；
22.步骤s3：所有模拟粒子随机游走结束后，统计每个体素网格内的经过的高能粒子数和低能粒子数，网格电场强度估计为ei nh/n
×
δe，其中ei为原始场强数值，δe为电场强度单位值；
23.步骤s4：根据焦耳热计算公式qj＝j
·
e计算射频电场在一个时间步长内的产生的热能，根据pennes方程
[0024][0025]
在当前网格模型中的差分形式，使用有限差分法求解传热方程，将结果反馈到温度场中。
[0026]
优选地，采用cuda对粒子游走过程进行了gpu加速。
[0027]
优选地，使用cuda对传热计算进行并行化硬件加速。
[0028]
优选地，所述细胞凋亡计算步骤中，根据温度场分布，依照阿累尼乌斯方程计算细胞累积收到的热损伤
[0029][0030]
其中a为消融损伤系数，δe为激活能量，r为通用气体常数，t为温度，当ω超过阈值时，将细胞设置为凋亡状态，计算由消融热效应造成的组织细胞凋亡结果，从而推算消融区域大小。
[0031]
优选地，使用cuda对损伤计算进行并行化硬件加速。
[0032]
根据本发明提供的一种基于蒙特卡洛法的射频消融快速仿真系统，包括如下模块：
[0033]
基础仿真模型构建模块：输入ct/mr形式的影像数据，根据影像数据生成体素网格，构建基础仿真模型；
[0034]
参数估计模块：输入基础仿真模型，基于影像数据的分割结果，细化模型边界，对边界体素进行细分，对于基础仿真模型计算域中的每一处，根据影像数据的数值，估算差异化参数；
[0035]
计算域约化模块：基于多维度条件对计算域进行聚类约化；
[0036]
电/热物理场耦合计算模块：使用蒙特卡洛法对电物理场耦合计算进行简化，使用有限差分法对热物理场进行计算，迭代进行快速的热/电物理场耦合计算；
[0037]
细胞凋亡计算模块：计算细胞凋亡结果，判断是否满足停止条件，若判断结果为是，则输出消融区域；若判断结果为否，则跳转至计算域约化模块。
[0038]
与现有技术相比，本发明具有如下的有益效果：
[0039]
1、本发明通过参数估计，实现了高低精度混合的体素网格建模，缓解了传统建模方式导致的计算量冗余的问题。
[0040]
2、本发明通过计算域约化，实现了对具有对称性、等价性位置的自动归约，自动化地将不需要计算的区域排除出计算域，减少了计算量。
[0041]
3、本发明通过蒙特卡洛法，对电场计算进行了合理简化，结合gpu硬件加速，实现了快速的电场计算。
[0042]
4、本发明通过cuda进行gpu加速，实现了电/热耦合场的加速运算，提供了实时的手术区域温度场计算。
附图说明
[0043]
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
[0044]
图1为基于蒙特卡洛法的射频消融快速仿真方法的算法流程图。
[0045]
图2为蒙特卡洛法示意图。
[0046]
图3为消融区域仿真结果图。
[0047]
图4为消融温度场仿真结果图。
具体实施方式
[0048]
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。
[0049]
如图1至4所示，本发明提供了一种基于蒙特卡洛法的射频消融快速仿真方法及系统，实现了对于射频消融仿真从模型构建、数值计算、硬件加速三个层面的全面加速，分别为：模型构建：根据影像信息对计算域进行聚类预处理，构建约化仿真模型；数值计算：使用蒙特卡洛法对电物理场耦合计算进行合理简化，使用速度更快的有限差分法对热物理场进行计算，迭代进行快速的热/电物理场耦合计算；硬件加速：基于cuda使用gpu硬件加速仿真计算过程中的模型构建和数值计算过程。具体包括：输入ct/mr形式的影像数据，根据影像
数据生成体素网格，构建基础仿真模型；参数估计步骤：输入基础仿真模型，基于影像数据的分割结果，细化模型边界，对边界体素进行细分，对于基础仿真模型计算域中的每一处，根据影像数据的数值，估算差异化的导热率、电阻率、血流灌注率等参数；计算域约化步骤：基于多维度条件对计算域进行聚类约化；电/热物理场耦合计算步骤：使用蒙特卡洛法对电物理场耦合计算进行简化，使用有限差分法对热物理场进行计算，迭代进行快速的热/电物理场耦合计算；细胞凋亡计算步骤：计算细胞凋亡结果，判断是否满足停止条件，若判断结果为是，则输出消融区域；若判断结果为否，则跳转至计算域约化步骤。
[0050]
进一步地，每个体素单元之间使用pennes生物导热方程的差分形式确定热流传导方向和关系。所述计算域约化步骤包括：基于影像数据的分割结果确定组织类型，根据组织边界调整计算域，将目标组织边界外的部位排除出计算域；根据消融探针位置，将超出距离阈值以及不在手术范围内的区域排除出计算域；依据仿真中间结果，对相变、凋亡等动态过程造成的参数变化进行响应；对已经变性、烧焦的区域调整模型参数，根据对称性条件，对模型计算域中电/热属性以及距离构成的计算条件相似的位置进行聚类约化。
[0051]
使用蒙特卡洛法对电物理场耦合计算包括：
[0052]
步骤s1：确定射频发射位置为起始位置，根据预设参数生成数量为n的模拟粒子，有比例为tb/th的粒子为高能粒子，其余为低能粒子，其中tb代表加热的目标温度，th代表组织细胞的最高温度；
[0053]
步骤s2：计算粒子的随机游走，每运行一个时间步长，模拟粒子根据模型的电传导可行方向游走m步，直到到达边界停止；采用cuda对粒子游走过程进行了gpu加速。
[0054]
步骤s3：所有模拟粒子随机游走结束后，统计每个体素网格内的经过的高能粒子数和低能粒子数，网格电场强度估计为ei nh/n
×
δe，其中ei为原始场强数值，δe为电场强度单位值；
[0055]
步骤s4：根据焦耳热计算公式qj＝j
·
e计算射频电场在一个时间步长内的产生的热能，根据pennes方程
[0056][0057]
在当前网格模型中的差分形式，使用有限差分法求解传热方程，将结果反馈到温度场中。使用cuda对传热计算进行并行化硬件加速。
[0058]
所述细胞凋亡计算步骤中，根据温度场分布，依照阿累尼乌斯方程计算细胞累积收到的热损伤
[0059][0060]
其中a为消融损伤系数，δe为激活能量，r为通用气体常数，t为温度，当ω超过阈值时，将细胞设置为凋亡状态，计算由消融热效应造成的组织细胞凋亡结果，从而推算消融区域大小。使用cuda对损伤计算进行并行化硬件加速。
[0061]
更为详细的，本发明提供的基于蒙特卡洛法的射频消融快速仿真方法，包括：
[0062]
s1：基础仿真模型构建。输入ct/mr形式的影像数据，根据影像数据生成体素网格，每个体素单元之间使用pennes生物导热方程的差分形式确定热流传导方向和关系，构建低精度的基础仿真模型。在构建了低精度的基础仿真模型后，对区域内各个位置的属性进行
初始化，初始化参数根据预设的组织基本参数确定。
[0063]
s2：参数估计。输入基础仿真模型，基于影像数据的分割结果，细化模型边界，对边界体素进行进一步细分以匹配脏器形状。对于仿真模型计算域中的每一处，根据影像数据的数值，估算差异化的导热率、电阻率、血流灌注率等参数。
[0064]
s3：计算域约化。基于影像数据的分割结果确定组织类型，根据组织边界调整计算域，将目标组织边界外的部位排除出计算域；根据消融探针位置，将超出距离阈值以及不在手术范围内的区域排除出计算域；依据仿真中间结果，对相变、凋亡等动态过程造成的参数变化进行响应，对已经变性、烧焦的区域调整模型参数，根据对称性条件，对模型计算域中电/热属性以及距离构成的计算条件相似的位置进行聚类约化。
[0065]
经过s2和s3，我们可以获得精细化后的仿真模型，接下来是数值计算部分。
[0066]
s4：电/热物理场耦合计算。使用有限差分法和蒙特卡洛法，进行电场和温度场的耦合计算。其中对于电场分布的计算使用了蒙特卡洛法的粒子派发算法进行加速计算。每次迭代会计算得到当前条件下射频产热和人体传热造成的温度场分布变化。
[0067]
s4中蒙特卡洛法计算射频产热的具体步骤包括：
[0068]
s4-mc1:确定射频发射位置(射频探针上的槽口等发射位置)为起始位置，根据预设参数生成数量为n的模拟粒子，有比例为tb/th的粒子为高能粒子，其余为低能粒子，其中tb代表加热的目标温度，th代表组织细胞的最高温度。
[0069]
s4-mc2:计算粒子的随机游走，每运行一个时间步长，模拟粒子根据模型的电传导可行方向游走m步，直到到达边界停止。由于此部分粒子间运动无具有相互独立的特点，使用了cuda对粒子游走过程进行了gpu加速。
[0070]
s4-mc3:所有模拟粒子随机游走结束后，统计每个体素网格内的经过的高能粒子数和低能粒子数，网格电场强度可以估计为ei nh/n
×
δe，其中ei为原始场强数值，δe为电场强度单位值。
[0071]
s4-mc4:根据焦耳热计算公式qj＝j
·
e计算射频电场在一个时间步长内的产生的热能，然后根据pennes方程
[0072][0073]
在当前网格模型中的差分形式，使用有限差分法求解传热方程，将结果反馈到温度场中。由于此部分各个位置之间传热方程求解相互独立，使用cuda对传热计算进行了并行化硬件加速。
[0074]
s5：细胞凋亡计算。根据温度场分布，依照阿累尼乌斯方程计算细胞累积收到的热损伤
[0075][0076]
其中a为消融损伤系数，δe为激活能量，r为通用气体常数，t为温度，当ω超过阈值(默认为1)时，将细胞设置为凋亡状态。计算由消融热效应造成的组织细胞凋亡结果，从而推算消融区域大小。由于此部分各位置损伤累积相互独立，使用cuda对损伤计算进行了并行化硬件加速。
[0077]
如果不满足停止条件(可以使用运算时间阈值、核心温度阈值、消融灶大小阈值
等)，则跳转步骤3，进行下一轮迭代；如果满足，则输出最终消融区域，包括消融区域的短径、长径以及具体形状。
[0078]
使用本算法进行数量为125000的体素网格中的不同加速方式下的仿真性能比较，可以得到如下结果：
[0079][0080]
其中cpu使用的是无加速的计算，cpu多线程使用c 的thread库实现，gpu加速使用cuda进行实现。
[0081]
根据本发明提供的一种基于蒙特卡洛法的射频消融快速仿真系统，包括基础仿真模型构建模块：输入ct/mr形式的影像数据，根据影像数据生成体素网格，构建基础仿真模型；参数估计模块：输入基础仿真模型，基于影像数据的分割结果，细化模型边界，对边界体素进行细分，对于基础仿真模型计算域中的每一处，根据影像数据的数值，估算差异化参数；计算域约化模块：基于多维度条件对计算域进行聚类约化；电/热物理场耦合计算模块：使用蒙特卡洛法对电物理场耦合计算进行简化，使用有限差分法对热物理场进行计算，迭代进行快速的热/电物理场耦合计算；细胞凋亡计算模块：计算细胞凋亡结果，判断是否满足停止条件，若判断结果为是，则输出消融区域；若判断结果为否，则跳转至计算域约化模块。
[0082]
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本技术的实施例和实施例中的特征可以任意相互组合。