基于信赖域算法的中子解谱方法研究

基于信赖域算法的中子解谱方法研究
林存宝;刘国福;杨俊
【摘要】提出了一种新的基于信赖域算法的中子能谱求解方法.使用O5S软件仿真了能量范围为0.25～16 MeV、能量间隔为0.25 MeV的BC501A液体闪烁体探测器响应函数.利用该响应函数仿真验证了入射中子能谱分别为单峰和多峰情况下算法的解谱效果,并通过D-T单能中子源、241 Am-Be和252 Cf连续中子源对解谱方法进行了实验验证.结果表明,使用信赖域算法求得的能谱与参考谱线具有较好的一致性,初步验证了解谱方法的有效性.%A novel method for unfolding neutron spectrum based on trust region algorithm was presented. The response functions of the BC501A liquid scintillation detector in the energy range from 0. 25 MeV to 16 MeV with an energy interval of 0. 25 MeV were generated by the code O5S. Both the single and multi-peak spectra were simulated to test the capability of the new method with the calculated response functions. Meanwhile, D-T monoenergetic neutron source,
241Arn-Be continuous neutron source and 232Cf continuous neutron source were also measured to verify the performance of this approach. The results are in good agreement with that of the reference spectra and the effectiveness of the above proposed method is demonstrated.
【期刊名称】《原子能科学技术》
【年(卷),期】2012(046)012
【总页数】6页(P1480-1485)
【关键词】中子解谱;信赖域算法;O5S
【作者】林存宝;刘国福;杨俊
【作者单位】国防科学技术大学仪器科学与技术系,湖南长沙410073;国防科学技
术大学仪器科学与技术系,湖南长沙410073;国防科学技术大学仪器科学与技术系,
湖南长沙410073
【正文语种】中文
【中图分类】TL812.2
在核物理研究、反应堆控制、加速器技术及空间探测、核辐射防护等领域，需精确了解辐射场内中子能谱分布。

液体闪烁体探测器由于具有较强的n／γ甄别能力、较高的探测效率及较好的光输出能力，被广泛应用于中子能谱的测量。

基于液体闪烁体探测器求解中子能谱本质上是核物理领域的一个反演问题，根据求解原理的不同，中子解谱的方法主要有最小二乘法［1］、蒙特卡罗法［2］、正则化法［3］、遗传算法［4］和神经网络法［5］等。

这些方法虽适用于一般情况下的解谱问题，但常存在算法参数不易选取、精度较低、求解耗时等不足。

因此，寻找新的更加有效的解谱方法一直是中子解谱领域的研究热点。

本文提出一种基于信赖域算法的中子解谱方法，该方法无需中子能谱的先验信息，对参数设置依赖较小，能较快实现对中子能谱的准确求解，本工作欲通过仿真及实验验证该方法的有效性。

1 中子解谱原理
液体闪烁体探测器测量得到的脉冲幅度谱Ni（H）是探测器固有响应函数Ri（H，E）与入射中子能量分布Φ（E）的卷积，则：
其中：Ri（H，E）dE定义为从E到E＋dE能量范围内的中子产生幅度为H的脉
冲的微分几率。

为方便使用计算机进行数据处理，实际应用中常利用多道分析器获得脉冲幅度谱的离散值，式（1）的离散形式为：
其中：Ni为多道分析器第i道的计数；Φj为第j个能量间隔的辐射强度；Rij为第
i个脉冲幅度间隔与第j个能量间隔耦合的响应矩阵。

为方便计算，常将式（2）表示为矩阵形式：
中子解谱的目的即为由探测器测量得到的脉冲幅度谱N和探测器固有响应矩阵R
求解出中子能谱Φ。

由于探测器固有响应矩阵R通常是一病态矩阵，R或N中的
微小误差常会导致求得的中子能谱Φ产生较大波动，从而与真实值相差甚远。

因此，如何准确求解式（3）所示的病态方程组是中子解谱研究所要解决的核心问题。

2 响应函数计算
探测器响应函数可通过蒙特卡罗方法仿真计算得到。

中子从中子源辐射出去到进入探测器内部，其与探测器内部及周围物质相互作用过程可通过蒙特卡罗方法仿真，这些反应产生的光输出同样可进行理论模拟。

常用的仿真液体闪烁体探测器响应函数的程序有O5S［6］、MCNP［7］、SCINFUL［8］、NRESP［9］和 GEANT ［10］等，本实验验证中所用BC501A型探测器与参考文献［11］中所用的是同一探测器，该探测器响应函数由O5S程序仿真得到。

在入射中子能量小于10MeV 时，1 H（n，n）1 H反应为主，入射中子能量大于10MeV时，12C（n，α）9Be和12 C（n，3α）反应变得越来越重要。

