放疗绝对剂量的数学算法模型*

谢天赐 张彬 贺泊 李昊鹏 秦壮 钱金钱石锲铭 Lewis Elfed 孙伟民†

1) (哈尔滨工程大学物理与光电工程学院,中国教育部纤维集成光学重点实验室,哈尔滨 150001)

2) (黑龙江大学电子工程学院,哈尔滨 150001)

3) (利莫瑞克大学,光纤传感器研究中心,利莫瑞克,爱尔兰)

本文提出一种通过物理模型计算放疗过程中每一个组织深度处绝对剂量的算法,它可代替蒙特卡罗仿真的部分工作且耗费时间更少.这个算法是基于对照射野内X射线产生电子的能量注量的积分运算,并考虑了射线的能谱及二次散射线,得到了后向散射对表面剂量的贡献比例,同时得到前向散射、后向散射及原射线剂量贡献的关系.比较了二次光子和二次电子的三维能谱,得出该能谱是粒子注量关于粒子能量和粒子运动方向的函数.为了得到每一深度处的光子注量,计算了有连续能谱的X射线的期望质量衰减系数.上述算法计算得到的绝对剂量与蒙特卡罗方式仿真的结果趋势一致,两者的差异在于算法未考虑高于二次的散射线.最后将算法应用到非均匀模体剂量计算,能准确反映其中剂量分布特点且具有较小的误差.

1 引 言

在放射治疗领域,蒙特卡罗 (Monte Carlo)仿真是目前公认的最精确的剂量计算方法,它可以精确地模拟加速器和患者组织的各种结构的真实辐射传输[1].蒙特卡罗仿真基于粒子运输计算和概率统计的原理,可以在计算机上进行癌症风险评估[2]、合理制定放疗计划[3]、直线加速器调试[4]等工作,但要想得到较高的仿真精度会耗费大量的时间.

为了更直观地反映辐射和组织之间的相互作用,减少计算时间,研究人员尝试提出一些理论模型和计算方法来代替蒙特卡罗仿真[5-7].Jette[8]计算了不同能量光子束从一个相互作用点产生的所有康普顿电子的剂量分布,但只计算了单能射线,且没有考虑射野中心轴以外其他地方对剂量分布的影响.Tillikainen等[9]给出了窄射线束条件下各深度辐射剂量的计算方法并考虑能谱和散射等相关因素,但其计算仍依赖于蒙特卡罗仿真的射线束.同时对于笔形束卷积法、各向异性解析法、筒串卷积叠加法在计算非均匀模体时存在较大误差[10-13],特别是对于笔形束卷积法,由于并未考虑电子的横向平衡条件,在低密度的肺模体中计算剂量明显偏高[12].

本文以康普顿电子的特性及其产生的能量沉积为基础来解决放疗中的剂量计算问题,主要工作是计算原射线和散射线在各个深度上产生的电子能量注量,并通过积分能量注量得到射野中心轴线上的绝对剂量与水模体深度的关系.这些研究都是基于康普顿效应,康普顿效应在高能X射线(如加速电压为6 MV的X射线)与物质相互作用中起着重要作用.本文以医用加速器的非单能射线为射线源,考虑了射线的二次散射影响,将其分为前向散射和后向散射,计算出前向、后向散射对表面剂量的贡献比例.比较了二次散射光子和二次反冲电子的三维能谱,该能谱表示粒子通量与能量和粒子方向的关系.通过对非均匀模体的剂量计算发现,本研究的算法能真实反映模体内轴向剂量和横向剂量分布.

2 理论模型与分析方法

2.1 绝对剂量表达式

在医用X射线照射水模体的物理模型中(如图1所示),为了计算射野中心轴(z轴)计算点P的剂量,建立了以水模体表面的射野中心点为坐标原点的三维坐标系,其中z轴的方向向下,2l为射野的边长,hv为入射光子的能量,O'为射野平面C的中心点,微元K的长、宽、高分别为dx,dy,Δz,P点所在深度为z.为了建立剂量形成的物理模型,水模体被分成了zmax/Δz个厚度为 Δz的薄层,其中zmax为水模体最大深度,Δz为薄层水模体厚度.

