鄧力 李瑞 王鑫 付元光

1) (北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所， 北京 100094)2) (中物院高性能數(shù)值模擬軟件中心， 北京 100088)(2020 年2 月24日收到; 2020 年3 月23日收到修改稿)

蒙特卡羅方法(MC)是模擬核探測(cè)問(wèn)題的理想方法， 用中子照射客體， 中子誘發(fā)產(chǎn)生非彈γ和俘獲γ， 通過(guò)特征γ射線能譜和時(shí)間譜分析， 確定客體核素組成和重量百分比. 本文基于非彈γ和俘獲γ時(shí)間門(mén)測(cè)量技術(shù)， 給出了脈沖源發(fā)射下探測(cè)器響應(yīng)計(jì)數(shù)公式. 在中子與核作用產(chǎn)生次級(jí)光子方面， 采用期望值估計(jì)(expect value estimator， EVE)產(chǎn)光. 為了避免大量小權(quán)光子模擬帶來(lái)的計(jì)算存儲(chǔ)量增加， 設(shè)計(jì)了EVE產(chǎn)光與直接估計(jì)(direct estimator， DE)產(chǎn)光耦合. 僅增加少量計(jì)算時(shí)間， 便實(shí)現(xiàn)了特征γ射線解譜. 數(shù)值模擬在自主MC軟件JMCT上開(kāi)展， 計(jì)算結(jié)果初步驗(yàn)證了方法的正確有效性.

1 引 言

蒙特卡羅方法(Monte Carlo， MC)具有模擬核探測(cè)問(wèn)題的天然優(yōu)勢(shì)， 近年已廣泛用于X-射線熒光分析， 在線中子俘獲瞬發(fā)γ射線分析，基于γ射線光譜的脈沖中子孔隙度測(cè)井， 隱藏爆炸物探測(cè)等. 相比X光常規(guī)探測(cè)， 中子誘發(fā)γ射線探測(cè)具有穿透能力強(qiáng)， 容易穿透包括鋼、常規(guī)/化學(xué)武器包殼、核武器內(nèi)部結(jié)構(gòu)， 可用于確定客體的化學(xué)變化等， 這是核探測(cè)的優(yōu)勢(shì)所在.

目前核探測(cè)的主要手段是利用高能中子照射客體， 產(chǎn)生次級(jí)光子， 通過(guò)特征γ射線能譜、時(shí)間譜測(cè)量， 確定客體核素組成及份額. 在中子作用下，絕大部分元素都可以發(fā)射可辨認(rèn)的特征γ譜線， 慢中子可以引起除氮以外所有元素的非彈性散射反應(yīng)， 并發(fā)射特征γ譜線. 利用特征γ射線原級(jí)線光子不隨入射中子能量變化這一特點(diǎn)， 通過(guò)特征γ射線的能峰特征， 確定客體內(nèi)所含核素， 通過(guò)特征γ射線直穿貢獻(xiàn)， 確定客體核素份額. 方法可用于化學(xué)/常規(guī)武器甄別、隱藏爆炸物/毒品探測(cè)、放射性石油測(cè)井和探測(cè)器靈敏度優(yōu)化設(shè)計(jì)等.

