Shielding Data for 400 MeV/u ¹²C Ions on Iron and Water Targets

重离子在物质中的能量损失集中在射程末端的布拉格峰处，且在布拉格峰区具有相对高的生物学效应。与传统的放疗方法相比，重离子治癌在物理特性及生物特性上具有独到的优势^[1-2]，该技术被认为是21世纪最理想的放疗技术之一。在我国，中国科学院近代物理研究所已在甘肃省兰州市和武威市建造了两台医用重离子加速器(HIMM)，其中武威HIMM装置已于2020年4月1日正式投入临床使用。

目前，世界范围内建成并投入运营的医用重离子加速器均为碳离子加速器。为满足不同深度的临床治疗需求，典型的¹²C能量为400 MeV/u，该能量的碳离子在生物体中的射程接近30 cm。

粒子加速器在束流注入、引出、传输及治疗等过程中均会出现束流损失，从而产生次级辐射。因此加速器在建造之前，必须进行详细的屏蔽设计来保护环境、公众及工作人员等的辐射安全。保守考虑，在实际的辐射屏蔽设计中，往往选用粒子的最大能量作为屏蔽设计的依据^[3-5]。

加速器装置通用的辐射屏蔽计算方法有两种，一种是蒙特卡罗模拟法，即采用FLUKA^[6]、MCNP^[7]等蒙特卡罗程序进行模拟计算；另一种是传统的经验公式法，如Moyer^[8]模型、Tesch公式^[9]等。早期，经验公式法是加速器装置屏蔽设计的最主要方式；21世纪以来，随着计算机技术的飞速发展，蒙特卡罗模拟法逐渐成为加速器辐射屏蔽设计的主要方法。蒙特卡罗模拟法可以对复杂模型进行直接的分析计算，但在计算时需要设计人员手动输入和设置几何结构、物理列表、反应模型等参数，这其中任一环节出现失误，均可能使计算结果产生较大的偏差。因此在实际的屏蔽设计中，往往需要通过其他方法对蒙特卡罗程序的模拟结果进行验证，而传统的经验公式法多针对GeV以上能区的质子加速器辐射屏蔽估算，用于重离子加速器屏蔽计算较为不便。

为了给医用重离子加速器提供一种专用的快速计算屏蔽厚度或对蒙特卡罗模拟结果进行验证的方法，本文采用FLUKA程序计算了400 MeV/u碳离子打厚靶在不同屏蔽材料中的屏蔽参数。

FLUKA是一款广泛应用于粒子加速器辐射源项和屏蔽计算的蒙特卡罗程序包，它包含了详细的强子和原子核相互作用模型、强子和电磁过程之间的耦合等物理作用原理，可以用来计算带电粒子输运和带电粒子与物质的相互作用^[6]。

屏蔽参数计算采用的模型如图1所示，400 MeV/u碳离子束流打圆柱形靶体，以打靶位置为球心在半径30至36 m处设置球壳状屏蔽体。靶体与屏蔽体内表面之间及屏蔽体外的空间均设置为真空。整个屏蔽体被划分为多层来减小模拟计算的统计方差^[10-11]，同时也便于采用USRBDX选项卡统计不同屏蔽层的粒子信息。如图1中红色部分所示，取自球心向外延伸的圆锥与屏蔽体球壳相交的范围作为统计区域来统计发射角0°~10°、20°~30°、40°~50°、60°~70°、80°~90°、120°~130°、170°~180°范围内的粒子信息。

在使用FLUKA程序计算时，调用DPMJET、IONTRANS和LOW-NEUT选项卡模拟重离子和低能中子的输运，调用PHYSICS选项卡对重离子核反应中的蒸发碎片、电磁离解等过程进行设置，中子及其他粒子的输运阈值设置为10^–3 eV。设置了AMB74选项卡来计算周围剂量当量，它采用的是ICRP74号出版物和Pelliccioni等^[12-13]的报告中关于高能中子部分的通量－剂量转换系数。