图1 绝对剂量的算法模型示意图Fig.1.Schematic diagram of absolute dose algorithm model.

射野中心轴线上P点的剂量D可以通过累加计算得到,该累加来自于对每一个深度处薄层内的射线与物质作用产生的电子到达P点时沉积的能量.在考虑射线在组织中的衰减后,P点处的剂量表示为

其中D′(m,z) 表示在深度为mΔz处的薄层C产生电子到达P点所造成的剂量,m为正整数.

根据吸收剂量的定义[14],某一点的吸收剂量在数值上等于能量注量.水模体某一深度处D′可以近似由原射线和二次散射线产生的能量注量表示：

2.2 原射线产生的剂量

在K点产生的电子的经过距离(图1中的红色虚线)后到达计算点P,φ1为电子运动方向与z轴的夹角,吸收剂量是薄层C内所有微元对P点的电子能量注量贡献的积分,由下列函数表示：

其中Ψ(m,φ1) 表示电子的能量注量,k(rkp) 表示电子的能量注量随距离的衰减规律,ne表示物质的每克电子数,单位为 g—1,对于水模体取值为 3.343 ×1023g—1.Ψ(m,φ1) 的表达式[15]为

其中 e-τmΔz表示非单能射线光子注量的衰减规律,τ表示期望质量衰减系数,Φ(φ1,hv) 表示水模体表面z= 0处光子与物质相互作用产生的电子注量,e-τmΔz与Φ(φ1,hv) 的乘积表示深度为mΔz的薄层C内的光子与物质相互作用产生的电子注量.E(φ1,hv)表示电子的能量,积分的上下限为电子的最大能量和最小能量.φ1的表达式为

其中x和y是微元K的坐标.Φ(φ1,hv) 表达式为

其中ρ(φ1,hv) 表示二次电子的微分截面(光子与物质相互作用的概率相对方向角的微分).n(hv)表示入射光子的注量,表达式为

其中A1,A2,t1和t2为一般参数,没有物理含义.电子的微分截面ρ(φ1,hv) 是一个重要的参数,决定了电子在不同方向的运动特性.ρ(φ1,hv) 的表达式为[16]

其中re表示电子的经典半径,β=hv′/hv[16]表示二次光子的能量hv′与入射光子能量hv的比值,θ1表示二次光子运动方向与入射光子方向的夹角.β是关于θ1的函数：

而θ1为关于φ1的函数：

其中α的表达式为

其中mec2为单个电子的静止能量.

E(φ1,hv)是基于康普顿效应的能量守恒原则的函数,其表达的电子能量可以通过入射光子能量与散射光子能量做差得到,函数表达式为

其中hv′表示二次光子的能量.

在(3)式中的k(rkp) 可通过研究窄电子束条件下能量注量在单一方向的变化规律得到.从电子束的百分深度剂量 (percent depth dose,PDD)曲线随射野的变化规律可以得出,当射野足够小(窄电子束)时,电子的能量注量将遵循下列指数衰减规律[17]：

根据图1中的几何关系,其中的rkp的表达式为

而µ2表示单能电子的衰减系数,可以通过射线产生的二次电子的最大穿透深度推导,该深度在数值上与最大剂量点对应的深度相等[17],电子的能量随着反冲角φ1增大而减小(见第3.3节分析),因此µ2是关于反冲角φ1的函数.不妨设两者关系为线性：

其中µ0表示反冲角φ1为0时电子的衰减系数.根据不同能量的电子对应不同射程,这里µ0取值4,kρ取值3.

2.3 二次射线产生的剂量

图2 (a)和 (b)算法模型示意图Fig.2.Schematic diagram of (a) and (b) algorithm model.

其中hv′′表示三级光子的能量,θ2表示二次光子散射角,φ2表示二次电子反冲角.

其中二次散射角θ2是一个关于φ1和θ1的函数：

二次散射角θ2是关于φ1和θ1的函数：

