摘要:蒙特卡罗方法(MC)是模拟核探测问题的理想方法, 用中子照射客体, 中子诱发产生非弹γ和俘获γ, 通过特征γ射线能谱和时间谱分析, 确定客体核素组成和重量百分比. 本文基于非弹γ和俘获γ时间门测量技术, 给出了脉冲源发射下探测器响应计数公式. 在中子与核作用产生次级光子方面, 采用期望值估计(expect value estimator, EVE)产光. 为了避免大量小权光子模拟带来的计算存储量增加, 设计了EVE产光与直接估计(direct estimator, DE)产光耦合. 仅增加少量计算时间, 便实现了特征γ射线解谱. 数值模拟在自主MC软件JMCT上开展, 计算结果初步验证了方法的正确有效性.
关键词: 蒙特卡罗方法/
核探测/
非弹γ射线/
俘获γ射线/
期望值估计

English Abstract

全文HTML

--> --> -->

蒙特卡罗方法(Monte Carlo, MC)具有模拟核探测问题的天然优势, 近年已广泛用于X-射线荧光分析, 在线中子俘获瞬发γ射线分析, 基于γ射线光谱的脉冲中子孔隙度测井, 隐藏爆炸物探测等. 相比X光常规探测, 中子诱发γ射线探测具有穿透能力强, 容易穿透包括钢、常规/化学武器包壳、核武器内部结构, 可用于确定客体的化学变化等, 这是核探测的优势所在.
目前核探测的主要手段是利用高能中子照射客体, 产生次级光子, 通过特征γ射线能谱、时间谱测量, 确定客体核素组成及份额. 在中子作用下, 绝大部分元素都可以发射可辨认的特征γ谱线, 慢中子可以引起除氮以外所有元素的非弹性散射反应, 并发射特征γ谱线. 利用特征γ射线原级线光子不随入射中子能量变化这一特点, 通过特征γ射线的能峰特征, 确定客体内所含核素, 通过特征γ射线直穿贡献, 确定客体核素份额. 方法可用于化学/常规武器甄别、隐藏爆炸物/毒品探测、放射性石油测井和探测器灵敏度优化设计等.
由于核探测环境复杂, 涉及空间r、能量E、方向Ω、时间t七维变量的Boltzmann方程求解, 一般采用连续能量MC求解. 目前包括著名的MCNP^[1]程序, 特征γ射线解谱还存在一定困难, 虽然F8计数理论上能够算出探测器响应-脉冲高度谱, 但代价大, 解很难收敛. 另外, 中子产生次级光子采用直接估计, 存在随机因素的影响, 无法保证不漏掉某些重要核素的贡献. 文献[2-6]在这方面开展了大量研究, 针对特征γ射线解谱, 发展了多项技巧. 我们参考了其部分研究成果, 近期在自主MC粒子输运软件JMCT^[7-12]上开发了特征γ射线分类标识计算及非弹γ射线和俘获γ射线时间箱计数方法. 特别在中子与核作用产生次级光子处理上, 先采用期望值估计(expect value estimator, EVE)产光, 分别统计原级线光子和原级连续光子对探测器的直穿贡献. 为了避免EVE产生的大量小权光子模拟带来的计算量和存储量的增加, 做完直穿贡献后, 光子历史结束, 立即回到原来的直接估计(direct estimator, DE)产光模式, 每次最多产生一个光子, 对其进行跟踪, 只记录散射贡献, 这样在不增加计算存储量下, 实现了特征γ射线解谱.

