CN101477203B - 一种解析蒙特卡罗剂量计算方法 - Google Patents

Abstract

本发明公开了一种解析蒙特卡罗剂量计算方法，在传统的DSA剂量计算模型的基础上，使用三维非均匀修正替代传统DSA方法的一维修正，提高了计算的精度；使用子束入射到介质表面的能量注量计算所有深度的比释总能，替代传统DSA方法单独计算各个深度能量注量，减小各深度计算的采样点数目，从而减小计算时间；根据不同的计算区域，分别采用不同的计算点密度：感兴趣区域的计算点密度高，其他区域计算点密度低，在保证用户关心区域计算精度的同时，减少计算时间。

Description

一种解析蒙特卡罗剂量计算方法

技术领域

本发明涉及辐射剂量测量和计算中一种解析蒙特卡罗剂量计算方法。

背景技术

剂量计算方法一般分为解析方法和蒙特卡罗方法。解析方法具有快速和在均匀区域精度较高的优势，但在组织非均匀区域存在较大误差；蒙特卡罗方法在所有区域均可保持较高精度，但由于较为耗时，限制了应用范围。

蒙特卡罗模拟得到的剂量分布阵列DSA(Dose Spread Array)可以表征在某一位置(Location)处由射线与介质之间的相互作用所产生的带电粒子的运动在介质中所造成的能量吸收(剂量)。DSA模型可作为独立的解析蒙特卡罗剂量计算方法使用。

传统DSA模型在计算非均匀介质剂量时，只考虑了作用栅元的电子密度和作用栅元与沉积栅元连线上各栅元电子密度平均值的比值来修正非均匀组织对剂量分布的影响，而没有充分考虑带电粒子在介质中运动的真实情况，即带电粒子并不是完全沿作用栅元到沉积栅元之间的直线运动。

传统DSA模型在计算介质剂量时，需单独计算各个深度的能量注量，较为耗时。此外，在均匀采样点设置的前提下，要在感兴趣区域获得较高的剂量，需要较长的时间。

发明内容

本发明的目的在于克服传统DSA计算模型在计算精度和计算速度上的不足，提供一种解析蒙特卡罗剂量计算方法，不仅能得到更好的精度，而且能得到更好的计算速度。

本发明的技术方案如下：

一种解析蒙特卡罗剂量计算方法，其特征在于包括以下步骤：

(1)获得计算参数，包括以下内容：

蒙特卡罗程序(例如EGSnrc的子程序edknrc)获得的能量沉积核h：h是一系列单能或一定能谱的光子入射到某一均匀介质中，对应于球坐标(θ，R)的辐射能量值矩阵；

一系列单能或一定能谱的光子在某一介质中的质量衰减系数μ/ρ；

一系列单能或一定能谱的电子在某一介质中的射程；

计算模型栅元划分信息；

计算模型栅元物理密度信息；

计算模型栅元电子密度信息；

计算模型中计算区域边界标示信息；

计算模型中感兴趣区域标示文件；

放射源的能谱信息，用我们发展的能谱反演软件获得；

射线在介质表面形成的形状，由用户给出；

(2)利用下列MDSA模型进行计算：

使用MDSA模型计算公式，如下：

$D\left(r\right)={\∫}_E\∫\∫{\∫}_V{T}_E\left(r^{\′}\right)\frac{\mathrm{\ρ}\left(r^{\′}\right)}{{\mathrm{\ρ}}_0}{c}^{2}h(E,c_{(r,r^{\′})}(r-r^{\′}))d^3{r}^{\′}\mathrm{dE}---\left(1\right)$

其中，T_E(r′)为微分比释总能；ρ(r′)是作用栅元的物理密度；h为预先计算好的物理密度为ρ₀的均匀介质中的能量沉积核；c_(r，r′)代表了三维非均匀修正系数，是以某个作用栅元为圆心，在垂直于作用栅元到沉积栅元连线方向，以带电粒子在介质中以射程为半径的圆面为底，以作用栅元到沉积栅元连线为高，所得的圆柱体内的栅元电子密度的平均值：

$c=c_{(r,r^{\′})}={\∫}_0^1{\mathrm{\ρ}}_{\mathrm{rel}}[r^{\′}-l(r^{\′}-r)]\mathrm{dl}---\left(2\right)$