来自所有方向的二次光子所产生的电子到达点P,所有电子在此过程中的运动轨迹形成“圆锥”(图2红色线和黑色虚线组成).由该“圆锥”造成的能量注量需要针对其体积作积分运算,但是根据(16)式的运算,需要对该能量注量计算“圆锥”体积积分的均值以得到前向散射的式子

(24)式与(25)式计算的是二维积分的均值,但在数值上与“圆锥”积分的均值相等.

2.4 射线的期望线性衰减系数τ

任意单能的线性衰减系数µ为[16]

其中ρ表示物质的密度,单位为 g/cm-3,ne表示物质的每克电子数,σe表示康普顿效应截面.σe是一个关于光子能量的表达式[16].

因此线性衰减系数µ是一个关于光子能量hv的函数µ(hv) ,在深度为z处的非单能射线强度表示为

其中hvmax表示能谱函数n(hv) 里射线能量的最大值.

将(7)式代入(28)式,得出光子强度关于深度的函数关系,发现其仍然遵循指数衰减规律.期望质量衰减系数τ用(29)式拟合I(z) 的函数曲线得到

其中y0表示组织深度为0时的光子强度.这个指数衰减表示了原射线的光强的变化规律.

3 结果与分析

3.1 蒙特卡罗对绝对剂量的仿真与计算

图3所示为真实加速器(型号为Varian IX 3937)测量 PDD实验数据 (experiment data)与基于蒙特卡罗仿真出来PDD数据(Monte Carlo simulation)的对比,纵坐标为以最大剂量值为单位的归一化剂量(normalized dose),横坐标为水模体深度(depth).实验条件为6 MV的X射线、照射野面积为 10 cm × 10 cm,源皮距为 100 cm.在每个深度处两组数据的差异小于 ± 3.6%.

图4为上述仿真的加速器结构下6 MV 的X射线能谱 (6 MV X-ray spectrum)及来自 (7)式的仿真函数n(hv) (fitting functionn(hv) ),纵坐标为光子通量 (photon flux),横坐标为粒子能量 (energyhv).两者的相关系数R2大于0.993,能谱数据采集于射野中心轴距离射线源100 cm处.

图3 蒙特卡罗 PDD 与真实测量数据的对比Fig.3.Comparison of PDD between Monte Carlo simulation and experimental data.

图4 6 MV 射线能谱及其仿真函数Fig.4.6 MV X-ray spectrum and its fitting function.

将该蒙特卡罗仿真出的加速器的射线能谱、期望质量衰减系数τ= 0.701 以及µ0=4 和kρ=3 分别代入(6)式、(4)式与(15)式,计算得到原射线和二次射线产生的绝对剂量,并将其与蒙特卡罗仿真剂量做对比来验证算法的正确性,图5为经过上述算法得到的绝对剂量D(z) (calculation ofD(z))与仿真的绝对剂量 (simulated by Monte Carlo,MC)的对比,其中算法耗时约为 5 min,蒙特卡罗软件Egsnrc耗时约为30 min,本方法的计算效率明显高于蒙特卡罗仿真方法.D(z) 的计算没有考虑以下因素： 1)原射线入射前的各种散射X光子和污染电子； 2)由原射线在水模体中产生的多次(大于二次)散射X光子.这两个因素是仿真与计算数据在剂量建成区和指数衰减区存在差异的原因,两者的差异曲线如图5中黑色虚线所示.由于原射线中含有的散射X光子和电子的能量低、数目少,因素1造成的剂量差异在建成区中快速衰减[18].随着深度z增大,射线被水模体的散射程度增加,由此造成的剂量差异在指数衰减区缓慢增加,此外,因素2造成的剂量差异还与散射校正因子(Scp)[19]相关.对差异曲线利用补偿函数fd进行拟合,其表达式如下

拟合线性度(R2)好于0.95,通过该函数补偿保证了上述算法计算的剂量与仿真的一致性.

图5 蒙特卡罗仿真的绝对剂量和计算绝对剂量 D (z) 的比较Fig.5.Comparison of absolute dose between Monte Carlo simulation and calculation of D (z) .