高能中子与核发生非弹(用(n, n' )表示)和俘获(用(n, γ)表示)反应时, 将产生次级光子, 用次级光子的特征γ射线来确定客体的核素组成和份额, 这是核探测相对X光探测的优势所在. 辐射俘获是吸收反应中最重要的反应之一, 其反应产物之一就是γ射线, 产生次级光子的主要反应道有(n, γ)和裂变(用(n, f )表示). 如果中子源的能量和强度较高, 则中子与核发生(n, n' )反应的概率增大, 并产生非弹γ射线. 表1给出H, C, N等11种核素发射非弹γ谱线和俘获γ谱线的能量^[13], 表2给出烈性炸药(TNT)和某些化学武器中所含元素的重量百分比^[13].

元素	反应类型	特征γ 谱线能量/MeV
H	辐射俘获	2.2233
C	非弹性散射	4.433
N	辐射俘获非弹性散射	1.8848, 5.2692, 5.5534, 6.3224, 7.2991, 10.8290, 2.3128, 4.4444, 5.1059, 7.0280
O	非弹性散射	2.7419, 3.6841, 6.1310, 6.9170, 7.1190
F	非弹性散射	0.1090, 0.1971, 1.2358, 1.3480, 1.3565
P	辐射俘获非弹性散射	2.1542, 3.5228, 3.9003, 4.6713, 6.7853, 1.2661, 2.2334
S	辐射俘获	0.8411, 2.3797, 2.9311, 3.2208, 4.4308, 4.8698, 5.4205,
Cl	辐射俘获	0.5167, 0.7884, 1.1647, 1.9509, 1.9591, 2.8639, 5.7153, 6.1109, 6.6195, 7.4138
As	辐射俘获非弹性散射	6.2941, 6.8094, 7.0192, 0.2646, 0.2795, 0.5725
Al	辐射俘获非弹性散射	0.9840, 2.9598, 4.1329, 4.2522, 7.7239, 0.8438, 1.0144, 2.2118
Fe	辐射俘获非弹性散射	0.3522, 1.7251, 5.9203, 6.0185, 7.6311, 7.6455, 8.8860, 9.2980, 0.8468, 1.2383, 1.4082, 1.8105, 2.1129, 2.5985

表1H, C, N, O等核素发射俘获γ谱线和非弹性散射γ谱线能量
Table1.Energy of spectrum line from inelastic γ and capture γ about H, C, N, O, etc.

元素	TNT	沙林/GB	神经性 毒气/VX	芥子气/HD	孁烂性 毒气/L
氢(H)	2.2	7.1	9.7	5.0	1.0
碳(C)	37.0	34.3	49.4	30.2	11.4
氮(N)	18.5		5.2
氧(O)	42.3	22.9	12.0
氟(F)		13.6
磷(P)		22.1	11.6
硫(S)			12.0	20.1
氯(Cl)				44.7	51.3
砷(As)					36.1

表2烈性炸药(TNT)和某些化学武器中所含元素的重量百分比
Table2.Weight percentage of elements in some spirited detonators (TNT) and chemical weapons.

2

-->

2.1.非弹γ和俘获γ时间门测量方法

对于某些问题, 需要通过发射脉冲源, 通过时间箱计数来区别非弹γ射线和俘获γ射线, 图1给出放射性石油测井中, 碳氧比(C/O)能谱测井的定时测量图^[14-16]. 用脉冲方式发射中子, 通常0—10 μs为脉冲门发射时间间隔, 10—20 μs为本底门时间间隔, 20—90 μs为晚俘获门时间间隔, 依据这种逻辑关系, 可以得到: 净非弹γ计数=脉冲门谱计数–本底谱计数; 俘获γ计数=晚俘获门计数. 这个时间门测量过程同样适合其他含时核探测问题的模拟.
图 1 碳氧比测井中子引发非弹性散射γ与俘获γ定时逻辑图
Figure1. The timing diagram of neutron induced inelastic γ and capture γ in C/O well-logging.

对应图1给出不同时间门下探测器脉冲高度谱计算公式. 设N_δ(E, t)为以δ(t)脉冲方式发射的E₀ = 14.1 MeV氘氚中子源探测器中测到的γ能谱的时间响应. 根据中子-光子输运方程的线性性质, 对应于任意时间分布S(t)的中子源, 探测仪中测到的γ能谱时间响应为

$N(E,t) = \int_0^t {S(t'){N_\delta }(E,t - t'){\rm{d}}t'} ,$

其中E表示能量, 单位MeV; t表示时间, 单位μs.
在任意测量门[t_a, t_b]内记录的γ能谱为

$\int_{{t_a}}^{{t_b}} N(E,t){\rm{d}}t \!=\! \int_{{t_a}}^{{t_b}} {{\rm{d}}t} \int_0^t S(t'){N_\delta }(E,t \!-\! t'){\rm{d}}t'.$

令中子脉冲时间分布函数S(t)为周期函数, 其周期为τ, 有

$S(t \pm nt) = S(t),~~ n = 1,2, \cdots ,$

其中t∈[0, τ]. 又假定S(t)为宽度为τ的函数

$S(t) =\begin{cases}{{S_0}{f_\delta }(t),}&{\;0 \leqslant t < \tau ,}\\{0,}&{\;\tau \leqslant t < T,}\end{cases}$

其中S₀为一个脉冲内释放的中子总数; T为终态时间; ${f_\delta }(t)$为脉冲时间分布函数, 满足归一条件$\int_0^\tau {{f_\delta }(t){\rm{d}}t} = 1$. 则有

$\begin{split}N(E,t) =\; & {S_0}\left[ \sum\limits_{n = 1}^\infty \int_{ - nT}^{ - nT + \tau } {f_\delta }(t'){N_\delta }(E,t - t'){\rm{d}}t' \right.\\ & \left.+ \int_0^{\min (t,\tau )} {{f_\delta }(t'){N_\delta }(E,t - t'){\rm{d}}t'} \right]\\ =\;& {S_0}\left[\sum\limits_{n = 1}^\infty \int_0^\tau {f_\delta }(t'){N_\delta }(E,t - t' + nT){\rm{d}}t' \right.\\ & \left.+ \int_0^{\min (t,\tau )} {{f_\delta }(t'){N_\delta }(E,t - t'){\rm{d}}t'} \right],\\ & \qquad t \in [0,T].\\[-10pt]\end{split}$