靶的材料选择铁和水，靶厚度略大于400 MeV/u碳离子在该靶材中的射程，以确保打靶可以最大限度地产生次级辐射。靶体的材料、尺寸以及利用SRIM软件^[14]算得的碳离子射程如表1所列。

混凝土造价便宜、屏蔽适应性广、与土建契合性高，是加速器屏蔽设计中最常用的材料。铁和铅等重材料也经常被用来与混凝土结合使用，以实现局部屏蔽的功能。本次计算选用了密度2.3 g/cm³的混凝土、7.874 g/cm³的铁及11.35 g/cm³的铅作为屏蔽材料。混凝土材料成分如表2所列。

图2是400 MeV/u碳离子打30 cm厚水靶产生的次级辐射剂量分布图，剂量单位是Sv·ion^–1。从图中可以看出，400 MeV/u碳离子打靶产生的次级辐射场具有很强的前冲性，与束流方向夹角越大，其周围剂量当量水平越低。

医用重离子加速器屏蔽外的周围剂量当量主要由次级中子贡献，其他次级粒子的贡献几乎可以忽略。从辐射源项的角度考虑，$\gamma $射线在未经任何屏蔽时带来的周围剂量当量比中子低两个量级以上；在屏蔽体中$\gamma $射线带来的周围剂量当量也要比中子低一个量级^[3-5]。因此本文对打靶产生的中子和质子周围剂量当量进行了统计，图3是部分发射角下中子和质子周围剂量当量占总周围剂量当量的比例随屏蔽厚度变化的曲线，图中横轴的起始点表示距离中心靶点30 m远的屏蔽球壳内表面处。从图中可以看出，400 MeV/u碳离子打厚铁靶和水靶产生的次级辐射未经屏蔽时，周围剂量当量几乎全部由次级中子与次级质子贡献；经一定厚度屏蔽后，周围剂量当量几乎全部由次级中子贡献。当混凝土屏蔽体的厚度大于50 cm时，次级中子占总周围剂量当量的比例基本达到90%以上；当铁屏蔽体厚度大于20 cm时，次级中子占总周围剂量当量的比例约为95%；当铅屏蔽体厚度大于10 cm时，次级中子占总周围剂量当量的比例约为99%。因此在医用重离子加速器的屏蔽设计中，要重点关注次级中子的屏蔽。

加速器打靶产生次级辐射场具有以下两个特征：是一个类点源场，剂量水平与关注点距打靶位置距离r的平方成反比；剂量随屏蔽厚度的增加呈指数衰减。在束流集中损失的屏蔽计算中，通常使用双参数的指数式(1)来描述周围剂量当量在屏蔽体中的衰减^[15]。

式中：H是每个粒子打靶产生的周围剂量当量，单位是Sv·ion^–1；E_p是粒子能量；r是粒子打靶位置与统计位置之间的距离，单位是m；θ是束流方向与统计位置的夹角；H₀是辐射源项，单位是Sv·m²·ion^–1；d是屏蔽层的厚度，单位是cm；ρ是屏蔽体材料密度，单位是g/cm³；λ(θ)是周围剂量当量在屏蔽材料中的衰减长度，单位是g/cm²；α是r方向与屏蔽面法线之间的夹角，g(α)取值为cosα，本次计算使用的模型屏蔽体为球状，因此cosα恒为1。

下文的剂量衰减曲线中，将以H·r²作为纵轴，其单位是Sv·m²·ion^–1，用来表征单一粒子带来的辐射剂量大小。以H·r²作为纵轴一方面便于直接观察一定角度范围内周围剂量当量随屏蔽厚度增加的衰减趋势，而不体现与距离r的平方成反比的规律；另一方面也便于对数据进行拟合处理从而得到H₀与λ的值。

图4是400 MeV/u碳离子打靶产生的周围剂量当量在屏蔽体中的衰减曲线。从图中可以看出，前向的衰减曲线在刚到达屏蔽体表面时会出现爬高或衰减缓慢的现象。这是由于前向产生的次级粒子中，中高能成分占绝对优势，刚到达屏蔽体表面时会迅速损耗能量产生大量低能中子，从而出现明显的二次中子增殖效应。而侧向及后向的衰减曲线刚开始衰减较快，随着屏蔽厚度增加衰减速度逐渐减慢。这是由于侧向和后向产生的次级粒子中，低能中子占绝对优势，随着屏蔽厚度的增加，较低能量中子所占的比例越来越少，能谱逐渐硬化，最终达到能谱平衡状态。

