1.   引　言

中子多重性(NMC)分析方法作为一种无损检测手段, 能够实现对密闭容器内放射性物质的定量分析, 在深度核裁军、核材料库房管理中发挥着重要作用^[1]. Bohnel^[2]于1985年提出的Bohnel方法是当前中子多重计数中最基本和最有用的方法之一. 为提高探测效率, Bohnel曾使用负载硼的塑料闪烁体(BC454)/锗酸铋(BGO)磷光体检测器阵列进行了研究^[3]. 随着美俄等国逐渐限制³He的对外出口, 造成供需不平衡, 导致价格的不断攀升. 为寻找有效的替代品, 液闪探测器逐渐登上舞台. 基于经典NMC分析方法发展而来的快中子多重性(FNMC)分析方法, 不需要经过高聚合物的慢化过程, 能够克服中子多重性测量过程中的不足, 有效保留中子的能量和时间信息, 在核材料衡算和属性认证中发挥着越来越重要的作用^[4]. 中子多重性测量根据有无外部中子源可以分为被动测量和主动测量, 如图1所示.

被动测量法不需要外部中子源的诱发, 主要依托于核材料的自发裂变反应, 多用于自发裂变率较高的材料, 例如²⁴⁰Pu, ²³⁸U. 针对此类材料的快中子多重性测量方法发展较为成熟, 建立了分析模型推导出了快中子多重性测量方程并进行检验修正^[5]. 在三阶方程的基础上还推导出了高阶测量方程, 但由于实际测量偏差较大, 应用效果不够理想^[6]. 对不同丰度、不同封装材料下的钚材料进行了研究, 并提出采用了$ M $, $ \alpha $系数修正的方法缩小测量偏差^[7,8]. 主动测量法通过外部中子源照射核材料进而促使其发生诱发裂变, 主要用于诱发裂变率较高的材料, 例如²³⁵U, ²³⁹Pu. 针对此类材料的快中子多重性测量方法起步较晚, 未见公开发表的测量方程. 美国密歇根大学Hua等^[9,10]采用有机闪烁体对²³⁷Np进行了快中子多重性测量, 但采用的却是热中子分析模型. 近年来, 研究机构加大了对铀材料的快中子多重性测量方法的重视, 对不同质量、密度、丰度的铀材料以及Am-Li中子源能谱对测量结果的影响进行了大量研究^[11]. 有研究通过建立铀材料的数学模型, 提出了$ R = F \cdot P \cdot Q \cdot T $的方程形式, 但依旧未能提出具体的快中子多重性测量方程, 不利用快中子多重性测量的深入研究^[12]. 经典的NMC方程实现了对样品参数的求解, 但却不适合于存在散射串扰的快中子多重性测量^[13], 推导铀材料的快中子多重性测量方程, 对核材料属性认证具有现实意义.

2.   原理模型

参考文献[14], 采用Geant4搭建一套井型测量系统进行模拟测量研究. 该测量系统由24个液闪探测器组成, 共3层, 每层均匀排列分布8个BC501A液闪探测器. 单个探测器的几何尺寸为$ \phi $120 mm × 60 mm, 探测器前端面距离探测系统的中心轴20 cm, 探测系统的几何中心设置球形样品(90%²³⁵U和10%²³⁸U), 空腔填充空气. 在距离样品正中心上方10 cm处放置外部中子源, 源强设置为1 × 10⁴ n/s, 能量为0.3 MeV, 系统布局结构如图2所示.

整个模拟测量过程共分为两部分, 一是利用Geant4搭建一套模拟探测系统, 通过对样品、探测器、物理过程、测量阈值进行设置, 实现对中子裂变及输运的全过程模拟, 生成list-mode数据文件, 记录下液闪探测器探测到的中子时间信息; 二是利用Matlab构建一个虚拟移位寄存器, 通过设置符合门宽、预延迟时间, 实现对中子多重计数率的模拟测量. 具体原理模型如图3所示.

3.   铀材料快中子多重性测量方程推导

3.1   基本假设

中子输运是一个极其复杂的过程, 在实际的探测中, 中子计数会受到一系列因素的影响^[15]. 为降低方程推导的复杂性和可行性, 通过对多重性计数影响因素的讨论引入基本假设, 简化方程的推导过程, 主要包括以下5点:

1)假设所有被俘获的中子都会引起诱发裂变;

2)假设中子衰退时间小于符合门宽;

3)假设所有裂变中子从同一点发射出, 即点模型假设, 保证中子探测器的效率和裂变概率在样品体积范围内是不变的;

4)假设每个中子在第1个探头中产生信号后都有相同的概率在第2个探头中产生信号, 但不会产生第3个信号, 即只考虑发生一重散射串扰的情况;