由于核探測(cè)環(huán)境復(fù)雜， 涉及空間r、能量E、方向Ω、時(shí)間t七維變量的Boltzmann方程求解， 一般采用連續(xù)能量MC求解. 目前包括著名的MCNP[1]程序， 特征γ射線解譜還存在一定困難，雖然F8計(jì)數(shù)理論上能夠算出探測(cè)器響應(yīng)-脈沖高度譜， 但代價(jià)大， 解很難收斂. 另外， 中子產(chǎn)生次級(jí)光子采用直接估計(jì)， 存在隨機(jī)因素的影響， 無(wú)法保證不漏掉某些重要核素的貢獻(xiàn). 文獻(xiàn)[2-6]在這方面開(kāi)展了大量研究， 針對(duì)特征γ射線解譜， 發(fā)展了多項(xiàng)技巧. 我們參考了其部分研究成果， 近期在自主MC粒子輸運(yùn)軟件JMCT[7-12]上開(kāi)發(fā)了特征γ射線分類(lèi)標(biāo)識(shí)計(jì)算及非彈γ射線和俘獲γ射線時(shí)間箱計(jì)數(shù)方法. 特別在中子與核作用產(chǎn)生次級(jí)光子處理上， 先采用期望值估計(jì)(expect value estimator， EVE)產(chǎn)光， 分別統(tǒng)計(jì)原級(jí)線光子和原級(jí)連續(xù)光子對(duì)探測(cè)器的直穿貢獻(xiàn). 為了避免EVE產(chǎn)生的大量小權(quán)光子模擬帶來(lái)的計(jì)算量和存儲(chǔ)量的增加， 做完直穿貢獻(xiàn)后， 光子歷史結(jié)束， 立即回到原來(lái)的直接估計(jì)(direct estimator， DE)產(chǎn)光模式， 每次最多產(chǎn)生一個(gè)光子， 對(duì)其進(jìn)行跟蹤，只記錄散射貢獻(xiàn)， 這樣在不增加計(jì)算存儲(chǔ)量下， 實(shí)現(xiàn)了特征γ射線解譜.

2 問(wèn)題描述

