一种核材料识别系统的Monte Carlo方法仿真研究
2012-09-23谯梁魏彪杨帆冯鹏周密
谯 梁 魏 彪 杨 帆 冯 鹏 周 密
（重庆大学光电技术及系统教育部重点实验室 重庆 400044）
核材料识别系统(Nuclear Material Identification System，NMIS)是源于主动噪声分析法测量原理的一种核查技术及方法[1]。它基于252Cf自发裂变中子源驱动，通过对核材料或核部件构成的核系统中的诱发中子和γ射线脉冲信号采集、分析与测量，获得未知核部件特定属性的时-频域特征标签，以此对核材料分析与识别[2−4]。NMIS主要目的是鉴定含有高富集度铀(Highly Enriched Uranium，HEU)的核武器及核材料，国外目前已应用于核军控、核裁军及各种高富集度核材料的识别工作中[5]，在国内该项研究工作尚处起步阶段[6]。
鉴于核材料使用的敏感性及强放射性辐射的危害性等因素，致使对核材料的分析与识别研究工作甚为棘手，仿真研究甚为重要。为此，本文基于Monte Carlo方法开展NMIS的理论仿真（模拟）研究，以此建立未知核材料的富集度、反应性与特征标签之间的联系，为核材料识别系统的研制提供保证，并对核军控核查技术的实验室仿真研究提供一种新途径。
1 系统结构及仿真模型
1.1 仿真模型的构造
Monte Carlo方法（又称随机抽样技巧或统计实验方法）是一种数值分析方法，它是数理统计与计算机相结合的产物[7]。它能有效地描述随机事件并模拟实际点的实验过程，特别适用于有随机性的粒子输运问题的模拟计算。MCNP (Monte Carlo Neutron and Photon Transport Code)是美国洛斯阿拉莫斯国家实验室基于Monte Carlo思想开发的一套模拟核实验的通用程序，用以解决中子、光子、电子辐射输运问题的仿真问题[8,9]。MCNP不直接求解输运方程，而是从物理原理出发，通过模拟大量粒子行为并记录统计学特征得到所需的结果。它有较强的几何建模能力和大量实时更新的反应截面数据库，广泛用于核物理特别是高能物理研究领域[10,11]。本文采用MCNP进行粒子输运过程模拟，得到一系列仿真数据。
仿真模型由自发裂变中子源、未知核部件和闪烁探测器组成。自发裂变中子源为252Cf（锎），裂变率为429800 Hz/mg，每次裂变平均放射出~4个中子和6个γ光子。仿真研究构造中采用的中子源质量为0.7 mg，故裂变率为3×105Hz。252Cf源距离铸件底部有7.62 cm，与铸件相连。待测铀部件高15.24 cm、外径12.7 cm、内径8.89 cm，放在高22.4 cm、外径15.6 cm、厚度0.15 cm的不锈钢罐里。钢罐外面为U型聚乙烯反射层，一侧嵌有两个BC430塑料闪烁探测器，用于探测中子诱发裂变粒子脉冲信号；另一侧放一个相同塑料闪烁探测器，用于探测中子源的裂变时间信号，其高和宽为7.62 cm、厚为10.16 cm。探测器周围布满铅屏蔽层，用于防止探测器间串扰。仿真模型构造如图1。
如图1，当252Cf中子源发生自发裂变后，产生的中子一部分进入核材料，产生散射、吸收反应和诱发裂变，裂变产物包含中子、γ光子等被探测器探测到，并在时域上形成一定分布[5]。该仿真模型为桶状铸件模型，包含235U、234U、236U和238U，各项同位素材料富集度如表1。
图1 NMIS仿真模型结构示意图 (a) 横剖面图；(b) 断面图1、3、5内为空气，2为铀材料铸件，4为钢罐，6和7为两路探测器通道，8为铅屏蔽层，9为聚乙烯反射层，10为系统外空间Fig.1 Schematic diagram of NMIS simulation model. (a) cross section; (b) horizontal section cell 1, 3 and 5 are void, cell 2 is uranium metal casting, cell 4 is steel can, cell 6 and 7 are two detector channels, cell 8 is lead,cell 9 is polyethylene reflector, cell 10 is outside world of the system.
表1 铀金属铸件的富集度一览表Table 1 Varying enrichment of uranium metal castings.
1.2 仿真流程
MCNP程序下对NMIS构造的模型进行仿真研究。为验证所需，在不影响仿真结论的前提下塑料闪烁探测器的探测效率被置为 100%。MCNP输入文件中把所有仿真模型组件和反射层描述成许多栅元，每个栅元由一个或多个曲面（或平面）围成，栅元内填充材料。所有栅元都在栅元卡(Cell Card)中列出，表面卡(Surface Card)则列出全部平面和曲面，信息卡(Information Card)列出所用全部材料、源信息及计数等。Monte Carlo方法的抽样与粒子密度成正比，因此计算中抽样问题不会构成障碍，可不予考虑。故堆芯中各栅元的重要性 IMP(Importance)均设为1，粒子进入系统外空间则自动湮灭。
1.3 仿真结果
仿真实验中分别将五种不同富集度的铀材料置于本文构造的NMIS模型中，通过外围两个闪烁探测器采集的诱发裂变中子时域分布如图2。
图2 五种不同富集度下探测器Ι (a)和探测器ІΙ (b)的中子随时间分布Fig.2 Time distribution of neutron in the detector Ι (a) and ІΙ (b) with five different degree of enrichment.
从图2看出，无论探测器I或II，10 ns以前中子随时间分布曲线几近重合，其原因是此时采集到的都是252Cf源发射的直射及散射中子，这部分粒子受材料富集度影响较小；10 ns以后曲线随235U富集度增加而有所不同，富集度越高待测核材料受诱发裂变概率就越大，从而产生更大量的诱发裂变中子。特别是15−40 ns内，曲线积分值与待测核材料的富集度直接相关，这也是NMIS进行核材料富集度识别的基础。
仿真系统采样频率为1 GHz，采样时间间隔为1 ns。因为直射γ射线以光速向前传播，γ峰出现在~1 ns，随后到达的是散射γ射线和未经“碰撞”的直射中子，如图 3。中子能量不同，飞行速度也就不同，因此中子峰在时域上存在较大展宽，最后到达探测器的是散射中子、诱发裂变中子和γ射线。富集度为93.15wt%的铀材料在NMIS中经252Cf源激发后，其时-频特征标签如图3所示。
2 NMIS时-频特征分析
用Monte Carlo方法模拟仿真NMIS工作流程，首先获得源通道（设置为1通道）和被探测的2个通道（分别设置为Ι通道和ІΙ通道）的脉冲信号，并由此得到自发裂变中子源（如252Cf中子源）及其诱发的粒子探测计数的时域分布。随后对3路探测器采集到的数据按一定长度N进行分块（实验中设置N=1024），分别自相关(Auto-Correlation, AC)、互相关(Cross-Correlation, CC)和自功率谱(Auto Power Spectral Densities, APSD)、互功率谱(Cross Power Spectral Densities, CPSD)等一系列时-频域分析运算[12−14]，进而得到能体现包括富集度、反应性等核材料特性的时-频特征标签。
2.1 自相关
在仿真模型中探测器的自相关函数可表示为：
式中，i=1时X1表示源探测器的时域计数，i=2、3时，X2或X3表示探测器Ι或探测器ІΙ的时域计数。t表示被积函数的时延。图4是三个探测器的自相关函数图谱。
图4 三路探测器信号自相关 (a) 源探测器；(b) 探测器Ι；(c) 探测器ІΙFig.4 The auto correlation of source detector (a), detector Ι (b) and detector ІΙ (c).
从图4看出，三路探测器在t=0 ns时均出现峰值。图4(a)中，252Cf源自相关的结果关于t=0轴对称，t=0 ns时的峰值表示252Cf源的裂变率，其值与实验设置相比同为3×105Hz；图4(b)中，t=0 ns时的值是探测器 Ι探测到的252Cf源的计数率；图4(c)中，t=0 ns时的值是探测器ІΙ探测到的252Cf源的计数率。仿真研究在理想环境下进行，因此图谱中扣除了本底辐照的影响。
2.2 互相关
探测器间的互相关函数表示为：
式中，当t<0时，源探测器和探测器Ι、ІΙ间的互相关计数被舍弃，是因为探测器事件依赖于裂变链与发生在其之前的源裂变事件，发生在源裂变之前的探测器事件并无实际意义。三个通道互相关运算所得结果如图5。
图5 三路探测器信号之间的互相关 (a) 源探测器与探测器Ι；(b) 源探测器与探测器ІΙ；(c) 探测器Ι与探测器ІΙFig.5 The cross correlation of three detectors. (a) detector source and Ι; (b) detector source and ІΙ; (c) detector Ι and ІΙ
从图5(a)、(b)中看出，两个探测器通道和252Cf中子源通道的互相关函数与脉冲中子测量的时域分布（见图 3）几近相同，均反映了源通道信号与探测器通道信号间关于两者时延t的共同信息。显然，由于环境噪声、系统噪声乃至电子学响应时间的影响，γ峰和中子峰往往存在时域交叠的情况。
由于探测器与待测核材料间的相对距离差异，互相关函数关于时延0点不对称，这在图5(c)中得到一定体现。但由于γ射线的存在，函数图像在t=0 ns时存在一个峰值。此外，不同于中子间明显的差异，γ光子的飞行速度导致飞行时间差无法在互相关函数上体现。
2.3 自功率谱密度
对自相关函数做快速傅立叶变换(Fast Fourier Transform, FFT)得自功率谱密度，其表达式为：
分别对三个通道的自相关函数做快速傅立叶变换，所得自功率谱如图6。
图6 三路探测器信号自功率谱 (a) 源探测器；(b) 探测器Ι；(c) 探测器ІΙFig.6 The auto power spectral density of source detector (a), detector Ι (b) and detector ІΙ (c).
图 6(a)是252Cf源的自功率谱密度随频率的分布，是一个常数，不随频率变化，只与252Cf裂变的计数率有关，它可监控源探测系统的运行状况；图6(b)和(c)可看出探测器І和ІΙ的自功率谱有些波动，由此可知反应堆的中子增殖因子和探测效率都很大，它们对探测器的计数率很敏感，故所有对计数率有影响的因素，如自发裂变、诱发裂变、本底、系统噪声乃至电子学响应时间等，都将对探测器的自功率谱有影响。
2.4 互功率谱密度
互功率谱密度是两个信号间共有信息的量度，互功率谱函数可表示为：
源-探测器的互功率谱函数：
探测器І-探测器ІІ的互功率谱函数：
式中，e2与e3分别表示探测器І和ІІ的探测效率，是每次源裂变发射的平均中子数，F0是平均裂变率，L是瞬发中子代时间，a是瞬发中子裂变链的衰变常数。
分别对三个通道相互之间的互相关做 FFT变换，得互功率谱函数的幅度谱如图7所示，它反映了源通道信号和被探测器І、ІІ测量通道信号间的共同信息。
图7 三路探测器信号互功率谱 (a) 源探测器与探测器І；(b) 源探测器与探测器ІІ；(c) 探测器І与探测器ІІFig.7 The cross power spectral density of three detectors. (a) detector source and І; (b) detector source and ІІ; (c) detector І and ІІ.
由图 7(a)、(b)和式(5)可知，当w<
由图 7(c)和式(6)可知，当w<
3 结语
用Monte Carlo方法对252Cf裂变中子源激发待测核材料的过程进行模拟仿真，获得五种不同富集度待测核材料的脉冲中子信号，并以富集度为93.15wt%的HEU为例，对其三路核信号进行时-频运算，得到不同种类的时-频特征标签，进而对其特征参数进行分析，探讨了NMIS系统的动态特性及物理意义。
仿真研究结果表明，本文构造的仿真模型、过程仿真及数据计算能较好地反映核裂变中粒子输运过程，可呈现NMIS系统中若干特征标签的有效性，为后续实际条件下相关实验的开展及特征谱线的解读和分析奠定了基础，并对核军控核查技术的实验室仿真研究有积极意义。
1 Valentine T E, Mattingly J K, Mihalczo J T. NMIS Measurements for Uranium Metal Annular Castings[D].1998
2 Mendel J M. Tutorial on higher-order statistics (spectra)in signal processing and system theory: theoretical results and some applications[J]. Proc IEEE, 1991, 79(3):278−305
3 Nomura T, Gotoh S, Yamaki K. Reactivity Measurements by Neutron Noise Analysis Using Two-Detector Correlation Method and Supercritical Reactor Noise Analysis[M]. Nippon Atomic Industry Group Co., Ltd
4 Valentine T E. MCNP-DSP: A Neutron and Gamma Ray Monte Carlo Calculation of Source-Driven Noise-Measured Parameters[M]. Ph.D. Dissertation, Univ of Tennessee, 1995
5 Mihalczo J T, Mullens J A, Mattigly J K, et al. Physical Description of Nuclear Materials Identification System(NMIS) Signatures[D]. 2000
6 魏彪, 任勇, 冯鹏, 等. 基于252Cf裂变中子源的核信息系统频谱分析[J]. 强激光与粒子束, 2009, 21(12):1898−1902 WEI Biao, REN Yong, FENG Peng, et al. Spectrum analysis of nuclear information system based on252Cf fission neutron source[J]. High Power Laser and Particle Beams, 2009, 21(12): 1898−1902
7 姜校亮.252Cf中子源慢化及屏蔽系统的蒙特卡罗模拟研究[D]. 2010 JIANG Xiaoliang. Study on the Structure of Moderator and Shield of252Cf Neutron Source with the Monte Carlo Method[D]. 2010
8 Briesmeister J F. MCNP General Monte Carlo N-Particle Transport Code Version 5[R]. Radiation Transport Group Los Alamos National Laboratory
9 吴冲, 张强, 孙志嘉, 等. 低能中子探测的GEANT4模拟研究[J]. 核技术, 2012, 29(2): 173−177 WU Chong, ZHANG Qiang, SUN Zhijia, et al.Simulation of low-energy neutron detection based on GEANT4[J]. Nuclear Techniques, 2012, 29(2): 173−177
10 裴鹿成, 张孝泽. 蒙特卡罗方法及其在粒子输运问题中的应用[M]. 北京: 科学出版社, 1980. 1 PEI Lucheng, ZHANG Xiaoze. Monte-Carlo Method and Applications in Transport Problems of Particles[M].Beijing: Science Press, 1980. 1
11 商静, 李为民. 电子束辐照装置屏蔽体预埋管道处的辐射场计算[J]. 核技术, 2011, 34(5): 372−376 SHANG Jing, LI Weimin. Radiation field calculation outside the duct-embedded shielding wall of an E-beam irradiator[J]. Nuclear Techniques, 2011, 34(5): 372−376
12 王永德, 王军. 随机信号分析基础[M]. 北京: 电子工业出版社, 2006 WANG Yongde, WANG Jun. Elements of Analysis for Random Signal[M]. Beijing: Publishing House of Electronics Industry, 2006
13 Perez R B, Munoz-Cobo J L, Valentine T E, et al.Incorporation of neutron and gamma multiplicity measurements in the ornl nuclear materials identification system (NMIS): Passive and Active Measurements[D].2008
14 Peña K, McConchie S, Mihalczo J. Prompt neutron time decay in single HEU and DU metal annular storage castings[D]. 2008
