譙 梁 魏 彪 楊 帆 馮 鵬 周 密

（重慶大學(xué)光電技術(shù)及系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室 重慶 400044）

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

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

1 系統(tǒng)結(jié)構(gòu)及仿真模型

1.1 仿真模型的構(gòu)造

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

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

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

圖1 NMIS仿真模型結(jié)構(gòu)示意圖 (a) 橫剖面圖；(b) 斷面圖1、3、5內(nèi)為空氣，2為鈾材料鑄件，4為鋼罐，6和7為兩路探測(cè)器通道，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程序下對(duì)NMIS構(gòu)造的模型進(jìn)行仿真研究。為驗(yàn)證所需，在不影響仿真結(jié)論的前提下塑料閃爍探測(cè)器的探測(cè)效率被置為 100%。MCNP輸入文件中把所有仿真模型組件和反射層描述成許多柵元，每個(gè)柵元由一個(gè)或多個(gè)曲面（或平面）圍成，柵元內(nèi)填充材料。所有柵元都在柵元卡(Cell Card)中列出，表面卡(Surface Card)則列出全部平面和曲面，信息卡(Information Card)列出所用全部材料、源信息及計(jì)數(shù)等。Monte Carlo方法的抽樣與粒子密度成正比，因此計(jì)算中抽樣問(wèn)題不會(huì)構(gòu)成障礙，可不予考慮。故堆芯中各柵元的重要性 IMP(Importance)均設(shè)為1，粒子進(jìn)入系統(tǒng)外空間則自動(dòng)湮滅。

1.3 仿真結(jié)果

仿真實(shí)驗(yàn)中分別將五種不同富集度的鈾材料置于本文構(gòu)造的NMIS模型中，通過(guò)外圍兩個(gè)閃爍探測(cè)器采集的誘發(fā)裂變中子時(shí)域分布如圖2。

圖2 五種不同富集度下探測(cè)器Ι (a)和探測(cè)器?Ι (b)的中子隨時(shí)間分布Fig.2 Time distribution of neutron in the detector Ι (a) and ?Ι (b) with five different degree of enrichment.

從圖2看出，無(wú)論探測(cè)器I或II，10 ns以前中子隨時(shí)間分布曲線(xiàn)幾近重合，其原因是此時(shí)采集到的都是252Cf源發(fā)射的直射及散射中子，這部分粒子受材料富集度影響較?。?0 ns以后曲線(xiàn)隨235U富集度增加而有所不同，富集度越高待測(cè)核材料受誘發(fā)裂變概率就越大，從而產(chǎn)生更大量的誘發(fā)裂變中子。特別是15?40 ns內(nèi)，曲線(xiàn)積分值與待測(cè)核材料的富集度直接相關(guān)，這也是NMIS進(jìn)行核材料富集度識(shí)別的基礎(chǔ)。

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

2 NMIS時(shí)-頻特征分析

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

2.1 自相關(guān)

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

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

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

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

2.2 互相關(guān)

探測(cè)器間的互相關(guān)函數(shù)表示為：

式中，當(dāng)t<0時(shí)，源探測(cè)器和探測(cè)器Ι、?Ι間的互相關(guān)計(jì)數(shù)被舍棄，是因?yàn)樘綔y(cè)器事件依賴(lài)于裂變鏈與發(fā)生在其之前的源裂變事件，發(fā)生在源裂變之前的探測(cè)器事件并無(wú)實(shí)際意義。三個(gè)通道互相關(guān)運(yùn)算所得結(jié)果如圖5。

圖5 三路探測(cè)器信號(hào)之間的互相關(guān) (a) 源探測(cè)器與探測(cè)器Ι；(b) 源探測(cè)器與探測(cè)器?Ι；(c) 探測(cè)器Ι與探測(cè)器?ΙFig.5 The cross correlation of three detectors. (a) detector source and Ι; (b) detector source and ?Ι; (c) detector Ι and ?Ι

從圖5(a)、(b)中看出，兩個(gè)探測(cè)器通道和252Cf中子源通道的互相關(guān)函數(shù)與脈沖中子測(cè)量的時(shí)域分布（見(jiàn)圖 3）幾近相同，均反映了源通道信號(hào)與探測(cè)器通道信號(hào)間關(guān)于兩者時(shí)延t的共同信息。顯然，由于環(huán)境噪聲、系統(tǒng)噪聲乃至電子學(xué)響應(yīng)時(shí)間的影響，γ峰和中子峰往往存在時(shí)域交疊的情況。

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

2.3 自功率譜密度

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

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

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

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

2.4 互功率譜密度

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

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

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

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

分別對(duì)三個(gè)通道相互之間的互相關(guān)做 FFT變換，得互功率譜函數(shù)的幅度譜如圖7所示，它反映了源通道信號(hào)和被探測(cè)器?、??測(cè)量通道信號(hào)間的共同信息。

圖7 三路探測(cè)器信號(hào)互功率譜 (a) 源探測(cè)器與探測(cè)器?；(b) 源探測(cè)器與探測(cè)器??；(c) 探測(cè)器?與探測(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)可知，當(dāng)w<>a時(shí)，成反比開(kāi)始衰減。

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

3 結(jié)語(yǔ)

用Monte Carlo方法對(duì)252Cf裂變中子源激發(fā)待測(cè)核材料的過(guò)程進(jìn)行模擬仿真，獲得五種不同富集度待測(cè)核材料的脈沖中子信號(hào)，并以富集度為93.15wt%的HEU為例，對(duì)其三路核信號(hào)進(jìn)行時(shí)-頻運(yùn)算，得到不同種類(lèi)的時(shí)-頻特征標(biāo)簽，進(jìn)而對(duì)其特征參數(shù)進(jìn)行分析，探討了NMIS系統(tǒng)的動(dòng)態(tài)特性及物理意義。

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

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]. 強(qiáng)激光與粒子束, 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 吳沖, 張強(qiáng), 孫志嘉, 等. 低能中子探測(cè)的GEANT4模擬研究[J]. 核技術(shù), 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 裴鹿成, 張孝澤. 蒙特卡羅方法及其在粒子輸運(yùn)問(wèn)題中的應(yīng)用[M]. 北京: 科學(xué)出版社, 1980. 1 PEI Lucheng, ZHANG Xiaoze. Monte-Carlo Method and Applications in Transport Problems of Particles[M].Beijing: Science Press, 1980. 1

11 商靜, 李為民. 電子束輻照裝置屏蔽體預(yù)埋管道處的輻射場(chǎng)計(jì)算[J]. 核技術(shù), 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 王永德, 王軍. 隨機(jī)信號(hào)分析基礎(chǔ)[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