田晨揚(yáng)，過惠平，歐陽曉平，姚青旭，張立波，白德平，焦東明，田慶云

(1.火箭軍工程大學(xué)，陜西 西安 710025；2.96607部隊(duì)，陜西 寶雞 721000；3.西北核技術(shù)研究所，陜西 西安 710000)

多年來，國際原子能機(jī)構(gòu)一直在逐步發(fā)展核材料核查監(jiān)控技術(shù)。目前國際核擴(kuò)散形勢不容樂觀，發(fā)展創(chuàng)新的核材料探測識別技術(shù)具有重大現(xiàn)實(shí)意義[1-3]。根據(jù)鈾材料的物理特性，測量其裂變產(chǎn)物的緩發(fā)γ能譜是一種有前景的技術(shù)方法，可為鈾材料中235U的豐度分析提供技術(shù)手段。國外研究者已將緩發(fā)γ能譜測量與分析成果應(yīng)用于許多領(lǐng)域，如核材料安全管控、海關(guān)貨物檢查、反應(yīng)堆工程及退役核設(shè)施的輻射監(jiān)測等[4-9]。

本文在前期研究獲得的衰變鏈模型[10]基礎(chǔ)上，進(jìn)行實(shí)驗(yàn)測量與驗(yàn)證，以加速器驅(qū)動的脈沖中子源為照射源、鈾樣品為照射對象、HPGe探測器為主要探測器，測量得到鈾材料緩發(fā)γ實(shí)驗(yàn)譜，結(jié)合衰變鏈模型得到緩發(fā)γ計算譜。通過實(shí)驗(yàn)譜來檢驗(yàn)理論計算的準(zhǔn)確性，為鈾材料豐度的實(shí)際分析提供技術(shù)依據(jù)。

1 理論計算

1.1 緩發(fā)γ射線計算公式

中子照射鈾材料后會產(chǎn)生上千種裂變產(chǎn)物，分布在百余條質(zhì)量鏈上，每條質(zhì)量鏈上的核素根據(jù)衰變模式的不同分為β-衰變鏈和EC衰變鏈。在衰變鏈模型[10]中，曾將衰變鏈簡化為基態(tài)線性鏈和激發(fā)態(tài)線性鏈，根據(jù)Bateman方程推導(dǎo)了零時前后各級核素數(shù)目的計算公式，建構(gòu)了計算緩發(fā)γ能譜的MCNPX程序。

計算公式中主要包括以下物理量：零時刻t0、測量開始時刻t1、測量結(jié)束時刻t2、衰變分支比rj、核素j衰變常量λj、核素j在零時的核素數(shù)目Nj(t0)、系數(shù)cjk(j>0，0≤k

在測量期間，衰變鏈任意核素j發(fā)生衰變的原子數(shù)目為nj(j=1，2，3，…)。對于母核有：

(1)

(2)

(3)

(4)

1.2 緩發(fā)γ發(fā)射譜蒙特卡羅計算程序

使用蒙特卡羅程序MCNPX建立緩發(fā)γ發(fā)射譜的探測器響應(yīng)模型，利用衰變數(shù)據(jù)和裂變產(chǎn)額數(shù)據(jù)[11-12]進(jìn)行具體計算。計算過程如下。

1) 輸入計算參數(shù)，讀取裂變產(chǎn)物質(zhì)量數(shù)，讀取每個核素的信息，包括質(zhì)子數(shù)、衰變模式、分支比和半衰期等。

2) 在1個衰變鏈中，如果某個核素的質(zhì)子數(shù)最小，則以它為母核建立1個β-基態(tài)鏈；如果某個核素的質(zhì)子數(shù)最大，則以它為母核建立1個EC基態(tài)鏈；如果核素以激發(fā)態(tài)存在，則以它為母核建立激發(fā)態(tài)鏈。

