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

摘要

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

说明书

技术领域

本发明涉及人体内放射剂量强度分布的数值计算，主要是用于模拟X射线在人体内部的传播以及剂量分布。

背景技术

放射治疗是一种重要的肿瘤治疗手段，70％的肿瘤患者在病程的不同时期，需要接受放射治疗。然而，在X射线照射肿瘤的时候，也对正常的器官和组织造成了影响。为了提高治疗效率，并降低对正常组织的损伤。在放射治疗之前，需要计算X射线在人体内部的分布，用于辅助医护人员确定肿瘤区域的X射线剂量是否达到了治疗剂量，正常组织的放射线剂量是否在安全阈值以下。

使用能量输运理论，可以准确的计算出X射线在人体内部的剂量分布。然而，人体是极度不均匀的介质，需要将人体离散为大量的体素，才能精确的构建人体的数字模型。这使得能量输运理论的数值求解过程中，需要大内存计算机，并且计算速度也很慢。蒙特卡洛是一种灵活的数值模拟方法，可以模拟X射线在复杂的人体内的传播过程。并且，根据概率统计的大数定律，随着模拟X射线的数目的增多，蒙特卡洛仿真的计算结果可以趋近能量输运理论的计算结果。

蒙特卡洛仿真的引入，解决了能量输运理论中产生的大内存使用的问题，但计算速度很慢。当前，人们引入显卡并行加速的方法提高了计算效率。同时，人们也在不断的提出新的改良算法，提高蒙特卡洛仿真的计算效率。在放射治疗的过程中，肿瘤的结构尺寸会随着治疗的进行而发生改变。人们需要一种高效的X射线的蒙特卡洛算法，可以在放射治疗的过程中，快速的计算体内剂量分布。

双步蒙特卡洛方法在博士论文《蒙特卡洛仿真的加速与电场矢量化及其在生物光学成像中的应用》中已有提及，但其应用点为可见光波段，本专利应用于X射线领域。针对于人体组织，可见光波段将受到极强的瑞利散射，所以并不适于人体组织，而X射线受到的瑞利散射较弱。另外，X射线在人体内部传播时，还将产生康普顿散射，这是可见光在人体内部传播时，不会发生的。另一方面，本发明的后期应用主要为肿瘤放射治疗中的X射线剂量分布，因此，将主要研究肿瘤区域的X射线剂量。这与可见光波段的双步蒙特卡洛方法是不同的。本发明会用到以下参考文献：

[1]SalvatF,Fernández-VareaJM,SempauJ.PENELOPE-2006:AcodesystemforMonteCarlosimulationofelectronandphotontransport[C]//WorkshopProceedings.2006,4:7.

[2]HubbellJH,VeigeleWJ,BriggsEA,etal.Atomicformfactors,incoherentscatteringfunctions,andphotonscatteringcrosssections[J].Journalofphysicalandchemicalreferencedata,1975,4(3):471-538.

发明内容

1、本发明的目的：

本发明提出一种双步蒙特卡洛仿真方法，用于在放射治疗的过程中，对于同一个模拟对象的体内的X射线剂量分布，进行高速的数值计算。本方法的核心思想是，因为X射线是一种穿透性能很强的电磁波，其在人体内的传播轨迹，因为人体肿瘤结构的改变而发生变化的概率极小，利用这一特点，我们提出了双步蒙特卡洛仿真，用于仿真X射线在人体内部的传播以及剂量分布。

2、本发明所采用的技术方案：

因为X射线也属于电磁波范畴，在蒙特卡洛仿真中，X射线被命名为光子，这是X射线的蒙特卡洛仿真的惯用命名方法。用于能量小于500kev的放射线剂量计算的双步蒙特卡洛模拟方法，包括如下步骤：

步骤1、进行蒙特卡洛仿真，获取有效种子数值数组；

步骤2、所述步骤一中在获得有效种子数值数组后，利用有效种子数值，进行二次蒙特卡洛仿真，高速的计算人体内部的放射线剂量分布。

更进一步的具体实施例中，所述的步骤1中包括如下步骤：

步骤1.1随机数的产生：在每一个光子的蒙特卡洛仿真开始之前，由当前时间的数值，组合为一个整数，定义为初始种子数值,将该数值输入到传统的伪随机数产生函数中，每产生一个伪随机数，种子数值将会由传统的伪随机数产生函数更新，并再次代入伪随机数产生函数，不断循环上述过程，将会得到一个伪随机数数组；

