田晨扬，过惠平，欧阳晓平，姚青旭，张立波，白德平，焦东明，田庆云

(1.火箭军工程大学，陕西 西安 710025；2.96607部队，陕西 宝鸡 721000；3.西北核技术研究所，陕西 西安 710000)

多年来，国际原子能机构一直在逐步发展核材料核查监控技术。目前国际核扩散形势不容乐观，发展创新的核材料探测识别技术具有重大现实意义[1-3]。根据铀材料的物理特性，测量其裂变产物的缓发γ能谱是一种有前景的技术方法，可为铀材料中235U的丰度分析提供技术手段。国外研究者已将缓发γ能谱测量与分析成果应用于许多领域，如核材料安全管控、海关货物检查、反应堆工程及退役核设施的辐射监测等[4-9]。

本文在前期研究获得的衰变链模型[10]基础上，进行实验测量与验证，以加速器驱动的脉冲中子源为照射源、铀样品为照射对象、HPGe探测器为主要探测器，测量得到铀材料缓发γ实验谱，结合衰变链模型得到缓发γ计算谱。通过实验谱来检验理论计算的准确性，为铀材料丰度的实际分析提供技术依据。

1 理论计算

1.1 缓发γ射线计算公式

中子照射铀材料后会产生上千种裂变产物，分布在百余条质量链上，每条质量链上的核素根据衰变模式的不同分为β-衰变链和EC衰变链。在衰变链模型[10]中，曾将衰变链简化为基态线性链和激发态线性链，根据Bateman方程推导了零时前后各级核素数目的计算公式，建构了计算缓发γ能谱的MCNPX程序。

计算公式中主要包括以下物理量：零时刻t0、测量开始时刻t1、测量结束时刻t2、衰变分支比rj、核素j衰变常量λj、核素j在零时的核素数目Nj(t0)、系数cjk(j>0，0≤k

在测量期间，衰变链任意核素j发生衰变的原子数目为nj(j=1，2，3，…)。对于母核有：

(1)

(2)

(3)

(4)

1.2 缓发γ发射谱蒙特卡罗计算程序

使用蒙特卡罗程序MCNPX建立缓发γ发射谱的探测器响应模型，利用衰变数据和裂变产额数据[11-12]进行具体计算。计算过程如下。

1) 输入计算参数，读取裂变产物质量数，读取每个核素的信息，包括质子数、衰变模式、分支比和半衰期等。

2) 在1个衰变链中，如果某个核素的质子数最小，则以它为母核建立1个β-基态链；如果某个核素的质子数最大，则以它为母核建立1个EC基态链；如果核素以激发态存在，则以它为母核建立激发态链。

3) 根据入射中子类型(热中子、裂变中子或高能中子)读取235U的裂变产额数据。

4) 读取γ能量和发射概率数据，利用1.1节推导的公式，计算对应衰变链上每个核素的衰变原子数目，获得不同能量的缓发γ射线数目。

5) 按照质量数A从66到172循环操作以上步骤，根据特定算法，对所有裂变产物的缓发γ射线数据进行排序，将γ发射谱数据输出到文件SPECT.txt中。

2 实验测量

2.1 可行性分析

由于理论计算的多参数和复杂性，亟需进行特定实验进行检验。通过实验，可达到：1) 比较铀材料缓发γ能谱的计算结果与实验结果，验证计算涉及的核数据库参数的准确性；2) 通过对铀材料的主动探测，生成缓发γ能谱实验结果数据库，得到用于分析铀丰度的特征峰参数，建立缓发γ能谱铀丰度分析方法。

文献[13]给出了高浓铀热中子辐照后的缓发γ能谱。在实验中，用慢化252Cf源对40 g的高浓铀样品照射100 s，照射后样品冷却1 050 s，测量350 s，测得的缓发γ能谱如图1a所示。

按照以上实验条件进行理论计算，得到的缓发γ能谱如图1b所示。根据衰变链计算，得到1 403个同位素，γ射线能量在0.8～1.55 MeV之间。

利用MCNPX程序对探测器响应进行模拟计算，得到缓发γ射线的模拟能谱，如图1c所示。在模拟时，假设探测器平均能量分辨率为1 keV，利用F8计数卡设置高斯能量展宽参数。因文献未给出具体的计数单位，只能从谱形上对计算谱和实验谱进行对比，可看出计算谱与实验谱吻合较好，表明可将理论模型用于缓发γ能谱的计算。

