CN105740573B - 一种用于放射线剂量计算的双步蒙特卡洛模拟方法 - Google Patents

一种用于放射线剂量计算的双步蒙特卡洛模拟方法 Download PDF

Info

Publication number
CN105740573B
CN105740573B CN201610115951.2A CN201610115951A CN105740573B CN 105740573 B CN105740573 B CN 105740573B CN 201610115951 A CN201610115951 A CN 201610115951A CN 105740573 B CN105740573 B CN 105740573B
Authority
CN
China
Prior art keywords
photon
numerical value
random number
pseudo random
seed
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201610115951.2A
Other languages
English (en)
Other versions
CN105740573A (zh
Inventor
蔡夫鸿
何赛灵
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Taizhou Anqiling Intelligent Technology Co ltd
Original Assignee
Suzhou Wanghao Information Technology Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Suzhou Wanghao Information Technology Co Ltd filed Critical Suzhou Wanghao Information Technology Co Ltd
Priority to CN201610115951.2A priority Critical patent/CN105740573B/zh
Publication of CN105740573A publication Critical patent/CN105740573A/zh
Application granted granted Critical
Publication of CN105740573B publication Critical patent/CN105740573B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
    • G16H50/50—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for simulation or modelling of medical disorders

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Medical Informatics (AREA)
  • Public Health (AREA)
  • Biomedical Technology (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Pathology (AREA)
  • Epidemiology (AREA)
  • General Health & Medical Sciences (AREA)
  • Primary Health Care (AREA)
  • Measurement Of Radiation (AREA)

Abstract

本发明公开了一种用于放射线剂量计算的双步蒙特卡洛模拟方法，包括如下步骤：步骤1、进行蒙特卡洛仿真，获取有效种子数值数组；步骤2、所述步骤一中在获得有效种子数值数组后，利用有效种子数值，进行二次蒙特卡洛仿真，高速的计算人体内部的放射线剂量分布。本发明提出双步蒙特卡洛算法，可根据有效种子数值，进行作用于肿瘤区域的光子的仿真，避免了进行没有通过肿瘤区域的光子的仿真，从而提高了计算效率。

Description

一种用于放射线剂量计算的双步蒙特卡洛模拟方法
技术领域
本发明涉及人体内放射剂量强度分布的数值计算，主要是用于模拟X射线在人体内部的传播以及剂量分布。
背景技术
放射治疗是一种重要的肿瘤治疗手段，70％的肿瘤患者在病程的不同时期，需要接受放射治疗。然而，在X射线照射肿瘤的时候，也对正常的器官和组织造成了影响。为了提高治疗效率，并降低对正常组织的损伤。在放射治疗之前，需要计算X射线在人体内部的分布，用于辅助医护人员确定肿瘤区域的X射线剂量是否达到了治疗剂量，正常组织的放射线剂量是否在安全阈值以下。
使用能量输运理论，可以准确的计算出X射线在人体内部的剂量分布。然而，人体是极度不均匀的介质，需要将人体离散为大量的体素，才能精确的构建人体的数字模型。这使得能量输运理论的数值求解过程中，需要大内存计算机，并且计算速度也很慢。蒙特卡洛是一种灵活的数值模拟方法，可以模拟X射线在复杂的人体内的传播过程。并且，根据概率统计的大数定律，随着模拟X射线的数目的增多，蒙特卡洛仿真的计算结果可以趋近能量输运理论的计算结果。
蒙特卡洛仿真的引入，解决了能量输运理论中产生的大内存使用的问题，但计算速度很慢。当前，人们引入显卡并行加速的方法提高了计算效率。同时，人们也在不断的提出新的改良算法，提高蒙特卡洛仿真的计算效率。在放射治疗的过程中，肿瘤的结构尺寸会随着治疗的进行而发生改变。人们需要一种高效的X射线的蒙特卡洛算法，可以在放射治疗的过程中，快速的计算体内剂量分布。
双步蒙特卡洛方法在博士论文《蒙特卡洛仿真的加速与电场矢量化及其在生物光学成像中的应用》中已有提及，但其应用点为可见光波段，本专利应用于X射线领域。针对于人体组织，可见光波段将受到极强的瑞利散射，所以并不适于人体组织，而X射线受到的瑞利散射较弱。另外，X射线在人体内部传播时，还将产生康普顿散射，这是可见光在人体内部传播时，不会发生的。另一方面，本发明的后期应用主要为肿瘤放射治疗中的X射线剂量分布，因此，将主要研究肿瘤区域的X射线剂量。这与可见光波段的双步蒙特卡洛方法是不同的。本发明会用到以下参考文献：
[1]Salvat F,Fernández-Varea J M,Sempau J.PENELOPE-2006:A code systemfor Monte Carlo simulation of electron and photon transport[C]//WorkshopProceedings.2006,4:7.
[2]Hubbell J H,Veigele W J,Briggs E A,et al.Atomic form factors,incoherent scattering functions,and photon scattering cross sections[J].Journal of physical and chemical reference data,1975,4(3):471-538.
发明内容
1、本发明的目的：
本发明提出一种双步蒙特卡洛仿真方法，用于在放射治疗的过程中，对于同一个模拟对象的体内的X射线剂量分布，进行高速的数值计算。本方法的核心思想是，因为X射线是一种穿透性能很强的电磁波，其在人体内的传播轨迹，因为人体肿瘤结构的改变而发生变化的概率极小，利用这一特点，我们提出了双步蒙特卡洛仿真，用于仿真X射线在人体内部的传播以及剂量分布。
2、本发明所采用的技术方案：
因为X射线也属于电磁波范畴，在蒙特卡洛仿真中，X射线被命名为光子，这是X射线的蒙特卡洛仿真的惯用命名方法。用于能量小于500kev的放射线剂量计算的双步蒙特卡洛模拟方法，包括如下步骤：
步骤1、进行蒙特卡洛仿真，获取有效种子数值数组；
步骤2、所述步骤一中在获得有效种子数值数组后，利用有效种子数值，进行二次蒙特卡洛仿真，高速的计算人体内部的放射线剂量分布。
更进一步的具体实施例中，所述的步骤1中包括如下步骤：
步骤1.1随机数的产生：在每一个光子的蒙特卡洛仿真开始之前，由当前时间的数值，组合为一个整数，定义为初始种子数值,将该数值输入到传统的伪随机数产生函数中，每产生一个伪随机数，种子数值将会由传统的伪随机数产生函数更新，并再次代入伪随机数产生函数，不断循环上述过程，将会得到一个伪随机数数组；
步骤1.2光子发射：通过入射光子初始坐标与光子发射方向以及放射线产生器与人体之间的源皮距，得到光子到达人体表面时的位置坐标：
Sxi+SSD×Dxi,Syi+SSD×Dyi,Szi+SSD×Dzi，
式中Sxi,Syi,Szi为入射光子初始坐标，Dxi,Dyi,Dzi为光子发射方向，SSD为源皮距，光子的初始能量定义为E；
步骤1.3光子迁移：对于能量小于500kev的X射线，光子在人体内部传播，将会产生瑞利散射、康普顿散射、光电吸收，并以此判断光子单次自由迁移的状态，根据光子所在的体素对应的人体组织类型，通过查表的方式获得当前组织的三种平均自由程，分别为瑞利散射自由程MFP_r、康普顿散射自由程MFP_c以及光电效应自由程MFP_p，上述三个自由程将决定光子单次自由迁移的距离，以及完成单次自由迁移后，发生的散射或吸收的概率，从而判断下一次自由迁移的初始状态；
步骤1.4有效种子数值的判定：在光子迁移的过程中，如果光子经过了用户预设的肿瘤区域，将会把步骤1.1中的初始种子数值定义为有效种子数值；反之，则不做此定义；
步骤1.5边界处理：判断光子是否处于皮肤层，并即将进入空气层，如果不是，将回到步骤1.3，再次进行光子迁移；如果是，则光子到达人体模型的边界，将从人体内部出射，该光子的蒙特卡洛模拟将结束，如果本次蒙特卡洛模拟的初始种子数值被定义为有效种子数值，就将该数值存储，供第二步蒙特卡洛仿真中使用；反之，该初始种子数值将不做任何处理；
1.6判断蒙特卡洛仿真是否结束：当已模拟的光子数目超出预设值时，蒙特卡洛仿真结束，停止模拟，反之，重新回到步骤1.1，开始新的光子的模拟，同时，已模拟的光子数目加一。
更进一步的具体实施例中，所述的步骤2中包括如下步骤：
步骤2.1有效随机数种子数值提取：从步骤1所存储的有效种子数值文件中，依次读取有效种子数值，将有效随机数种子数值依次输入到传统的伪随机数产生函数中，将可以产生一个伪随机数，每产生一个伪随机数，种子数值将会由传统的伪随机数产生函数更新，并再次代入伪随机数产生函数，用于产生下一个伪随机数，接下来的每一个步骤，都将使用伪随机数进行随机抽样，从而判断光子的传播方式，当有效种子数值读取结束时，进入步骤2.5；
步骤2.2光子发射：通过入射光子初始坐标与光子发射方向以及放射线产生器与人体之间的源皮距，得到此时光子的位置坐标：
Sxi+SSD×Dxi,Syi+SSD×Dyi,Szi+SSD×Dzi，
式中Sxi,Syi,Szi为入射光子坐标，Dxi,Dyi,Dzi为光子发射方向，SSD为源皮距，光子的初始能量定义为E
步骤2.3光子迁移：对于能量小于500kev的X射线，光子在人体内部传播，将会产生瑞利散射、康普顿散射、光电吸收，并以此判断光子单次自由迁移的状态，在算法中，根据光子所在的体素对应的人体组织类型，通过查表的方式获得当前组织的三种平均自由程，分别为瑞利散射自由程MFP_r、康普顿散射自由程MFP_c以及光电效应自由程MFP_p，上述三个自由程将决定光子单次自由迁移的距离，以及完成单次自由迁移后，发生的散射或吸收的概率，从而判断下一次自由迁移的初始状态；当光子在体素[i,j,k]处发射光电吸收时，将使用式子Dose[i,j,k]＝Dose[i,j,k]+E将该体素吸收的光子能量记录在数组Dose[i,j,k]中；
步骤2.4边界处理：判断光子是否处于皮肤层，并即将进入空气层，如果不是，将回到步骤2.3，再次进行光子迁移；如果是，则光子到达人体模型的边界，将从人体内部出射，该光子的蒙特卡洛模拟将结束；然后，重新回到步骤2.1，开始新的光子的模拟；
步骤2.5数值计算结束，得到的人体内的剂量分布。
更进一步的具体实施例中，所述的步骤1.1随机数的产生具体产生方法如下：由当前时间组合为一个种子数值Seed，该数值的前四位为当前年份，中间四位为当前日期，最后四位为当前时间，精确到分钟；将该数值输入到传统的伪随机数产生函数中，将可以产生一个平均分布在(0,1)中的伪随机数ε，每产生一个伪随机数，种子数值将会由传统的伪随机数产生函数更新，并再次代入伪随机数产生函数，用于产生下一个平均分布在(0,1)中的伪随机数ε；上述过程不断循环，将会得到一个伪随机数数组，其统计特性为平均分布在(0,1)中；接下来的每一个步骤，都将使用伪随机数进行随机抽样，从而判断光子的传播方式；对于伪随机数列而言，当初始的种子数值确定了，其数列所产生的伪随机数数值也随之确定。
更进一步的具体实施例中，所述的步骤2.1有效随机数种子数值提取具体方法如下：
根据当前模拟的光子的序号，导入有效种子数值Seed_array[i]，其中i是一个自然数，其取值上限为Seed_array的数目，代表当前模拟的光子的序号；将该数值输入到传统的伪随机数产生函数中，将可以产生一个平均分布在(0,1)中的伪随机数ε，每产生一个伪随机数，种子数值将会由传统的伪随机数产生函数更新，并再次代入伪随机数产生函数，用于产生下一个平均分布在(0,1)中的伪随机数ε；上述过程不断循环，将会得到一个伪随机数数组，其统计特性为平均分布在(0,1)中；接下来的每一个步骤，都将使用伪随机数进行随机抽样，从而判断光子的传播方式；对于伪随机数列而言，当初始的种子数值确定了，其数列所产生的伪随机数数值也随之确定，因此，由有效种子数值驱动的光子，都将经过肿瘤区域。
更进一步的具体实施例中，所述的步骤1.1或步骤2.1中，将人体三维数值模型定义为三维数组Density，该三维数组是利用计算机断层扫描仪器获得的人体内部的影像信息；Density[i,j,k]表示在三维网格空间中，体素[i,j,k]区域对应CT密度值；Density[i,j,k]中的数值分别表示人体组织；另外，Density[i,j,k]中有部分区域是仪器与皮肤之间的空气区域；X射线机穿过空气，到达皮肤的距离，定义为源皮距SSD，SSD为大于0的自然数，表示光源到人体表面的距离，定义三维数组Dose[i,j,k]，用于存储人体内的剂量分布，定义蒙特卡洛仿中，需要模拟的光子数量为N，定义已模拟的光子数为变量Ns，Ns的初始值为0，再每一次进行步骤1.1时，Ns＝Ns+1。
更进一步的具体实施例中，所述的步骤1.3或步骤2.3中光子迁移步骤，入射光子的能量以及光子到达区域的人体组织类型，会影响平均自由程；在光子光子传播过程中，其发生瑞利散射、康普顿散射、光电吸收的概率，由MFP_r、MFP_c以及MFP_p的数值决定；瑞利散射不改变光子的能量E，只改变光子的运行方向Dxi,Dyi,Dzi，康普顿散射既改变光子的能量E，也改变光子的运行方向Dxi,Dyi,Dzi，光电吸收发生后，光子的能量将会被吸收，该光子的模拟结束。
更进一步的具体实施例中，所述的光子迁移步骤中的光子迁移采用如下方法：光子迁移自由程记为MFP_l，满足公式1/MFP_l＝1/MFP_r+1/MFP_c+1/MFP_p；自由程的步长的抽样表达式如下：l＝-Ln(ε)/MFP_l；光子行进了当前的自由程后，其空间位置为：
Sxi＝Sxi+l×Dxi,
Syi＝Syi+l×Dyi,
Szi＝Szi+l×Dzi
再次根据当前的位置Sxi,Syi,Szi，获取当前的瑞利散射自由程、康普顿散射自由程以及光电效应自由程；在完成了一次自由程的迁移后，将会依概率发生瑞利散射、康普顿散射或者光电吸收。
更进一步的具体实施例中，所述的发生瑞利散射、康普顿散射或者光电吸收过程中：
定义1/MFP_l＝ul，1/MFP_r＝ur，1/MFP_c＝uc，1/MFP_p＝up再次生成伪随机数ε，当ε<ur/ul时，发生瑞利散射，光子的迁移方向将会发生改变，其中，迁移方向的改变由球坐标下的偏转角θ与方位角决定，偏转角θ的数值由瑞利散射公式进行随机抽样获取；当ur/ul<ε<(ur+uc)/ul，发生康普顿散射，光子能量E将会被衰减，光子迁移方向将会发生改变，其中，迁移方向的改变由球坐标下的偏转角θ与方位角决定，偏转θ的数值由康普顿散射公式进行随机抽样获取；当ε>(ur+uc)/ul发生光电吸收效应，光子能量被位置Sxi,Syi,Szi对应的体素[i,j,k]所吸收，将使用式子Dose[i,j,k]＝Dose[i,j,k]+E将该体素吸收的光子能量记录在数组Dose[i,j,k]中；
所述的产生方式如下：产生一个随机数ε；，那么ψ＝2πε；
所述的光子迁移方向的改变，在发生偏转角θ与方位角的改变后，光子迁移方向为：
在完成上面的计算后，更新光子迁移方向Dx＝Dx',Dy＝Dy',Dz＝Dz'。3、本发明所产生的有益效果：
本发明提出双步蒙特卡洛算法，在第一步中，记录仿真过程中作用于肿瘤区域的光子的有效种子数值。在未来需要再使用蒙特卡洛算法对同一个患者的肿瘤区域做剂量分布计算时，可根据有效种子数值，进行作用于肿瘤区域的光子的仿真，避免了进行没有通过肿瘤区域的光子的仿真，从而提高了计算效率。计算效率的提高与有效种子数值的数目有关。当第一步蒙特卡洛仿真中，总的光子仿真数目是N(N>1000000)，而有效种子数值的数目为M(0<M<N)，计算效率提高系数为N/M。
相比于可见光波段的双步蒙特卡洛方法，本发明处理X射线在人体内部的传播过程，系统的考虑了瑞利散射与康普顿散射的物理过程。而在可见光波段的蒙特卡洛中，只考虑瑞利散射。此外，本发明中，将通过肿瘤区域的X射线记为有效模拟，可在放射治疗过程中，提高肿瘤放射治疗的软件验证的效率。这也是可见光波段的双步蒙特卡洛仿真所不具备的。
附图说明
图1.第一步蒙特卡洛计算的流程图。
图2.第二步蒙特卡洛计算的流程图。
具体实施方式
以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还可以有各种变化和改进。
本发明的生物探测中的高速蒙特卡洛模拟方法，包括如下步骤：
步骤1、进行蒙特卡洛仿真，获取有效种子数值数组；因为X射线也属于电磁波范畴，在蒙特卡洛仿真中，X射线被命名为光子，这是X射线的蒙特卡洛仿真的惯用命名方法。
步骤1.1随机数的产生：由当前时间组合为一个种子数值Seed，该数值的前四位为当前年份，中间四位为当前日期，最后四位为当前时间，精确到分钟；将该数值输入到传统的伪随机数产生函数中，将可以产生一个平均分布在(0,1)中的伪随机数ε，每产生一个伪随机数，种子数值将会由传统的伪随机数产生函数更新，并再次代入伪随机数产生函数，用于产生下一个平均分布在(0,1)中的伪随机数ε；上述过程不断循环，将会得到一个伪随机数数组，其统计特性为平均分布在(0,1)中；接下来的每一个步骤，都将使用伪随机数进行随机抽样，从而判断光子的传播方式；对于伪随机数列而言，当初始的种子数值确定了，其数列所产生的伪随机数数值也随之确定。我们可以将人体三维数值模型定义为三维数组Density，该三维数组是利用计算机断层扫描仪器(Computed Tomography,CT)获得的人体内部的影像信息。Density[i,j,k]表示在三维网格空间中，体素[i,j,k]区域对应CT密度值。Density[i,j,k]中的数值可分别表示血液、肌肉、皮肤等人体组织。另外，Density[i,j,k]中有部分区域是仪器与皮肤之间的空气区域。X射线机穿过空气，到达皮肤的距离，定义为源皮距(Source Skin Distance，SSD,SSD为大于0的自然数，表示光源到人体表面的距离)。在算法中，还将定义三维数组Dose[i,j,k]，用于记录人体内部的剂量分布，定义蒙特卡洛仿中，需要模拟的光子数量为N，定义已模拟的光子数为变量Ns，Ns的初始值为0，再每一次进行步骤1.1时，Ns＝Ns+1；
步骤1.2光子发射：通过用户输入的入射光子坐标(Sxi,Syi,Szi)与光子发射方向(Dxi,Dyi,Dzi)以及源皮距SSD，光子沿着发射方向行进SSD，移动到人体的表面，此时光子的位置坐标为，Sxi+SSD×Dxi,Syi+SSD×Dyi,Szi+SSD×Dzi，光子的初始能量定义为E；
步骤1.3光子迁移：对于能量小于500kev的X射线，光子在人体内部传播，将会产生瑞利散射、康普顿散射、光电吸收，并以此判断光子单次自由迁移的状态，在算法中，根据光子所在的体素对应的人体组织类型，通过查表的方式获得当前组织的三种平均自由程(mean free path，MFP)，分别为瑞利散射自由程MFP_r、康普顿散射自由程MFP_c以及光电效应自由程MFP_p。入射光子的能量以及光子到达区域的人体组织类型，都会影响平均自由程；在光子传播过程中，完成每一次自由迁移后，其发生瑞利散射、康普顿散射、光电吸收的概率，由MFP_r、MFP_c以及MFP_p的数值决定；瑞利散射不改变光子的能量E，只改变光子的运行方向Dxi,Dyi,Dzi，康普顿散射既改变光子的能量E，也改变光子的运行方向Dxi,Dyi,Dzi，能量的改变由传统的康普顿散射理论计算；当光子在体素[i,j,k]处发射光电吸收时，将使用式子Dose[i,j,k]＝Dose[i,j,k]+E将该体素吸收的光子能量记录在数组Dose[i,j,k]中；
步骤1.4有效种子数值的判定：在光子迁移的过程中，如果光子经过了用户预设的肿瘤区域，将会把步骤1.1中的初始种子数值定义为有效种子数值。反之，则不做此定义；
步骤1.5边界处理：判断光子是否处于皮肤层，并即将进入空气层，如果不是，将回到步骤1.3，再次进行光子迁移；如果是，则光子到达人体模型的边界，将从人体内部出射，该光子的蒙特卡洛模拟将结束，如果本次蒙特卡洛模拟的初始种子数值被定义为有效种子数值，就将该数值存储在硬盘中，供第二步蒙特卡洛仿真中使用；反之，该初始种子数值将不做任何处理。
1.6判断蒙特卡洛仿真是否结束：当Ns≥N时，蒙特卡洛仿真结束，停止模拟，反之，重新回到步骤1.1，开始新的光子的模拟。
步骤1主要用于存储有效种子数值，供步骤2使用。步骤2与步骤1的主要不同点在于，步骤1使用随机的初始种子数值，步骤2使用有效种子数值。因此，步骤2中模拟的光子，都将通过肿瘤区域。
步骤2、所述步骤1中在获得有效种子数值数组后，在治疗过程中，无论人体内部的肿瘤形状如何变化，均可利用有效种子数值，进行高速的二次蒙特卡洛仿真。利用有效种子数值驱动的光子，其传播路径中，必然将经过肿瘤区域。
步骤2.1有效随机数种子数值提取：从步骤1所存储的有效种子数值文件中，依次读取有效种子数值。当文件读取结束时，执行步骤2.5。反之，根据当前模拟的光子的序号，导入有效种子数值Seed_array[i]，其中i是一个自然数，其取值上限为Seed_array的数目，代表当前模拟的光子的序号；将该数值输入到传统的伪随机数产生函数中，将可以产生一个平均分布在(0,1)中的伪随机数ε，每产生一个伪随机数，种子数值将会由传统的伪随机数产生函数更新，并再次代入伪随机数产生函数，用于产生下一个平均分布在(0,1)中的伪随机数ε；上述过程不断循环，将会得到一个伪随机数数组，其统计特性为平均分布在(0,1)中；接下来的每一个步骤，都将使用伪随机数进行随机抽样，从而判断光子的传播方式；对于伪随机数列而言，当初始的种子数值确定了，其数列所产生的伪随机数数值也随之确定。我们可以将人体三维数值模型定义为三维数组Density，该三维数组是利用计算机断层扫描仪器(Computed Tomography,CT)获得的人体内部的影像信息。Density[i,j,k]表示在三维网格空间中，体素[i,j,k]区域对应CT密度值。Density[i,j,k]中的数值可分别表示血液、肌肉、皮肤等人体组织。另外，Density[i,j,k]中有部分区域是仪器与皮肤之间的空气区域。X射线机穿过空气，到达皮肤的距离，定义为源皮距(Source Skin Distance，SSD,SSD为大于0的自然数，表示光源到人体表面的距离)。在算法中，还将定义三维数组Dose[i,j,k]，用于记录人体内部的剂量分布。
步骤2.2光子发射：通过用户输入的入射光子坐标(Sxi,Syi,Szi)与光子发射方向(Dxi,Dyi,Dzi)以及源皮距SSD，光子沿着发射方向行进SSD，移动到人体的表面，此时光子的位置坐标为Sxi+SSD×Dxi,Syi+SSD×Dyi,Szi+SSD×Dzi，光子的初始能量定义为E；
步骤2.3光子迁移：对于能量小于500kev的X射线，光子在人体内部传播，在完成单次自由迁移后，将会产生瑞利散射、康普顿散射、光电吸收，在算法中，根据光子所在的体素对应的人体组织类型，通过查表的方式获得当前组织的三种平均自由程(mean free path，MFP)，分别为瑞利散射自由程MFP_r、康普顿散射自由程MFP_c以及光电效应自由程MFP_p。入射光子的能量以及光子到达区域的人体组织类型，都会影响平均自由程；在光子传播过程中，其发生瑞利散射、康普顿散射、光电吸收的概率，由MFP_r、MFP_c以及MFP_p的数值决定；瑞利散射不改变光子的能量E，只改变光子的运行方向Dxi,Dyi,Dzi，康普顿散射既改变光子的能量E，也改变光子的运行方向Dxi,Dyi,Dzi，能量的改变由传统的康普顿散射理论计算；当光子在体素[i,j,k]处发射光电吸收时，将使用式子Dose[i,j,k]＝Dose[i,j,k]+E将该体素吸收的光子能量记录在数组Dose[i,j,k]中；
步骤2.4边界处理：判断光子是否处于皮肤层，并即将进入空气层，如果不是，将回到步骤2.3，再次进行光子迁移；如果是，则光子到达人体模型的边界，将从人体内部出射，该光子的蒙特卡洛模拟将结束。然后，重新回到步骤2.1，开始新的光子的模拟；
步骤2.5数值计算结束。
在完成了2.1-2.5后，得到的人体内的剂量分布。
所述的步骤1.3和步骤2.3中的光子迁移采用如下方法：光子迁移自由程记为MFP_l，满足公式1/MFP_l＝1/MFP_r+1/MFP_c+1/MFP_p。自由程的步长的抽样表达式如下：l＝-Ln(ε)/MFP_l。光子行进了当前的自由程后，其空间位置为
Sxi＝Sxi+l×Dxi,
Syi＝Syi+l×Dyi,。
Szi＝Szi+l×Dzi
再次根据当前的位置Sxi,Syi,Szi。获取当前的瑞利散射自由程、康普顿散射自由程以及光电效应自由程。在完成了一次自由程的迁移后，将会依概率发生瑞利散射、康普顿散射或者光电吸收。定义1/MFP_l＝ul，1/MFP_r＝ur，1/MFP_c＝uc，1/MFP_p＝up再次生成伪随机数ε，当ε<ur/ul时，发生瑞利散射，光子的迁移方向将会发生改变，其中，迁移方向的改变由球坐标下的偏转角θ与方位角决定，偏转角θ的数值由瑞利散射公式进行随机抽样获取；当ur/ul<ε<(ur+uc)/ul，发生什么康普顿散射，光子能量E将会被衰减，光子迁移方向将会发生改变，其中，迁移方向的改变由球坐标下的偏转角θ与方位角决定，偏转θ的数值由康普顿散射公式进行随机抽样获取；当ε>(ur+uc)/ul发生光电吸收效应，光子能量被位置Sxi,Syi,Szi对应的体素[i,j,k]所吸收，将使用式子Dose[i,j,k]＝Dose[i,j,k]+E将该体素吸收的光子能量记录在数组Dose[i,j,k]中。
所述的产生方式如下：产生一个随机数ε；，那么ψ＝2πε；
所述的光子迁移方向的改变，在发生偏转角θ与方位角的改变后，光子迁移方向为：
在完成上面的计算后，更新光子迁移方向Dx＝Dx',Dy＝Dy',Dz＝Dz'。
所述的有效种子数值数组Seed_array，是一个整数数列，其中的每一个元素都为一个整数，代表了有效种子数值。
所述的瑞利散射偏转角，其概率抽样方法描述如下：
i.使用RITA算法[参考文献1]，抽样出满足概率分布函数[F(x,Z)]2的随机数x2。其中
另外，a＝α(Z-5/16)。α为物理学中的fine-structure常数，Z＝1to 99，由发生瑞利散射的原子种类所决定[参考文献2]。
ii.由x2获得其中me为电子质量常数，c为光速常数。
iii.产生随机数ε
iv.当回到步骤i
v.获得cosθ。
所述的康普顿散射偏转角，其概率抽样方法描述如下：
i.产生随机数ε，当ε>a1/(a1+a2)时，另k＝1，反之，k＝2；其中a1＝Ln(1+2κ)，其中me为电子质量常数，c为光速常数。
此时，我们将会把1.1中的种子数值定义为有效种子数值，并将其保存在有效种子数值数组Seed_array中，并在第二步蒙特卡洛仿真中使用；
ii.产生一个随机数ε，其中
iii.由τ，计算获得
iv.计算其中
在中，Ui表示第i层电子壳的电离能，fi表示i层电子壳的电子数，另外
在的表达式中，Ji；0＝Ji(0)。
其中，Ji(pz)＝∫∫ρi(P)dpxdpy，ρi(P)为原子第i层电子壳的动量概率分布[参考文献1]。
v.产生随机数ε，当ε>T(cosθ)时，回到步骤i
vi.获得cosθ。

Claims (6)

1.一种用于放射线剂量计算的双步蒙特卡洛模拟方法，其特征在于包括如下步骤：
步骤1、进行蒙特卡洛仿真，获取有效种子数值数组；
步骤1.1随机数的产生：在每一个光子的蒙特卡洛仿真开始之前，由当前时间的数值，组合为一个整数，定义为初始种子数值,将该数值输入到传统的伪随机数产生函数中，每产生一个伪随机数，种子数值将会由传统的伪随机数产生函数更新，并再次代入伪随机数产生函数，不断循环上述过程，将会得到一个伪随机数数组；
步骤1.2光子发射：通过入射光子初始坐标与光子发射方向以及放射线产生器与人体之间的源皮距，得到光子到达人体表面时的位置坐标：
Sxi+SSD×Dxi,Syi+SSD×Dyi,Szi+SSD×Dzi，
式中Sxi,Syi,Szi为入射光子初始坐标，Dxi,Dyi,Dzi为光子发射方向，SSD为源皮距，光子的初始能量定义为E；
步骤1.3光子迁移：对于能量小于500kev的X射线，光子在人体内部传播，将会产生瑞利散射、康普顿散射、光电吸收，并以此判断光子单次自由迁移的状态，根据光子所在的体素对应的人体组织类型，通过查表的方式获得当前组织的三种平均自由程，分别为瑞利散射自由程MFP_r、康普顿散射自由程MFP_c以及光电效应自由程MFP_p，上述三个自由程将决定光子单次自由迁移的距离，以及完成单次自由迁移后，发生的散射或吸收的概率，从而判断下一次自由迁移的初始状态；
步骤1.4有效种子数值的判定：在光子迁移的过程中，如果光子经过了用户预设的肿瘤区域，将会把步骤1.1中的初始种子数值定义为有效种子数值；反之，则不做此定义；
步骤1.5边界处理：判断光子是否处于皮肤层，并即将进入空气层，如果不是，将回到步骤1.3，再次进行光子迁移；如果是，则光子到达人体模型的边界，将从人体内部出射，该光子的蒙特卡洛模拟将结束，如果本次蒙特卡洛模拟的初始种子数值被定义为有效种子数值，就将该数值存储，供第二步蒙特卡洛仿真中使用；反之，该初始种子数值将不做任何处理；
1.6判断蒙特卡洛仿真是否结束：当已模拟的光子数目超出预设值时，蒙特卡洛仿真结束，停止模拟，反之，重新回到步骤1.1，开始新的光子的模拟，同时，已模拟的光子数目加一；
步骤2、所述步骤1.5中在获得有效种子数值数组后，利用有效种子数值，进行二次蒙特卡洛仿真，高速的计算人体内部的放射线剂量分布；
步骤2.1有效随机数种子数值提取：从步骤1所存储的有效种子数值文件中，依次读取有效种子数值，将有效随机数种子数值依次输入到传统的伪随机数产生函数中，将产生一个伪随机数，每产生一个伪随机数，种子数值将会由传统的伪随机数产生函数更新，并再次代入伪随机数产生函数，用于产生下一个伪随机数，接下来的每一个步骤，都将使用伪随机数进行随机抽样，从而判断光子的传播方式，当有效种子数值读取结束时，进入步骤2.5；
步骤2.2光子发射：通过入射光子初始坐标与光子发射方向以及放射线产生器与人体之间的源皮距，得到此时光子的位置坐标：
Sxi+SSD×Dxi,Syi+SSD×Dyi,Szi+SSD×Dzi，
式中Sxi,Syi,Szi为入射光子坐标，Dxi,Dyi,Dzi为光子发射方向，SSD为源皮距，光子的初始能量定义为E
步骤2.3光子迁移：对于能量小于500kev的X射线，光子在人体内部传播，将会产生瑞利散射、康普顿散射、光电吸收，并以此判断光子单次自由迁移的状态，在算法中，根据光子所在的体素对应的人体组织类型，通过查表的方式获得当前组织的三种平均自由程，分别为瑞利散射自由程MFP_r、康普顿散射自由程MFP_c以及光电效应自由程MFP_p，上述三个自由程将决定光子单次自由迁移的距离，以及完成单次自由迁移后，发生的散射或吸收的概率，从而判断下一次自由迁移的初始状态；当光子在体素[i,j,k]处发生光电吸收时，将使用式子Dose[i,j,k]＝Dose[i,j,k]+E将该体素吸收的光子能量记录在数组Dose[i,j,k]中；
步骤2.4边界处理：判断光子是否处于皮肤层，并即将进入空气层，如果不是，将回到步骤2.3，再次进行光子迁移；如果是，则光子到达人体模型的边界，将从人体内部出射，该光子的蒙特卡洛模拟将结束；然后，重新回到步骤2.1，开始新的光子的模拟；
步骤2.5数值计算结束，得到的人体内的剂量分布；
所述的步骤1.3或步骤2.3中光子迁移步骤，入射光子的能量以及光子到达区域的人体组织类型，会影响平均自由程；在光子传播过程中，其发生瑞利散射、康普顿散射、光电吸收的概率，由MFP_r、MFP_c以及MFP_p的数值决定；瑞利散射不改变光子的能量E，只改变光子的运行方向Dxi,Dyi,Dzi，康普顿散射既改变光子的能量E，也改变光子的运行方向Dxi,Dyi,Dzi，光电吸收发生后，光子的能量将会被吸收，该光子的模拟结束。
2.根据权利要求1所述的用于放射线剂量计算的双步蒙特卡洛模拟方法，其特征在于：所述的步骤1.1随机数的产生具体产生方法如下：由当前时间组合为一个种子数值Seed，该数值的前四位为当前年份，中间四位为当前日期，最后四位为当前时间，精确到分钟；将该数值输入到传统的伪随机数产生函数中，将产生一个平均分布在(0,1)中的伪随机数ε，每产生一个伪随机数，种子数值将会由传统的伪随机数产生函数更新，并再次代入伪随机数产生函数，用于产生下一个平均分布在(0,1)中的伪随机数ε；上述过程不断循环，将会得到一个伪随机数数组，其统计特性为平均分布在(0,1)中；接下来的每一个步骤，都将使用伪随机数进行随机抽样，从而判断光子的传播方式；对于伪随机数列而言，当初始的种子数值确定了，其数列所产生的伪随机数数值也随之确定。
3.根据权利要求1所述的用于放射线剂量计算的双步蒙特卡洛模拟方法，其特征在于：所述的步骤2.1有效随机数种子数值提取具体方法如下：
根据当前模拟的光子的序号，导入有效种子数值Seed_array[i]，其中i是一个自然数，其取值上限为Seed_array的数目，代表当前模拟的光子的序号；将该数值输入到传统的伪随机数产生函数中，将产生一个平均分布在(0,1)中的伪随机数ε，每产生一个伪随机数，种子数值将会由传统的伪随机数产生函数更新，并再次代入伪随机数产生函数，用于产生下一个平均分布在(0,1)中的伪随机数ε；上述过程不断循环，将会得到一个伪随机数数组，其统计特性为平均分布在(0,1)中；接下来的每一个步骤，都将使用伪随机数进行随机抽样，从而判断光子的传播方式；对于伪随机数列而言，当初始的种子数值确定了，其数列所产生的伪随机数数值也随之确定，因此，由有效种子数值驱动的光子，都将经过肿瘤区域。
4.根据权利要求1所述的用于放射线剂量计算的双步蒙特卡洛模拟方法，其特征在于:所述的步骤1.1或步骤2.1中，将人体三维数值模型定义为三维数组Density，该三维数组是利用计算机断层扫描仪器获得的人体内部的影像信息；Density[i,j,k]表示在三维网格空间中，体素[i,j,k]区域对应CT密度值；Density[i,j,k]中的数值分别表示人体组织；另外，Density[i,j,k]中有部分区域是仪器与皮肤之间的空气区域；X射线机穿过空气，到达皮肤的距离，定义为源皮距SSD，SSD为大于0的自然数，表示光源到人体表面的距离，定义三维数组Dose[i,j,k]，用于存储人体内的剂量分布，定义蒙特卡洛仿中，需要模拟的光子数量为N，定义已模拟的光子数为变量Ns，Ns的初始值为0，在模拟方法的循环中，每一次运行步骤1.1时，Ns＝Ns+1。
5.根据权利要求4所述的用于放射线剂量计算的双步蒙特卡洛模拟方法，其特征在于：所述的光子迁移步骤中的光子迁移采用如下方法：光子迁移自由程记为MFP_l，满足公式1/MFP_l＝1/MFP_r+1/MFP_c+1/MFP_p；自由程的步长的抽样表达式如下：l＝-Ln(ε)/MFP_l；光子行进了当前的自由程后，其空间位置为：
Sxi＝Sxi+l×Dxi,
Syi＝Syi+l×Dyi,
Szi＝Szi+l×Dzi
再次根据当前的位置Sxi,Syi,Szi，获取当前的瑞利散射自由程、康普顿散射自由程以及光电效应自由程；在完成了一次自由程的迁移后，将会依概率发生瑞利散射、康普顿散射或者光电吸收。
6.根据权利要求5所述的用于放射线剂量计算的双步蒙特卡洛模拟方法，其特征在于：所述的发生瑞利散射、康普顿散射或者光电吸收过程中：
定义1/MFP_l＝ul，1/MFP_r＝ur，1/MFP_c＝uc，1/MFP_p＝up再次生成伪随机数ε，当ε＜ur/ul时，发生瑞利散射，光子的迁移方向将会发生改变，其中，迁移方向的改变由球坐标下的偏转角θ与方位角决定，偏转角θ的数值由瑞利散射公式进行随机抽样获取；当ur/ul＜ε＜(ur+uc)/ul，发生康普顿散射，光子能量E将会被衰减，光子迁移方向将会发生改变，其中，迁移方向的改变由球坐标下的偏转角θ与方位角决定，偏转θ的数值由康普顿散射公式进行随机抽样获取；当ε＞(ur+uc)/ul发生光电吸收效应，光子能量被位置Sxi,Syi,Szi对应的体素[i,j,k]所吸收，将使用式子Dose[i,j,k]＝Dose[i,j,k]+E将该体素吸收的光子能量记录在数组Dose[i,j,k]中；
所述的即方位角产生方式如下：产生一个随机数ε；，那么
所述的光子迁移方向的改变，在发生偏转角θ与方位角的改变后，光子迁移方向为：
在完成上面的计算后，更新光子迁移方向Dx＝Dx′,Dy＝Dy′,Dz＝Dz′。
CN201610115951.2A 2016-03-02 2016-03-02 一种用于放射线剂量计算的双步蒙特卡洛模拟方法 Active CN105740573B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610115951.2A CN105740573B (zh) 2016-03-02 2016-03-02 一种用于放射线剂量计算的双步蒙特卡洛模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610115951.2A CN105740573B (zh) 2016-03-02 2016-03-02 一种用于放射线剂量计算的双步蒙特卡洛模拟方法

Publications (2)

Publication Number Publication Date
CN105740573A CN105740573A (zh) 2016-07-06
CN105740573B true CN105740573B (zh) 2019-10-11

Family

ID=56248870

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610115951.2A Active CN105740573B (zh) 2016-03-02 2016-03-02 一种用于放射线剂量计算的双步蒙特卡洛模拟方法

Country Status (1)

Country Link
CN (1) CN105740573B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107292075B (zh) * 2016-04-06 2020-06-12 南京中硼联康医疗科技有限公司 增进放射治疗系统计算效益的方法
CN106291650A (zh) * 2016-08-31 2017-01-04 广州市岱尼欣贸易有限公司 基于蒙特卡罗的剂量测量方法
EP3338859B1 (en) * 2016-12-22 2019-06-26 RaySearch Laboratories AB System and method for modelling of scattering in ion radiotherapy treatment planning
CN106932810B (zh) * 2017-04-01 2018-02-23 西安一体医疗科技有限公司 一种伽玛射线剂量的卷积计算方法
CN108618796A (zh) * 2018-02-09 2018-10-09 南方医科大学 一种基于任务驱动的蒙特卡洛散射光子模拟方法
CN116973966A (zh) * 2019-06-21 2023-10-31 清华大学 辐射分析方法、装置和计算机可读存储介质
CN110702706B (zh) * 2019-09-20 2022-05-20 天津大学 一种能谱ct系统输出数据的模拟方法
CN110911007B (zh) * 2019-12-29 2023-07-25 杭州科洛码光电科技有限公司 基于成像光谱仪的生物组织光学参数重构方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103295188A (zh) * 2012-02-28 2013-09-11 上海联影医疗科技有限公司 X射线蒙特卡洛模拟的路径积分方法
CN103699816A (zh) * 2013-12-12 2014-04-02 深圳先进技术研究院 基于蒙特卡洛模拟的蛋白质热力学分析方法
CN105335215A (zh) * 2015-12-05 2016-02-17 中国科学院苏州生物医学工程技术研究所 一种基于云计算的蒙特卡洛仿真加速方法及系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103295188A (zh) * 2012-02-28 2013-09-11 上海联影医疗科技有限公司 X射线蒙特卡洛模拟的路径积分方法
CN103699816A (zh) * 2013-12-12 2014-04-02 深圳先进技术研究院 基于蒙特卡洛模拟的蛋白质热力学分析方法
CN105335215A (zh) * 2015-12-05 2016-02-17 中国科学院苏州生物医学工程技术研究所 一种基于云计算的蒙特卡洛仿真加速方法及系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
应用Geant4模拟放疗、诊断和防护线质下的X射线能谱;李明生 等;《中国辐射卫生》;20120331;第21卷(第1期);摘要，引言，第1节 *
蒙特卡洛仿真的加速与电场矢量化及其在生物光学成像中的应用;蔡夫鸿;《中国博士学位论文全文数据库 基础科学辑》;20130815(第2013年第08期);正文第1章，第2章，第5章 *

Also Published As

Publication number Publication date
CN105740573A (zh) 2016-07-06

Similar Documents

Publication Publication Date Title
CN105740573B (zh) 一种用于放射线剂量计算的双步蒙特卡洛模拟方法
US10692283B2 (en) Geometric model establishment method based on medical image data
Doolan et al. Patient-specific stopping power calibration for proton therapy planning based on single-detector proton radiography
CN104066479B (zh) 使用异质性补偿迭加进行用于放射疗法的剂量计算
CN107392977B (zh) 单视图切伦科夫发光断层成像重建方法
US10223480B2 (en) Method for modeling and accounting for cascade gammas in images
RU2736917C1 (ru) Способ анализа элементов и отношений масс элементов ткани и способ построения геометрической модели на основе медицинского изображения
CN102945328A (zh) 基于gpu并行运算的x射线造影图像仿真方法
Chao et al. Organ dose conversion coefficients for 0.1–10 MeV electrons calculated for the VIP-Man tomographic model
Lemaréchal et al. GGEMS-Brachy: GPU GEant4-based Monte Carlo simulation for brachytherapy applications
CN108310677A (zh) 基于医学影像数据的平滑几何模型建立方法
CN108451508A (zh) 基于多层感知机的生物自发荧光三维成像方法
CN110097580A (zh) 一种超声图像标志物运动追踪方法
Penfold Image reconstruction and Monte Carlo simulations in the development of proton computed tomography for applications in proton radiation therapy
Lazos et al. A software data generator for radiographic imaging investigations
CN104258505A (zh) 肿瘤放射治疗剂量个体化验证仿真体模及其建立和应用
Jones et al. Characterization of Compton-scatter imaging with an analytical simulation method
Jiang et al. 3D in vivo dose verification in prostate proton therapy with deep learning-based proton-acoustic imaging
CN106408543A (zh) 一种用于锥束ct图像散射修正的阻挡光栅优化方法及装置
CN107480445A (zh) 基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法
Zhang et al. Advantages of MCNPX-based lattice tally over mesh tally in high-speed Monte Carlo dose reconstruction for proton radiotherapy
CN104605823A (zh) 一种基于pfMC模型的荧光扩散光学断层图像重建方法
Zheng et al. GPU accelerated stochastic origin ensemble method with list-mode data for compton camera imaging in proton therapy
CN105073006A (zh) 针对计算机断层摄影检查的x射线剂量分布计算
Gilling A GATE Monte Carlo dose analysis from Varian XI conebeam computed tomography

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20230727

Address after: 318050 room 313, building 1, No. 199, Zhuguang street, Lubei street, Luqiao District, Taizhou City, Zhejiang Province

Patentee after: Taizhou anqiling Intelligent Technology Co.,Ltd.

Address before: 215500 Building 1, Dongnan Hui, Xianshi Road, High tech Zone, Changshu, Suzhou, Jiangsu

Patentee before: SUZHOU WANGHAO INFORMATION TECHNOLOGY Co.,Ltd.