步骤1.2光子发射：通过入射光子初始坐标与光子发射方向以及放射线产生器与人体之间的源皮距，得到光子到达人体表面时的位置坐标：

Sx_i+SSD×Dx_i,Sy_i+SSD×Dy_i,Sz_i+SSD×Dz_i，

式中Sx_i,Sy_i,Sz_i为入射光子初始坐标，Dx_i,Dy_i,Dz_i为光子发射方向，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光子发射：通过入射光子初始坐标与光子发射方向以及放射线产生器与人体之间的源皮距，得到此时光子的位置坐标：

Sx_i+SSD×Dx_i,Sy_i+SSD×Dy_i,Sz_i+SSD×Dz_i，

式中Sx_i,Sy_i,Sz_i为入射光子坐标，Dx_i,Dy_i,Dz_i为光子发射方向，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，只改变光子的运行方向Dx_i,Dy_i,Dz_i，康普顿散射既改变光子的能量E，也改变光子的运行方向Dx_i,Dy_i,Dz_i，光电吸收发生后，光子的能量将会被吸收，该光子的模拟结束。

更进一步的具体实施例中，所述的光子迁移步骤中的光子迁移采用如下方法：光子迁移自由程记为MFP_l，满足公式1/MFP_l＝1/MFP_r+1/MFP_c+1/MFP_p；自由程的步长的抽样表达式如下：l＝-Ln(ε)/MFP_l；光子行进了当前的自由程后，其空间位置为：

Sx_i＝Sx_i+l×Dx_i,

Sy_i＝Sy_i+l×Dy_i,

Sz_i＝Sz_i+l×Dz_i

再次根据当前的位置Sx_i,Sy_i,Sz_i，获取当前的瑞利散射自由程、康普顿散射自由程以及光电效应自由程；在完成了一次自由程的迁移后，将会依概率发生瑞利散射、康普顿散射或者光电吸收。

更进一步的具体实施例中，所述的发生瑞利散射、康普顿散射或者光电吸收过程中：