a——实验结果；b——发射谱计算结果；c——MCNPX模拟结果图1 高浓铀缓发γ能谱Fig.1 High-enriched uranium delayed γ spectrum

2.2 实验设置

在实验中，使用慢化的中子源对铀材料诱发裂变，使用高精度探测器获取缓发γ能谱。直线加速器驱动的脉冲中子源示于图2，在该装置中，被加速到5 MeV的电子束轰击钨靶(4.2 g/cm2)产生轫致辐射光子。光子穿过铝层(用来过滤电子)后，被铍转换器截获产生中子。9Be(γ,n)8Be反应的阈能约为1.7 MeV，反应截面与光子能量的关系如图3所示，图中数据源自ENDF/B-Ⅶ核数据库[14]。轫致辐射光子主要分布在低能区，因此反应产生的中子主要是热中子和超热中子，如图4所示。因光子最高能量低于其裂变阈能，不会发生光致裂变，所以被照射样品的缓发γ射线只由低能中子诱发裂变产生。加速器在20 μA电流下运行时，中子发射率为5×109s-1。设置中子源几何参数，并采用MCNPX程序进行计算。验证几何参数时，在照射过程中对金箔进行活化，具体方法为：在铍圆柱体不同位置放置已知质量的薄片金箔，照射结束后，通过计算中子活化产生的198Au 411 keV γ射线峰面积，得到入射热中子通量密度。

图2 直线加速器驱动的光子-中子源Fig.2 Photon-neutron source driven by linear accelerator

图3 根据ENDF/B-Ⅶ核数据库得到的9Be(γ,n)8Be反应截面分布Fig.3 9Be(γ,n)8Be cross section distribution based on ENDF/B-Ⅶ nuclear data library

实验测量采用长照射运行模式，具体方法为：1) 将辐照的样品放在照射区，用中子照射；2) 关闭电子束后，将辐照后的样品传送到采集区，期间有1～2 min的延迟；3) 将辐照后的样品放在能获取本底谱的采集区，用探测器测量一段时间。

图4 Be转换器中单个电子产生的中子通量密度Fig.4 Neutron flux density of single electron in Be converter

将得到的缓发γ能谱归一化到计数率单位，将本底扣除，再将实验谱与计算预测谱进行比较。实验参数可作为后续模拟计算的输入参数。后续计算时，由于直线加速器的高重复率，可忽略中子源的脉冲特性，并将中子通量密度归一化为积分电流。利用理论模型可查看中间计算结果，如中子通量、裂变率、光子源项、个别同位素贡献等。数据采集区使用探测效率为40%的HPGe探测器，距离样品23 cm。进行计数率死时间控制时，在样品和探测器之间放置铅板。在测量前后，使用文献[11]中的标准源对探测器进行能量刻度。分析刻度结果，并用最小二乘法拟合成二次多项式。

使用标准源的独立γ射线峰数据对探测器进行能量刻度，对独立γ射线峰进行高斯拟合得到标准偏差，对HPGe探测器的测得峰进行能量展宽。对于高斯能量展宽，γ射线能量与峰半高宽(FWHM)之间的关系[15]如下：

(5)

式中：E为γ射线能量；a、b、c为拟合系数。

HPGe探测器的能量分辨率刻度结果列于表1，其高斯能量展宽拟合示于图5。由图5可知，HPGe探测器的能量分辨率与γ射线能量近似呈正比，线性较好。根据图5可确定拟合系数。利用MCNPX程序得到计算能谱时，再输入这些拟合系数。

表1 HPGe探测器能量分辨率刻度结果Table 1 Energy resolution calibration result of HPGe detector

图5 探测器高斯能量展宽拟合Fig.5 Gaussian energy broadening fitting for detector

2.3 测量结果

图6 铀样品测量实验布局(a)和MCNPX 3D视图(b)Fig.6 Experimental arrangement of uranium sample measurement (a) and MCNPX 3D view (b)