将(5)式代入(2)式, 对t积分便得任意时间门内的γ能谱强度.
脉冲门内的测量值为

${N_{\rm{I}}} = \int_0^\tau {N(E,t){\rm{d}}t} .$

本底门内测量值为

${N_{{\rm{II}}}} = \int_\tau ^{2\tau } {N(E,t){\rm{d}}t} .$

根据中子进入客体诱发的各种γ射线的时间特点, 将(5)式中的N_δ(E, t)分解为4个部分, 即

${N_\delta }(E,t) = \sum\limits_{i = 1}^4 {N_i^\delta (E){f_i}(t)} ,$

其中$N_i^\delta (E), i = 1, 2, 3, 4$分别为单位强度δ(t)脉冲中子源引起的非弹性γ射线、慢化过程中的俘获γ射线、热中子俘获γ射线和活化反应γ射线的γ能谱强度; f_i(t), i = 1, 2, 3, 4分别为上述4种γ射线的时间谱, 满足归一条件$\displaystyle\int_0^\infty {f{}_i(t){\rm{d}}t = 1} $, (i = 1, 2, 3, 4).
经过一系列公式推导, 得到脉冲门能谱

$\begin{split}&{N_{\rm{I}}}(E) - {N_{{\rm{II}}}}(E)\\\approx \; & {S_0}\left\{N_1^\delta (E) + N_2^\delta (E)\left[\int_0^\tau \left( {1 - \frac{{2t}}{\tau }} \right){f_2}(t){\rm{d}}t \right.\right.\\ & \left.\left.- \int_\tau ^{2\tau } {\left( {2 - \frac{t}{\tau }} \right){f_2}(t){\rm{d}}t}\right] - \frac{\tau }{2}N_3^\delta (E){f_3}(\tau ) \right\}.\end{split}$

可见脉冲门与本底门测量值之差除了非弹性γ能谱, 还受到少量慢化过程俘获γ射线和少量未扣除干净的热中子俘获γ谱的“污染”. 考虑到活化γ射线的发射时刻很晚, 加之活化γ射线的发射寿命$\gg T $. 因此, 在脉冲门与本底门谱之差中, 活化γ射线的污染是很小的. 另外, 由于慢化过程中的中子俘获数目远远小于热中子俘获的数目, 所以, $N_3^\delta (E) \gg N_2^\delta (E)$, 可以认为污染源主要来自热中子俘获γ射线, 其中污染量大小除与中子源的特征及客体物理性质有关, 探测仪器材料的适当选择也起一定的作用.
要准确算出探测器响应的测量能谱N(E, t), 首先应设法算出N_δ(E, t), 然后与中子源脉冲时间谱形f_δ(t)卷积, 并分别对两个测量门的时间间隔积分, 相减后得到. 因此, MC模拟的关键是算出t时刻单位强度δ-脉冲源相应的探测器计数N_δ(E, t), 然后与给定的中子脉冲谱形函数f_δ(t)卷积, 再对指定的时间门积分, 得到任意时间门内的测量值.
2 -->

2.2.N(E, t)的MC模拟

快中子非弹性散射, 它所要测量的是非弹性γ射线. 理论上通过解中子-光子-电子耦合输运求出N(E, t). 由于探测器相对整个问题系统很小, 尽管从源发出了大量的中子, 但能够进入探测器的次级光子数还是非常有限的, 依据少量的次级光子要算准探测器响应几乎是不可能的. 因此, 探测器响应-脉冲高度谱计算通常分三步进行.
第一步: 解中子-光子耦合输运方程, 求出进入探测器表面的次级γ流,

$\begin{split}{J_\gamma }({E_0},t) =\; & \int_0^t \int\nolimits_{{r} \in S} \int\nolimits_{{\varOmega } \cdot {n} < 0}\int_0^{{E_0}}{n}\\ & \times {\varOmega }\phi ({r},E,{\varOmega },t){\rm{d}}t{\rm{d}}{r}{\rm{d}}{\varOmega }{\rm{d}}E ,\end{split}$

