一種核材料識別系統(tǒng)的Monte Carlo方法仿真研究

譙 梁 魏 彪 楊 帆 馮 鵬 周 密

（重慶大學光電技術及系統(tǒng)教育部重點實驗室 重慶 400044）

核材料識別系統(tǒng)(Nuclear Material Identification System，NMIS)是源于主動噪聲分析法測量原理的一種核查技術及方法[1]。它基于252Cf自發(fā)裂變中子源驅(qū)動，通過對核材料或核部件構成的核系統(tǒng)中的誘發(fā)中子和γ射線脈沖信號采集、分析與測量，獲得未知核部件特定屬性的時-頻域特征標簽，以此對核材料分析與識別[2?4]。NMIS主要目的是鑒定含有高富集度鈾(Highly Enriched Uranium，HEU)的核武器及核材料，國外目前已應用于核軍控、核裁軍及各種高富集度核材料的識別工作中[5]，在國內(nèi)該項研究工作尚處起步階段[6]。

鑒于核材料使用的敏感性及強放射性輻射的危害性等因素，致使對核材料的分析與識別研究工作甚為棘手，仿真研究甚為重要。為此，本文基于Monte Carlo方法開展NMIS的理論仿真（模擬）研究，以此建立未知核材料的富集度、反應性與特征標簽之間的聯(lián)系，為核材料識別系統(tǒng)的研制提供保證，并對核軍控核查技術的實驗室仿真研究提供一種新途徑。

1 系統(tǒng)結構及仿真模型

1.1 仿真模型的構造

Monte Carlo方法（又稱隨機抽樣技巧或統(tǒng)計實驗方法）是一種數(shù)值分析方法，它是數(shù)理統(tǒng)計與計算機相結合的產(chǎn)物[7]。它能有效地描述隨機事件并模擬實際點的實驗過程，特別適用于有隨機性的粒子輸運問題的模擬計算。MCNP (Monte Carlo Neutron and Photon Transport Code)是美國洛斯阿拉莫斯國家實驗室基于Monte Carlo思想開發(fā)的一套模擬核實驗的通用程序，用以解決中子、光子、電子輻射輸運問題的仿真問題[8,9]。MCNP不直接求解輸運方程，而是從物理原理出發(fā)，通過模擬大量粒子行為并記錄統(tǒng)計學特征得到所需的結果。它有較強的幾何建模能力和大量實時更新的反應截面數(shù)據(jù)庫，廣泛用于核物理特別是高能物理研究領域[10,11]。本文采用MCNP進行粒子輸運過程模擬，得到一系列仿真數(shù)據(jù)。

仿真模型由自發(fā)裂變中子源、未知核部件和閃爍探測器組成。自發(fā)裂變中子源為252Cf（锎），裂變率為429800 Hz/mg，每次裂變平均放射出~4個中子和6個γ光子。仿真研究構造中采用的中子源質(zhì)量為0.7 mg，故裂變率為3×105Hz。252Cf源距離鑄件底部有7.62 cm，與鑄件相連。待測鈾部件高15.24 cm、外徑12.7 cm、內(nèi)徑8.89 cm，放在高22.4 cm、外徑15.6 cm、厚度0.15 cm的不銹鋼罐里。鋼罐外面為U型聚乙烯反射層，一側嵌有兩個BC430塑料閃爍探測器，用于探測中子誘發(fā)裂變粒子脈沖信號；另一側放一個相同塑料閃爍探測器，用于探測中子源的裂變時間信號，其高和寬為7.62 cm、厚為10.16 cm。探測器周圍布滿鉛屏蔽層，用于防止探測器間串擾。仿真模型構造如圖1。

如圖1，當252Cf中子源發(fā)生自發(fā)裂變后，產(chǎn)生的中子一部分進入核材料，產(chǎn)生散射、吸收反應和誘發(fā)裂變，裂變產(chǎn)物包含中子、γ光子等被探測器探測到，并在時域上形成一定分布[5]。該仿真模型為桶狀鑄件模型，包含235U、234U、236U和238U，各項同位素材料富集度如表1。