本文所用铀样品规格为9.1 cm×9.1 cm×0.635 cm，质量为1.2 kg，235U丰度为90%。实验中，在铀样品前放置聚乙烯层以慢化多余的中子。实验各时间段分布如下：在20 mA直线加速器电流下辐照30 min，辐照后用30 s将样品传输到采集区，测量60 min，死时间小于2.5%。实验布局和MCNPX 3D视图如图6所示，铀样品的缓发γ实验谱与计算谱示于图7a，实验结果较计算结果稍大，原因在于理论模型未考虑辐射俘获及本底等产生的γ计数。MCNPX计算得到的缓发γ发射谱示于图7b，包括12 000条谱线。采用ORTEC交互峰拟合程序[16]分析缓发γ实验谱与计算谱，计算得到各自的独立峰面积。实验谱与计算谱的峰数据列于表2。由表2可知，对同一γ峰面积，计算值与实验值偏差较小，在2.72%～9.7%之间。

由图7可知，探测器得到的缓发γ能谱中，较强的γ能峰集中在0.8～1.55 MeV区域，可在此区域选择特征峰来分析样品的铀丰度。发射谱由于谱线太密集，容易混淆，不便利用其进行铀丰度分析。

2.4 降低方差与实验准确性检验

由于电子、中子、光子的理论输运模型较理想化，会带来结果的偏差。利用MCNPX进行计算时，建立探测器响应较费时，为提高MCNPX计算效率，需减小方差，同时还需进行被动法γ源测量，检验实验的准确性。

图7 探测器采集的缓发γ实验谱与计算谱(a)和理论计算的缓发γ发射谱(b)Fig.7 Delayed γ spectra of experiment and calculation for detector (a) and delayed γ emission spectrum of theoretical calculation (b)

表2 缓发γ射线单个峰面积的实验和计算结果Table 2 Experiment and calculation results of single peak area for delayed γ-ray

为提高MCNPX对探测器响应的计算效率，采取部分确定论方法减小方差。在每个相互作用点，采用随机数模拟粒子输运，确定粒子从探测器散射或逃逸的概率，减少粒子输运时间。在MCNPX中，使用DXTRAN球和F5计数卡可达到该目的。

减小方差的方法如下。

1) 光子规模控制。由于对缓发γ能谱的分析主要是对能量较高的γ射线进行分析，可将能量较低的光子从输运中去除，节省计算时间。

2) 源方向偏倚。其他各向同性的光子源射入探测器，快速产生大量光子输运事件。对每种几何设置，分析探测器对方向偏倚的灵敏度，并将结果归一化到均匀源的发射概率上。

3) 源概率偏倚。缓发γ源包括几千条γ射线，各自的发射概率可相差几个数量级。偏倚函数对所有γ射线提供一个均匀抽样概率。

采用上述降低方差措施，计算已知同位素源的探测器响应。采用与缓发γ能谱相同的探测系统测量同位素γ源实验谱，将结果归一化到平均计数率，并扣除本底。同时根据源项参数进行MCNPX模拟，得到计算谱。60Co、152Eu、88Y的实验谱和计算谱示于图8，实验谱和计算谱的对应峰强度列于表3。由图8可知，计算谱与实验谱吻合较好。

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

同位素源峰能量/keV峰强度/s-1计算值实验值计算值与实验值之比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 γ峰为特征峰，测量不同铀丰度下的特征峰计数，峰计数与235U丰度的关系如图9所示。由图9可知，特征峰计数与235U丰度近似呈正比，因此可在已知特征峰计数情况下，计算出样品的235U丰度。

图9 94Y、89Rb、138Cs的特征峰计数与铀丰度的关系Fig.9 Relationship between characteristic peak count and uranium enrichment for 94Y, 89Rb and 138Cs

3 结论

使用线性加速器驱动的脉冲中子源辐照铀材料，采用能量线性良好的HPGe探测器进行数据采集，得到缓发γ能谱的实验结果。对比缓发γ计算谱与实验谱发现，两者吻合良好，表明理论计算模型具有可靠性。为检验实验结果的准确性，对60Co、152Eu、88Y同位素源采用相同探测器和MCNPX计算参数进行谱测定，得到的计算谱与实验谱也较一致，说明实验结果是准确的。最后，在0.8～1.55 MeV能区选择94Y 918.8 keV γ峰、89Rb 1 248 keV γ峰和138Cs 1 436 keV γ峰作为特征峰，测量了不同铀丰度下的特征峰计数，根据特征峰计数与铀丰度的关系可计算出235U丰度，达到分析铀材料的目的。研究表明，基于缓发γ能谱可对未知含量的铀样品进行235U丰度分析，为现场核查提供技术参考。