核爆炸地面感生放射性辐射环境的蒙特卡罗模拟

商 鹏，黄流兴，牛胜利，左应红，朱金辉

(西北核技术研究所，西安 710024)

核爆炸释放的中子经大气输运后，一部分入射到地面的中子被地面介质中的某些元素俘获，发生活化反应产生放射性核素。这些放射性核素衰变放出β射线和γ射线，称为感生放射性。土壤的感生放射性是地面放射性沾染的来源之一，在爆高较高的条件下，土壤感生放射性是地面放射性沾染的主要来源[1]，对周围人员和环境的影响不容忽视。

土壤的感生放射性主要由中子与土壤中的多种元素与中子发生俘获反应引起的。原子核俘获中子后被活化，形成激发态的复合核，该复合核会在很短的时间内(约10-14s)退激发向外辐射出一条或多条γ射线，称为中子俘获反应，即 (n,γ) 反应[2]。土壤中能被活化的元素很多，一般而言，活化后的感生放射性主要来自28Al，56Mn，24Na，59Fe，对应的靶核分别是27Al，55Mn，23Na，58Fe。

一般用距地面一定高度处的γ剂量率来表征感生放射性形成的辐射环境。文献[1, 3-4]中采用点源近似法[5]，结合经验累积系数修正推导了探测点处的γ放射性活度公式，并通过半无限空间积分，进一步给出了土壤感生放射性在探测点处产生的γ剂量率表达式。上述方法在计算过程中进行了较多的简化和近似，如中子注量采用能量小于0.5 eV的等效热中子[6]近似及中子活化截面采用等效热中子平均活化截面近似等，只能给出粗略的结果。

为合理而相对准确地计算感生放射性，武祯等[7]针对电子直线加速器感生放射性的计算，对比分析了点源近似法、径迹长度法、功率吸收法和蒙特卡罗模拟法，研究表明，蒙特卡罗模拟法考虑反应全面，可同时计算直接和间接感生放射性，计算步骤简单，结果准确，是计算感生放射性的“最好”方法。本文细致分析了感生放射性的物理过程，基于蒙特卡罗方法建立了核爆炸地面感生放射性两步计算方法。

1 核爆炸地面感生放射性两步计算方法

蒙特卡罗模拟方法能较逼真地描述大量粒子在介质中输运的物理过程，且收敛速度与问题维度无关，程序结构简单，易于实现，因此，蒙特卡罗方法也被喻为“最后的方法”。模块化粒子输运程序包PHEN[8]由西北核技术研究所开发，能模拟中子、光子、电子及质子与物质的相互作用过程，如中子的弹性散射、非弹性散射和辐射俘获等，光子的光电吸收、康普顿散射和电子对产生等。它还可针对特定物理问题进行进一步开发，统计多种物理量。目前，PHEN程序并不具备直接模拟感生放射性的功能，需通过构建中子输运分层统计模型和γ辐射体源模型，并采用模型参数转换技术使二者结合，进而开展模拟计算，即地面感生放射性两步计算方法，计算流程如图1所示。

1.1 中子输运分层统计模型

核爆炸后输运到地面的有直射和散射中子，本文建立了中子输运分层统计模型，如图2所示。

中子源[6]各向同性发射中子，在均匀大气中输运，最终入射到地面土壤中，与土壤中的靶核发生相互作用，产生放射性核素28Al，56Mn，24Na，59Fe，记录不同深度处放射性核素的数密度。考虑到不同深度的土壤活化程度并不相同，传统方法中采用的均匀化假设与实际情况相差较大，所以采用中子输运分层统计方法。中子与土壤介质中的原子核发生相互作用，可能发生(n,γ)，(n,α)，(n,p)，(n,2n)，(n,3n)等多种反应。其中，感生放射性γ剂量率计算主要关注中子在土壤中输运并发射γ射线的过程，即(n,γ)反应过程。因此，需记录中子输运过程中发生的(n,γ)反应的次数。在PHEN程序包中，将中子在输运过程中与土壤某核素发生的所有(n,γ) 反应进行统计，得到土壤中该核素被活化后生成的放射性核素个数；对不同的核素重复上述过程，得到所有放射性核素的数密度。

上述计算过程可表示为

(1)

其中：Φ(E)为分能群中子注量，cm-2；σc(E)为(n,γ)反应的截面，b；常数η为任意可用于归一化处理的标量，也称为归一化因子。当η的单位为b-1·cm-1时，式(1)的计算结果即为单位体积内发生(n,γ)反应的总数nc，它在数值上与放射性核素的数密度nrad相等。