仿真过
程中12 C（n，α）9Be和12 C（n，3α）的反应截面和角分布由ENDF／B4和
ENDF／B5评价核数据库得到，通过与实验结果的对比进行了一些修正。

图1所示为使用O5S程序仿真φ5．08cm×5．08cm的BC501A型液体闪烁体探测器响应函数，其中，中子能量范围为1～16MeV，能量间隔为1MeV。

图1 O5S程序仿真BC501A型液体闪烁体探测器对1～16MeV入射中子的响应函数Fig．1 Response function generated by O5S for BC501Aliquid scintillation detector in energy range of 1－16MeV
为检验O5S程序仿真得到的探测器响应函数的准确性，利用Cockcroft－Walton 加速器产生的约为15MeV的单能中子与仿真结果进行了比较，图2所示为实验测量的脉冲幅度谱与仿真结果的对比。

从图2可看出，二者在整个中子测量能段内符合得较好。

在低能段，实验结果与仿真结果相比稍高一些，主要是由于实验测量结果中未扣除由于大厅散射所引起的散射本底的影响。

3 信赖域算法原理
对于基于脉冲幅度谱和探测器固有响应函数求解中子能谱，由于响应函数通常是一病态矩阵，因此上述问题是一典型的不适定问题。

处理该类问题最有效的方法是正则化方法，本文采用的信赖域算法已被充分证明是一种正则化方法［12－13］，且已在大气图像重建、遥感参数反演等地球物理反问题中获得成功的应用［14－15］。

信赖域算法的基本思想是在每次迭代时强制性地要求新的迭代点与当前迭代点之间的距离不超过某一控制量，即在以当前迭代点为中心的1个邻域内对1个近似于原问题的简单模型进行求解，从而有效防止因步长过大导致算法失败的情况发生。

利用信赖域方法反演中子能谱正是充分利用了算法的这一特点，有效地实现了对中子能谱的准确求解。

图2 中子源脉冲幅度谱实验结果与仿真结果对比Fig．2 Comparison of measured and calculated pulse height distribution of neutron source
对于求解病态线性方程组RΦ＝N，可转化为求解如下的优化问题：
其中，F［Φ］为将原问题转化为优化问题后的目标函数。

将该优化问题转化为信
赖域子问题形式如下：
其中：J［s］为转化为信赖域子问题后的目标函数；gk＝grad（F［Φk］）＝RTRΦk－RTN为目标函数F［Φ］的梯度；Bk＝Hess（F［Φk］）＝RTR为目标函数F［Φ］的Hessian阵近似矩阵；s为每次迭代时的最优迭代值；Δ为信赖域
半径。

使用信赖域算法求解问题步骤如下。

1）给定初始值Φ1非负，Δ1≥0，0＜η1≤η2＜1，0≤γ1≤1≤γ2≤γ3，并计算F ［Φ1］，令k∶＝1。

2）若满足某种停机准则，则停止迭代；否则，求解子问题（式（5））得到sk。

3）计算F，其中，Aredk＝F［Φk＋sk］－F［Φk］，Predk ＝Jk［0］－Jk ［sk］，更新方程解：
同时更新信赖域半径：
4）计算gk＋1，令k∶＝k＋1，转步骤2。

本文采用文献［16］中提出的带信赖域技巧的截断共轭梯度法求解步骤2中信赖
域子问题，子问题的求解步骤如下。

1）给定初始点s1＝0，计算g1＝g（s1），d1 ＝－g1 ＝－grad（F［Φ1］）。

2）若‖R（Φk＋sk）－N‖ ≤τ‖R（Φk）－N‖，0＜τ＜1为常数，则停机；否则，计算，转步骤4；否则，计算
3）若‖sk＋akdk‖ ≥Δk，转步骤4；否则，计算sk＋1 ＝sk＋akdk、gk＋1 ＝
gk＋akBkdk、βk ＝，令k∶＝k＋1，转步骤2。

4）计算ak≥0，使其满足‖sk＋akdk‖＝Δk，令＝sk＋akdk，停机。

步骤1中需用户指定参数ηi（i＝1，2）及γj（j＝1，2，3），一种典型的取法为η1 ＝0．6，η2 ＝0．6，γ1 ＝0．8，γ2 ＝1．5，γ3 ＝2［17］。

4 解谱方法的仿真与实验验证
4．1 仿真验证
核物理实验耗资巨大、操作复杂，且带有一定的危险性，因此在对解谱方法有效性进行实验验证之前，需先对其进行仿真验证，以有效地降低各种风险。