从图4可以看出，三种材料的屏蔽效果具有很大差异，这是主要是由其中子反应截面的差异导致的。中高能中子在进入屏蔽材料后，主要通过非弹性散射反应迅速降低能量；当能量低于非弹性散射阈能时，主要通过弹性散射来降低能量；当能量足够低时，主要通过俘获反应被吸收^[16]。图5给出了混凝土、铁和铅的中子反应总截面、非弹反应截面、弹性散射截面的对比图。图中纵轴为σ·ρ·N_A·M ^–1，单位是cm^–1，用来表征单位厚度的屏蔽材料与中子发生反应的概率大小。其中σ是单一原子核的中子反应截面，单位cm²，作图所用的俘获截面数据均取自ENDF数据库^[17]；ρ是屏蔽材料的密度，单位g/cm³；N_A是阿伏伽德罗常数，单位mol^–1；M是屏蔽材料的摩尔质量，单位是g·mol^–1。混凝土的σ·ρ·N_A·M^–1计算方法如式(2)所示，

式中：n_i是第i种元素归一化的原子数目占比；σ_i是第i种元素的中子反应截面；M_i是第i种元素的摩尔质量。

从图5可以看出，对于非弹性散射反应(主要针对几MeV以上能区的中子)，铁和铅的反应截面差距不大，混凝土的反应截面要比铁和铅的小一个量级左右。对于弹性散射反应，铁的反应截面整体上明显大于铅的，除过中子能量约10^–3 eV以下时铁的反应截面略小于混凝土的，其余能量段均大于混凝土的；在10^–1 eV以下能区铅的反应截面明显小于混凝土的，在10^–1 eV至10⁵ eV能区铅和混凝土的反应截面基本一致，10⁵ eV以上时，铅的反应截面比混凝土的大。对于俘获反应(主要针对低能中子)，铁的反应截面比铅的高一个量级左右，铅的反应截面比混凝土的高一倍左右。整体来看，铁在各个能量段的中子反应截面都要大于铅和混凝土的；铅的非弹性散射反应截面明显大于混凝土的，其低能中子俘获反应截面也稍大于混凝土的，但其弹性散射截面比混凝土的小；另外铅、铁等重材料对俘获反应伴生的$\gamma $射线的屏蔽能力也要比混凝土的更强^[13]。

在混凝土、铁屏蔽体中，前向的周围剂量当量在屏蔽厚度0到100 cm时衰减较为缓慢，屏蔽厚度达到100 cm以后衰减速度会增快。结合Agosteo等^[18-19]之前的研究中使用的拟合方法，选用屏蔽厚度100 cm以上的范围对混凝土和铁屏蔽体的前项数据进行拟合，将式(1)简化为拟合式(3)：

在混凝土、铁屏蔽体中，侧向和后向周围剂量当量在屏蔽厚度0到100 cm时衰减较快，100 cm以上时衰减变慢。如果仅利用式(3)进行拟合，会造成辐射源项的严重低估^[10]。因此采用双指数式(4)对混凝土、铁屏蔽体的侧、后向数据进行分段拟合，拟合区间分别选取0到100 cm以及100到600 cm。

铅屏蔽体厚度达到400 cm以上时，周围剂量当量的衰减趋势加快，但考虑到厚度400 cm以上时周围剂量当量非常低，实际屏蔽设计工作中几乎不可能采用厚度400 cm以上的铅屏蔽体，故选取100 到400 cm的屏蔽厚度，采用拟合式(5)对铅屏蔽体中的衰减数据进行拟合。