其中S为探测器表面, n为探测器表面外法向矢量, E₀为入射光子能量.
第二步: 解光子-电子耦合输运方程, 求出探测器响应函数R(E₀, h). 以碘化钠闪烁探测器为例, 其响应函数形式为

$\begin{split}R({E_0},h) ={}& \eta ({E_0})\int_0^{{E_0}} D({E_0},E)G(E,h){\rm{d}}E, \\& 0 \leqslant h < {E_0},\end{split}$

其中$G(E, h)$为高斯函数; $D({E_0}, E)$为能量沉积谱; $\eta ({E_0})$为探测器效率; h为能量道, 单位MeV. 由于响应函数中含有高斯函数G, 为了保证求积精度, 分点要足够多, 通常在全能区[0, E_max]分256道.
根据Berger等^[17]和Jin等^[18]的研究, 不同入射方向的光子对探测器响应影响很小, 因此, 方向可近似为各向同性, 当探测器形状及材料一定后, 响应函数R(E₀, h)(矩阵形式)只需计算一次.
第三步: 卷积积分得到探测器响应脉冲高度谱

$N(h,t) = \int_0^{{E_{\max }}} {R({E_0},h){J_\gamma }({E_0},t){\rm{d}}{E_0}} .$

$N(h, t)$求解涉及瞬态中子-光子-电子耦合输运计算, 计数包括时间-能量联合谱, 是MC粒子输运计算中难度最大, 模拟最复杂的过程.
2 -->

2.3.特征γ射线标识计数

中子发生非弹(n, n')或俘获(n, γ), (n, f )反应时, 将产生次级光子, 这些光子统称为特征γ射线. 根据次级光子能量特征, 特征γ射线分为原级线光子和原级连续光子, 其定义如下:

${E_\gamma }(i,j) =\begin{cases}{E_G^{(i)},}&{{\rm{LP}} \ne 2,{\text{原级线光子}},}\\{E_G^{(i)} + \dfrac{{{A_i}}}{{{A_i} + 1}}{E_{\rm{n}}},}& {\rm{LP}} = 2,{E_n} > {E_{{\rm{line}}}},\\ {}&{\text{原级连续光子}},\end{cases}$

其中i为核; j为反应道; E_n为中子发生非弹或俘获反应时的能量; 通常取E_line = 0.001 MeV作为原级线光子和原级连续光子的分界能量; LP为ENDF数据库反应律序号；A_i为碰撞核i的原子序数；E_G⁽ⁱ⁾为特征γ能量.
LP ≠ 2对应的γ原级线光子, 它不随入射中子能量变化, 是甄别客体核素组成的核心关键. 在表1中给出了主要核素γ谱线能量峰位置, 表2给出烈性炸药(TNT)和某些化学武器中所含元素的重量百分比^[13], 依据表1和表2所列参考值, 可以快速识别出客体是否为危禁品. MCNP程序模拟次级光子时, 没有把这两类光子区别开来, 因此, 次级光子能谱中的一些特征γ峰会被散射能谱磨平, 无法确定客体的核素组成.
当前ENDF最新评价核数据库提供了两种中子产光子模式: 1) 30 × 20产光模式, 中子从10^–5 eV到20 MeV分30个能群, 每个中子群对应20个等高度谱线, 次级光子的发射方向按各项同性处理; 2)按碰撞核对应反应道的产光概率产光子. 模拟采用第2)种产光模式.
2 -->

2.4.EVE产光

设中子与i核发生j种反应的产光概率为p_{i, j} (j = 1, 2, ···, J ), 满足归一条件$\displaystyle\sum\limits_{j = 1}^J {{p_{i, j}}} = 1$, 这里J为i核对应的反应道总数.
DE方法　已知i, 抽随机数ξ, 求出满足不等式

$\sum\limits_{k = 1}^{j - 1} {{p_{i,k}}} \leqslant \xi < \sum\limits_{k = 1}^j {{p_{i,k}}} $

的j, 则j反应道产光, 光子权重为w_γ.
EVE方法　已知i, 按概率p_{i, j} (j = 1, 2, ···, J )全部产光, 相应光子权重为

${w_{i,j}} = {p_{i,j}}{w_\gamma },~j = 1, \cdots ,J.$