高能中子與核發(fā)生非彈(用(n， n' )表示)和俘獲(用(n， γ)表示)反應(yīng)時(shí)， 將產(chǎn)生次級(jí)光子， 用次級(jí)光子的特征γ射線來(lái)確定客體的核素組成和份額， 這是核探測(cè)相對(duì)X光探測(cè)的優(yōu)勢(shì)所在. 輻射俘獲是吸收反應(yīng)中最重要的反應(yīng)之一， 其反應(yīng)產(chǎn)物之一就是γ射線， 產(chǎn)生次級(jí)光子的主要反應(yīng)道有(n， γ)和裂變(用(n， f )表示). 如果中子源的能量和強(qiáng)度較高， 則中子與核發(fā)生(n， n' )反應(yīng)的概率增大， 并產(chǎn)生非彈γ射線. 表1給出H， C， N等11種核素發(fā)射非彈γ譜線和俘獲γ譜線的能量[13]，表2給出烈性炸藥(TNT)和某些化學(xué)武器中所含元素的重量百分比[13].

表 1 H， C， N， O等核素發(fā)射俘獲γ譜線和非彈性散射γ譜線能量Table 1. Energy of spectrum line from inelastic γ and capture γ about H， C， N， O， etc.

表 2 烈性炸藥(TNT)和某些化學(xué)武器中所含元素的重量百分比Table 2. Weight percentage of elements in some spirited detonators (TNT) and chemical weapons.

2.1 非彈γ和俘獲γ時(shí)間門(mén)測(cè)量方法

對(duì)于某些問(wèn)題， 需要通過(guò)發(fā)射脈沖源， 通過(guò)時(shí)間箱計(jì)數(shù)來(lái)區(qū)別非彈γ射線和俘獲γ射線， 圖1給出放射性石油測(cè)井中， 碳氧比(C/O)能譜測(cè)井的定時(shí)測(cè)量圖[14-16]. 用脈沖方式發(fā)射中子， 通常0—10 μs為脈沖門(mén)發(fā)射時(shí)間間隔， 10—20 μs為本底門(mén)時(shí)間間隔， 20—90 μs為晚俘獲門(mén)時(shí)間間隔，依據(jù)這種邏輯關(guān)系， 可以得到: 凈非彈γ計(jì)數(shù)=脈沖門(mén)譜計(jì)數(shù)—本底譜計(jì)數(shù); 俘獲γ計(jì)數(shù)=晚俘獲門(mén)計(jì)數(shù). 這個(gè)時(shí)間門(mén)測(cè)量過(guò)程同樣適合其他含時(shí)核探測(cè)問(wèn)題的模擬.

對(duì)應(yīng)圖1給出不同時(shí)間門(mén)下探測(cè)器脈沖高度譜計(jì)算公式. 設(shè)Nδ(E， t)為以δ(t)脈沖方式發(fā)射的E0= 14.1 MeV氘氚中子源探測(cè)器中測(cè)到的γ能譜的時(shí)間響應(yīng). 根據(jù)中子-光子輸運(yùn)方程的線性性質(zhì)， 對(duì)應(yīng)于任意時(shí)間分布S(t)的中子源， 探測(cè)儀中測(cè)到的γ能譜時(shí)間響應(yīng)為

其中E表示能量， 單位MeV; t表示時(shí)間， 單位μs.

圖 1 碳氧比測(cè)井中子引發(fā)非彈性散射γ與俘獲γ定時(shí)邏輯圖Fig. 1. The timing diagram of neutron induced inelastic γ and capture γ in C/O well-logging.

在任意測(cè)量門(mén)[ta， tb]內(nèi)記錄的γ能譜為

令中子脈沖時(shí)間分布函數(shù)S(t)為周期函數(shù)，其周期為τ， 有

其中t∈[0， τ]. 又假定S(t)為寬度為τ的函數(shù)

其中S0為一個(gè)脈沖內(nèi)釋放的中子總數(shù); T為終態(tài)時(shí)間; fδ(t) 為脈沖時(shí)間分布函數(shù)， 滿足歸一條件. 則有

將(5)式代入(2)式， 對(duì)t積分便得任意時(shí)間門(mén)內(nèi)的γ能譜強(qiáng)度.

脈沖門(mén)內(nèi)的測(cè)量值為

本底門(mén)內(nèi)測(cè)量值為

根據(jù)中子進(jìn)入客體誘發(fā)的各種γ射線的時(shí)間特點(diǎn)， 將(5)式中的Nδ(E， t)分解為4個(gè)部分， 即

其中 Niδ(E)，i=1，2，3，4 分別為單位強(qiáng)度δ(t)脈沖中子源引起的非彈性γ射線、慢化過(guò)程中的俘獲γ射線、熱中子俘獲γ射線和活化反應(yīng)γ射線的γ能譜強(qiáng)度; fi(t)， i = 1， 2， 3， 4分別為上述4種γ射線的時(shí)間譜， 滿足歸一條件1， 2， 3， 4).

經(jīng)過(guò)一系列公式推導(dǎo)， 得到脈沖門(mén)能譜

可見(jiàn)脈沖門(mén)與本底門(mén)測(cè)量值之差除了非彈性γ能譜， 還受到少量慢化過(guò)程俘獲γ射線和少量未扣除干凈的熱中子俘獲γ譜的“污染”. 考慮到活化γ射線的發(fā)射時(shí)刻很晚， 加之活化γ射線的發(fā)射壽命?T . 因此， 在脈沖門(mén)與本底門(mén)譜之差中， 活化γ射線的污染是很小的. 另外， 由于慢化過(guò)程中的中子俘獲數(shù)目遠(yuǎn)遠(yuǎn)小于熱中子俘獲的數(shù)目， 所以，N3δ(E)?N2δ(E) ， 可以認(rèn)為污染源主要來(lái)自熱中子俘獲γ射線， 其中污染量大小除與中子源的特征及客體物理性質(zhì)有關(guān)， 探測(cè)儀器材料的適當(dāng)選擇也起一定的作用.