5)假设次级诱发裂变中子与初级裂变中子具有相同能谱, 在这个假设下, 两个中子源具有相同的探测效率、裂变概率和诱发裂变多重性.

3.2   推导过程

铀材料快中子多重性测量方程是在NMC分析方程的推导过程基础上, 根据物理过程的差异做出相应改变推导而出的, NMC分析方程的具体推导过程见文献[16, 17]. 在钚材料的快中子多重性测量方程的推导中, 通过建立一系列的基本假设, 将未知参数简化为3个, 即自发裂变率$ F $、样品增殖系数$ M $和$ (\alpha, n) $反应比例系数$ \alpha $. 而在铀材料多重性分析中通常不考虑$ (\alpha, n) $反应的影响. 因此在铀材料的快中子多重性测量中, 样品特征参数可以简化为两个: 自发裂变率$ F $、样品增殖系数$ M $. 本文引入概率母函数、阶乘矩和总体矩等数学工具推导测量方程, 通过对物理过程的分析建立起样品参数与多重计数率之间的函数关系.

方程从中子产生、输运、探测、计数等物理过程出发进行推导, 根据诱发裂变中子重数分布和中子在样品内输运过程中的数目变化, 求得一次诱发裂变的样品发射中子重数分布(即中子出射率), 再根据探测器的探测和计数方式推出真符合重数的总体矩$ {r_k}\left( {F, \alpha, M} \right) $; 另一方面, 根据前景计数分布f、背景计数分布b与真符合重数分布r之间的关系, 推出由测量计数表示的真符合重数样本矩$ {r_k}\left( {f, b} \right) $. 令总体矩与样本矩相等即可得到参数估计方程, 进而得到测量方程.

1)源事件的概率母函数

首先对增殖系数M作如下定义:

