由于放射性核素会对人体健康造成极大威胁，不当的使用会造成严重的社会恐慌。核素识别是有效鉴定危险放射性物质的手段，在当今社会有着极为重要的应用，是机场、海关、铁路以及公共区域等场所安检过程中不可或缺的环节。利用放射性核素的γ能谱对核素进行识别是一种常见的手段，传统的识别方法有γ能谱的全能峰面积法、逐次差引法以及逆矩阵法^[1]。然而，传统的核素识别方法存在计算复杂、耗费时间较长等问题。

近年来，随着人工神经网络的发展，越来越多的学者应用基于神经网络的方法进行核素鉴定。弟宇鸣等^[2]通过改变测量时间和测量距离获得¹⁵²Eu、¹³³Ba、⁶⁰Co、¹³⁷Cs的γ能谱训练数据，针对不同核素特征采取不同数据预处理方法，再使用反向传播(back-propagation, BP)神经网络进行识别。史东生等^[3]利用粒子群算法优化BP神经网络，识别效果有较大的提升。王崇杰等^[4]利用径向基神经网络对¹⁵²Eu、²²Na、⁶⁰Co、¹³⁷Cs进行识别与活度定量分析，误差控制在3%内。刘议聪等^[5]收集了7种放射性核素γ能谱数据，利用主成分分析算法对γ能谱进行降维后，作为神经网络的输入数据，对比了BP网络和径向基网络训练结果，发现径向基网络效果更佳。

神经网络在γ能谱识别中的应用降低了传统γ能谱识别工作的复杂性。然而，神经网络需要大量的训练数据，若要识别更多种类的核素，γ能谱数据收集和预处理时间将会非常漫长。针对此问题，另一类广泛应用在核素识别领域中的方法更具优势，此类方法首先建立核素γ能谱模板库，再通过线性回归拟合未知核素，分析每个核素的相对贡献率和误差，判断某种核素是否存在。Roemer等^[6]利用Geant4模拟γ能谱用于手提式探测器识别核素，认为模拟γ能谱用于识别未知核素是有效可行的。Hogan等^[7]采用多元线性回归识别γ能谱，并提出如何自动地选取候选集结果。Carlevaro等^[8]利用遗传算法优化线性拟合系数，其计算结果可以较好地拟合实际核素活度。Alamaniotis等^[9]等利用进化算法多目标优化实现对未知γ能谱的识别，采用Pareto最优原则，返回最有可能的核素结果。然而，上述这些求解方法在高维情况下，解向量容易陷入局部最优，识别出错的概率较大。

本文尝试建立由26种常见的放射性核素组成的γ能谱模板数据库，在此基础上，利用最小二乘法和带约束条件的非线性规划算法，在各自的作用域分别寻找最优解，从而充分利用两种算法所得结果的信息，综合研判分析γ能谱，进而识别出未知γ能谱可能存在的核素。

1 建立核素γ能谱模板数据库 1.1 数据获取

国际原子能机构《核安保丛书》第5版^[10]列出了4种类型常见的放射性核素，共26种(见表 1)，本文尝试建立由这26种核素组成的γ能谱模板数据库。

26种放射性核素的γ能谱基于Geant4^[11]平台，通过Monte Carlo模拟生成。模拟的γ探测器采用敏感体积为4 cm×4 cm×4 cm的NaI晶体材料，其中心位于坐标系原点，放射源位置为距探测器1 m的z轴负方向，模拟的空间分布示意图如图 1所示。放射源使用Geant4 Particle Gun自定义核素，模拟衰变粒子总数为10⁷，相当于活度0.99×10⁷ Bq的放射源，γ能谱模拟后添加Gauss分布展宽N(0, 0.03)以模拟探测器的能量分辨率响应。

模拟获得每一个放射性核素γ能谱，横坐标为1~2 000 keV的能量分布，纵坐标为单位能量下的粒子计数，并对能谱进行标准化，用于后续数据分析。标准化的方法如下：

$\int_1^{2\;000} {C\left( E \right){\rm{d}}E = 1} . $

(1)

其中：E代表能量，C(E)代表对应能量下的粒子计数。通过式(1)可以将粒子总数标准化为1。

1.2 RBF神经网络拟合γ能谱

径向基函数(radial basis function, RBF)有着较好的学习能力和归纳能力，RBF神经网络结构简洁，收敛速度快，逼近性能好，在很多方面性能优于BP神经网络^[12]，因此RBF神经网络常用来逼近非线性连续函数。RBF神经网络通常有3层，即输入层、隐含层和输出层，其基本的网络结构如图 2所示。

RBF神经网络隐含层的激活函数采用Gauss核函数，其形式为

${\rm{g}}\left( {E, {\rm{ }}c, {\rm{ }}\sigma } \right) = {{\rm{e}}^{-\frac{{ {{\left\| {E - c} \right\|}^2}}}{{2{\sigma ^2}}}}}. $