要準(zhǔn)確算出探測(cè)器響應(yīng)的測(cè)量能譜N(E， t)，首先應(yīng)設(shè)法算出Nδ(E， t)， 然后與中子源脈沖時(shí)間譜形fδ(t)卷積， 并分別對(duì)兩個(gè)測(cè)量門(mén)的時(shí)間間隔積分， 相減后得到. 因此， MC模擬的關(guān)鍵是算出t時(shí)刻單位強(qiáng)度δ-脈沖源相應(yīng)的探測(cè)器計(jì)數(shù)Nδ(E， t)， 然后與給定的中子脈沖譜形函數(shù)fδ(t)卷積，再對(duì)指定的時(shí)間門(mén)積分， 得到任意時(shí)間門(mén)內(nèi)的測(cè)量 值.

2.2 N(E， t)的MC模擬

快中子非彈性散射， 它所要測(cè)量的是非彈性γ射線. 理論上通過(guò)解中子-光子-電子耦合輸運(yùn)求出N(E， t). 由于探測(cè)器相對(duì)整個(gè)問(wèn)題系統(tǒng)很小， 盡管從源發(fā)出了大量的中子， 但能夠進(jìn)入探測(cè)器的次級(jí)光子數(shù)還是非常有限的， 依據(jù)少量的次級(jí)光子要算準(zhǔn)探測(cè)器響應(yīng)幾乎是不可能的. 因此， 探測(cè)器響應(yīng)-脈沖高度譜計(jì)算通常分三步進(jìn)行.

第一步: 解中子-光子耦合輸運(yùn)方程， 求出進(jìn)入探測(cè)器表面的次級(jí)γ流，

其中S為探測(cè)器表面， n為探測(cè)器表面外法向矢量， E0為入射光子能量.

第二步: 解光子-電子耦合輸運(yùn)方程， 求出探測(cè)器響應(yīng)函數(shù)R(E0， h). 以碘化鈉閃爍探測(cè)器為例，其響應(yīng)函數(shù)形式為

其中 G (E，h) 為 高斯函數(shù); D (E0，E) 為能量沉積譜;η(E0) 為探測(cè)器效率; h為能量道， 單位MeV. 由于響應(yīng)函數(shù)中含有高斯函數(shù)G， 為了保證求積精度，分點(diǎn)要足夠多， 通常在全能區(qū)[0， Emax]分256道.根據(jù)Berger等[17]和Jin等[18]的研究， 不同入射方向的光子對(duì)探測(cè)器響應(yīng)影響很小， 因此， 方向可近似為各向同性， 當(dāng)探測(cè)器形狀及材料一定后，響應(yīng)函數(shù)R(E0， h)(矩陣形式)只需計(jì)算一次.

第三步: 卷積積分得到探測(cè)器響應(yīng)脈沖高度譜

N(h，t) 求解涉及瞬態(tài)中子-光子-電子耦合輸運(yùn)計(jì)算， 計(jì)數(shù)包括時(shí)間-能量聯(lián)合譜， 是MC粒子輸運(yùn)計(jì)算 中難度最大， 模擬最復(fù)雜的過(guò)程.

2.3 特征γ射線標(biāo)識(shí)計(jì)數(shù)