实验装置为φ5．08cm×5．08cm的BC501A型液体闪烁体探测器，通过O5S程序成功仿真了能量范围为0．25～16MeV、能量间隔为0．25MeV的探测器响应函数。

利用计算出的探测器固有响应函数，分别对单峰和多峰情况下的单能入射中子进行仿真验证。

仿真过程中验证中子计数率归一化后单峰情况下的单能入射中子，图3所示为入射中子能量为2．5MeV时解谱值与理论值的对比。

由于实际应用中的能谱测量常涉及较宽的能量范围，因此，仿真过程同样对中子计数率归一化后多峰情况下的单能入射中子进行了验证，图4所示为入射中子能量分别为1．25、3、5、7．5和11．25MeV时解谱值与理论值的对比结果。

图3 单峰入射中子情况下解谱值与理论值的对比Fig．3 Comparison of unfolded and theoretic spectra for incident neutron with single－peak
图4 多峰入射中子情况下解谱值与理论值的对比Fig．4 Comparison of unfolded and theoretic spectra for incident neutron with multi－peak
在空间粒子辐射环境、反应堆内部及天然放射性物质辐射粒子中，中子能量分布不仅具有较宽的范围，且各能量处通常具有不同的中子计数率。

针对这种情况，仿真过程中对中子计数率与能量成反比、同时含有多个峰的单能入射中子情况进行了验证。

图5所示为入射中子能量分别为1、3．5、6、7．75和12．5MeV，中子计
数率与入射中子能量成1／E关系时解谱值与理论值的对比结果。

图5 入射中子计数率与能量成反比时解谱值与理论值的对比Fig．5 Comparison of unfolded and theoretic spectra for incident neutron with count rate inversely proportional to neutron energy
从上述仿真验证的结果中可看出，使用信赖域算法求得的结果与理论值符合得很好，这为解谱方法的实验验证提供了可靠的保证。

4．2 实验验证
解谱方法实验验证的数据测量在中国工程物理研究院核物理与化学研究所完成。

实验所用探测器为φ5．08cm×5．08cm 的BC501A型液体闪烁体探测器，
通过硅油与英国EMI公司生产的9807B型光电倍增管直接耦合，实验测量硬件电路图如图6所示，所用器件除探测器及光电倍增管外均来自于ORTEC公司。

图6 实验测量硬件电路图Fig．6 Hardware circuit diagram of experimental measurement
实验过程中分别对D－T单能中子源、241 Am－Be连续中子源和252 Cf连续中子源进行验证。

D－T中子源是由Cockcroft－Walton加速器产生的能量约为
15MeV的单能中子，出射中子经过准直器后，被屏蔽体包围的BC501A探测器探测到。

由于屏蔽体的屏蔽作用，使得环境本底与中子信号相比贡献很小，通常忽略不计。

图7为对D－T中子源测量得到的脉冲幅度谱解谱得到的结果。

对其峰值进行高斯拟合，得到高斯分布的均值为14．984 1MeV、方差为0．369 2MeV、半高宽为0．801 2MeV，与预期的15MeV 十分吻合，且能量展宽不明显，解谱结果中的微小波动主要是受环境本底的影响。

图7 使用信赖域算法反演D－T中子源得到的解谱值Fig．7 Unfolded result of
D－T neutron source by trust region algorithm
为减小环境本底的影响，两种连续中子源均采取悬空测量方式，其中241 Am－
Be中子源放置在距探测器中心115cm、距地面175cm的位置，252Cf中子源放置在距探测器中心100cm、距地面150cm的位置。

环境本底的信号是通过在探测器与中子源之间放置金属棒测量得到的。

图8所示为使用扣除本底后的脉冲幅度谱反演得到的241 Am－Be中子源能谱与ISO标准241 Am－Be中子源能谱的对比结果。

由图8可看出，尽管实验所用中子源与ISO标准中子源在制作工艺和使用时间长短方面不同，但使用信赖域算法求得的结果与标准241 Am－Be中子源能谱具有位置相同的峰值分布，且整体上趋势一致。

图9所示为使用信赖域算法反演得到的252Cf能谱与参考值［18］的对比结果。

标准252 Cf能谱曲线近似于麦克斯韦裂变谱，满足表达式：
其中：C为归一化因子；T为一常数。

参考文献［18］中参考值T＝1．42，采用本文提出的信赖域算法拟合后求得的结果为T＝1．414 7，两者符合较好。

从图9中可看出，本文提出的方法求得的结果比参考值略偏低，原因可能是受到响应函数求解精度的限制，导致低能段求解不准，归一化后影响了整体趋势及拟合效果。