圖1 NMIS仿真模型結構示意圖 (a) 橫剖面圖；(b) 斷面圖1、3、5內(nèi)為空氣，2為鈾材料鑄件，4為鋼罐，6和7為兩路探測器通道，8為鉛屏蔽層，9為聚乙烯反射層，10為系統(tǒng)外空間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輸入文件中把所有仿真模型組件和反射層描述成許多柵元，每個柵元由一個或多個曲面（或平面）圍成，柵元內(nèi)填充材料。所有柵元都在柵元卡(Cell Card)中列出，表面卡(Surface Card)則列出全部平面和曲面，信息卡(Information Card)列出所用全部材料、源信息及計數(shù)等。Monte Carlo方法的抽樣與粒子密度成正比，因此計算中抽樣問題不會構成障礙，可不予考慮。故堆芯中各柵元的重要性 IMP(Importance)均設為1，粒子進入系統(tǒng)外空間則自動湮滅。

1.3 仿真結果

仿真實驗中分別將五種不同富集度的鈾材料置于本文構造的NMIS模型中，通過外圍兩個閃爍探測器采集的誘發(fā)裂變中子時域分布如圖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源發(fā)射的直射及散射中子，這部分粒子受材料富集度影響較?。?0 ns以后曲線隨235U富集度增加而有所不同，富集度越高待測核材料受誘發(fā)裂變概率就越大，從而產(chǎn)生更大量的誘發(fā)裂變中子。特別是15?40 ns內(nèi)，曲線積分值與待測核材料的富集度直接相關，這也是NMIS進行核材料富集度識別的基礎。

仿真系統(tǒng)采樣頻率為1 GHz，采樣時間間隔為1 ns。因為直射γ射線以光速向前傳播，γ峰出現(xiàn)在~1 ns，隨后到達的是散射γ射線和未經(jīng)“碰撞”的直射中子，如圖 3。中子能量不同，飛行速度也就不同，因此中子峰在時域上存在較大展寬，最后到達探測器的是散射中子、誘發(fā)裂變中子和γ射線。富集度為93.15wt%的鈾材料在NMIS中經(jīng)252Cf源激發(fā)后，其時-頻特征標簽如圖3所示。

2 NMIS時-頻特征分析

用Monte Carlo方法模擬仿真NMIS工作流程，首先獲得源通道（設置為1通道）和被探測的2個通道（分別設置為Ι通道和?Ι通道）的脈沖信號，并由此得到自發(fā)裂變中子源（如252Cf中子源）及其誘發(fā)的粒子探測計數(shù)的時域分布。隨后對3路探測器采集到的數(shù)據(jù)按一定長度N進行分塊（實驗中設置N=1024），分別自相關(Auto-Correlation, AC)、互相關(Cross-Correlation, CC)和自功率譜(Auto Power Spectral Densities, APSD)、互功率譜(Cross Power Spectral Densities, CPSD)等一系列時-頻域分析運算[12?14]，進而得到能體現(xiàn)包括富集度、反應性等核材料特性的時-頻特征標簽。

2.1 自相關

在仿真模型中探測器的自相關函數(shù)可表示為：

式中，i=1時X1表示源探測器的時域計數(shù)，i=2、3時，X2或X3表示探測器Ι或探測器?Ι的時域計數(shù)。t表示被積函數(shù)的時延。圖4是三個探測器的自相關函數(shù)圖譜。

圖4 三路探測器信號自相關 (a) 源探測器；(b) 探測器Ι；(c) 探測器?ΙFig.4 The auto correlation of source detector (a), detector Ι (b) and detector ?Ι (c).

從圖4看出，三路探測器在t=0 ns時均出現(xiàn)峰值。圖4(a)中，252Cf源自相關的結果關于t=0軸對稱，t=0 ns時的峰值表示252Cf源的裂變率，其值與實驗設置相比同為3×105Hz；圖4(b)中，t=0 ns時的值是探測器 Ι探測到的252Cf源的計數(shù)率；圖4(c)中，t=0 ns時的值是探測器?Ι探測到的252Cf源的計數(shù)率。仿真研究在理想環(huán)境下進行，因此圖譜中扣除了本底輻照的影響。