3.2 原射线与二次射线

在上述实验条件下计算射野薄层m= 1的能量注量,得到该射野层原射线剂量贡献(primary ray)、前向散射剂量贡献(forward ray)和后向散射剂量贡献(backward ray)三者随深度z的关系,如图6所示.是由原射线与水模体表面作用产生的电子造成的,可类比成以深度为0处的零射野为电子源的电子束造成的吸收剂量,所以不存在表面剂量,剂量建成区从剂量为0开始.也可以类比成电子束造成的剂量,它与的不同在于电子是由来自任意散射角的二次光子产生的,这意味着一定数量的电子可以到达深度为0处的位置,因此其存在表面剂量.代表的电子运动方向与和的相反,为了方便比较,将的图像也放在z轴的位置.

图6 比较由原射线产生的剂量 、前向散射产生的剂量 、后向散射产生的剂量 随深度的变化Fig.6.Comparison of dose caused by primary ray ,forward scatter and backscatter with depth.

图7 向前和向后的二次射线产生的剂量随深度的变化Fig.7.Dose caused by scattered ray in forward and back directions with depth.

3.3 二次电子与二次光子的能谱

利用(9)式—(12)式之间的转换可得到hv关于φ1和E的函数关系hv(φ1,E),并将其代入(6)式中,得到图4所示的6 MV X射线产生的二次电子的能谱Φ(φ1,E) ,如图8所示.

图8中横坐标为电子反冲角 (recoil angle,Ra)φ1,其范围为 0 ＜φ1＜ π/2,纵坐标为电子的能量E,竖坐标为粒子注量 (particle flux,Pf),图中的二维图为三维能谱图的三视图.当反冲角接近0时,电子的数目达到最大值,此时,电子的能量为0.44 MeV.从图8中还可以看出二次电子分布在φ1与E的数量乘积小于1的范围内.

二次射线的能谱不同于原射线的能谱,且其在不同散射角 (scattered angle)θ1上也不尽相同.根据(7)式以及二次光子的微分截面,可以得到二次光子关于每个散射角下的能谱,其表达式如下

图8 二次电子的粒子注量随能量E和反冲角 φ1 的变化Fig.8.Particle flux of secondary electron with energy E and recoil angle φ1 .

图9 二次光子的粒子注量随能量 h v′ 和散射角 θ1 的变化Fig.9.Particle flux of secondary photon with energy h v′ and scattered angle θ1 .

根据 (9)式,(32)式中的hv是关于hv′的函数.

图9为图4所示的6 MV X射线产生的二次光子的三维能谱横坐标为光子的散射角 (scattered angle,Sa)θ1,其范围为 0 ＜θ1＜ π,图中二维图为三维能谱的三视图.当散射角为0时,光子的数目达到最大值,此时光子的能量为0.46 MeV.二次光子分布在θ1与hv′的乘积小于1.1的范围内.

对比两种粒子的能谱可以发现,当电子的反冲角φ1= 0 时,光子的散射角θ1= π,此时,入射的光子与电子发生对心碰撞,光子的能量范围在0—0.25 MeV,光子注量峰值约为0.55 × 10—6cm2·MeV—1,电子的能量范围在0—5.75 MeV,电子注量峰值约为 6 × 10—6cm2·MeV—1.当电子的反冲角φ1=π/2 时,光子的散射角θ1= 0,此时,入射的光子从电子旁边掠过,光子没有能量损失,电子注量为0.当散射角θ1= π/2 时,光子三维能谱图中“山脊”中有一个最低点,其代表的粒子注量为 0.5 × 10—6cm2·MeV—1,该点对应的光子能量为 0.17 MeV,而此时的电子反冲角φ1= 0.75,对应电子三维能谱图“山脊”曲线斜率绝对值的最小值.

3.4 PDD的计算与仿真