3) 根據(jù)入射中子類型(熱中子、裂變中子或高能中子)讀取235U的裂變產(chǎn)額數(shù)據(jù)。

4) 讀取γ能量和發(fā)射概率數(shù)據(jù)，利用1.1節(jié)推導(dǎo)的公式，計算對應(yīng)衰變鏈上每個核素的衰變原子數(shù)目，獲得不同能量的緩發(fā)γ射線數(shù)目。

5) 按照質(zhì)量數(shù)A從66到172循環(huán)操作以上步驟，根據(jù)特定算法，對所有裂變產(chǎn)物的緩發(fā)γ射線數(shù)據(jù)進(jìn)行排序，將γ發(fā)射譜數(shù)據(jù)輸出到文件SPECT.txt中。

2 實(shí)驗(yàn)測量

2.1 可行性分析

由于理論計算的多參數(shù)和復(fù)雜性，亟需進(jìn)行特定實(shí)驗(yàn)進(jìn)行檢驗(yàn)。通過實(shí)驗(yàn)，可達(dá)到：1) 比較鈾材料緩發(fā)γ能譜的計算結(jié)果與實(shí)驗(yàn)結(jié)果，驗(yàn)證計算涉及的核數(shù)據(jù)庫參數(shù)的準(zhǔn)確性；2) 通過對鈾材料的主動探測，生成緩發(fā)γ能譜實(shí)驗(yàn)結(jié)果數(shù)據(jù)庫，得到用于分析鈾豐度的特征峰參數(shù)，建立緩發(fā)γ能譜鈾豐度分析方法。

文獻(xiàn)[13]給出了高濃鈾熱中子輻照后的緩發(fā)γ能譜。在實(shí)驗(yàn)中，用慢化252Cf源對40 g的高濃鈾樣品照射100 s，照射后樣品冷卻1 050 s，測量350 s，測得的緩發(fā)γ能譜如圖1a所示。

按照以上實(shí)驗(yàn)條件進(jìn)行理論計算，得到的緩發(fā)γ能譜如圖1b所示。根據(jù)衰變鏈計算，得到1 403個同位素，γ射線能量在0.8～1.55 MeV之間。

利用MCNPX程序?qū)μ綔y器響應(yīng)進(jìn)行模擬計算，得到緩發(fā)γ射線的模擬能譜，如圖1c所示。在模擬時，假設(shè)探測器平均能量分辨率為1 keV，利用F8計數(shù)卡設(shè)置高斯能量展寬參數(shù)。因文獻(xiàn)未給出具體的計數(shù)單位，只能從譜形上對計算譜和實(shí)驗(yàn)譜進(jìn)行對比，可看出計算譜與實(shí)驗(yàn)譜吻合較好，表明可將理論模型用于緩發(fā)γ能譜的計算。

a——實(shí)驗(yàn)結(jié)果；b——發(fā)射譜計算結(jié)果；c——MCNPX模擬結(jié)果圖1 高濃鈾緩發(fā)γ能譜Fig.1 High-enriched uranium delayed γ spectrum

2.2 實(shí)驗(yàn)設(shè)置

在實(shí)驗(yàn)中，使用慢化的中子源對鈾材料誘發(fā)裂變，使用高精度探測器獲取緩發(fā)γ能譜。直線加速器驅(qū)動的脈沖中子源示于圖2，在該裝置中，被加速到5 MeV的電子束轟擊鎢靶(4.2 g/cm2)產(chǎn)生軔致輻射光子。光子穿過鋁層(用來過濾電子)后，被鈹轉(zhuǎn)換器截獲產(chǎn)生中子。9Be(γ,n)8Be反應(yīng)的閾能約為1.7 MeV，反應(yīng)截面與光子能量的關(guān)系如圖3所示，圖中數(shù)據(jù)源自ENDF/B-Ⅶ核數(shù)據(jù)庫[14]。軔致輻射光子主要分布在低能區(qū)，因此反應(yīng)產(chǎn)生的中子主要是熱中子和超熱中子，如圖4所示。因光子最高能量低于其裂變閾能，不會發(fā)生光致裂變，所以被照射樣品的緩發(fā)γ射線只由低能中子誘發(fā)裂變產(chǎn)生。加速器在20 μA電流下運(yùn)行時，中子發(fā)射率為5×109s-1。設(shè)置中子源幾何參數(shù)，并采用MCNPX程序進(jìn)行計算。驗(yàn)證幾何參數(shù)時，在照射過程中對金箔進(jìn)行活化，具體方法為：在鈹圓柱體不同位置放置已知質(zhì)量的薄片金箔，照射結(jié)束后，通過計算中子活化產(chǎn)生的198Au 411 keV γ射線峰面積，得到入射熱中子通量密度。