T_E(r′)使用子束入射到介质表面的能量注量计算所有深度的微分比释总能，即根据各个小子束在介质表面的划分，随着深度的增加，每个小野的尺寸按比例扩大，保持小野数目不变，相对于传统DSA方法，减小各深度计算的采样点数目，达到减小计算时间的目的；

在使用MDSA模型计算时，根据不同的计算区域，分别采用不同的计算点密度：感兴趣区域的计算点密度高，其他区域计算点密度低。

所述的三维非均匀修正系数，其特征在于以某个作用栅元为圆心，在垂直于作用栅元到沉积栅元连线方向，以带电粒子在介质中以射程为半径的圆面为底，以作用栅元到沉积栅元连线为高，所得的圆柱体内的栅元电子密度的平均值替代传统DSA模型一维修正，提高了计算的精度。

所述的微分比释总能，其特征在于使用各个小子束在计算模型介质表面的划分，随着深度的增加，每个小野的尺寸按比例扩大，保持小野数目不变，相对于传统的DSA方法，减小各深度计算的采样点数目，达到减小计算时间的目的。

所述的在不同的计算区域采用不同的计算点密度，其特征在于感兴趣区域的计算点密度高，其他区域计算点密度低。这种方法可以在保证用户关心区域计算精度的同时，减少计算时间。

加速器的能谱信息；可利用我们发展的反演能谱软件获得；

一系列单能或一定能谱的光子在某一介质中的质量衰减系数μ/ρ；可在美国标准化实验室(NIST)官方网站获得；

一系列单能或一定能谱的电子在某一介质中的射程；可在美国标准化实验室(NIST)官方网站获得；

本发明的剂量计算方法，不仅能得到更好的精度，而且能得到更好的计算速度。

附图说明

图1是本发明算法实现示意图。

图2是本发明的计算模型三维非均匀修正示意图。

图3是本发明计算实例图。

具体实施方式

实现本发明的计算过程，包括以下内容：

1.精确多算法光子加速器源反演方法，通过如下步骤实现：

第一步，获得水模中单能深度剂量曲线PPD，单能深度剂量曲线是通过国际公开使用的蒙特卡罗程序EGSnrc模拟单向垂直入射的光子在10cm×10cm或者更大射野，射野的能量沉积获得，能量从0～30MeV或者更大能量；

第二步，利用三维水箱测量获取百分深度剂量曲线；利用自动水箱测量加速器PDD曲线：把加速器头和机架角度调整到0度，水箱放置在水平地面，水箱中心点与加速器头中心点连线垂直于水平地面，加速器开野大小设置为10cm×10cm或者更大射野；利用加速器的电离室探头按照一定的深度间隔扫描；

第三步，加速器源反演；

(1)加速器源的光子能谱采用如下的数学模型实现：

$\mathrm{Min\σ}=\sqrt{\frac{1}{m-(N+4)}\underset{j}{\overset{m}{\mathrm{\Σ}}}{[D^{\′}\left(z_j\right)-D\left(z_j\right)]}^2}---\left(1\right)$

其中，

D′(z_j)为根据反演出来的能谱重建出来的PDD曲线；