图10为通过蒙特卡罗仿真得到的不同射野PDD与算法计算得到PDDD(z)%(以D(z) 最大剂量为单位)的对比,其中源皮距为 100 cm,射线源为图4 所示的 6 MV X 射线.在射野 3 cm × 3 cm范围内,计算与仿真PDD最大误差剂量约为2%(表面剂量例外).两者在剂量建成区及指数衰减区的差异原因在第3.1节中已指出.在大多数文献中[20-22],都认为表面剂量来自原射线中的散射线及污染电子,部分来自于组织中的后向散射[21].而后向散射贡献的这一部分可由(25)式计算得到,为相对最大剂量的5%,如图10中D(z)%与y轴的交点.在指数衰减区,射野 3 cm × 3 cm 下的蒙特卡罗仿真数据与D(z) 的计算数据的差异和4 cm ×4 cm 射野与 3 cm × 3 cm 射野蒙特卡罗仿真数据差异相似,后者主要是因为当射野面积增加时,散射体积的增加导致了更多的散射线,这从侧面验证了计算与仿真的差异在于仿真的剂量会有更多来自散射线的贡献.不难得出,当照射野接近0 × 0 时,D(z)%的计算值将和真实的 PDD 非常接近.

图10 比较随不同射野 (射野大小 3 cm × 3 cm 和 4 cm ×4 cm)变化的蒙特卡罗仿真数据 D MC 与计算的 D (z)% (射野大小 3 cm × 3 cm)Fig.10.Comparison of Monte Carlo simulation PDD with different fields (field size 3 cm × 3 cm and 4 cm × 4 cm)and calculation PDD D (z)% (field size 3 cm × 3 cm).

3.5 非均匀模体剂量计算

非均匀模体由三层物质组成,分别为水、肺、水,厚度均为 5 cm.利用图4 所示的 6 MV 射线,射野为 3 cm × 3 cm,模体中的肺密度ρ设置为0.1,0.2,0.4 和 1 g/cm3.不同密度物质的剂量D(z)的计算差异在于光子注量的衰减和电子能量注量的衰减的不同,对应参数分别为τ和µ2,其中τ可以根据2.4节计算.µ2的取值过程则取决于横向电子平衡.

如果射野的尺寸不能满足横向电子平衡,电子会将部分能量带到射野外,同时物质的密度变化会改变电子的射程,也会导致横向电子不平衡.6 MV射线在射野为3 cm × 3 cm刚好达到横向电子平衡时的肺密度为ρT,当ρ≥ρT时,肺模体内达到横向电子平衡,(15)式 (µ2=kρ ·φ1+µ0)中kρ与µ0仍保持不变,此时肺密度的改变只改变光子强度的衰减程度.当肺模体的密度ρ＜ρT时,不满足横向电子平衡,在φ1接近 π/2 时,肺模体密度越低,从射野边缘逃逸的电子射程越大,此时µ2越小； 在φ1接近0时,电子的能量注量都在模体内沉积,当模体密度变小时,能量注量的沉积变少,此时的µ2越大.不同肺模体密度ρ对应的参数kρ与µ0关系如表1 所列,kρ与µ0均关于ρ呈等差变化,公差数与射野大小和射线能量相关.

表1 不同肺模体密度 ρ 对应参数 kρ 和µ0Table 1.Parameters kρ and µ0 of lung phantom densities ρ.

图11所示为4种不同肺密度的PDD曲线,对比肺模体密度为0.1,0.2与1 g/cm3时的结果发现,横向电子不平衡带来的剂量跌落随着密度的降低而增大.肺模体密度0.1 g/cm3的PDD在第一个水肺交界处(如图11中深度为5 cm处曲线所示)相对深度为10 cm的剂量跌幅约为17%,密度0.2 g/cm3的剂量跌幅约为 5%,肺模体密度 0.4 g/cm3的PDD曲线在10 cm后剂量增长约为4.5%,这一结果与文献[23]一致.但在深度6—10 cm处的计算剂量偏低,与文献[23]差别最大为3%,这主要是因为算法未考虑到从第一层水模体中进入肺模体的多次散射带来的剂量贡献,而更小射野散射带来的剂量贡献会更小,因此图11的计算结果与更小射野的非均匀模体PDD更加接近[24].