定义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发生光电吸收效应，光子能量被位置Sx_i,Sy_i,Sz_i对应的体素[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，该三维数组是利用计算机断层扫描仪器(ComputedTomography,CT)获得的人体内部的影像信息。Density[i,j,k]表示在三维网格空间中，体素[i,j,k]区域对应CT密度值。Density[i,j,k]中的数值可分别表示血液、肌肉、皮肤等人体组织。另外，Density[i,j,k]中有部分区域是仪器与皮肤之间的空气区域。X射线机穿过空气，到达皮肤的距离，定义为源皮距(SourceSkinDistance，SSD,SSD为大于0的自然数，表示光源到人体表面的距离)。在算法中，还将定义三维数组Dose[i,j,k]，用于记录人体内部的剂量分布，定义蒙特卡洛仿中，需要模拟的光子数量为N，定义已模拟的光子数为变量Ns，Ns的初始值为0，再每一次进行步骤1.1时，Ns＝Ns+1；

步骤1.2光子发射：通过用户输入的入射光子坐标(Sx_i,Sy_i,Sz_i)与光子发射方向(Dx_i,Dy_i,Dz_i)以及源皮距SSD，光子沿着发射方向行进SSD，移动到人体的表面，此时光子的位置坐标为，Sx_i+SSD×Dx_i,Sy_i+SSD×Dy_i,Sz_i+SSD×Dz_i，光子的初始能量定义为E；

步骤1.3光子迁移：对于能量小于500kev的X射线，光子在人体内部传播，将会产生瑞利散射、康普顿散射、光电吸收，并以此判断光子单次自由迁移的状态，在算法中，根据光子所在的体素对应的人体组织类型，通过查表的方式获得当前组织的三种平均自由程(meanfreepath，MFP)，分别为瑞利散射自由程MFP_r、康普顿散射自由程MFP_c以及光电效应自由程MFP_p。入射光子的能量以及光子到达区域的人体组织类型，都会影响平均自由程；在光子传播过程中，完成每一次自由迁移后，其发生瑞利散射、康普顿散射、光电吸收的概率，由MFP_r、MFP_c以及MFP_p的数值决定；瑞利散射不改变光子的能量E，只改变光子的运行方向Dx_i,Dy_i,Dz_i，康普顿散射既改变光子的能量E，也改变光子的运行方向Dx_i,Dy_i,Dz_i，能量的改变由传统的康普顿散射理论计算；当光子在体素[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，该三维数组是利用计算机断层扫描仪器(ComputedTomography,CT)获得的人体内部的影像信息。Density[i,j,k]表示在三维网格空间中，体素[i,j,k]区域对应CT密度值。Density[i,j,k]中的数值可分别表示血液、肌肉、皮肤等人体组织。另外，Density[i,j,k]中有部分区域是仪器与皮肤之间的空气区域。X射线机穿过空气，到达皮肤的距离，定义为源皮距(SourceSkinDistance，SSD,SSD为大于0的自然数，表示光源到人体表面的距离)。在算法中，还将定义三维数组Dose[i,j,k]，用于记录人体内部的剂量分布。

步骤2.2光子发射：通过用户输入的入射光子坐标(Sx_i,Sy_i,Sz_i)与光子发射方向(Dx_i,Dy_i,Dz_i)以及源皮距SSD，光子沿着发射方向行进SSD，移动到人体的表面，此时光子的位置坐标为Sx_i+SSD×Dx_i,Sy_i+SSD×Dy_i,Sz_i+SSD×Dz_i，光子的初始能量定义为E；

步骤2.3光子迁移：对于能量小于500kev的X射线，光子在人体内部传播，在完成单次自由迁移后，将会产生瑞利散射、康普顿散射、光电吸收，在算法中，根据光子所在的体素对应的人体组织类型，通过查表的方式获得当前组织的三种平均自由程(meanfreepath，MFP)，分别为瑞利散射自由程MFP_r、康普顿散射自由程MFP_c以及光电效应自由程MFP_p。入射光子的能量以及光子到达区域的人体组织类型，都会影响平均自由程；在光子传播过程中，其发生瑞利散射、康普顿散射、光电吸收的概率，由MFP_r、MFP_c以及MFP_p的数值决定；瑞利散射不改变光子的能量E，只改变光子的运行方向Dx_i,Dy_i,Dz_i，康普顿散射既改变光子的能量E，也改变光子的运行方向Dx_i,Dy_i,Dz_i，能量的改变由传统的康普顿散射理论计算；当光子在体素[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。光子行进了当前的自由程后，其空间位置为

Sx_i＝Sx_i+l×Dx_i,

Sy_i＝Sy_i+l×Dy_i,。

Sz_i＝Sz_i+l×Dz_i

再次根据当前的位置Sx_i,Sy_i,Sz_i。获取当前的瑞利散射自由程、康普顿散射自由程以及光电效应自由程。在完成了一次自由程的迁移后，将会依概率发生瑞利散射、康普顿散射或者光电吸收。定义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发生光电吸收效应，光子能量被位置Sx_i,Sy_i,Sz_i对应的体素[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)]²的随机数x²。其中

另外，a＝α(Z-5/16)。α为物理学中的fine-structure常数，Z＝1to99，由发生瑞利散射的原子种类所决定[参考文献2]。

ii.由x²获得其中m_e为电子质量常数，c为光速常数。

iii.产生随机数ε

iv.当回到步骤i

v.获得cosθ。

所述的康普顿散射偏转角，其概率抽样方法描述如下：

i.产生随机数ε，当ε>a₁/(a₁+a₂)时，另k＝1，反之，k＝2；其中a₁＝Ln(1+2κ)，

其中m_e为电子质量常数，c为光速常数。

此时，我们将会把1.1中的种子数值定义为有效种子数值，并将其保存在有效种子数值数组Seed_array中，并在第二步蒙特卡洛仿真中使用；

ii.产生一个随机数ε，其中

iii.由τ，计算获得

iv.计算其中

在中，U_i表示第i层电子壳的电离能，f_i表示i层电子壳的电子数，另外

在的表达式中，J_i；0＝J_i(0)。

其中，J_i(p_z)＝∫∫ρ_i(P)dp_xdp_y，ρ_i(P)为原子第i层电子壳的动量概率分布[参考文献1]。

v.产生随机数ε，当ε>T(cosθ)时，回到步骤i

vi.获得cosθ。