中子發(fā)生非彈(n， n')或俘獲(n， γ)， (n， f )反應(yīng)時(shí)， 將產(chǎn)生次級(jí)光子， 這些光子統(tǒng)稱(chēng)為特征γ射線. 根據(jù)次級(jí)光子能量特征， 特征γ射線分為原級(jí)線光子和原級(jí)連續(xù)光子， 其定義如下:

其中i為核; j為反應(yīng)道; En為中子發(fā)生非彈或俘獲反應(yīng)時(shí)的能量; 通常取Eline= 0.001 MeV作為原級(jí)線光子和原級(jí)連續(xù)光子的分界能量; LP為ENDF數(shù)據(jù)庫(kù)反應(yīng)律序號(hào)；Ai為碰撞核i的原子序數(shù)；EG(i)為特征γ能量.

LP ≠ 2對(duì)應(yīng)的γ原級(jí)線光子， 它不隨入射中子能量變化， 是甄別客體核素組成的核心關(guān)鍵. 在表1中給出了主要核素γ譜線能量峰位置， 表2給出烈性炸藥(TNT)和某些化學(xué)武器中所含元素的重量百分比[13]， 依據(jù)表1和表2所列參考值， 可以快速識(shí)別出客體是否為危禁品. MCNP程序模擬次級(jí)光子時(shí)， 沒(méi)有把這兩類(lèi)光子區(qū)別開(kāi)來(lái)， 因此，次級(jí)光子能譜中的一些特征γ峰會(huì)被散射能譜磨平， 無(wú)法確定客體的核素組成.

當(dāng)前ENDF最新評(píng)價(jià)核數(shù)據(jù)庫(kù)提供了兩種中子產(chǎn)光子模式: 1) 30 × 20產(chǎn)光模式， 中子從10—5eV到20 MeV分30個(gè)能群， 每個(gè)中子群對(duì)應(yīng)20個(gè)等高度譜線， 次級(jí)光子的發(fā)射方向按各項(xiàng)同性處理;2)按碰撞核對(duì)應(yīng)反應(yīng)道的產(chǎn)光概率產(chǎn)光子. 模擬采用第2)種產(chǎn)光模式.

2.4 EVE產(chǎn)光

設(shè)中子與i核發(fā)生j種反應(yīng)的產(chǎn)光概率為pi，j(j = 1， 2， ···， J )， 滿足歸一條件， 這里J為i核對(duì)應(yīng)的反應(yīng)道總數(shù).

DE方法已知i， 抽隨機(jī)數(shù)ξ， 求出滿足不等式的j， 則j反應(yīng)道產(chǎn)光， 光子權(quán)重為wγ.

EVE方法已知i， 按概率pi，j(j = 1， 2， ···，J )全部產(chǎn)光， 相應(yīng)光子權(quán)重為

不難證明兩種估計(jì)方法的數(shù)學(xué)期望是一致的.DE存在隨機(jī)因素， 某些小概率大貢獻(xiàn)事件會(huì)因?yàn)殡S機(jī)因素少抽或漏抽， 這對(duì)隱藏爆炸物探測(cè)這類(lèi)問(wèn)題是不允許的. 故采用EVE是必要的. 我們?cè)缙陂_(kāi)發(fā)研制的MCCO程序， 對(duì)EVE產(chǎn)出的光子進(jìn)行全部跟蹤， 用統(tǒng)計(jì)估計(jì)計(jì)數(shù)[19]， 增加了大量計(jì)算存儲(chǔ)量. 考慮到特征γ射線探測(cè)主要關(guān)心的是次級(jí)光子的直穿貢獻(xiàn)， 實(shí)際上跟蹤所有小權(quán)次級(jí)光子的散射過(guò)程沒(méi)必要. 為此， 我們?cè)O(shè)計(jì)了現(xiàn)在的組合產(chǎn)光模式.

組合產(chǎn)光模式對(duì)EVE產(chǎn)生的大量小權(quán)光子， 僅做直穿估計(jì)， 之后回到原來(lái)的DE產(chǎn)光， 每次最多產(chǎn)生一個(gè)次級(jí)光子， 對(duì)該光子進(jìn)行全程跟蹤， 只統(tǒng)計(jì)散射貢獻(xiàn)， 最后的數(shù)值模擬解由兩部分組成， 即

由于直穿估計(jì)可以解析計(jì)算， 所花時(shí)間可以忽略不計(jì). 因此， 采用組合產(chǎn)光模式后， 實(shí)現(xiàn)了特征γ射線能量-時(shí)間解譜， 而總計(jì)算時(shí)間和存儲(chǔ)量增加很少.

3 數(shù)值實(shí)驗(yàn)