图11 不同肺密度的水肺水模体的射野中心轴百分深度剂量Fig.11.Percentage depth dose of the central axis of the field of radiation in water-lung-water phantom with different lung densities.

在第二个肺水的交界面(如图11中深度为10 cm处曲线所示),低密度肺模体PDD相对于1 g/cm3肺模体PDD出现明显的剂量增长,这主要是因为光子强度在低密度肺模体中衰减较少,在达到第二个肺水交界处后光子强度相对较高.此处与射线从空气入射到水模体中类似,会在水下出现第二个最大剂量点,而这个最大剂量点相对深度低于第一层水模体的最大剂量点深度1.5 cm,这是由于6 MV射线在第二个最大剂量点处由于射线的衰减能谱发生变化,低能射线的成分增加,次级电子的射程变短.

图12所示为在深度为7和11 cm处采集的离轴比曲线,所有离轴比剂量均以肺密度为1 g/cm3的水肺水模体在离轴0 cm处剂量为单位.在深度7 cm 处,密度为 1 g/cm3肺模体对应的离轴比曲线在射野边缘剂量跌落幅度最大,跌落幅度随着密度的减小而减小,与高密度的离轴剂量相比,较低密度的离轴剂量在射野边缘会更快进入剂量跌落状态,因为较低密度模体的射野内电子会因为横向电子不平衡将部分能量带离射野.同时还可以发现,在射野外,低密度的肺模体剂量明显高于1 g/cm3肺模体对应的剂量,这是由于低密度肺模体的电子能量注量流失到射野外,流失的这部分能量注量随着密度的降低而增加.在深度为11 cm水模体处,来自较低密度肺模体的剂量因为光子强度衰减较少明显高于较高密度的剂量,此时电子横向平衡被恢复,来自各个密度肺模体的射野外剂量变化规律保持一致.

图12 不同模体深度处的离轴比曲线 (a) 7 cm； (b) 11 cmFig.12.Off-axis ratio curves at different phantom depths：(a) 7 cm； (b) 11 cm.

4 结 论

本文通过构建放疗过程中剂量形成的物理模型来得到吸收剂量的数学算法.加入补偿函数后,该算法计算得到的吸收剂量与蒙特卡罗仿真的结果相符,利用该算法计算了在一定射野范围内不同组织深度处的吸收剂量,在计算中考虑了非单能射线的能谱及二次散射线.其中补偿函数在剂量建成区呈现快速衰减,在指数衰减区呈缓慢上升趋势,前者是一个与原射线中的散射线和电子相关的变量,后者是一个与康普顿效应产生多次(大于二次)散射射线相关的变量.因此该算法能保证与仿真结果的一致性,且相比蒙特卡罗仿真耗时更短.

讨论了每一深度处由原射线及二次散射线对射野中心轴产生的剂量贡献,并分别讨论了前向散射和后向散射.其中后向散射的剂量贡献呈指数衰减规律,前向散射的剂量贡献是一个先上升后下降的函数,发现模体表面剂量部分来自组织中的后向散射,并计算出这部分所占比例.发现了二次散射线的剂量贡献在偏离作用点处存在峰值.在此研究中,得到了医用6 MV X射线产生的二次光子和二次电子的三维能谱,该能谱是粒子注量与粒子能量和方向相关的函数.对比两种粒子的能谱,发现两者的粒子注量的分布在较小方向角和较小的能量范围内更为集中,并分析了三种散射角或反冲角下光子与电子的特征差异.

对非均匀模体的剂量计算中发现,电子在低密度肺模体有更大的射程,因此电子的能量注量衰减系数相比密度较高的更小,本研究在计算中考虑到了该系数随着密度和电子运动方向角变化,同时调整不同密度模体中光子强度衰减系数,准确计算了由于不同密度肺模体带来横向电子不平衡的PDD曲线和离轴比曲线.通过对比MC仿真PDD数据和对比他人研究非均匀模体剂量分布发现,本算法在小射野PDD的计算上有更小的误差.