2.2 互相關

探測器間的互相關函數(shù)表示為：

式中，當t<0時，源探測器和探測器Ι、?Ι間的互相關計數(shù)被舍棄，是因為探測器事件依賴于裂變鏈與發(fā)生在其之前的源裂變事件，發(fā)生在源裂變之前的探測器事件并無實際意義。三個通道互相關運算所得結果如圖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中子源通道的互相關函數(shù)與脈沖中子測量的時域分布（見圖 3）幾近相同，均反映了源通道信號與探測器通道信號間關于兩者時延t的共同信息。顯然，由于環(huán)境噪聲、系統(tǒng)噪聲乃至電子學響應時間的影響，γ峰和中子峰往往存在時域交疊的情況。

由于探測器與待測核材料間的相對距離差異，互相關函數(shù)關于時延0點不對稱，這在圖5(c)中得到一定體現(xiàn)。但由于γ射線的存在，函數(shù)圖像在t=0 ns時存在一個峰值。此外，不同于中子間明顯的差異，γ光子的飛行速度導致飛行時間差無法在互相關函數(shù)上體現(xiàn)。

2.3 自功率譜密度

對自相關函數(shù)做快速傅立葉變換(Fast Fourier Transform, FFT)得自功率譜密度，其表達式為：

分別對三個通道的自相關函數(shù)做快速傅立葉變換，所得自功率譜如圖6。

圖6 三路探測器信號自功率譜 (a) 源探測器；(b) 探測器Ι；(c) 探測器?ΙFig.6 The auto power spectral density of source detector (a), detector Ι (b) and detector ?Ι (c).

圖 6(a)是252Cf源的自功率譜密度隨頻率的分布，是一個常數(shù)，不隨頻率變化，只與252Cf裂變的計數(shù)率有關，它可監(jiān)控源探測系統(tǒng)的運行狀況；圖6(b)和(c)可看出探測器?和?Ι的自功率譜有些波動，由此可知反應堆的中子增殖因子和探測效率都很大，它們對探測器的計數(shù)率很敏感，故所有對計數(shù)率有影響的因素，如自發(fā)裂變、誘發(fā)裂變、本底、系統(tǒng)噪聲乃至電子學響應時間等，都將對探測器的自功率譜有影響。

2.4 互功率譜密度

互功率譜密度是兩個信號間共有信息的量度，互功率譜函數(shù)可表示為：

源-探測器的互功率譜函數(shù)：

探測器?-探測器??的互功率譜函數(shù)：

式中，e2與e3分別表示探測器?和??的探測效率，是每次源裂變發(fā)射的平均中子數(shù)，F(xiàn)0是平均裂變率，L是瞬發(fā)中子代時間，a是瞬發(fā)中子裂變鏈的衰變常數(shù)。

分別對三個通道相互之間的互相關做 FFT變換，得互功率譜函數(shù)的幅度譜如圖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裂變中子源激發(fā)待測核材料的過程進行模擬仿真，獲得五種不同富集度待測核材料的脈沖中子信號，并以富集度為93.15wt%的HEU為例，對其三路核信號進行時-頻運算，得到不同種類的時-頻特征標簽，進而對其特征參數(shù)進行分析，探討了NMIS系統(tǒng)的動態(tài)特性及物理意義。

仿真研究結果表明，本文構造的仿真模型、過程仿真及數(shù)據(jù)計算能較好地反映核裂變中粒子輸運過程，可呈現(xiàn)NMIS系統(tǒng)中若干特征標簽的有效性，為后續(xù)實際條件下相關實驗的開展及特征譜線的解讀和分析奠定了基礎，并對核軍控核查技術的實驗室仿真研究有積極意義。

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裂變中子源的核信息系統(tǒng)頻譜分析[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中子源慢化及屏蔽系統(tǒng)的蒙特卡羅模擬研究[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]. 北京: 電子工業(yè)出版社, 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