炸藥球模型幾何單位為cm， 對(duì)行李箱進(jìn)行安檢. 如圖2所示， 設(shè)行李箱長(zhǎng)、寬、高分別為80， 30，50， 內(nèi)放一半徑為3.51， 密度為ρ = 1.654 g/cm3的RDX炸藥小球(分子式為C2.63H4.69N0.658O0.85);探測(cè)器位于行李箱上方， 為3'' × 3''碘化鈉閃爍晶體正圓柱探測(cè)器， 密度為ρ = 3.67 g/cm3， 柱中心坐標(biāo)為(0， 0， 38.75); 采用14.1 MeV各項(xiàng)同性氘氚中子點(diǎn)源， 源位置為(0， 45， 0)， 計(jì)算進(jìn)入探測(cè)器的次級(jí)γ流Jγ(E， t). 采用源方向偏倚發(fā)射， 偏倚立體張角覆蓋行李箱. 對(duì)比程序選擇MCNP[1]， 采用相同的數(shù)據(jù)庫(kù)、樣本數(shù)和源方向偏倚技巧.

圖 2 行李箱模型示意圖Fig. 2. Sketch of luggage model.

表3給出JMCT與MCNP探測(cè)器次級(jí)光子流Jγ結(jié)果比較; 表4給出JMCT統(tǒng)計(jì)的箱子內(nèi)的核素組成及份額， 與表1給出的TNT炸藥H， C，N， O核素百分比基本相符(炸藥類(lèi)型不同); 圖3(a)和圖3(b)分別給出原級(jí)線光子和康普頓散射光子能譜; 圖3(c)和圖3(d)給出JMCT與MCNP次級(jí)γ流的能譜比較， 可以看出次級(jí)γ總流JMCT與MCNP結(jié)果符合良好， 能譜差異很小; 圖4給出次級(jí)γ流時(shí)間譜， JMCT與MCNP時(shí)間譜總計(jì)數(shù)相符， 時(shí)間譜細(xì)節(jié)上有些差異， 分析原因是由于時(shí)間箱分得過(guò)細(xì)， 某些時(shí)間箱計(jì)數(shù)的統(tǒng)計(jì)誤差偏大.

圖 3 次級(jí)γ射線能譜計(jì)算結(jié)果比較 (a)次級(jí)γ原級(jí)線光子能譜; (b)原級(jí)連續(xù)光子與Compton散射能譜; (c) JMCT次級(jí)γ總能譜; (d) MCNP次級(jí)γ總能譜Fig. 3. Comparison of calculated result about energy spectra of secondary γ: (a) JMCT primary line γ; (b) JMCT Compton γ;(c) JMCT total γ; (d) MCNP total γ.

表 3 JMCT與MCNP次級(jí)γ流計(jì)算結(jié)果比較Table 3. Comparison of calculated results about secondary γ between JMCT and MCNP.

表 4 H， C， N， O瞬發(fā)γ計(jì)數(shù)及份額Table 4. Count and percentage of prompt γ from H， C， N and O.

圖 4 JMCT與MCNP次級(jí)γ流時(shí)間譜比較Fig. 4. Comparison of secondary γ-fluent tine spectrum between JMCT and MCNP.

4 結(jié) 論

基于JMCT程序開(kāi)發(fā)了特征γ射線能量-時(shí)間聯(lián)合譜計(jì)算功能， 數(shù)值實(shí)驗(yàn)初步驗(yàn)證了軟件計(jì)算的正確高效性. 目前已開(kāi)展了數(shù)十種炸藥模型的計(jì)算， 構(gòu)建了相應(yīng)的數(shù)據(jù)庫(kù)， 為反演計(jì)算做準(zhǔn)備. 后續(xù)將配合實(shí)驗(yàn)樣機(jī)的研制， 對(duì)測(cè)量?jī)x器進(jìn)行刻度工作， 形成能快速做出判斷的計(jì)算機(jī)軟件正演/反演系統(tǒng)， 為在線檢測(cè)提供理論技術(shù)支持.