圖2 直線加速器驅(qū)動的光子-中子源Fig.2 Photon-neutron source driven by linear accelerator

圖3 根據(jù)ENDF/B-Ⅶ核數(shù)據(jù)庫得到的9Be(γ,n)8Be反應(yīng)截面分布Fig.3 9Be(γ,n)8Be cross section distribution based on ENDF/B-Ⅶ nuclear data library

實(shí)驗(yàn)測量采用長照射運(yùn)行模式，具體方法為：1) 將輻照的樣品放在照射區(qū)，用中子照射；2) 關(guān)閉電子束后，將輻照后的樣品傳送到采集區(qū)，期間有1～2 min的延遲；3) 將輻照后的樣品放在能獲取本底譜的采集區(qū)，用探測器測量一段時間。

圖4 Be轉(zhuǎn)換器中單個電子產(chǎn)生的中子通量密度Fig.4 Neutron flux density of single electron in Be converter

將得到的緩發(fā)γ能譜歸一化到計數(shù)率單位，將本底扣除，再將實(shí)驗(yàn)譜與計算預(yù)測譜進(jìn)行比較。實(shí)驗(yàn)參數(shù)可作為后續(xù)模擬計算的輸入?yún)?shù)。后續(xù)計算時，由于直線加速器的高重復(fù)率，可忽略中子源的脈沖特性，并將中子通量密度歸一化為積分電流。利用理論模型可查看中間計算結(jié)果，如中子通量、裂變率、光子源項(xiàng)、個別同位素貢獻(xiàn)等。數(shù)據(jù)采集區(qū)使用探測效率為40%的HPGe探測器，距離樣品23 cm。進(jìn)行計數(shù)率死時間控制時，在樣品和探測器之間放置鉛板。在測量前后，使用文獻(xiàn)[11]中的標(biāo)準(zhǔn)源對探測器進(jìn)行能量刻度。分析刻度結(jié)果，并用最小二乘法擬合成二次多項(xiàng)式。

使用標(biāo)準(zhǔn)源的獨(dú)立γ射線峰數(shù)據(jù)對探測器進(jìn)行能量刻度，對獨(dú)立γ射線峰進(jìn)行高斯擬合得到標(biāo)準(zhǔn)偏差，對HPGe探測器的測得峰進(jìn)行能量展寬。對于高斯能量展寬，γ射線能量與峰半高寬(FWHM)之間的關(guān)系[15]如下：

(5)

式中：E為γ射線能量；a、b、c為擬合系數(shù)。

HPGe探測器的能量分辨率刻度結(jié)果列于表1，其高斯能量展寬擬合示于圖5。由圖5可知，HPGe探測器的能量分辨率與γ射線能量近似呈正比，線性較好。根據(jù)圖5可確定擬合系數(shù)。利用MCNPX程序得到計算能譜時，再輸入這些擬合系數(shù)。

表1 HPGe探測器能量分辨率刻度結(jié)果Table 1 Energy resolution calibration result of HPGe detector

圖5 探測器高斯能量展寬擬合Fig.5 Gaussian energy broadening fitting for detector