400 MeV/u碳离子打靶产生的次级辐射场具有很强的前冲性，在束流正前向(0°方向)屏蔽厚度估算时，如果采用0°~10°数据的拟合参数，会对辐射源项值H₀产生约2到5倍的低估，因此补充了发射角0°~1°的相关数据。根据前文提到的拟合方法，使用Origin软件对衰减曲线进行了拟合，拟合结果如表3至表8所列。表中，式(4)的拟合结果列入H₁、λ₁、H₀、λ₀四栏；式(3)和式(5)的拟合结果列入H₀、λ₀两栏；$R_0^2 $、$R_1^2 $分别表示两次拟合数据的相关系数。

福建莆田市即将建造一台医用重离子加速器，下面以该装置水平治疗室的屏蔽设计为例，介绍本文给出的屏蔽参数的使用方法。该装置加速粒子为¹²C，最大能量为400 MeV/u，治疗室内最大束流损失参数为1×10⁸ pps。治疗室的FLUKA计算模型如图6所示，屏蔽体材料选用密度2.3 g/cm³的混凝土，靶体为厚度30 cm的水靶。

FLUKA模拟计算得到的治疗室周围剂量当量率分布如图7所示。蒙特卡罗计算结果显示，当屏蔽体外周围剂量当量率限值设定为2.5 μSv/h时，图6中A位置所在方向需屏蔽体厚度438.21 cm；B位置所在方向(约45°发射角处)是侧向场强最大的位置，所需屏蔽体厚度为168.31 cm；C位置所在方向需屏蔽体厚度128.91 cm；D位置所在方向所需屏蔽体厚度70.6 cm。

下面采用本文得到的屏蔽参数计算所需的厚度。周围剂量当量率与单个粒子带来的周围剂量当量的换算关系如式(6)所示。式中，H^*是周围剂量当量率，单位是μSv/h；H是单个粒子产生的周围剂量当量，单位是Sv·ion^–1；I是平均束流损失，设计值是1×10⁸ pps。当H^*=2.5 μSv/h时，计算得到H=6.944 4×10^–18 Sv·ion^–1。

选取图5中A、B、C、D四个位置进行计算。A位置选用0°~1°的屏蔽参数，将下列数据代入式(3)中：

计算得到A位置所需屏蔽体厚度d为447.13 cm。

与A位置同样的方法，选用40°~50°的屏蔽参数进行计算，计算得到B位置所需的屏蔽体厚度d为170.29 cm。

C、D位置的计算方法与A位置稍有不同，选用80°~90°、170°~180°的数据进行计算。计算时先通过H₁、λ₁(厚度100 cm以下的拟合数据)计算得到厚度d₁，若d₁<100 cm，则d₁即为所求厚度；若d₁>100 cm，则使用H₀、λ₀(厚度100 cm以上的拟合数据)进行计算，得到的厚度d₀即为所求厚度。经计算，C位置所需的屏蔽厚度为126.25 cm，D位置所需屏蔽厚度为70.36 cm。

使用FLUKA直接模拟计算和通过本文给出的屏蔽参数得到的结果如表9所列。可以看出两种方法的计算结果符合很好，相对误差小于3%，说明本文给出的屏蔽参数是可靠的，能够直接应用于医用重离子加速器的屏蔽设计工作。

为了给医用重离子加速器提供一种专用的快速计算屏蔽厚度或对蒙特卡罗模拟结果进行验证的方法，本文使用FLUKA程序完成了400 MeV/u碳离子打厚铁靶、水靶的屏蔽参数计算。计算结果显示打靶产生的次级辐射场有很强的前冲性。未经任何屏蔽时，打靶产生的周围剂量当量几乎全部由次级中子与次级质子贡献；经一定厚度屏蔽后，周围剂量当量几乎全部由次级中子贡献。接着给出了打靶产生的次级辐射在混凝土、铁、铅屏蔽中的周围剂量当量衰减曲线，拟合得出了相应的辐射源项值和衰减长度。计算得到的屏蔽参数在莆田市医用重离子加速器的治疗室屏蔽设计中得到了应用，屏蔽参数的计算结果与FLUKA直接建模的计算结果符合很好。本文给出的屏蔽参数可以为同类医用重离子加速器的屏蔽设计工作提供可靠的参考数据。