设距爆心为R，距地面深度为x，能量为E的中子注量为Φ(R,x,E)，则距爆心R处，不同深度x，第i种的放射性核素的数密度nrad(R,x，i)可表示为

(2)

其中：ni为第i种靶核的原子数密度，cm-3，与土壤中该元素的含量ci及天然丰度αi有关，可表示为

(3)

其中：NA为阿伏加德罗常数，mol-1；ρg为模拟区域土壤介质的密度，取1.4 g·cm-3；Mi为第i种靶核的摩尔质量，g·mol-1，数值与相对原子质量Ai相等。对比式(1)和式(2)可知，归一化因子ηi=ni×10-24b-1·cm-1。

计算中使用的几种元素的计算参数及由式(3)计算得到的归一化因子如表1所列。

表1 27Al，55Mn，23Na，58Fe的相关计算参数Tab.1 Computational parameters of 27Al，55Mn，23Na and 58Fe

1.2 γ辐射分层体源模型

核爆炸发生后，土壤活化产生的放射性核素发生衰变，向外各向同性地放出γ射线，这个过程相对较慢，导致γ剂量率随时间变化。而放射性核素单次衰变放出γ射线并被探测器探测的过程可视为瞬态过程。

28Al，56Mn，24Na，59Fe衰变放出的γ射线参数如表2所列。由表2可知，衰变γ具有离散的能量值，且分支比rB也有所差别，这是由放射性核素自身性质决定的。

表2 28Al，56Mn，24Na，59Fe的衰变γ射线参数[9]Tab.2 Parameters of decay γ rays for 28Al，56Mn，24Na and 59Fe

为模拟放射性核素出射γ射线对探测点探测剂量的贡献，需将中子输运分层统计模型计算得到的放射性核素数密度分布转换为γ体源模型中不同深度处的γ射线的发射强度。同时，γ辐射源还应当体现表2列出的能量、分支比及产额等特性。放射性核素发射γ射线数与放射性核素数密度及分层栅元体积相关，因此，构建γ辐射分层体源模型，模型特征参数由中子输运分层统计模型进行参数转换得到，如图3所示。

如1次衰变的γ射线产额为Y，能量离散分布，分布概率即为对应能量的γ射线分支比。将源权重赋值为Y，则模拟结果的物理含义为体源中发生1次衰变放出的Y个具有不同能量的γ射线对探测点探测剂量的贡献。

由上述分析可知，γ辐射分层体源模型中的辐射源为一个能量离散分布、位置具有特定分布及各向同性的复杂体源。离散能量抽样和方向各向同性抽样算法可参考文献[13]。γ射线在分层体源模型中的位置抽样采取封闭圆柱体包含多个体源栅元的舍选方法(rejection method,RM)，如图3所示。由于地面不同深度处的放射性核素数密度不同，衰变发射γ射线的强度也不同，即每个栅元发射γ射线的概率与栅元内的放射性核素个数成正比。设栅元j经归一化后的γ射线发射概率为Pj，在模拟程序中对γ射线发射位置进行抽样时，设置一个虚拟圆柱体将栅元体源完全包围，并在虚拟圆柱体内进行均匀抽样；若抽样得到的粒子位置位于栅元j中，则以Pj的概率接受该位置抽样，否则舍去本次抽样，重新抽取位置信息；重复上述舍选过程，直到一次抽样的位置信息被接受。虚拟圆柱体需完全包围发射γ射线的多个栅元，但不应设置得过大，否则将降低位置抽样效率。

利用PHEN程序包分别计算28Al，56Mn，24Na，59Fe单次放射性核素衰变放射的γ射线对探测点的平均剂量贡献Di。对于整个多层γ辐射体源而言，栅元j中第i种放射性核素的数密度ni,j与体积Vj的乘积之和即为第i种放射性核素总数Ni，表示为

(4)

如果放射性核素全部衰变，则衰变产生的总剂量D(t)为所有放射性核素的剂量贡献之和，表示为

(5)

(6)

其中：λi为第i种放射性核素的衰变常量；t为衰变时间。由式(6)可知，γ衰变剂量率受到放射性活度的影响，同样具备指数衰减规律。

2 两步计算方法应用实例

