一種核材料識別系統的Monte Carlo方法仿真研究

譙 梁 魏 彪 楊 帆 馮 鵬 周 密

（重慶大學光電技術及系統教育部重點實驗室 重慶 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<>a時，成反比開始衰減。

由圖 7(c)和式(6)可知，當w<>a時，反比衰減。

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