高春宇，汤秀章，陈欣南，范澄军，陈雁南，李雨芃，吕建友
(中国原子能科学研究院，北京 102413)
近年来，核材料的非法转移、核扩散等问题时刻威胁着国土安全，宇宙射线μ子散射成像作为一种新兴的核材料检测技术逐渐引起世界各国的重视。μ子是一种天然辐射源，没有引入额外辐照危害，同时对高原子序数物质敏感，穿透能力强，在核材料检测技术应用中具备天然优势。在集装箱、货物材料走私检测等要求时效性的应用场景，需要在较短时间内(几分钟)快速准确地检查出是否藏匿特殊核材料，从而发展了一种不需要重建图像而只依据少量μ子数据特征快速实现材料识别的方法。英国核武器研究所(AWE)与布里斯托尔大学合作开发了基于分箱聚类理论的散射成像材料识别算法，1 min可以区分10 cm×10 cm×10 cm铀块与铅块[1]。美国洛斯阿拉莫斯国家实验室(LANL)与决策公司开发了多模态被动式成像算法，可以在1 min内检测出40英尺集装箱内藏匿的核材料[2]。印度萨哈核物理研究所采用了模式识别方法，识别具有相似模式的材料，拒绝不同模式的材料[3]。目前国内的研究主要基于聚类算法。清华大学研究了基于支持向量机的聚类分析和分类器，根据散射密度的平均值和标准偏差进行物质识别[4]。中国科学技术大学通过计算待识别目标物与参考物之间散射角的灰色关联度快速识别材料[5]。
本文提出一种将卷积神经网络(CNN)应用于μ子成像材料识别的方法，以满足快速准确检测出藏匿特殊核材料的实际应用需求。利用宇宙射线μ子穿过待测材料后的散射角统计信息，根据其与材料原子序数的正相关性实现材料的分类识别。材料识别模型是主要由卷积层、池化层和全连接层构成的卷积神经网络结构，实验测量的散射角数据或蒙特卡罗模拟Geant4数据通过卷积层进行数据特征提取，池化层降低数据维度，最后在全连接层给出模型分类结果。利用Geant4搭建μ子散射成像实验模型，并在中国原子能科学研究院的CIAE-MNT宇宙射线μ子成像检测实验装置上进行验证实验，分别测试铝、铁、钨3种不同原子序数的10 cm×10 cm×10 cm样品，实验材料识别模型结构上增加卷积层数以提取更深层次的数据特征，进一步引入残差结构、特征矩阵改进网络结构来提高模型的准确度。并分别测试模拟模型与实验模型在1、3、5和10 min不同测量时间下的识别准确度。
宇宙射线μ子穿过待测样品时，会发生多次库仑散射现象，依据其散射角分布与材料原子序数呈正比的规律，实现对材料的检测识别。根据莫里埃理论[6-7]，投影到平面上的μ子的散射角θ近似服从期望为0的高斯分布，高斯分布的宽度为：
(1)
辐射长度Lrad由经验公式[8]计算：
(2)
其中：A为相对原子质量；
ρ为材料密度；
Z为材料原子序数。
定义散射密度λLrad为标准μ子(βc=1，动量为p0)穿过单位厚度物质后散射角分布的均方根，则由式(3)可知[9]：
(3)
则有：
(4)
即散射密度λLrad与材料原子序数呈正相关，可将散射密度作为材料识别的特征量[10]。PoCA算法对μ子的多次库仑散射做出了简化，假设为一次散射的结果，对经过探测区域的有效μ子重建入射、出射径迹，并计算散射角[11]。将探测区域均匀划分为若干网格，获取每个网格内所有μ子的散射角。根据散射密度的定义，得到第i个网格散射密度的估计值λi为：
(5)
2.1 实验装置
样品测试实验在中国原子能科学研究院研制的一套基于漂移管探测器的CIAE-MNT宇宙射线μ子成像检测实验装置上进行。该装置主要由漂移管μ子探测器、信号电子学、数据采集系统、成像处理软件等构成。漂移管探测器阵列共6层，由上千支漂移管互成90°垂直交叉排列构成，上部分3层探测器组用于测量入射μ子的位置信息，下部分3层用于测量μ子经过待测样品后的散射角度和位置，中间100 cm高摆放待测样品。漂移管是长250 cm、直径3 cm的铝管，中间阳极丝引出测量信号，管内充Ar/CO2混合气体，探测面积最大为220 cm×220 cm。本实验样品较小，利用100 cm×100 cm的探测面积，将10 cm×10 cm×10 cm的样品放入探测区域中心。
μ子穿过待测样品前后的位置和角度信息通过上下两组漂移管探测器组分别测量得到，并通过测量μ子穿过漂移管电离产生电子的漂移时间得到其精确位置信息。μ子径迹重建包括径迹寻找和径迹拟合两步，径迹寻找是搜索一定时间窗口内漂移管阵列记录的击中信号的集合，通过一定的规则进行判断，并挑选出μ子径迹信息；
径迹拟合采用基于最小距离平方和拟合二维投影面上的径迹，使径迹到各层信号位置的距离最小。根据重建的μ子入射、出射径迹，采用PoCA算法计算得到μ子的散射角分布，依据散射角均方根与材料原子序数的对应关系得到待测样品的空间分布。
2.2 仿真模拟
采用Geant4程序进行宇宙射线μ子散射成像的蒙特卡罗模拟[12]，根据实验装置的实际尺寸构建物理模型，模拟高能粒子与探测器和被探测材料的相互作用，跟踪粒子的输运过程，显示粒子径迹，如图1所示。模型由上下两部分共6层探测器阵列构成，上部分和下部分各3组探测器分别用来记录入射和出射μ子径迹，每层探测器高度间隔30 cm，探测面积100 cm×100 cm，上下两部分探测器之间为探测区域，间距高度100 cm。选择铝、铁、钨作为低中高原子序数的代表材料，将10 cm×10 cm×10 cm的铝、铁、钨块样品放置在探测区域中心，与实验时的摆放位置一致。
图1 Geant4模拟实验装置结构Fig.1 Experimental facility structure simulated by Geant4
CRY软件库由美国劳伦斯利弗莫尔国家实验室(LLNL)研发，可链接到Geant4程序，根据模拟的需要设置生成不同时间、不同纬度与不同海拔高度上的宇宙射线粒子簇射分布。使用CRY软件库引入海平面宇宙线μ子能谱及角分布源项，以更真实地模拟天然μ子源[13]。考虑到地球海平面处μ子的平均注量率为10 000 m-2·min-1[14]，模拟显示在100 cm×100 cm的探测区域上，1 min约有1 600个μ子能穿过6层探测器，用于重建有效径迹。因此，在本模型中，整个探测区域上1 600条μ子散射角数据等效于1 min的测量时间，依次模拟不同的测量时间，输出各层探测器记录的μ子位置信息。
3.1 用于模拟的卷积神经网络模型
模拟模型定义一个包含两层卷积层、两层池化层和1层全连接层的卷积神经网络结构，如图2所示。网络的输入数据由PoCA算法对Geant4模拟数据进行计算得到，为50×50的散射密度矩阵。
图2 模拟条件下基于卷积神经网络的材料识别模型Fig.2 Material discrimination model based on convolutional neural network under simulation condition
第1层卷积层定义64个卷积核，第2层卷积层定义32个卷积核，卷积核大小均为3×3，卷积层通过卷积核与输入数据进行卷积计算，以扫描的方式对输入数据进行特征提取。使用非线性激活函数ReLU函数，对卷积层输出结果做非线性映射，定义如下：
f(x)=max(0,x)
(6)
ReLU函数保留大于等于0的值，其余所有小于0的值直接改写为0，以修正线性单元。池化层使用3×3的池化窗口扫描特征图，采用最大值池化函数选取池化窗口中所有值的最大值，进一步降低数据维度，输入数据中微小的变化、冗余将不会改变池化的输出，使网络更加具有健壮性。经过多层卷积层和池化层的处理后，50×50的散射密度矩阵被抽象成高阶特征，由全连接层对提取的特征进行非线性组合以得到预测输出，给出分类结果。由于线性模型的表达能力不够，本文加入了非线性因素将卷积层输出结果做非线性映射，计算得到概率矩阵[P铝,P铁,P钨]，并进一步根据预测值与真实值计算模型误差，使用交叉熵损失函数进行反向传播，更新网络权重，优化模型[15]。
3.2 实验条件下的卷积神经网络模型
散射密度是材料识别的判据，根据宇宙射线μ子成像检测实验装置测得的μ子位置及角度信息拟合μ子入射、出射径迹，进一步计算统计信息得到散射密度矩阵。为增强散射数据本身的特征，更多地保留原始数据的信息，扩大数据维度，将100 cm×100 cm探测区域划分为100×100的网格，每个网格大小为1 cm×1 cm。对整个探测区域计算散射密度得到100×100的散射密度矩阵。为了与模拟模型进行定量对比，对于探测区域，同样以1 min的测量时间等效于1 600条有效μ子位置信息。分别基于不同测量时间下实验样品铁和钨立方体材料区域的散射密度矩阵，计算散射密度大于0的像素数量的概率密度，并绘制正态分布曲线，如图3所示。可看出，测量时间为1 min时，非零数据点较少，矩阵提供的有效数据有限，这给分类模型的特征提取工作带来了难度。测量时间为10 min时，μ子的散射角数据较多，铁块、钨块的散射密度矩阵内大于0的像素数量增加，呈现出较好的分布。同时，相对于中等原子序数材料铁，重核材料钨的非零数据点更多。
a——铁；
b——钨图3 不同测量时间下材料区域内散射密度大于0的像素数量的概率密度Fig.3 Probability density of number of pixel whose scattering density is greater than zero in material area at different measurement time
神经网络的训练过程往往需要大量样本数据，而实验数据的总样本数量有限，为此进一步改进网络结构弥补数据量少的不足，本文从多角度提取了μ子的散射数据特征，包括最大值、最小值、均值、标准差，在传统结构上加入了残差结构及特征矩阵，充分利用散射角数据。同时设置测试集尽可能大，样本量是训练集的25%，以测试模型泛化能力，使得测试结果可信。
实验模型经过优化后为包含4层卷积层、4层池化层、1层全连接层和残差结构的卷积神经网络，其中卷积层从输入的散射密度矩阵中提取数据特征，获得局部感知信息，池化层减少网络参数，控制过拟合现象，相对于模拟模型结构更加复杂。网络结构如图4所示，散射密度矩阵输入到卷积神经网络中，要依次经过4层卷积-池化结构，其中第1层卷积层有8个卷积核，卷积核大小为3×3，分别与散射密度矩阵进行卷积运算以提取不同特征，输出100×100×8的特征图，并对输出结果做非线性映射。在池化层使用最大值池化函数，对特征图降维到50×50×8。同样地，在第2、3、4层进行卷积池化运算，最终输出维度为25×6×6的特征图。
图4 应用于实验环境的卷积神经网络材料识别模型Fig.4 Material discrimination model based on convolutional neural network applied to experiment
其中，经过第2层的卷积-池化操作后，考虑到可能会丢失原始散射密度矩阵X的部分数据特征，引入第2条线，加入残差结构，第2层最终输出的结果H(X)为：
H(X)=F(X)+X
(7)
其中，F(X)为第2层卷积-池化结构的运算结果。与传统卷积神经网络的输出H(X)=F(X)相比，残差结构将F(X)与原始数据X相加，相当于对输入X加入一个微小变化，使得残差在网络的反向传播过程中梯度信息更易传播。
为进一步加强数据特征，提高网络的分类准确度，引入第3条线，将提取的原始数据的特征矩阵C=[最小值，最大值，均值，标准差]与输出的特征图相加后进入全连接层。全连接层对所有得到的特征进行加权组合，计算得到概率矩阵[P铝,P铁,P钨]，其中概率最大的类别为模型分类结果。对测试集数据输入到神经网络后的识别分类过程进行追踪，以更加清晰地了解模型中每步过程的效果。图5为测试集中的1个散射密度矩阵从输入模型到给出识别结果所经过的所有操作。可看出，浅层的卷积池化结构更多的是对图像边缘提取特征，提取到的信息全面而概括，同时也会有部分关键信息提取，随着层数的加深，特征图越来越抽象，相当于对特征图的提纯过程。经过4层的卷积-池化后，100×100的散射密度矩阵被抽象成6×6的高阶特征，在全连接层输出概率矩阵[P铝,P铁,P钨]=[0,0.01,0.99]，模型预测该样本是铝的概率为0，铁的概率为0.01，钨的概率为0.99，取其中的最大值0.99，即该样本识别为钨。
图5 散射密度矩阵输入到卷积神经网络后的识别分类过程Fig.5 Classification process when scattering density matrix is input into convolutional neural network
当成像区域存在多个、置于不同位置的待测物块时，对整个探测区域的散射密度矩阵提取异常数据位置，进而输入材料识别模型进行判别。
定义材料xi被正确识别的准确度taccuracy为：
(8)
其中：xi∈X，X={铝，铁，钨}；
n为待测材料的数量，n=3。定义总体准确度ttotal为：
(9)
利用Geant4对10 cm×10 cm×10 cm的铝、铁和钨3种材料分别进行1、3、5和10 min测量时间的模拟。将每种材料识别的准确度以及总体准确度作为评价模型的指标，得到的结果列于表1。
表1 模拟条件下卷积神经网络模型的材料识别准确度Table 1 Discrimination accuracy by convolutional neural network model under simulation condition
测量时间为1 min时，μ子的散射事例数据量有限，散射密度矩阵中非零数据比例较大，导致准确度不高。由于铁的原子序数介于铝、钨之间，相对于二者来说较易混淆，在3种材料中识别准确度最低，仅为82.1%。3 min时，3种材料的识别准确度均可达到90%以上。5 min时，μ子穿过待测物体的散射角数据较多，用来训练网络的散射密度矩阵中非零网格的比例有所提高，3种材料的准确度均在98%之上，总体准确度达到98.5%。测量时间为10 min时，铝与钨块的识别准确度均已达到100.0%，铁块的误报率也仅为0.1%，总体材料识别准确度达到99.9%。统计10 min时模拟与实验条件下μ子穿过钨块的投影角分布，实验结果与模拟结果吻合良好(图6a)。进一步计算不同测量时间下模拟与实验模型对钨块识别的准确度如图6b所示，由于在实验中探测器具有一定的位置分辨率，同时网络的训练数据集较少，导致钨块的识别准确度略低于模拟结果。但是随着测量时间的增加，实验模型的准确度逐渐接近于模拟模型，差距逐渐缩小，在10 min时两个模型对钨块的准确度都达到了100.0%，由此说明了模拟与实验平台是相对可靠的。
a——投影角分布；
b——材料识别准确度图6 模拟与实验结果的对比Fig.6 Comparison between simulation and experiment results
在实验装置上，同样对10 cm×10 cm×10 cm的铝块、铁块和钨块进行1、3、5和10 min的测量实验，结果列于表2。当测量时间为3 min时，实验模型对重核材料钨的准确度达到了98.1%，由于低、中原子序数材料的散射密度较小，探测器的位置分辨率对其散射角的测量影响更大，导致铝、铁样品的识别准确度相对较低。测量时间为5 min时，3种材料的准确度均达到90%以上，总体准确度为94.6%。测量时间为10 min时，该模型对钨的准确度达到100.0%。结果表明，基于卷积神经网络的材料识别方法能在较短的测量时间内实现对低、中、高原子序数材料的检测，并有效识别重核材料钨。
表2 实验条件下卷积神经网络模型的材料识别准确度Table 2 Discrimination accuracy by convolutional neural network model under experiment condition
为研究不同层数的卷积-池化结构对神经网络准确度的影响，以实验条件下测量时间10 min获得的散射密度矩阵作为训练集，分别建立包含2层、3层、4层以及5层卷积-池化结构的网络模型，测试结果如图7所示，可看出，4层、5层卷积-池化结构能更加有效地提取到散射密度矩阵的特征，提升材料识别的准确度。此外，由于5层结构相比于4层结构准确度提升效果并不显著，且复杂的网络结构可能会过度学习训练集数据，易出现过拟合现象。因此实验模型选择使用4层卷积-池化结构。
图7 不同数量卷积结构的卷积神经网络模型的材料识别准确度Fig.7 Material discrimination accuracy of convolutional neural network models with different numbers of convolutional structures
使用机器学习中的随机森林、经典BP网络分类算法，与本文提出的卷积神经网络算法进行材料识别效能的比较(表3)，结果表明，卷积神经网络材料识别模型在μ子散射的材料识别上具有一定的优势。
表3 不同算法模型材料识别准确度的对比Table 3 Comparison of material discrimination accuracy of different algorithm models
根据宇宙射线μ子穿过待测物体时散射密度与材料原子序数的相关性，将卷积神经网络的方法应用于μ子散射成像中的材料识别。通过卷积进行特征提取，通过池化降低数据维度，不断迭代调整权重来优化神经网络模型，并分别构建了不同结构的卷积神经网络，分析网络的测试结果，选择其中最优的网络结构。本文构建了不同的卷积神经网络材料识别模型以分别应用于模拟数据及实验测量数据，当测量时间为1、3、5和10 min时，模拟模型的总体识别准确度分别为86.1%、95.6%、98.5%和99.9%；
实验模型的总体识别准确度分别为68.4%、81.3%、94.6%和97.1%。本文采用的卷积神经网络不失为一种有效的μ子成像材料识别方法，在实验装置1 m3的探测区域上对10 cm×10 cm×10 cm的钨块5 min即可实现99.1%的准确度，但由于实际应用需要大量学习训练，还要考虑到更复杂的探测情况，要克服这些限制，还需要更深入地进一步研究，优化模型结构。