(2)

式(2)中：‖E-c‖²可看作输入值能量E与c的Euclid距离的平方；c为Gauss核函数中心，决定核函数峰值的位置；σ控制径向基的作用范围。以一维Gauss核函数为例，图 3所示展示了不同c、σ所作用的区域。

输出层为隐含层的线性叠加求和，即 $f(x) = \sum\limits_{i = 1}^m {{w_i} \cdot {g_i}(x)} $，m为隐含层神经元个数。RBF网络通过反向传播机制，不断缩小输出层值和真实值之间的差距，损失函数如式(3)所给出，通过降低损失值，不断优化参数c、w、σ，以实现函数的拟合。

$f(E) = \sum\limits_{i = 1}^m {w \cdot g(E, c, \sigma )} , $

(3)

$L(c, w, \sigma ) = \frac{1}{{2\;000}}\sum\limits_{j = 1}^{2\;000} {{{[f({E_j}) - {C_j}]}^2}} . $

(4)

式中：m为隐含层神经元个数，w为隐含层到输出层权重集合，C_j为待拟合γ能谱的第j道标准化后的粒子数。

利用上述方法，本文针对每一种核素，分别训练一组RBF神经网络，共获得26组RBF神经网络。本文训练RBF神经网络基于MATLAB平台完成，采用newrb函数，以γ能谱的能量作为神经网络的输入，以与之对应的标准化后的粒子计数作为目标值，网络输出拟合的γ能谱计数值，并计算输出值与目标值的均方差。训练过程中通过不断迭代增加隐含层神经元个数，使RBF神经网络能更好地拟合真实数据，减小均方差。

经多次实验后，将神经网络训练误差目标设定为10^-8，使得RBF网络拟合值不会偏离实际值过多，可保证特征峰形状与位置的准确性。另外，实验发现，设定隐含层神经元个数最大值为100时，可满足大多数核素误差要求，同时也能起到很好的压缩数据的作用。若神经网络隐含层神经元个数达到设定上限时，仍未能满足误差目标，则需调整径向基函数分布范围σ，再重新训练。RBF神经网络训练流程图如图 4所示。

将26组RBF仿真结果汇集于表 2，核素γ能谱拟合误差几乎都在10^-8数量级及以下，可用于本文核素识别方法。

通过RBF神经网络拟合核素γ能谱建立核素模板库，可压缩数据，大大减少数据本身所占用的储存空间。就手提式探测器而言，该类探测器对存储空间有着严苛的要求，难以实现对多种核素γ能谱模板数据的录入，使用RBF神经网络拟合可有效解决此问题，并且具有良好的扩展性。图 5列出4种常见核素γ能谱的RBF神经网络的拟合结果。

2 核素识别模型与算法原理 2.1 模型与评价

本文通过上述方法建立了26种核素的γ能谱模板库。对待分析的目标γ能谱T，采用线性叠加的方法进行识别^[8]，其数学表达式为

$\mathit{\boldsymbol{T}} = \sum\limits_{i = 1}^m {{p_i}{\mathit{\boldsymbol{R}}_i}} . $

(5)

其中：m为核素的种类，R_i为第i个核素γ能谱模板，p_i为第i个核素γ能谱的贡献率，其值越大表示越可能存在该核素。对此，本文采用下述两种方法求解p_i，并分析判断目标γ能谱T可能存在的核素种类。

在统计学中，常用决定系数^[13]度量模型拟合的好坏，如式(6)所示。该值越接近与1，表明拟合回归效果越好。针对本文问题，采用此指标判断模型的解释力。

${R^2} = 1 - \frac{{\sum\limits_i {{{({y_i} - {f_i})}^2}} }}{{\sum\limits_i {{{({y_i} - \bar y)}^2}} }}. $

(6)

其中：y_i代表真实值，y代表真实值的均值，f_i代表拟合值。

2.2 最小二乘法

最小二乘法是一种常见的数学优化模型^[14]，通过最小化误差平方，寻找出与实际数据最佳的数据匹配组合，其数学表达式为

$\min\sum\limits_{i = 1}^n {{{({{\hat y}_i} - {y_i})}^2}} . $

(7)

根据最小二乘法的数学思想，推导出式(8)，用于求解向量p。

$\mathit{\boldsymbol{p}} = {({\mathit{\boldsymbol{R}}^{\rm{T}}}\cdot\mathit{\boldsymbol{R}})^{ - 1}}\cdot{\mathit{\boldsymbol{R}}^{\rm{T}}}\cdot\mathit{\boldsymbol{T}}. $

(8)

其中：矩阵R_{2 000×26}为已建立的γ能谱模板数据库，2 000为γ能谱的能量维度，26为核素种类个数；T_{2 000×1}为待分析的目标γ能谱；向量p_26×1为最终求解结果，即各核素的拟合系数。