$D\left(z\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\left|a_i\right|D(E_i,z)+\left|D_c\right|(e^{-\left|{\mathrm{\μ}}_p\right|z}-|v\left|e^{-\left|{\mathrm{\μ}}_e\right|z}\right);$

σ为均方根误差(Root Mean Squared Error)；

m为用于拟合时的测量PDD数据个数；

N+4为需要拟合的系数个数；N为能群的个数；

D(E_i，z)为能量为E_i的深度剂量曲线；

能谱为Ф(E)＝|a_i|；其中a_i可以是常系数也可以是带有待定系数的经验公式；

D_c为归一化系数；

μ_p为伴随射线入射的高能散射光平均衰减系数；

μ_e为在轫致辐射光子和次级光子衰减系数；

v为与入射表面剂量有关的系数。

(2)加速器源信息反演：通过利用传统成熟的非线性反演算法包括：

C1：Levenberg-Marquardt；C2：Quasi-Newton；C3：Gradient；C4：Conjugate-Gradient；C5：Newton；C6：Principal-Axis；C7：NMinimize算法，根据已知的测量PDD曲线D(z)和单能PDD曲线D(E，z)，进行数据拟合求得光子能谱Φ(E)＝|a_i|；

第四步，多算法结果比较；根据以上提到的反演算法的计算结果：光子能谱以及对应的均方根误差，综合评价最优的计算结果；综合评价方法如下：

测量PDD与根据反演所得能谱反演计算D(z)对比；选取计算D(z)与测量PDD在“建成”区附近的均方根误差比较小，采用公式(1)计算误差，并且总的均方根误差比较小的计算结果；

第五步，最优结果输出与显示；将以上比较的最优结果同时以图像以及数据的形式输出。

根据反演得到的需要计算的光子能谱，通过介质上的作用点(x₁，y₁，z₁)对沉积点(x₂，y₂，z₂)的剂量贡献(如图1所示)：

根据作用点R₁(x₁，y₁，z₁)的坐标，可计算出到达作用点R₁处的原射线与介质表面的交点R₀坐标，进而求得放射源(x_s，y_s，z_s)到R₀处的距离r₀，根据(x_s，y_s，z_s)的相对微分注量及距离平方反比公式得到R₀处相对微分注量。

根据作用点R₁(x₁，y₁，z₁)的坐标和放射源(x_s，y_s，z_s)，求得放射源(x_s，y_s，z_s)到R₁处的距离r₁。

R₁处的相对微分注量为：

${\mathrm{\ψ}}_E\left(R_1\right)={\mathrm{\ψ}}_E\left(R_0\right)\exp (-{\∫}_{r_0}^{r_1}\frac{\mathrm{\μ}}{\mathrm{\ρ}}(E,l)\mathrm{\ρ}\left(l\right)\mathrm{dl})---\left(2\right)$

注意：是沿作用点r₁与r₀间原射线径迹的积分，可离散为原射线径迹上的栅元的求和。

上述操作已得到微分比释总能的信息，以下操作将获得点核的信息：

根据作用点R₁(x₁，y₁，z₁)的坐标与沉积点R₂(x₂，y₂，z₂)的坐标，可得到两点间的物理长度L，以及两点连线关于作用点r与r₀间原射线径迹方向的夹角θ：

$L=\sqrt{{(x_2-x_1)}^2+{(y_2-y_1)}^2+{(z_2-z_1)}^2}---\left(3\right)$

$\mathrm{\θ}=\mathrm{\π}-\arccos \frac{{r_{01}}^2+L^2-{r_{02}}^2}{{2r}_{01}L}---\left(4\right)$

其中，r₀₁和r₀₂分别R₀到R₁的距离和R₀到R₂的距离，可用其坐标求得：

$r_{01}=\sqrt{{(x_0-x_1)}^2+{(y_0-y_1)}^2+{(z_0-z_1)}^2}---\left(5\right)$

$r_{02}=\sqrt{{(x_2-x_0)}^2+{(y_2-y_0)}^2+{(z_2-z_0)}^2}---\left(6\right)$

考虑到作用点R₁(x₁，y₁，z₁)与沉积点R₂(x₂，y₂，z₂)之间的非均匀组织修正，需要得到用以作用点R₁为圆心，在垂直于作用点到沉积点连线方向，以带电粒子在介质中的射程为半径的圆面为底面，以作用栅元到沉积栅元连线为高，所得的圆柱体内的栅元电子密度的平均值代替连线上各栅元电子密度的平均值c_(r，r′)(如图2所示)，调用密物理度为ρ₀的均匀介质的点核矩阵h(E，c_(r，r′)(r-r′))，根据(θ，c·L)选出对应的点核具体值，若角度或半径不是矩阵对应的数值，可使用线性插值的方法，得到任意角度、半径的点核值。将先前得到的微分比释总能乘以对应的点核具体值，就得到了能量为E的光子通过作用点R₁(x₁，y₁，z₁)对沉积点R₂(x₂，y₂，z₂)的剂量贡献。将所有作用点的剂量贡献相加就得到沉积点R₂(x₂，y₂，z₂)的剂量。

4.MDSA的输出结果，包含以下内容：

全空间剂量分布的文本文件。

5.计算实例，包含以下内容：

需要计算能量为6MV的光子，通过作用点(0，0，20)对沉积点(0，20，20)的剂量贡献(如图3所示)：

根据作用点R₁的坐标(0，0，20)，可计算出到达作用点R₁处的原射线与介质表面的交点R₀坐标(0，0，0)，进而求得放射源(0，0，-100)到R₀处的距离r₀＝100cm，根据(0，0，0)的相对微分注量及距离平方反比公式得到R₀处相对微分注量ψ_E(R₀)＝1。

根据作用点R₁(x₁，y₁，z₁)的坐标和放射源(0，0，-100)，求得放射源(0，0，-100)到R₁处的距离r₁＝120cm。

R₁处的相对微分注量为： ${\mathrm{\ψ}}_E\left(R_1\right)={\mathrm{\ψ}}_E\left(R_0\right)\exp (-{\∫}_{r_0}^{r_1}\frac{\mathrm{\μ}}{\mathrm{\ρ}}(E,l)\mathrm{\ρ}\left(l\right)\mathrm{dl})---\left(7\right)$

注意：是沿作用点r₁与r₀间原射线径迹的积分，可离散为原射线径迹上的栅元的求和。

上述操作已得到微分比释总能的信息，以下操作将获得点核的信息：

根据作用点R₁的坐标(0，0，20)与沉积点R₂的坐标(0，20，20)，可得到两点间的物理长度L＝20cm，以及两点连线关于作用点r与r₀间原射线径迹方向的夹角θ＝90°：

考虑到作用点与沉积点之间的非均匀组织修正，需要得到用以作用点R₁为圆心，在垂直于作用点到沉积点连线方向，以带电粒子在介质中的射程为半径的圆面为底，以作用栅元到沉积栅元连线为高，所得的圆柱体内的栅元电子密度的平均值代替连线上各栅元电子密度的平均值c_(r，r′)，调用物理密度为1.0kg/cm³的均匀介质的点核矩阵h(E，c_(r，r′)(r-r′))，根据(90°，c·20)选出对应的点核具体值，若角度或半径不是矩阵对应的数值，可使用线性插值的方法，得到任意角度、半径的点核值。将先前得到的微分比释总能乘以对应的点核具体值，就得到了能量为6MV的光子通过作用点对沉积点的剂量贡献。将所有作用点的剂量贡献相加就得到沉积点的剂量。

Claims (1)

1.一种解析蒙特卡罗剂量计算方法，其特征在于包括以下步骤：

(1)获得计算参数，包括以下内容：

蒙特卡罗程序获得的能量沉积核h：h是一系列单能或一定能谱的光子入射到某一均匀介质中，对应于球坐标(θ，R)的辐射能量值矩阵；

一系列单能或一定能谱的光子在某一介质中的质量衰减系数μ/ρ；

一系列单能或一定能谱的电子在某一介质中的射程；

计算模型栅元划分信息；

计算模型栅元物理密度信息；

计算模型栅元电子密度信息；

计算模型中计算区域边界标示信息；

计算模型中感兴趣区域标示文件；

获得放射源的能谱信息；

射线在介质表面形成的形状，由用户给出；

(2)利用下列MDSA模型进行计算：

使用MDSA模型计算公式，如下：

$D\left(r\right)={\∫}_E\∫\∫{\∫}_V{T}_E\left(r^{\′}\right)\frac{\mathrm{\ρ}\left(r^{\′}\right)}{{\mathrm{\ρ}}_0}{c}^{2}h(E,c_{(r,r^{\′})}(r-r^{\′}))d^3{r}^{\′}\mathrm{dE}---\left(1\right)$

其中，T_E(r′)为微分比释总能；ρ(r′)是作用栅元的物理密度；h为预先计算好的物理密度为ρ₀的均匀介质中的能量沉积核；c_(r，r′)代表了三维非均匀修正系数，是以某个作用栅元为圆心，在垂直于作用栅元到沉积栅元连线方向，以带电粒子在介质中以射程为半径的圆面为底，以作用栅元到沉积栅元连线为高，所得的圆柱体内的栅元电子密度的平均值：

$c=c_{(r,r^{\′})}={\∫}_0^1{\mathrm{\ρ}}_{\mathrm{rel}}[r^{\′}-l(r^{\′}-r)]\mathrm{dl}---\left(2\right)$

T_E(r′)使用子束入射到介质表面的能量注量计算所有深度的微分比释总能，即根据各个小子束在介质表面的划分，随着深度的增加，每个小野的尺寸按比例扩大，保持小野数目不变，相对于传统DSA方法，减小各深度计算的采样点数目，达到减小计算时间的目的；

在使用MDSA模型计算时，根据不同的计算区域，分别采用不同的计算点密度：感兴趣区域的计算点密度高，其他区域计算点密度低。

Images