$$ M = {h'}(1), $$

(1)

其中$ {h'}(1) $是每个初级中子所致样品发射中子数的一阶矩. 每次发生诱发裂变时发射的中子重数是独立同分布的非负整值随机变量, 用概率母函数表示为

$$ {S_s}(z) = \sum\limits_{i = 0}^\infty {{p_i}{z^i}} \text{, } $$

(2)

其中p_i表示发射i个中子时的概率. 一次源事件发射中子数的概率母函数为

$$ S(z) = {S_s}(z) . $$

(3)

根据概率母函数的性质, 一次源事件发射中子数的一、二、三、四阶矩分别为

$$ {S'}(1) = {v_{s1}}, $$

(4)

$$ {S{''}}(1) = {v_{s2}}, $$

(5)

$$ {S{'''}}(1) = {v_{s3}}, $$

(6)

$$ {S{''''}}(1) = {v_{s4}} . $$

(7)

2)样品发射中子数的概率母函数

源事件产生1个中子, 从样品中发射出来的中子数的概率母函数记为$ {h_1}(z) $. 源事件产生n个中子, 则根据概率母函数的性质, 发射中子数的概率母函数为

$$ {h_n}\left( z \right) = {\left[ {{h_1}\left( z \right)} \right]^n}. $$

(8)

将样品中1个中子发生诱发裂变的概率记为$ {P_{if}} $, 诱发裂变产生i个中子的概率记为${P_{if\nu }}\left( i \right)$, 这i个中子所致样品发射中子数的概率母函数也为${h_i}\left( z \right)$.

一次诱发裂变产生中子数的概率母函数为

$$ U\left( z \right) = \sum\limits_{i = 0}^\infty {{z^i}{P_{if\nu }}} \left( i \right). $$

(9)

中子在样品中有可能被吸收, 也可能未被吸收, 未被吸收概率为1－$ {P_{if}} $. 假设中子在样品中只要被吸收都会诱发裂变, 忽略中子在样品中被吸收而没有诱发裂变的可能性, 根据概率母函数的性质有

$$ \begin{split} {h_1}\left( z \right) =\;& \left( {1 - {P_{if}}} \right)z + {P_{if}}\sum\limits_{i = 0}^\infty {{P_{if\nu }}\left( i \right){h_i}\left( z \right)} \\ = \;&\left( {1 - {P_{if}}} \right)z + {P_{if}}\sum\limits_{i = 0}^\infty {{P_{if\nu }}\left( i \right)h_1^i\left( z \right)} \\ = \;&\left( {1 - {P_{if}}} \right)z + {P_{if}}U\left[ {{h_1}\left( z \right)} \right]. \end{split} $$

(10)

为书写方便, 后面统一用$h\left( z \right)$代替${h_1}\left( z \right)$, 即

$$ h\left( z \right) \equiv {h_1}\left( z \right) = \left( {1 - {P_{if}}} \right)z + {P_{if}}U\left[ {h\left( z \right)} \right]. $$

(11)

(11)式给出了源事件发射(产生) 1个中子时, 样品发射中子数概率母函数$h\left( z \right)$与诱发裂变概率${P_{if}}$和诱发裂变发射中子数$U\left( z \right)$之间的关系. 考虑源事件发射(产生)中子数分布, 根据概率母函数的性质, 一次源事件中样品发射的中子数概率母函数为

$$ H\left( z \right) = s\left[ {h\left( z \right)} \right]. $$

(12)

3)源事件产生的泄漏中子数的阶乘矩

根据概率母函数的性质, 一次源事件对应的样品发射中子数的一、二、三阶阶乘矩为

$$ H'\left( 1 \right) = s'\left( 1 \right)h'\left( 1 \right), $$

(13)

$$ H''\left( 1 \right) = s''(1){\left[ {h'(1)} \right]^2} + s'(1)h''(1), $$

(14)

$$ \begin{split}H'''\left( 1 \right) = \;&s'''(1){\left[ {h'(1)} \right]^3} + 3s''(1)h''(1)h'(1) \\ &+ s'(1)h'''(1) ,\end{split} $$

(15)

$$\begin{split}H''''\left( 1 \right) =\;& s''''(1){\left[ {h'(1)} \right]^4} + 6s'''(1)h''(1){[h'(1)]^2} \\ &+ 4s''(1)h'''(1)h'(1) + 3s''(1){[h''(1)]^2}\\ & + s'(1)h''''(1).\\[-10pt] \end{split} $$

(16)

对(11)式求z = 1处的导数, 得

$$\begin{split} h'\left( 1 \right) = \;&\frac{{1 - {P_{if}}}}{{1 - {P_{if}}U'(h(1))}} = \frac{{1 - {P_{if}}}}{{1 - {P_{if}}U'(1)}} \\ =\;& \frac{{1 - {P_{if}}}}{{1 - {P_{if}}{v_{i1}}}} . \end{split}$$

(17)

考虑到二次诱发裂变中子阶乘矩可由核数据计算得到, 用$ {v_{i1}} $代替$ U'\left( 1 \right) $, 后面使用$ {v_{i2}}, {v_{i3}}, {v_{i4}} $代替${U{''}}(1), {U{'''}}(1), {U{''''}}(1)$. $ h'\left( 1 \right) $是每个源中子在输运过程中的增殖期望, 记为M, 则M为

$$ M = h'\left( 1 \right) = \frac{{1 - {P_{if}}}}{{1 - {P_{if}}{\nu _{i1}}}} . $$

(18)

求z = 1处的二阶导数, 可获得中子增殖的二阶矩、三阶矩和四阶矩:

$$ {h{''}}(1) = {M^2}{v_{i,2}}\frac{{M - 1}}{{{v_{i1}} - 1}} \text{; } $$

(19)

$$ \begin{split} h'''(z)=\;&\frac{{P}_{if}(1-{P}_{if})U'''(h(z)){(h'(z))}^{2}}{{\left[1-{P}_{if}U'(h(z))\right]}^{2}}\\ &+\frac{{P}_{if}(1-{P}_{if})U''(h(z))h''(z)}{{\left[1-{P}_{if}U'(h(z))\right]}^{2}}\\ &+\frac{2{({P}_{if})}^{2}(1-{P}_{if}){\left[U''(h(z))\right]}^{2}{(h'(z))}^{2}}{{\left[1-{P}_{if}U'(h(z))\right]}^{3}},\end{split} $$

$$ h'''(1) = {M^3}\frac{{M - 1}}{{{\nu _{i1}} - 1}}\left( {{\nu _{i3}} + 3\nu ^2_{i2}\frac{{M - 1}}{{{\nu _{i1}} - 1}}} \right) \text{; } $$

(20)

$$ \begin{split} h''''(z)=\;& \frac{{P}_{if}(1-{P}_{if})U''''(h(z)){(h'(z))}^{3}}{{\left[1-{P}_{if}U'(h(z))\right]}^{2}}+\frac{2{P}_{if}(1-{P}_{if})U'''(h(z))h'(z)h''(z)}{{\left[1-{P}_{if}U'(h(z))\right]}^{2}}\\ &+\frac{2{({P}_{if})}^{2}(1-{P}_{if})U'''(h(z)){(h'(z))}^{3}U''(h(z))}{{\left[1-{P}_{if}U'(h(z))\right]}^{3}} +\frac{{P}_{if}(1-{P}_{if})U'''(h(z))h'(z)h''(z)}{{\left[1-{P}_{if}U'(h(z))\right]}^{2}}\\ &+\frac{{P}_{if}(1-{P}_{if})U'''(h(z))h'''(z)}{{\left[1-{P}_{if}U'(h(z))\right]}^{2}}+\frac{2{({P}_{if})}^{2}(1-{P}_{if}){(U''(h(z)))}^{2}h'(z)h''(z)}{{\left[1-{P}_{if}U'(h(z))\right]}^{3}}\\ &+\frac{4{({P}_{if})}^{2}(1-{P}_{if})U''(h(z))U'''(h(z)){(h'(z))}^{3}}{{\left[1-{P}_{if}U'(h(z))\right]}^{3}} +\frac{4{({P}_{if})}^{2}(1-{P}_{if}){(U''(h(z)))}^{2}h'(z)h''(z)}{{\left[1-{P}_{if}U'(h(z))\right]}^{3}}\\ &+\frac{6{({P}_{if})}^{3}(1-{P}_{if}){(U''(h(z)))}^{3}{(h'(z))}^{3}}{{\left[1-{P}_{if}U'(h(z))\right]}^{4}},\end{split} $$

$$ \begin{split} h''''(1) = \;&{M^4}\frac{{M - 1}}{{{\nu _{i1}} - 1}}\left[ {\nu _{i4}} + 10{\nu _{i2}}{\nu _{i3}}\frac{{M - 1}}{{{\nu _{i1}} - 1}} + 15\nu ^3_{i2}{{\left(\frac{{M - 1}}{{{\nu _{i1}} - 1}}\right)}^2} \right] ; \end{split}$$

(21)

$$ H'\left( 1 \right) = \frac{{F{\nu _{sf,1}}}}{{F + {F_\alpha }}}\left( {1 + \alpha } \right)M, $$

(22)

$$ \begin{split}&H''\left( 1 \right) = {M^2}\frac{F}{{F + {F_\alpha }}}\left[ {{\nu _{sf,2}} + {\nu _{sf,1}}\left( {1 + \alpha } \right){\nu _{i2}}\frac{{M - 1}}{{{\nu _{i1}} - 1}}} \right], \end{split}$$

(23)

$$ H'''(1) = {M^3}\frac{F}{{F + {F_\alpha }}}\left[ {{\nu _{sf,3}} + 3{\nu _{sf,2}}{\nu _{i2}}\frac{{M - 1}}{{{\nu _{i1}} - 1}}} + {\nu _{sf,1}}(1 + \alpha){\nu _{i3}}\frac{{M - 1}}{{{\nu _{i1}} - 1}} + 3{\nu _{sf,1}}({1 + \alpha }){\nu _{i2}^2}{{\left( {\frac{{M - 1}}{{{\nu _{i1}} - 1}}} \right)}^2} \right], $$

(24)

$$ \begin{split} H''''\left( 1 \right) =\;& {M^4}\frac{F}{{F + {F_\alpha }}}\Bigg\{ {\nu _{sf,4}} + 6{\nu _{sf,3}}{\nu _{i2}}\frac{{M - 1}}{{{\nu _{i1}} - 1}} + 4{\nu _{sf,2}}{\nu _{i3}}\frac{{M - 1}}{{{\nu _{i1}} - 1}} + 15{\nu _{sf,2}}\nu ^2_{i2}{{\left(\frac{{M - 1}}{{{\nu _{i1}} - 1}}\right)}^2} \\ &+ {\nu _{sf,1}}\left( {1 + \alpha } \right)\frac{{M - 1}}{{{\nu _{i1}} - 1}}\left[ {\nu _{i4}} + 10{\nu _{i2}}{\nu _{i3}}\frac{{M - 1}}{{{\nu _{i1}} - 1}} + 15 \nu ^3 _{i2}{{\left(\frac{{M - 1}}{{{\nu _{i1}} - 1}}\right)}^2} \right] \Bigg\}. \\[-15pt] \end{split} $$

(25)

4)真符合重数的总体矩

假定同属一个源事件的中子中有n个被探测到, 在其中随机选取一个作为触发信号, 下面推导这个触发信号的真符合计数重数的概率分布$p\left( {i|n} \right)$.

假设源事件发生在$t = 0$时刻, 探测到的中子信号时间分布为$f\left( t \right)$, 计数门宽为G, 预延迟为PD. 在$\left( {t, t + {\text{d}}t} \right)$内能够探测到1个中子的概率为

$$ p{(t|n)_{|n = 1}} = p(t|1) = f\left( t \right){\text{d}}t, $$

(26)

在$\left( {t, t + {\text{d}}t} \right)$内能够探测到n个中子的概率为

$$ p(t|n) = nf\left( t \right){\text{d}}t. $$

(27)

以其作为触发信号的计数门为从$t + PD$时刻到$t + PD + G$时刻, 其余$n - 1$个中子信号中的任意1个处于此计数门中(记为事件A, 对此次触发的计数重数产生贡献)的概率为

$$ p\left( {A|n,t} \right) = p\left( {A|t} \right) = \int_{t + PD}^{t + PD + G} {f\left( s \right){\text{d}}s} . $$

(28)

所以$\left( {t, t + {\text{d}}t} \right)$内探测到1个中子, 并且余下的$n - 1$个中子中在计数门中有i个的可能性为

$$ p\left( {i,t|n} \right) = p\left( {t|1} \right)\left( {\begin{array}{*{20}{c}} {n - 1} \\ i \end{array}} \right){p^i}\left( {A|t} \right){\left[ {1 - p\left( {A|t} \right)} \right]^{n - 1 - i}}. $$

(29)

在整个时间范围对t做积分, 即得到随机选取n个信号中的一个作为触发信号, 所记重数为i (计数门中有i个中子信号)的概率为

$$ \begin{split} p\left( {i|n} \right)= \;&\int_0^\infty f(t)\left( {\begin{array}{*{20}{c}} {n - 1} \\ i \end{array}} \right){p^i}( {A|t}) \\ & \times {{\left[ {1 - p( {A|t})} \right]}^{n - 1 - i}}{\text{d}}t . \end{split}$$

(30)

每个中子被探测器探测到的概率称为探测器的探测效率, 记为ε, 则一个样品发射中子产生信号数的概率母函数为(考虑散射串扰, $ k $为散射串扰因子):

$$ e(z) = (1 - \varepsilon ) + \varepsilon (1 - k)z + \varepsilon k{z^2} . $$

(31)

一个源事件产生信号数的概率母函数为

$$ d(z) = H\left( {e(z)} \right) . $$

(32)

一个初级时间探测中子重数的前三阶阶乘矩为

$$ {d'}(1) = \varepsilon (1 + k){H'}(1), $$

(33)

$$ {d{''}}(1) = {\varepsilon ^2}{(1 + k)^2}{H{''}}(1) + 2\varepsilon k{H'}(1), $$

(34)

$$ {d{'''}}(1) = {\varepsilon ^3}{(1 + k)^3}{H{'''}}(1) + 6{\varepsilon ^2}k(1 + k){H{''}}(1), $$

(35)

$$ \begin{split} d'''(1) =\;& H''''(1){\varepsilon ^4}{\left( {1 + \kappa } \right)^4} + 12H'''(1){\varepsilon ^3}\kappa {\left( {1 + \kappa } \right)^2} \\ & +12H''(1){\varepsilon ^2}{\kappa ^2} . \\[-10pt]\end{split}$$

(36)

记一次源事件探测到的中子数为n的可能性为$D\left( n \right)$, 则对于每个探测到的中子都会触发一次计数, 则某一个触发中子属于探测到n个中子的源事件的可能性为

$$ \tilde D\left( n \right) = \frac{{nD\left( n \right)}}{{d'(1)}}, $$

(37)

所以, 某次计数重数为i的概率为

$$ r\left( i \right) = \sum\limits_{n = i}^\infty {\tilde D(n)p\left( {i|n} \right)}, $$

(38)

一、二、三阶矩为

$$ {r_1} = \frac{{d''(1)}}{{d'(1)}}\int_0^\infty {\left[ {f\left( t \right)\int_{t + PD}^{t + PD + G} {f\left( s \right){\text{d}}s} } \right]{\text{d}}t} = \frac{{\dfrac{1}{2}{\varepsilon ^2}{{\left( {1 + \kappa } \right)}^2}H''(1) + \varepsilon \kappa H'(1)}}{{d'(1)}}{f_{\rm d}}, $$

(39)

$$ {r_2} = \frac{{d'''(1)}}{{d'(1)}}\int_0^\infty {f\left( t \right){{\left[ {\int_{t + PD}^{t + PD + G} {f\left( s \right){\text{d}}s} } \right]}^2}{\text{d}}t} = \dfrac{{\dfrac{1}{3}{\varepsilon ^3}{{\left( {1 + \kappa } \right)}^3}H'''(1) + 2{\varepsilon ^2}\kappa (1 + \kappa )H''(1)}}{{d'(1)}}{f_{\rm t}}, $$

(40)

$$ {r_3} = \frac{{d''''(1)}}{{d'(1)}}\int_0^\infty {f\left( t \right){{\left[ {\int_{t + PD}^{t + PD + G} {f\left( s \right){\text{d}}s} } \right]}^3}{\text{d}}t} =\dfrac{{\dfrac{1}{4}{\varepsilon ^4}{{\left( {1 + \kappa } \right)}^4}H''''(1) + 3{\varepsilon ^3}\kappa {{(1 + \kappa )}^2}H'''(1) + 3{\varepsilon ^2}{\kappa ^2}H''(1)}}{{d'(1)}}{f_{\rm q}}. $$

(41)

$ {f_{\text{d}}}, {f_{\text{t}}}, {f_{\text{q}}} $分别是探测器的二重、三重和四重符合门因子:

$$ {f_{\text{d}}} = 2\int_0^\infty {\left[ {f\left( t \right)\int_{t + PD}^{t + PD + G} {f\left( s \right){\text{d}}s} } \right]{\text{d}}t}, $$

(42)

$$ {f_{\text{t}}} = 3\int_0^\infty {f\left( t \right){{\left[ {\int_{t + PD}^{t + PD + G} {f\left( s \right){\text{d}}s} } \right]}^2}{\text{d}}t}, $$

(43)

$$ {f_{\text{q}}} = 4\int_0^\infty {f\left( t \right){{\left[ {\int_{t + PD}^{t + PD + G} {f\left( s \right){\text{d}}s} } \right]}^3}{\text{d}}t} . $$

(44)

根据基本假设的第5条, 探测到的中子信号时间分布为

$$ f\left( t \right) = \frac{1}{\tau }{{\text{e}}^{ - t/\tau }},\quad t > 0, $$

(45)

所以

$$ {f_{\text{d}}} = {{\text{e}}^{ - PD/\tau }}\big( {1 - {{\text{e}}^{ - G/\tau }}} \big), $$

(46)

$$ {f_{\text{t}}} = {\big[ {{{\text{e}}^{ - PD/\tau }}\big( {1 - {{\text{e}}^{ - G/\tau }}} \big)} \big]^2}, $$

(47)

$$ {f_{\text{q}}} = {\big[ {{{\text{e}}^{ - PD/\tau }}\big( {1 - {{\text{e}}^{ - G/\tau }}} \big)} \big]^3}, $$

(48)

在快中子探测中由于中子衰减时间不满足指数关系, 所以$ {f_{\text{d}}} = {f_{\text{t}}} = {f_{\text{q}}} = 1 $.

5)测量计数与真符合重数的关系

探测到的中子既包含源事件诱发裂变产生的中子, 又包含与源事件无关的本底中子. 前景计数器以探测到的中子为触发, 总触发率为

$$ S = F\varepsilon {H'}(1) + {S_{bkg}}, $$

(49)

$$ {S_{bkg}} = {S_0} + B + {S_s} . $$

(50)

前景重数的阶乘矩为

$$ {f_k} = \frac{{F\varepsilon {H'}(1)}}{S}{g_k} + \left[ {1 - \frac{{F\varepsilon {H'}(1)}}{S}} \right]{b_k}, $$

(51)

其中b_k为本底中子触发重数, 即偶然符合重数, 与背景重数分布相同. g_k为源事件中子触发重数, 包含真符合重数r_k与偶然符合重数b_k:

$$ {g_k} = \sum\limits_{j = 0}^k {\left( {\begin{array}{*{20}{c}} k \\ j \end{array}} \right){r_j}{b_{k - j}}} \text{, } $$

(52)

$$ {f_k} - {b_k} = \frac{{F\varepsilon {H'}(1)}}{S}\bigg[ {\sum\limits_{j = 0}^k {\left( \begin{array}{l} k \hfill \\ j \hfill \\ \end{array} \right){r_j}{b_{k - j}} - {b_k}} } \bigg] . $$

(53)

对于$ k = 1, 2, 3 $使用${f_k}, {b_k}$的显式表达, 可以获得真符合重数的一、二阶矩:

$$ {r_1} = \frac{{S({f_1} - {b_1})}}{{F{d'}(1)}}, $$

(54)

$$ {r_2} = \frac{S}{{F{d'}(1)}}\left[ {({f_2} - {b_2}) - 2{b_1}({f_1} - {b_1})} \right], $$

(55)

$$ \begin{split} r =\;& \frac{S}{{F{d'}(1)}}[({f_3} - {b_3}) - 3{b_1}({f_2} - {b_2}) \\ &- 3{b_2}({f_1} - {b_1}) + 6{b_1}{({f_1} - {b_1})^2}] . \end{split}$$

(56)

令真符合重数的总体矩((39)—(41)式)与样本矩((54)—(56)式)相等得到

$$ S - {S_{bkg}} = F\varepsilon (1 + k){H'}(1) \text{, } $$

(57)

$$ S({f_1} - {b_1}) = F{f_d}\Big[\frac{{{\varepsilon ^2}{{(1 + k)}^2}{H{''}}(1)}}{2} + \varepsilon k{H'}(1)\Big] \text{, } $$

(58)

$$\begin{split} &\frac{{S[({f_2} - {b_2}) - 2{b_1}({f_1} - {b_1})]}}{2} \\ =\;& F{f_t}\Big[\frac{{{\varepsilon ^3}{{(1 + k)}^3}{H{'''}}(1)}}{6} + {\varepsilon ^2}k(1 + k){H{''}}(1)\Big], \end{split}$$

(59)

$$ \begin{split} &\frac{{S[({f_3} - {b_3}) - 3{b_1}({f_2} - {b_2}) - 3{b_2}({f_1} - {b_1}) + 6{b_1}{{({f_1} - {b_1})}^2}]}}{6} \\=\;& F{f_q}\left[\frac{{{\varepsilon ^4}{{(1 + k)}^4}{H{''''}}(1)}}{{24}} + \frac{{{\varepsilon ^3}k{{(1 + k)}^2}{H{'''}}(1)}}{2} + \frac{{{\varepsilon ^2}{k^2}{H{''}}(1)}}{2}\right]. \end{split} $$

(60)

(57)—(60)式的左端就是需要的一、二、三、四重记数Single, Double, Tribple和Quarts, 为了书写方便, 简单记为S, D, T, Q. 结合(22)—(25)式可得

$$ S = {S_0} + B + {S_s} + F\varepsilon \left( {1 + \kappa } \right){\nu _{s1}}M\text{, } $$

(61)

$$ D = \frac{{F{\varepsilon ^2}{{\left( {1 + \kappa } \right)}^2}{f_{\rm d}}{M^2}}}{2}\left[ {{\nu _{s2}} + \left( {\frac{{M - 1}}{{{\nu _{i1}} - 1}}} \right){\nu _{s1}}{\nu _{i2}}} \right] + F\varepsilon \kappa {f_{\rm d}}{\nu _{s1}}M , $$

(62)

$$ \begin{split} T =\;& \frac{{F{\varepsilon ^3}{{\left( {1 + \kappa } \right)}^3}{f_{\rm t}}{M^3}}}{6}\left[ {{\nu _{s3}} + \left( {\frac{{M - 1}}{{{\nu _{i1}} - 1}}} \right)3{\nu _{s2}}{\nu _{i2}}} \right. + \left( {\frac{{M - 1}}{{{\nu _{i1}} - 1}}} \right){\nu _{s1}}{\nu _{i3}} \\ & \left. { + 3{{\left( {\frac{{M - 1}}{{{\nu _{i1}} - 1}}} \right)}^2}{\nu _{s1}}\nu _{i2}^2} \right] + F{\varepsilon ^2}\kappa \left( {1 + \kappa } \right){f_{\rm t}}{M^2}\left[ {{\nu _{s2}} + \left( {\frac{{M - 1}}{{{\nu _{i1}} - 1}}} \right){\nu _{s1}}{\nu _{i2}}} \right], \end{split} $$

(63)

$$ \begin{split} Q = \;&\frac{{F{\varepsilon ^4}{{\left( {1 + \kappa } \right)}^4}{f_{\rm q}}{M^4}}}{{24}}\bigg\{ {{\nu _{s4}} + 6{\nu _{s3}}{\nu _{i2}}\frac{{M - 1}}{{{\nu _{i1}} - 1}} + 4{\nu _{s2}}{\nu _{i3}}\frac{{M - 1}}{{{\nu _{i1}} - 1}} } +15{\nu _{s2}} \nu ^2_{i2}{\left(\frac{{M - 1}}{{{\nu _{i1}} - 1}}\right)^2} \\ & + {{\nu _{s1}}\frac{{M - 1}}{{{\nu _{i1}} - 1}}\left[ {{\nu _{i4}} + 10{\nu _{i2}}{\nu _{i3}}\frac{{M - 1}}{{{\nu _{i1}} - 1}} + 15 \nu ^3_{i2}{{\left(\frac{{M - 1}}{{{\nu _{i1}} - 1}}\right)}^2}} \right]} \bigg\} \frac{{F{\varepsilon ^3}{{\left( {1 + \kappa } \right)}^2}k{f_{\rm q}}{M^3}}}{2}\bigg[ {{\nu _{s3}} + \left( {\frac{{M - 1}}{{{\nu _{i1}} - 1}}} \right)3{\nu _{s2}}{\nu _{i2}}} \\ & { + \left( {\frac{{M - 1}}{{{\nu _{i1}} - 1}}} \right){\nu _{s1}}{\nu _{i3}} + 3{{\left( {\frac{{M - 1}}{{{\nu _{i1}} - 1}}} \right)}^2}{\nu _{s1}}\nu _{i2}^2} \bigg] + \frac{{F{\varepsilon ^2}{k^2}{f_{\rm q}}{M^2}}}{2}\left[ {{\nu _{s2}} + \left( {\frac{{M - 1}}{{{\nu _{i1}} - 1}}} \right){\nu _{s1}}{\nu _{i2}}} \right], \\[-15pt] \end{split} $$

(64)

式中, $ {v_{s1}}, {v_{s2}}, {v_{s3}}, {v_{s4}} $分别是外部中子源诱发裂变发射中子数分布的一、二、三、四阶阶乘矩, ${v_{i1}}, $$ {v_{i2}}, {v_{i3}}, {v_{i4}}$分别是次级诱发裂变发射中子数分布的一、二、三、四阶阶乘矩, ε是探测器的探测效率, $ {f_{\text{d}}}, {f_{\text{t}}}, {f_{\text{q}}} $分别是探测器二重、三重、四重符合门因子. 由于主动中子测量是利用外部中子诱发铀材料裂变, 所以测得的诱发裂变率与铀材料质量、外部中子源强度以及耦合系数有关, 如下式所示:

$$ F = CmY . $$

(65)

耦合系数$ C $是有源多重性测量过程中的关键参数, 与样品的几何结构、密度、材料组成等有关, 难以通过理论计算的方法求出. Krick研究团队^[3]揭示了耦合系数和增殖系数之间的经验关系式:

$$ C = a - \frac{{b(M - 1)}}{{1 + c(M - 1)}} . $$

(66)

方程的求解过程为利用测得的二重计数率$ D $和三重计数率$ T $迭代消元求得增殖系数$ M $, 再将$ M $代入方程求解出诱发裂变率$ F $, 将已知的参数$ m, F, Y $代入公式求出耦合系数$ C $, 通过一系列的$ C $和$ M $拟合得到刻度系数$ a, b, c $, 进而得到耦合系数与增殖系数的具体表达式. 根据确定好的耦合-增殖关系曲线计算出耦合系数的值, 再由(65)式解出²³⁵U的有效质量.

4.   验证及分析

为对上述推导铀材料快中子多重性测量方程进行验证, 利用上述搭建的探测系统进行模拟测量. 设置一系列质量的样品, 测量结果如图4所示.

图4中子多重计数率的测量结果与文献[14]相似, 基本维持在同一量级内. 增殖系数分布见图5, 增殖系数会随着样品质量的增大而增大, 与文献[11]中的拟合曲线基本吻合, 说明了测量结果的准确性.

通过测得的裂变率$ F $、模拟测量时预设的样品质量以及中子源强度$ Y $可以发解出对应的耦合系数$ C $. 将耦合系数和增殖系数进行拟合(见图6), 得到拟合方程为$ C = a/(M + b) $, 系数$ a $ = 0.6173, $ b $ = –1.196, 拟合优度达0.9923, 反映了拟合曲线的高度吻合. 对耦合系数和增殖系数之间的经验公式(66)进行化简, 得到如下公式:

$$ C = a - \frac{b}{c} + \frac{{b/{c^2}}}{{M + (1 - c)/c}} . $$

(67)

通过对化简方程与拟合方程进行比较, 发现二者的函数形式高度相似, 说明了拟合方程的可靠性. 将部分测量数据代入拟合方程进行检验, 通过多重计数率的迭代消元, 得到增殖系数$ M $, 代入拟合方程求得耦合系数$ C $, 利用$ m, F, Y, C $的函数关系解得样品质量$ m $, 与理论值对比, 得到测量偏差如图7所示.

在2 kg范围内, 样品的质量求解偏差绝大多数保持在20%以内, 且小质量样品测量偏差偏小, 随着样品质量的增大, 求解偏差呈现扩大趋势. 分析原因主要有两部分: 一是多重性测量方程是建立在点模型的基础上的, 而在实际测量中放射源往往是以体源的形式呈现的. 二是耦合系数与增殖系数拟合曲线的偏差也会对测量结果造成一定的影响.

5.   总　结

本文在NMC分析方程的推导基础上, 采用概率母函数完成对铀材料快中子多重性测量方程进行推导. 根据铀钚材料物理过程的差异, 在铀材料的快中子多重性测量中忽略了$ (\alpha, n) $反应系数$ \alpha $, 并引入了常数$ k $修正快中子的散射串扰影响, 推导出了铀材料快中子多重性测量方程.

为检验方程的准确性, 分别利用Geant4和Matlab实现中子探测和中子计数的模拟过程, 得到不同质量情况下样品特征参数. 对耦合系数与增殖系数进行拟合, 得到方程$ C = a/(M + b) $, 与经验方程$C = a - \dfrac{{b(M - 1)}}{{1 + c(M - 1)}}$高度吻合. 将拟合方程代入计算, 得到2 kg范围内质量求解偏差小于20%, 证实了测量方程的可靠性. 此项研究建立了铀材料快中子多重性测量方程, 实现了铀材料的质量属性认证, 拓展了快中子多重性的应用范围, 对军控核查技术的发展具有重要意义.