2.3 非线性规划

非线性规划求解的是由一系列未知实函数组成的方程组和不等式定义的最优化问题^[15]。由于本文求解问题的目标函数是非线性函数，故采用带约束条件的非线性规划算法优化向量p，其数学表达式见式(9)。但是，高维空间下非线性规划在寻优过程往往会陷入局部最优解，导致无法获取全局最优解。因此，本文首先采用粒子群优化算法进行全局寻优，将粒子群算法的最优解作为非线性规划的初始解，再进行非线性规划寻优，该方法可减少陷入局部最优解的概率。

$\begin{array}{c} \min\int_1^{2\;000} {[\mathit{\boldsymbol{T}}(E) - \sum\limits_{i = 1}^m {{p_i}{\mathit{\boldsymbol{R}}_i}(E){]^2}} {\rm{d}}E} , \\ {\rm{subject}}\;{\rm{to}}\;\sum\limits_{i = 1}^m {{p_i}} = 1, {p_i} \ge 0. \end{array} $

(9)

粒子群优化算法最早由Kennedy^[16]提出，该算法利用粒子群之间的信息分享使粒子找到最优区域，从而避免了因初始位置选择不当使算法陷入局部最优解。

算法步骤为：1) 设定适应度函数为式(9)，并初始化一批粒子，每个粒子都有速度v_i和位置x_i两个属性，针对本文所要求解问题，每个属性维度均为26维，位置x_i的每个维度的意义表示对应核素的贡献率。2) 计算每个粒子的适应度函数值，找到初始个体最优适应值pbest_i和群体最优适应值gbest，再通过式(10)和(11)进行迭代，完成设定的迭代次数后，以粒子群的gbest值所对应的粒子位置作为非线性规划的初始解，再进行非线性规划求解。

$\begin{array}{c} v_{id}^k = w \cdot v_{id}^{k - 1} + {c_1}{\eta _1}({\rm{pbes}}{{\rm{t}}_{id}} - x_{id}^{k - 1}) + \\ {c_2}{\eta _2}({\rm{gbes}}{{\rm{t}}_d} - x_{id}^{k - 1}), \end{array} $

(10)

$x_{id}^k = x_{id}^{k - 1} + v_{id}^{k - 1}. $

(11)

其中：v_id^k表示第k次迭代粒子i速度矢量第d维分量，x_id^k表示第k次迭代粒子i位置矢量第d维分量。c₁、c₂为加速度常数；w为惯性权重，其值越大表示全局搜索能力越强；η₁、η₂为(0, 1)范围内的随机数。具体算法流程如图 6所示。

2.4 模型集成

多组实验测试结果表明，最小二乘法对特征峰较为敏感，而非线性规划对Compton平台较为敏感。因此，本文采用集成学习的思想，利用加权平均法将线性模型(8)和非线性模型(9)结合起来，综合二者的优点，从而有效地提升模型的泛化能力。

利用模型的决定系数确定二者的权重，最小二乘法所占权重由式(12)给出，同理计算非线性规划结果权重。

${w_1} = \frac{{R_{\rm{L}}^2}}{{R_{\rm{L}}^2 + R_{\rm{N}}^2}}. $

(12)

式(12)中：R_L²表示最小二乘法的决定系数，R_N²表示非线性规划的决定系数。集成模型的计算结果可表示为 $\sum\limits_{i = 1}^2 {{w_i} \cdot {p_i}} $，并通过分析多组核素γ能谱的计算结果，确定出核素相关系数阈值，利用此阈值判定未知核素γ能谱可能存在的核素种类。

3 核素识别结果分析

为验证上述方法可行性，首先通过Geant4模拟生成不同条件下单一放射性核素或多种放射性核素混合的5组γ能谱，再利用上述算法分别对其求解识别，并通过绘制求解结果的γ能谱图，与模拟核素γ能谱比较，进一步分析求解结果，确定出识别核素阈值；最后通过多组实验验证本文方法的准确性。

3.1 单核素的γ能谱识别分析

第1组核素γ能谱是利用Monte Carlo方法模拟¹³³Ba放射得到的。先对能谱进行标准化，然后利用本文算法在建立起的核素γ能谱模板库中分别寻找最优解，将计算结果按数值由大到小顺序记入表 3。

根据表 3结果，有充分理由相信核素γ能谱1是由核素¹³³Ba衰变得到的。进一步将3种模型的计算结果进行绘图验证，如图 7所示。

根据图 7拟合结果可进一步确定算法的识别结果与实际结果相吻合，识别正确。

3.2 双核素叠加的γ能谱识别分析

第2组γ能谱是模拟核素⁴⁰K和¹³³Ba叠加放射的结果，其活度比为5∶1。同样利用本文方法对γ能谱2进行求解识别，得到表 4所示的计算结果。