2.3 測量結(jié)果

圖6 鈾樣品測量實(shí)驗(yàn)布局(a)和MCNPX 3D視圖(b)Fig.6 Experimental arrangement of uranium sample measurement (a) and MCNPX 3D view (b)

本文所用鈾樣品規(guī)格為9.1 cm×9.1 cm×0.635 cm，質(zhì)量為1.2 kg，235U豐度為90%。實(shí)驗(yàn)中，在鈾樣品前放置聚乙烯層以慢化多余的中子。實(shí)驗(yàn)各時間段分布如下：在20 mA直線加速器電流下輻照30 min，輻照后用30 s將樣品傳輸?shù)讲杉瘏^(qū)，測量60 min，死時間小于2.5%。實(shí)驗(yàn)布局和MCNPX 3D視圖如圖6所示，鈾樣品的緩發(fā)γ實(shí)驗(yàn)譜與計算譜示于圖7a，實(shí)驗(yàn)結(jié)果較計算結(jié)果稍大，原因在于理論模型未考慮輻射俘獲及本底等產(chǎn)生的γ計數(shù)。MCNPX計算得到的緩發(fā)γ發(fā)射譜示于圖7b，包括12 000條譜線。采用ORTEC交互峰擬合程序[16]分析緩發(fā)γ實(shí)驗(yàn)譜與計算譜，計算得到各自的獨(dú)立峰面積。實(shí)驗(yàn)譜與計算譜的峰數(shù)據(jù)列于表2。由表2可知，對同一γ峰面積，計算值與實(shí)驗(yàn)值偏差較小，在2.72%～9.7%之間。

由圖7可知，探測器得到的緩發(fā)γ能譜中，較強(qiáng)的γ能峰集中在0.8～1.55 MeV區(qū)域，可在此區(qū)域選擇特征峰來分析樣品的鈾豐度。發(fā)射譜由于譜線太密集，容易混淆，不便利用其進(jìn)行鈾豐度分析。

2.4 降低方差與實(shí)驗(yàn)準(zhǔn)確性檢驗(yàn)

由于電子、中子、光子的理論輸運(yùn)模型較理想化，會帶來結(jié)果的偏差。利用MCNPX進(jìn)行計算時，建立探測器響應(yīng)較費(fèi)時，為提高M(jìn)CNPX計算效率，需減小方差，同時還需進(jìn)行被動法γ源測量，檢驗(yàn)實(shí)驗(yàn)的準(zhǔn)確性。

圖7 探測器采集的緩發(fā)γ實(shí)驗(yàn)譜與計算譜(a)和理論計算的緩發(fā)γ發(fā)射譜(b)Fig.7 Delayed γ spectra of experiment and calculation for detector (a) and delayed γ emission spectrum of theoretical calculation (b)

表2 緩發(fā)γ射線單個峰面積的實(shí)驗(yàn)和計算結(jié)果Table 2 Experiment and calculation results of single peak area for delayed γ-ray

為提高M(jìn)CNPX對探測器響應(yīng)的計算效率，采取部分確定論方法減小方差。在每個相互作用點(diǎn)，采用隨機(jī)數(shù)模擬粒子輸運(yùn)，確定粒子從探測器散射或逃逸的概率，減少粒子輸運(yùn)時間。在MCNPX中，使用DXTRAN球和F5計數(shù)卡可達(dá)到該目的。

減小方差的方法如下。

1) 光子規(guī)模控制。由于對緩發(fā)γ能譜的分析主要是對能量較高的γ射線進(jìn)行分析，可將能量較低的光子從輸運(yùn)中去除，節(jié)省計算時間。

2) 源方向偏倚。其他各向同性的光子源射入探測器，快速產(chǎn)生大量光子輸運(yùn)事件。對每種幾何設(shè)置，分析探測器對方向偏倚的靈敏度，并將結(jié)果歸一化到均勻源的發(fā)射概率上。