不难证明两种估计方法的数学期望是一致的. DE存在随机因素, 某些小概率大贡献事件会因为随机因素少抽或漏抽, 这对隐藏爆炸物探测这类问题是不允许的. 故采用EVE是必要的. 我们早期开发研制的MCCO程序, 对EVE产出的光子进行全部跟踪, 用统计估计计数^[19], 增加了大量计算存储量. 考虑到特征γ射线探测主要关心的是次级光子的直穿贡献, 实际上跟踪所有小权次级光子的散射过程没必要. 为此, 我们设计了现在的组合产光模式.
组合产光模式　对EVE产生的大量小权光子, 仅做直穿估计, 之后回到原来的DE产光, 每次最多产生一个次级光子, 对该光子进行全程跟踪, 只统计散射贡献, 最后的数值模拟解由两部分组成, 即

${\text{光子总计数}} = {\rm{EVE{\text{光子直穿计数}} + DE{\text{光子散射计数}}}}{\rm{.}}$

由于直穿估计可以解析计算, 所花时间可以忽略不计. 因此, 采用组合产光模式后, 实现了特征γ射线能量-时间解谱, 而总计算时间和存储量增加很少.

炸药球模型 几何单位为cm, 对行李箱进行安检. 如图2所示, 设行李箱长、宽、高分别为80, 30, 50, 内放一半径为3.51, 密度为ρ = 1.654 g/cm³的RDX炸药小球(分子式为C_2.63H_4.69N_0.658O_0.85); 探测器位于行李箱上方, 为3'' × 3''碘化钠闪烁晶体正圆柱探测器, 密度为ρ = 3.67 g/cm³, 柱中心坐标为(0, 0, 38.75); 采用14.1 MeV各项同性氘氚中子点源, 源位置为(0, 45, 0), 计算进入探测器的次级γ流J_γ(E, t). 采用源方向偏倚发射, 偏倚立体张角覆盖行李箱. 对比程序选择MCNP^[1], 采用相同的数据库、样本数和源方向偏倚技巧.
图 2 行李箱模型示意图
Figure2. Sketch of luggage model.

表3给出JMCT与MCNP探测器次级光子流J_γ结果比较; 表4给出JMCT统计的箱子内的核素组成及份额, 与表1给出的TNT炸药H, C, N, O核素百分比基本相符(炸药类型不同); 图3(a)和图3(b)分别给出原级线光子和康普顿散射光子能谱; 图3(c)和图3(d)给出JMCT与MCNP次级γ流的能谱比较, 可以看出次级γ总流JMCT与MCNP结果符合良好, 能谱差异很小; 图4给出次级γ流时间谱, JMCT与MCNP时间谱总计数相符, 时间谱细节上有些差异, 分析原因是由于时间箱分得过细, 某些时间箱计数的统计误差偏大.

程序	原级线 光子	原级连续 光子	散射光子	Jγ 总 光子流	偏差/%
JMCT	4.92519-7	0	4.34947-8	5.36014-7	0.4406
MCNP	无	无	无	5.38386-7	标准解
注: 偏差 = [Jγ(JMCT) – Jγ(MCNP)]/ Jγ(MCNP).

表3JMCT与MCNP次级γ流计算结果比较
Table3.Comparison of calculated results about secondary γ between JMCT and MCNP.

图 3 次级γ射线能谱计算结果比较　(a)次级γ原级线光子能谱; (b)原级连续光子与Compton散射能谱; (c) JMCT次级γ总能谱; (d) MCNP次级γ总能谱
Figure3. Comparison of calculated result about energy spectra of secondary γ: (a) JMCT primary line γ; (b) JMCT Compton γ; (c) JMCT total γ; (d) MCNP total γ.

元素	计数	份额比/%	统计误差/%
H	3.02643 × 10–11	0	0.56
C	1.77077 × 10–7	36	0.49
N	1.11146 × 10–7	23	0.12
O	2.03254 × 10–7	41	0.18
注: 偏差 = [Jγ(JMCT) – Jγ(MCNP)]/ Jγ(MCNP).

表4H, C, N, O瞬发γ计数及份额
Table4.Count and percentage of prompt γ from H, C, N and O.

图 4 JMCT与MCNP次级γ流时间谱比较
Figure4. Comparison of secondary γ-fluent tine spectrum between JMCT and MCNP.

基于JMCT程序开发了特征γ射线能量-时间联合谱计算功能, 数值实验初步验证了软件计算的正确高效性. 目前已开展了数十种炸药模型的计算, 构建了相应的数据库, 为反演计算做准备. 后续将配合实验样机的研制, 对测量仪器进行刻度工作, 形成能快速做出判断的计算机软件正演/反演系统, 为在线检测提供理论技术支持.