以爆高为200 m，典型氢弹(归一化到1 kt当量)的核爆炸场景为例，利用第1节所述地面感生放射性两步计算方法，研究了γ辐射分层体源的半径和厚度对感生放射性的影响，选取了一组较为合理的模型参数，并采用优化参数计算了地面上方1 m处的γ总剂量率及放射性核数剂量贡献比。

2.1 土壤层厚度及半径参数

文献[1]指出，地面上方1 m处的γ剂量率主要来自半径为10～20 m和深度为10～20 cm土壤薄层。本文通过建立中子输运分层统计模型，模拟计算了半径为20 m，土壤层厚度为2 cm和4 cm时，放射性核素28Al的数密度随深度的变化关系，如图4所示。

由图4可见，放射性核素数密度最大值均出现在10～20 cm深度处，与文献[1]所述一致。可用图5所示27Al的中子辐射俘获截面随能量的变化关系来解释：入射中子中能量较高的快中子成分在土壤中输运并在被慢化过程中产生更多的热中子，使土壤一定深度处的热中子水平高于地面，而热中子的辐射俘获截面比快中子高的多，导致由辐射俘获反应产生的放射性核素的数密度随深度呈先上升后下降的趋势。放射性核素的数密度随深度的变化关系如图6所示。

改变γ辐射体源模型土壤层的厚度，计算得到不同厚度γ辐射分层体源模型中28Al 的总数及衰变γ总剂量，如表3所列。由表3可知，采用蒙特卡罗方法模拟计算地面γ剂量，模型厚度小于60 cm时，有一部分衰变γ穿透土壤层，被终止输运，导致测点有效记数偏低；模型厚度大于60 cm后，记数基本稳定。因此，后续计算中，土壤层厚度取60 cm。

确定土壤层厚度后，改变γ辐射体源模型的半径，得到衰变γ总剂量随半径的变化关系，如图7所示。由图7可见，模型半径为40 m与模型半径为100 m时，衰变γ总剂量相差不超过5%。为节省粒子输运时间提高计算效率，取模型半径为40 m。

表3 不同厚度γ辐射分层体源模型中28Al 的总数及衰变γ总剂量Tab.3 N of 28Al and D for γ layered-radiator source model with different thickness

2.2 地面衰变γ剂量率计算

2.1节研究表明，采用蒙特卡罗方法计算土壤中放射性核素出射γ对地面上方探测点探测剂量的贡献时，取模型半径为40 m，厚度为60 cm可满足计算要求。

按照两步计算法计算流程，利用中子输运分层统计模型计算给出感生放射性核素的数密度分布，利用放射性核素数密度分布及放射性核素衰变特性构建γ辐射分层体源模型，模拟样本数为1×107个光子，分别计算28Al，56Mn，24Na，59Fe单次衰变的平均剂量贡献与放射性核素总数，如表4所列。

表4 28Al， 56Mn，24Na，59Fe的单次衰变平均剂量贡献与放射性核素总数Tab.4 Average dose contributions of single decay and total radionuclide amounts for 28Al， 56Mn，24Na and 59Fe

将放射性核素单次衰变平均剂量贡献与总数代入式(6)，可得土壤感生放射性总剂量率随时间的变化关系，如图8所示。由图8可见，总剂量率整体呈下降趋势，但下降速率发生了几次较大的变化，使总剂量率变化曲线表现出“台阶”状变化的形式。

这一时间特性同样可从图9所示放射性核素的剂量贡献比随时间的变化关系分析得到。由图9还可见，4种放射性核素分别在不同的时间区间上发挥主要作用，具体而言，在爆炸后的数分钟内，地面感生放射性主要来自28Al，强度比其他放射性核素高2个量级左右；在数小时范围区间，主要来自56Mn和24Na；而在200 h之后，基本只来自59Fe。

3 结论及展望

本文基于蒙特卡罗方法建立了核爆炸地面感生放射性两步计算方法。通过敏感性分析给出了构建中子输运分层统计模型和γ辐射分层体源模型的典型参数，修正了传统半经验公式方法中的等效热中子近似和等效热中子平均活化截面近似，考虑了土壤活化后放射性核素数密度随深度的变化及γ射线在空气中输运时的散射和衰减，并更细致地描述了感生放射性辐射环境形成的物理过程。本文方法还能应用于人工放射性辐射源活化其他材料后的γ放射性辐射环境计算，为高精度辐射环境计算提供了有效技术手段。