3) 源概率偏倚。緩發(fā)γ源包括幾千條γ射線，各自的發(fā)射概率可相差幾個數(shù)量級。偏倚函數(shù)對所有γ射線提供一個均勻抽樣概率。

采用上述降低方差措施，計算已知同位素源的探測器響應(yīng)。采用與緩發(fā)γ能譜相同的探測系統(tǒng)測量同位素γ源實(shí)驗(yàn)譜，將結(jié)果歸一化到平均計數(shù)率，并扣除本底。同時根據(jù)源項(xiàng)參數(shù)進(jìn)行MCNPX模擬，得到計算譜。60Co、152Eu、88Y的實(shí)驗(yàn)譜和計算譜示于圖8，實(shí)驗(yàn)譜和計算譜的對應(yīng)峰強(qiáng)度列于表3。由圖8可知，計算譜與實(shí)驗(yàn)譜吻合較好。

圖8 60Co、152Eu和88Y的γ能譜Fig.8 γ spectra of 60Co, 152Eu and 88Y

同位素源峰能量/keV峰強(qiáng)度/s-1計算值實(shí)驗(yàn)值計算值與實(shí)驗(yàn)值之比60Co1.13723.61±0.1825.66±0.190.92±0.011.33320.59±0.1722.60±0.180.91±0.01152Eu0.7784.99±0.055.87±0.050.85±0.020.8672.19±0.042.46±0.040.89±0.040.9644.82±0.045.63±0.050.86±0.011.0853.80±0.044.41±0.040.86±0.021.1124.24±0.044.88±0.040.87±0.021.4084.98±0.045.89±0.040.85±0.0188Y0.8982.01±0.022.17±0.020.93±0.031.8361.62±0.021.68±0.020.96±0.04

2.5 特征峰與鈾豐度

選擇94Y的918.8 keV γ峰、89Rb的1 248 keV γ峰、138Cs的1 436 keV γ峰為特征峰，測量不同鈾豐度下的特征峰計數(shù)，峰計數(shù)與235U豐度的關(guān)系如圖9所示。由圖9可知，特征峰計數(shù)與235U豐度近似呈正比，因此可在已知特征峰計數(shù)情況下，計算出樣品的235U豐度。

圖9 94Y、89Rb、138Cs的特征峰計數(shù)與鈾豐度的關(guān)系Fig.9 Relationship between characteristic peak count and uranium enrichment for 94Y, 89Rb and 138Cs

3 結(jié)論

使用線性加速器驅(qū)動的脈沖中子源輻照鈾材料，采用能量線性良好的HPGe探測器進(jìn)行數(shù)據(jù)采集，得到緩發(fā)γ能譜的實(shí)驗(yàn)結(jié)果。對比緩發(fā)γ計算譜與實(shí)驗(yàn)譜發(fā)現(xiàn)，兩者吻合良好，表明理論計算模型具有可靠性。為檢驗(yàn)實(shí)驗(yàn)結(jié)果的準(zhǔn)確性，對60Co、152Eu、88Y同位素源采用相同探測器和MCNPX計算參數(shù)進(jìn)行譜測定，得到的計算譜與實(shí)驗(yàn)譜也較一致，說明實(shí)驗(yàn)結(jié)果是準(zhǔn)確的。最后，在0.8～1.55 MeV能區(qū)選擇94Y 918.8 keV γ峰、89Rb 1 248 keV γ峰和138Cs 1 436 keV γ峰作為特征峰，測量了不同鈾豐度下的特征峰計數(shù)，根據(jù)特征峰計數(shù)與鈾豐度的關(guān)系可計算出235U豐度，達(dá)到分析鈾材料的目的。研究表明，基于緩發(fā)γ能譜可對未知含量的鈾樣品進(jìn)行235U豐度分析，為現(xiàn)場核查提供技術(shù)參考。