计算结果显示，该核素γ能谱最有可能为⁴⁰K和¹³³Ba两种核素的叠加，其余核素存在的可能性较小。进一步绘图验证，如图 8所示。

通过分析图 8中特征峰位置和形状的拟合程度，可判定该γ能谱2识别的结果为核素⁴⁰K和¹³³Ba。

第3组γ能谱模拟的是核素¹³⁷Cs和¹⁹²Ir叠加放射的结果，其活度比为2∶1。下面对核素γ能谱3利用本文方法求解分析，计算结果如表 5所示。

根据表 5结果，有理由认为该γ能谱是¹³⁷Cs和¹⁹²Ir两种核素衰变的结果，通过进一步绘图分析其特征峰拟合情况，验证计算结果。图 9表明，随着核素¹³⁷Cs比例增加，计算结果可更好地拟合662 keV左右的特征峰。低能段的偏差可能是由RBF网络仿真过程造成的，但从计算结果和特征峰位置的吻合程度来看，依旧可正确地识别出核素。

3.3 3种核素叠加的γ能谱识别分析

第4组γ能谱模拟的是¹³⁷Cs、¹⁹²Ir和¹⁸ F 3种核素叠加放射的结果，其活度比为2∶1∶1。对核素γ能谱4进行求解分析，计算结果如表 6所示。

进一步绘制出3种计算结果的拟合曲线，分析该γ能谱的特征，如图 10所示。可见，非线性规划模型未能识别出在670 keV附近的特征峰，而最小二乘对平台段拟合不敏感。因此，综合二者优势，集成模型有较好的鲁棒性，可正确识别出γ能谱4所对应核素。

第5组核素γ能谱模拟的是¹⁸F、⁶⁰Co和¹⁵²Eu 3种核素的放射结果，其活度比为1∶2∶1。利用本文方法对γ能谱5进行求解分析，结果见表 7。

根据非线性规划计算结果，该γ能谱可能是¹⁸F、¹⁵²Eu和⁴⁰K 3种核素的衰变结果；根据最小二乘法和集成模型计算结果，该γ能谱是¹⁸F、⁶⁰Co和¹⁵²Eu 3种核素的衰变结果。进一步绘图分析两种可能性，如图 11所示。

图 11结果显示，在1 170 keV和1 340 keV附近有两个微弱的特征峰，为核素⁶⁰Co放射产生的特征峰，而非线性规划计算结果未能与之匹配。因此，可判定核素γ能谱5的识别结果为¹⁸F、⁶⁰Co和¹⁵²Eu 3种核素。

通过比较上述求解分析所识别出的核素结果与实际模拟的核素结果可以发现，在对不同种类核素叠加放射下的γ能谱分析识别时，利用本文方法均可较为准确地识别出核素种类。由于在模型计算前需要对待检测γ能谱进行归一化处理，因此模型的计算结果与核素活度比例系数无明显联系，模型计算结果不能代表核素活度的相对大小。

3.4 模型识别结果验证

集成模型综合了线性模型和非线性模型的计算结果，具有一定的鲁棒性和泛化性。因此，采用集成模型的计算结果进行模型验证，设定标准的识别过程。综合分析上述5组γ能谱的识别结果，将拟合系数的阈值设定为0.15，当计算结果超过该阈值时，则认为存在此核素。

为验证本文模型的准确性，另模拟了多组核素的γ能谱，将模型的识别结果汇总于表 8。

由表 8可见，利用本文模型识别未知γ能谱准确率高，尤其是对单核素和双核素γ能谱的识别准确率较高。

4 结束语

本文通过建立放射性核素γ能谱模板数据库，采用最小二乘法、非线性规划两种算法对未知核素γ能谱进行识别，两种算法在寻找最优解过程中各有优势，计算结果可相互验证、相互补充，并通过集成模型综合二者优点，增强了模型的鲁棒性与泛化性，防止漏掉最佳识别结果。

根据本文方法识别核素种类的机理，此方法对γ能谱存在偏移的情况较为敏感，因此使用此方法前应对γ能谱进行偏移修正，或对探测器进行校准，防止偏移过大造成识别错误。

利用径向基神经网络模型识别核素γ能谱会有一定的误差，但不会影响特征峰的判断。从识别结果来看，本文模型具有一定的鲁棒性。因此，此模型允许待检测的γ能谱有一定噪声存在，不必要求探测器具有很高的分辨率。

由于环境本底的影响可能会使求解结果识别成错误的核素，因此对于未知核素的γ能谱，最好先进行全谱基底扣除处理，再利用此方法识别，会提升识别的准确率。

本文所建立的核素γ能谱模板数据库针对的是较为常见的核素，如果需要识别未在模板库中的核素，需先对模板库进行扩充，再进行识别。