粒子辐射问题计算通常有输运方程法、蒙特卡罗法(MC法)、实验测量法以及经验法等几种方法。蒙特卡罗计算法又称随机抽样法或统计试验法, 是基于计算机模拟的思想, 抓住物理过程的数量和几何特征, 进行数字模拟试验, 该方法是求解辐射输运问题的一种相当成熟和有效的方法, 而且它对于各种复杂问题, 具有良好的通用性, 实用性相当广泛, 几乎涉及核科学的各个领域。
2 蒙特卡罗方法的解题能力和特点粒子输运过程可以用玻耳兹曼方程加以描述, 然而, 以此基础上发展起来的近似数值方法如扩散近似法、离散坐标方法在处理截面与能量相关以及散射各向异性介质、复杂几何条件问题时碰到了较大困难。而蒙特卡罗方法在处理这类问题时得心应手, 有很强的解题能力, 并且近似较少, 接近于真实情况。
辐射与物质相互作用的物理过程服从概率统计规律, 描述和处理这些现象的数学工具是概率论与数理统计, 这正是蒙特卡罗方法的数学基础。对于粒子与介质相互作用问题的模拟, 最基本的实质是描述粒子的状态矢量S(r, Ω, ω, t)。
首先, 在粒子所在位置r处, 根据介质的核素组成, 抽样得到将发生作用的核素, 继而查询核反应截面数据库, 得到粒子与该核素可能发生的各种作用截面σi, 则σi/∑σi即为发生第i种作用的概率, 然后处理每种反应的相应结果。若粒子被吸收, 则该历史结束; 若被散射, 则可通过抽样方法计算出散射后的能量E、运动方向Ω, 并根据其平均自由程抽样, 将会得到下一次输运的步长, 从而确定粒子在空间新的位置。粒子输运终止条件通常是粒子逃离所关心的介质范围或在介质中被吸收等。
在满足统计无偏的条件下, 通常需要添加一些数学技巧, 如粒子权重(ω)、隐式模拟、分裂与轮盘技巧等, 以提高计算效率。对于结果的记录, 可分为控测器通量、面通量与注量、体通量、剂量等, 这些结果是通过对大量粒子进行输运后, 作统计平均后得到的, 故它们都是具有一定误差的随机变量。
当然, 蒙特卡罗方法在解粒子输运问题中依然存在收敛速度慢以及概率性质误差等问题。除此之外, 经验证明, 只有当系统的大小与粒子平均自由程可以相比较时, 即十多个平均自由程时, 结果较为满意。而对于大系统, 算出结果偏低, 但是数值方法适应于大系统的计算, 得到结果较好, 现在已经有人将两者结合起来使用, 取得了一定效果。
3 MC方法在核科学领域中的应用中子按能量大小可分为热中子、中能中子、快中子和高能中子。具有以上能量的中子与物质元素原子核的微观相互作用类型有:弹性碰撞、非弹性碰撞、辐射俘获。宏观过程可概括为减速、扩散等行为。中子与原子核轻的物质相互作用可使中子能量一次损失很多而最终降至热运动水平, 如果这时物质吸收截面又很小, 则最终会形成平衡辐射场。中子辐射场是指中子源以外的空间, 描述中子场的参数主要有注量、能谱以及角分布。另外, 在体积较大的中子源(例如反应堆)内部, 各部分发出的中子也会有自己的各种平衡分布情况。
MC模拟计算软件可以应用在原子能研究, 辐射防护, 剂量学和辐射生物物理学等核科学领域。用MC法计算粒子输运的大型通用软件包主要有MCNP[1]、EGS4、ITS、GEANT、SIMSET等。以下基于中子辐射的上述各种特性, 列出MC方法用于中子辐射研究领域的应用实例。
3.1 计算筛选并优化中子辐射照射条件调整中子/γ混合场的中子/γ比通常采用屏蔽法, 对γ射线做适当的过滤, 而利用MC法可对屏蔽材料进行模拟计算, 从而得出最优的屏蔽条件。Endo S[2]等进行了裂变中子辐射生物效应的研究, 由于所用的裂变中子源252Cf发出的γ辐射成份占33%, 为减小这一份额使辐射生物效应实验能够在中子剂量较高的辐射场中进行, 使用了MCNP中子和光子输运软件对Al、Fe、Pb和混凝土几种过滤材料进行MC模拟计算, 组织等效TE材料的元素组成取自ICRU报告的软组织数据, 利用Maxwell-Boltzmann公式得出中子能谱, 用于MCNP模拟计算的源谱输入, 计算的结果得到的最佳铅过滤板厚度为4 cm, 在此条件下, γ成份已减小到总剂量的6.7%, 而中子剂量只减小10%。设置好过滤装置后, 使用TE-TE电离室和7LiF(Mg, Cu, P)TLD进行了剂量测量。MCNP计算与实验测量在6 %以内相符合, 并得出了最优的细胞生物辐射条件、细胞存活和RBE值。研究证明, 减小γ成份后再观察辐射生物效应的方法是可行的, 这要比不经过过滤的照射而估算出γ成份效应的方法更准确、可靠。
3.2 计算硼中子俘获利用在组织深部产生的热中子与肿瘤细胞内的10B发生中子俘获反应, 这种方法称之为BNCE, 即硼中子俘获增强。增强的剂量与10B集中的浓度以及组织某深部热中子通量有关。研究表明, 使用MC方法模拟计算BNCE剂量时, 通过对准直系统进行调整, 可改进快中子辐射的BNCE效果。从剂量学角度说, 剂量增强与组织某深部处所测量出的热中子注量有关, 而且热中子注量与射野密切相关, 射野减小时, 热中子注量就会相应减小。
Pignol JP等[3]使用计算带电粒子输运的FLUKA软件包计算来自质子与铍靶反应的初始中子场, 在此基础上, 利用MCNP软件包模拟中子在准直器中的几何空间输运过程, 得出体模内某深部快中子谱、能量沉积、热中子注量和热中子谱, 在10 cm ×10 cm和20 cm ×20 cm射野下, 剂量分别增大4.6%和10.4%。该结果与TLD探测器的测量实验结果进行了比较, 二者符合得很好。FLUKA软件包和MCNP联合运用, 可以对整个物理过程进行较好的描述, 并具有很好的灵活性和有效性。
3.3 描述中子治疗设备放射源用于间接电离辐射的辐射治疗方案软件通常会使用吸收剂量的卷积计算或在一个体模内测量的吸收剂量数据来估计吸收剂量的分布。然而, 对于卷积计算会因为介质的非均匀性导致吸收剂量估计上的误差, MC法可提供更为精确的粒子输运和吸收剂量分布。MC计算的精确度取决于病人和设备的特殊数据资料, 准确和完整的放射源资料是进行MC治疗方案的第一步, 也是最重要的一步。
Bohm TD[4]等建立了入射带电粒子、Be靶、前准直器、过滤材料和准直器模型, LAHET软件包使用MC技术, 以内核层叠方式描述了核的相互作用和粒子反应截面物理机制, 进行入射带电粒子在Be靶中的输运计算, 得出产出中子的能谱和角分布, 然后再以相同的三维几何模型, 用MCNP软件模拟计算从靶中出来的中子输运过程, 之后穿过准直器装置。计算结果与在体模内测量的实验结果符合得很好。
临床中子照射恶性肿瘤使用一种常用的7Li(p, n)7Be加速器中子源, 其截面在几十keV至1.8 MeV能量段增加得很快, 如果入射质子的能量在阈值附近, 则产出的中子将有相对低的能量, 只需较少的慢化剂即可达到慢化的目的, 满足BNCT临床治疗的要求。因为其相对低的发射加速成本, 中子源和靶可以装成轻便式的结构, 这种基于能量阈值附近的BNCT技术可作为医院广泛采用的治疗方法。
C.L.Lee[5]等利用MC方法针对7Li(p, n)7Be加速器中子源, 研究了质子入射能量在阈值附近时中子束的剂量特性, 使用了指示中子束的穿透能力、周围健康组织吸收剂量与肿瘤吸收剂量以及治疗时间等几个物理量, 建立了计算模型, 使用MCNP软件对慢化剂效应、热中子衰减、光子衰减和靶冷却装置进行优化计算, 利用中子总注量和注量-比释动能因子, 计算出中子和光子剂量, 并与其它BNCT设备进行了比较。研究的结果表明质子的入射能量介于1.93~1.99 MeV之间, 使用5 cm的水慢化剂和薄6Li和Pb屏蔽层, 可提供治疗时间小于1小时、加速电流小于5 mA的临床所用的中子射线。
3.4 计算反应堆附近的中子、γ谱和相应的剂量学量利用MC软件可模拟裂变反应堆周围的中子辐射场, 计算出反应堆周围的中子能谱和剂量学参数。Pesic MP[6]等为一座实验反应堆(RB)建立了一个三维几何计算模型, 该模型描述的堆芯和中子/γ混合辐射场与真实情况比较接近。用MCNP计算中子和γ辐射获得了稳定、可靠的结果, 得出了RB反应堆周围不同辐射点处的中子和γ能谱, 利用注量-剂量转换因子, 将中子和γ注量率转换成空气中的吸收剂量率。对于反应堆全部能量范围内的中子和能量范围为30 keV~8.0 MeV的γ射线, 其统计误差小于10 %, 中子注量在轴向空间分布的不确定度为±1.5%, 并将计算值与国际比对的实验测量值进行了比较, 实际计算得出的辐射点位置上空气中的中子吸收剂量比测量值高24.6%。通过对这些差别进行分析, 得出的结论是RB反应堆运行功率高于所报告的数值。在对反应堆功率进行修正之后, 将6组中子能谱和实验测量值和计算值进行比较, 两者符合得非常好。同样, 实验的测量值与MCNP计算的能谱值对每个能量组除以该修正系数, 中子和γ谱以及相应剂量学物理量的实验测量值与计算值符合得很好。
4 结语利用MC方法进行中子辐射的输运计算近年多集中于中子辐射治疗[7~9], 其它领域的研究文献报道较少, 目前仍需解决的问题包括:针对所要模拟的设备, 建立更为精确的输运计算模型; 对MC计算软件应做进一步的改进, 减小对中子能谱计算误差; 增加对MC计算的实验验证; MCNP软件还缺少高能中子反应截面数据库, 对于微剂量学计算, 还需要次级带电粒子的径迹。
| [1] |
Briesmeister J F.MCNP, A general Monte Carlo N-particle transport code, Version 4B-Manual, Los Alamos[CP/CD]. NM: LANKLLA-12625-M; 1997.
|
| [2] |
Endo S, Stevens D L, Bonner P, et al. Reduction of the gammaray component from fission neutron source-optimization for biological irradiations and comparison with MCNP code[J]. Phys Med Biol, 1999, 44: 1207. DOI:10.1088/0031-9155/44/5/009 |
| [3] |
Pignol J P, Cuendet P, Brassart N, et al. Combined use of FLUKA and MCNP-4A for the Monte Carlo simulation of the dosimetry of 10B neutron capture enhancement of fast neutron irradiations[J]. Med Phys, 1998, 25(6): 885. DOI:10.1118/1.598264 |
| [4] |
Bohm T D, Deluca P M, Cox L J, et al. Monte Carlo calculations to characterize the source for neutron therapy facilities[J]. Med Phys, 1999, 26(5): 783. DOI:10.1118/1.598596 |
| [5] |
Lee C L, Zhou XL, Kudchadker RJ, et al. A Monte Carlo dosimetry-based evaluation of the 7 Li (p, n)7Be reaction near threshold for accelerator boron neutron capture therapy[J]. Med Phys, 2000, 27(1): 192. DOI:10.1118/1.598884 |
| [6] |
Pesic M P, Ninkovic M M. Comparison of the MCNP calculated and measured radiation field quantities near the RB reactor[J]. Health Phys, 1999, 77(3): 276. DOI:10.1097/00004032-199909000-00005 |
| [7] |
DeMarco J J, Solberg T D, Smathers J B. A CT-base Monte Carlo simulation tool for dosimetry panning and analysis[J]. Med Phys, 1998, 25: 1. DOI:10.1118/1.598167 |
| [8] |
Zamenhof R, Redmond E II, Solares G, et al. Monte Carlo- based treatment planning for boron neutron capture therapy using custom designed models automatically generated from CT data[J]. Int J Radiat Oncol Biol Phys, 1996, 35(2): 383. DOI:10.1016/0360-3016(96)00084-3 |
| [9] |
Wallace S A, Allen B J, Mathur J N. Monte Carlo calculations of epithermal boron neutron capture therapy with heavy water[J]. Phys Med Biol, 1995, 40(10): 1599. DOI:10.1088/0031-9155/40/10/003 |