图8 使用信赖域算法反演241 Am－Be中子源得到的解谱值与参考值的对比Fig．8 Comparison of 241 Am－Be neutron source energy spectrum unfolded by trust region algorithm and ISO reference spectrum of 241 Am －Be neutron source
图9 使用信赖域算法反演252 Cf中子源得到的解谱值与参考值的对比Fig．9 Comparison of 252 Cf neutron source energy spectrum unfolded by trust region algorithm and reference spectrum of 252 Cf neutron source
5 结论
本文基于信赖域算法开展了中子能谱测量相关研究，使用O5S程序仿真了
BC501A型液体闪烁体探测器响应函数，并利用该响应函数仿真验证了各种情况下的中子解谱问题。

通过实验验证对 D－T单能中子源、241 Am－Be和252 Cf 连续中子源进行了求解，与参考值的对比结果表明了该方法的有效性。

感谢中国工程物理研究院核物理与化学研究所陈渊教授和言杰博士的帮助，同时感谢该所中子加速器组工作人员给予的积极配合。

［1］ PEREY F G．Least－squares dosimetry unfolding：The program STAY’SL，ORNL／TM－6062［R］．Oak Ridge，USA：Oak Ridge National Laboratory，1977．
［2］ SANNA R，OBRIEN K．Monte－Carlo unfolding of neutron spectra ［J］．Nuclear Instruments and Methods A，1971，91（4）：573－576．［3］ DEHIMI F Z，SEGHOUR A，ABAIDIA S H．Unfolding of neutron spectra with fisher regularization［J］．Nuclear Science，2010，57（2）：768－774．
［4］ FREEMAN D W，EDWARDS D R，BOLON A E．Genetic algorithms：A new technique for solving the neutron spectrum unfolding problem ［J］．Nuclear Instruments and Methods A，1999，425（3）：549－576．［5］ SHARGHIIDO A，BONYADI M R，ETAATI G R，et al．Unfolding the neutron spectrum of a NE213scintillator using artificial neural networks ［J］． Applied Radiation and Isotopes，2009，67（10）：1 912－1 918．［6］ TEXTOR R E，VERBINSKI V V．O5S：A Monte Carlo code for calculating pulse height distributions due to monoenergetic neutron incident on organic scintillators，ORNL－4160［R］．Oak Ridge，USA：Oak Ridge National Laboratory，1968．
［7］ BRIESMEISTER J F． MCNPTM： A general Monte Carlo N－particle
transport code，LA－13709－M［R］．Los Alamos，USA：Los Alamos National Laboratory，2000．
［8］ DICKEN J K．SCINFUL：A Monte Carlo based computer program to determine a scintillator full energy response to neutron detection for Enbetween 0．1and 80MeV，ORNL－6462［R］．Oak Ridge，USA：Oak Ridge National Laboratory，1988．
［9］ DIETZE G，KLEIN H．O5S4and NEFF4：Monte Carlo codes for the calculation of neutron response functions and detection efficiencies for NE213scintillation detectors，PTB－ND－22［R］．Braunschweig，Germany： Physikalisch Technische Bundesanstalt，1982．
［10］PATRONIS N，KOKKORIS M，GIANTSOUNDI D，et al．Aspects of GEANT4Monte Carlo calculations of the BC501Aneutron detector ［J］．Nuclear Instruments and Methods A，2007，578（1）：351－355．［11］YAN J，LIU R，LI C，et al．Application of artificial neural networks for unfolding neutron spectra by using a scintillation detector ［J］．Science China： Physics， Mechanics ＆ Astronomy，2011，54（3）：465－469．
［12］WANG Y F，YUAN Y X．Convergence and regularity of trust region methods for nonlinear illposed inverse problems［J］．Inverse Problems，2005，21（3）：821－838．
［13］WANG Y F．On the regularity of trust region－cg algorithm，with application to image deconvolution problem［J］．Science in China A，2003，46（3）：312－325．
［14］WANG Y F，YUAN Y X，ZHANG H C，et al．A trust region－CG
algorithm for deblurring problem in atmospheric image reconstruction ［J］．Science in China A，2002，45（6）：731－740．
［15］WANG Y W，WANG Y F，XUE Y，et al．A new algorithm for remotely sensed image texture classification and segmentation ［J］．International Journal of Remote Sensing，2004，25（19）：4 043－4 050．
［16］王彦飞．反演问题的计算方法及其应用［M］．北京：高等教育出版社，2007：62－66．
［17］WEN Z W，WANG Y F．A new trust region algorithm for image restoration［J］．Science in China A，2005，48（2）：169－184．［18］HAMBSCH F J，TUDORAB A，VLADUCAB G，et al．Prompt fission neutron spectrum evaluation for 252 Cf（SF）in the frame of the multimodal fission model［J］．Annals of Nuclear Energy，2005，32（10）：1 032－1 046．。