基于主成分分析和Mahalanobis距离的异常γ能谱识别

删除或更新信息，请邮件至freekaoyan#163.com(#换成@)

本站小编 Free考研考试/2020-04-15

赵日 ^1,2 , 刘立业 ² , 李君利 ¹
1. 清华大学工程物理系, 北京 100084;
2. 中国辐射防护研究院, 太原 030006

收稿日期：2017-02-23
基金项目：国家自然科学基金项目（11605162）
作者简介：赵日(1988-), 男, 博士研究生
通信作者：刘立业, 研究员, E-mail:liuliye@cirp.org.cn

摘要：为了提高货物或人体放射性筛查的可靠性，提出了一种基于主成分分析和Mahalanobis距离的异常γ能谱识别方法。该方法首先对大量不含异常放射性的测量对象产生的正常γ能谱进行主成分分析，提取出其所有主成分，并按从大到小的顺序，选取前若干主成分构成子空间；将正常及待识别γ能谱在此子空间上投影，得到它们的Mahalanobis距离，通过比较这些距离的相对大小实现对异常γ能谱的识别。Monte Carlo模拟实验和实际测试实验表明，在子空间信息量占原始信息比例大于99%时该方法可靠有效。
关键词：γ能谱异常识别人工智能主成分分析Mahalanobis距离
Anomaly gamma spectra detection based on principal component analysis and the Mahalanobis distance
ZHAO Ri^1,2, LIU Liye², LI Junli¹
1.Department of Engineering Physics, Tsinghua University, Beijing 100084, China;
2.China Institute of Radiation Protection, Taiyuan 030006, China

Abstract: A principal component analysis (PCA) and Mahalanobis distance (MD) based anomaly gamma spectra detection method is developed to improve the reliability of radiation scanning of goods and human bodies. This method first extracts all the principal components (PCs) of large numbers of benign gamma spectra by PCA and selects several largest PCs to form a subspace. The algorithm then projects the benign, unknown gamma spectra to this subspace, calculates their MDs, and completes the anomaly detection by comparing these MDs. Monte Carlo simulations and actual tests show that the method is reliable and effective when the subspace has more than 99% of the original information.
Key words: gamma spectraanomaly detectionartificial intelligenceprincipal component analysisMahalanobis distance
异常γ能谱识别是指在特定的测量环境下识别出由含异常放射性核素的测量对象形成的γ能谱，从而完成对测量对象的异常放射性筛查。本文中，由含异常放射性核素的测量对象形成的γ能谱被称为异常γ能谱，此外则为正常γ能谱。在许多场景下，异常γ能谱识别具有重要的实际意义，例如港口、边境的货物放射性筛查，核工业从业人员的体表和体内放射性污染监测等。目前，异常γ能谱识别依赖于传统的基于全能峰的能谱分析技术^[1-3]，然而这一方法在实际分析中容易出现误识别与漏识别，特别是在能谱计数率较低、道址与能量关系不确定以及能谱中存在类似全能峰的假峰结构的情况下。
近年来，在将人工智能思想应用到γ能谱分析方面，国外已经开展了初步的理论和应用研究：Boardman等及Runkle等将机器学习应用到边境车辆的放射性检测^[3-4]，Dragovic等、Heggemanna等及Varley等则借助人工智能算法根据γ能谱实现了土壤放射性特征分析^[5-7]，Wang等综合多种子空间学习技术提高了中子探测器的中子分辨能力^[8]，而Hata等则借助支持向量机完成了乏燃料容器放射性强度分类^[9]。国内目前这方面的研究很少，杨杨在对地层岩性测井识别的研究中使用了基于自然γ能谱的模糊聚类进行辅助识别^[10]。
本文提出了一种基于主成分(principal component, PC)分析和Mahalanobis距离的异常γ能谱识别方法。该方法通过对大量正常γ能谱的主成分分析提取其共有特征，并利用Mahalanobis距离这一指标进行特征符合度的比较，最终实现未知γ能谱的异常识别。
1 方法原理1.1 主成分分析主成分分析方法已成为了统计分析、人工智能领域最常用的无监督线性降维方法^[11]。考虑一个由m个正常γ能谱构成的样本集Ψ={x_i|x_i∈Rⁿ; i=1, 2, …, m}，每个γ能谱由n道数据构成，即n维。各道间的协方差矩阵为

${\mathit{\boldsymbol{C}}_{jk}} = \sum\limits_{i = 1}^m {\frac{{({\mathit{\boldsymbol{x}}_{ij}} - {{\mathit{\boldsymbol{\tilde x}}}_j}){{({\mathit{\boldsymbol{x}}_{ik}} - {{\mathit{\boldsymbol{\tilde x}}}_k})}^{\rm{T}}}}}{{m - 1}}} .$

(1)

式中：i=1, 2, …, m; j=1, 2, …, n; k=1, 2, …, n。C_jk是协方差矩阵的第j行第k列元素，x_ij是第i个能谱的第j道计数，${{{\mathit{\boldsymbol{\tilde x}}}_j}}$是第j道的平均值。
对C进行特征值分解，

$\mathit{\boldsymbol{C}} = \mathit{\boldsymbol{\zeta \lambda }}{\mathit{\boldsymbol{\zeta }}^{\rm{T}}}.$

(2)

式中：ζ的各列是C的特征向量；λ是对角矩阵，各对角元是C的特征值。
C的特征向量即集合Ψ的主成分^[12]，依特征值从大至小排列，可得Ψ的第1, 2, …，n个主成分PC₁, PC₂, …，PC_n，这n个主成分最终形成对Ψ的完整描述。
1.2 Mahalanobis距离对于上述原始的或经主成分分析降维后的γ能谱集合，某γ能谱y(属于或不属于该集合)距该集合的Mahalanobis距离^[11]d为

$d = \sqrt {{{\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{\tilde x}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{C}}^{ - 1}}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{\tilde x}}} \right)} .$

(3)

其中：C与1.1节定义相同，即集合中能谱各道间的协方差矩阵；${\mathit{\boldsymbol{\tilde x}}}$是平均能谱。
Mahalanobis距离是一种广义距离，因为考虑了各变量间的协方差，与普通Euclid距离相比，它能消除量纲及各变量间相关性的影响，本质上可以认为它是某样本点属于某集合的概率度量。
2 材料与方法2.1 Monte Carlo模拟实验本文使用MCNP(Monte Carlo N particle transport code)模拟3英寸×5英寸×16英寸(7.62 cm×12.7 cm×40.64 cm)的NaI探测器测量谱，测量对象为高60 cm、底面半径约为15 cm的充水圆柱体。
正常γ能谱的模拟：源能量为1 460 keV，以模拟天然核素⁴⁰K。为了仿真现实场景，每一次模拟计算时，圆柱体的半径r、距探测器的距离l、源抽样数n以及源的能量E均在一定范围取随机值，即0≤r≤20, 10≤l≤15, 1×10⁶≤n≤2×10⁶，1 430 keV ≤E≤1 490 keV。能谱展宽按典型的NaI探测器实测展宽函数计算。记录能谱时，从0到1 800 keV等间距取512道。通过模拟，最终得到正常能谱625个。
异常γ能谱的模拟：在正常γ能谱场景中的源抽样中加入661 keV能量的光子以模拟异常核素¹³⁷Cs，通过选取661 keV与1 460 keV γ光子的抽样比例实现不同程度的异常能谱。共模拟得到3个异常γ能谱，对应661 keV光子与1 460 keV光子的抽样比例分别为0.1:1、0.5: 1及1:1。下文将此3个模拟谱分别称为MC-1、MC-2、MC-3，MC-1异常程度最轻微，MC-3异常程度最严重。图 1中显示了MC-1与MC-3的形态。

图 1 模拟的异常γ能谱形态

图选项

2.2 核电厂全身计数器实测γ能谱的描述正常γ能谱：本文收集了某在运行核电厂内照射监测实验室全身计数器(whole body counter，WBC)实测的2万余个正常γ能谱。全身计数器是直接测量由人体内放射性核素发出的γ射线的装置，用于判断被测人员是否因摄入放射性核素而受到内照射。本文收集的γ能谱由一台Canberra公司生产的型号为FASTSCAN的WBC测得。这些能谱的采集时间为2010—2014年间，每个能谱的测量有效时间(live time)均为60 s，均由1 024道计数构成，并均进行了相应的能量和效率刻度。
异常γ能谱：由于核电厂极少发生摄入放射性核素而产生内照射的情形，因此收集的由该电厂WBC记录到的异常能谱仅2个。这2个能谱的文件格式、测量时间、道数等情况均与上述正常能谱相同。为了便于描述，本文将这两个能谱记为WBC-1、WBC-2，其形态如图 2所示。

图 2 WBC记录到的异常γ能谱形态

图选项

2.3 能谱数据预处理2.3.1 道址的重新划分各个γ能谱的道与能量的关系并不完全相同，需要将各能谱重新按能量划分道址，并将原始γ能谱数据通过一定方法归并到新的道址中。
假设将从E_min至E_max能量范围划分为t个区间，即t道：[E₁, E₂), [E₂, E₃), …, [E_t, E_t+1]。其中：E₁=E_min，E_t+1=E_max。能谱x_i的道址-能量刻度曲线即函数关系f_i，则新道k的计数为：

$\begin{array}{l}\mathit{\boldsymbol{x}}_{ik}^{{\rm{rebin}}} = \sum\limits_{j = 1}^n {{\mathit{\boldsymbol{x}}_{ij}}} \frac{{p - q}}{{{f_i}\left( {j + 1} \right) - {f_i}\left( j \right)}},\\i = 1,2, \cdots ,m,k = 1,2, \cdots ,t.\end{array}$

(4)

式中：$\mathit{\boldsymbol{x}}_{ik}^{{\rm{rebin}}}$表示重新并道后的第i个能谱的第k道计数，p为f_i(j+1)、E_k、E_k+1三者的中间值，q为f_i(j)、E_k、E_k+1三者的中间值。式(4) 实质就是将原能谱中每道的计数按该道落入新道的能量比例进行重新分配。
为了降低统计涨落的影响，新道应当设置得宽一些，可按道宽为该处全能峰半高宽3倍的准则确定各道的能量区间。
2.3.2 归一化为了突出γ能谱的形态特征，削弱计数对形态的影响，重新划分道址后还需要对能谱数据进行归一化处理。对于重新划分道址后的γ能谱$\mathit{\boldsymbol{x}}_i^{{\rm{rebin}}}$，其归一化为：

$\begin{array}{l}\mathit{\boldsymbol{x}}_{ij}^{{\rm{norm}}} = \mathit{\boldsymbol{x}}_{ij}^{{\rm{rebin}}}/\mathop {{\rm{max}}}\limits_{1 \le j \le t} \mathit{\boldsymbol{x}}_{ij}^{{\rm{rebin}}},\\i = 1,2, \cdots ,m,k = 1,2, \cdots ,t.\end{array}$

(5)

式中$\mathit{\boldsymbol{x}}_{ij}^{{\rm{norm}}}$表示归一化后的第i个能谱的第j道计数。
2.3.3 中心化在主成分分析前，一般还要进行数据中心化，即将各个能谱分别减去平均能谱：

$\begin{array}{l}\mathit{\boldsymbol{x}}_{ij}^c = \mathit{\boldsymbol{x}}_{ij}^{{\rm{norm}}} - \mathit{\boldsymbol{\tilde x}}_j^{{\rm{norm}}} = \mathit{\boldsymbol{x}}_{ij}^{{\rm{norm}}} - \frac{1}{m}\sum\limits_{i = 1}^m {\mathit{\boldsymbol{\tilde x}}_{ij}^{{\rm{norm}}}} ,\\i = 1,{\rm{ }}2, \cdots ,m,k = 1,2, \cdots ,t.\end{array}$

(6)

式中$\mathit{\boldsymbol{x}}_{ij}^c$表示中心化后的第i个能谱的第j道计数。
2.4 算法流程本文方法包含两个主要步骤，即正常γ能谱的信息提取和未知γ能谱的异常识别。
在正常γ能谱样本的信息提取中，首先需要对正常γ能谱样本集Ψ进行数据预处理，得到预处理后的样本集Ψ′；然后对Ψ′进行主成分分析，得到其全部t个主成分(重新确定道址后能谱变为t道)，并按对应特征值从大到小的排序，选取前k个主成分PC₁, PC₂, …, PC_k构成k维子空间Ω；将Ψ′在Ω上投影，得到降维后的新集合Ψ″；计算Ψ″中各点到原点的Mahalanobis距离，并得到此距离的95%分位值d^*。在正常γ能谱Mahalanobis距离的具体概率分布函数未知的情况下，d^*可近似视为此分布的95%概率分位数。
在此基础上进行未知样本的异常识别时，同样先将未知能谱y进行数据预处理，然后将y′投影在子空间Ω上得到y″，计算y″与原点的Mahalanobis距离d，并与d^*比较，如果d大于d^*则表明该γ能谱异常，否则为正常。根据d^*的定义，此时的识别置信度为95%。
本文基于MATLAB实现了各模块的算法，算法流程如图 3所示。

图 3 算法流程框架示意图

图选项

3 结果与讨论3.1 模拟实验的结果γ能谱预处理时，在50~1 800 keV的能量范围内，按道宽为该处全能峰半高宽3倍的准则进行道址的重新划分，最终划分为13道。划分前后能谱的形态如图 4所示。

图 4 道址重新划分前后的γ能谱

图选项

为了揭示主成分与原始数据在信息含量上的关系，本文研究了随着子空间Ω维数k的增加，Ω上投影Ψ″的方差占投影前原始方差即Ψ′方差的百分比，即var(Ψ″)/var(Ψ′)，其中var表示方差。结果如图 5所示，在曲线的最初阶段，随着k的增加，var(Ψ″)/var(Ψ′)快速上升，很快达到并超过90%的比例，此后曲线缓慢上升，并不断缓慢接近100%，最终在k=13时完全达到100%。这一结果说明了排序靠前的PC包含了Ψ′的大部分信息，随着PC的下标“i”的增大，后面的PC包含的信息越来越少，但最终所有的PC完整地包含Ψ′全部的信息。另外，随着序数的增加，PC信息含量的变化趋势可从图 6中更清晰地显示出来。图 6中的横坐标为PC的序号i，纵坐标为数据集Ψ′在该PC_i方向投影的方差占Ψ′原始方差的比例，即var(Ψ′_i)/var(Ψ′)，其中Ψ′_i为Ψ′在PC_i方向上的投影。Ψ′在PC₁方向投影的方差为Ψ′原始方差的86.6%，PC₂方向的这一比例为30.5%，而对于PC₁₀这一比例为7.7%。

图 5 子空间信息含量随维数k的变化趋势

图选项

图 6 主成分PC_i信息含量随i的变化趋势

图选项

在异常γ能谱识别测试阶段，首先选取了k=2即选取PC₁和PC₂作为子空间Ω的基，图 7显示了此时的识别效果。可见，正常γ能谱集合Ψ′在此二维Ω空间(平面)上聚集在一起，MC-2、MC-3的投影点则远离此聚团，显示了一定的异常识别能力。然而此时，MC-1的投影点却依然处在聚团中，无法得到正确识别。MC-1、MC-2、MC-3 γ能谱对应的Mahalanobis距离d₁、d₂、d₃分别为0.63、3.86、4.92，而Ψ″中所有样本点的Mahalanobis距离的95%分位值d^*为2.33，即d₁＜d^*＜d₂＜d₃，这一计算结果与图 8显示的结果相符合。

图 7 模拟实验中k=2时的识别效果

图选项

图 8 d₁、d₂、d₃、d^*随子空间维数k的变化

图选项

图 8显示了随着k的增加，d^*与d₁、d₂、d₃关系的变化。当k＜6时，d₁＜d^*＜d₂＜d₃，MC-1无法正确识别；当6≤k≤11时，d^*＜d₁＜d₂＜d₃满足识别要求。当k＞11时，Ψ′在Ω上的投影Ψ″的协方差矩阵高度病态，无法对其求逆，即此时无法完成异常识别。
k=6时，var(Ψ″)/var(Ψ′)=99.1%，而k=5时，var(Ψ″)/var(Ψ′)=97.7%，可见在选择k时，如果以var(Ψ″)/var(Ψ′) ≥99%为准则则能保证最终的识别效果，但同时还要避免k取值过大而造成计算过程中奇异矩阵的出现。
3.2 核电厂WBC实测γ能谱的实验结果在γ能谱预处理阶段，选取能谱的50~3 000 keV能量范围重新划分道址，划分的准则与模拟实验相同。重新划分后，最终得到17道。划分前后的能谱形态如图 9所示。对2万余个预处理后的γ能谱进行主成分分析，得到共17个主成分。选取前k个主成分构成子空间，并计算投影前后的数据方差比例，即var(Ψ″)/var(Ψ′)。当k≤11时，var(Ψ″)/var(Ψ′)＜99%，k=12时，var(Ψ″)/var(Ψ′)=99.2%，因此按模拟实验时得出的k选择准则，异常识别时k取12即可。计算k=12时各能谱的Mahalanobis距离并得到其95%分位值d^*。

图 9 道址重新划分前后的WBC的γ能谱

图选项

对待识别能谱WBC-1、WBC-2用同样的方法进行预处理，并在前12个主成分构成的子空间上进行投影，计算它们的Mahalanobis距离d₁、d₂。
WBC-1：d^*=4.7, d₁=16.56，即d₁＞d^*，可知此能谱在95%置信度下为异常能谱；
WBC-2：d^*=4.7, d₂=17.58，即d₂＞d^*，可知此能谱在95%置信度下也为异常能谱。
3.3 讨论以上的结果表明，在合理设置相关参数的前提下，本文提出的基于主成分分析和Mahalanobis距离的异常γ能谱识别技术有效可靠；同时，具体实现流程也显示出了该方法几个明显的优点：1) 该方法以大量样本为识别依据，这就使得最终的识别能够建立在对样本特征更全面更充分掌握的基础上，避免了只着眼于单个能谱而造成误判；2) 能谱的计数信息被弱化，而与全谱相关的形态信息被提取和凸显出来，使得用于甄别的信息量比传统技术关注的峰结构所包含的信息量更具特异性，这也提高了低计数率情形下识别的准确性；3) 由于其依据大量样本共同特征来进行识别，该方法抵抗假峰干扰的能力较强。
4 总结为了弥补传统异常γ能谱识别技术在能谱计数率低、道址与能量关系不确定以及存在假峰时不够可靠的弊端，本文提出了基于主成分分析和Mahalanobis距离的异常γ能谱人工智能识别方法。通过原理研究和实验测试，建立了该技术的基础算法流程，包括γ能谱数据预处理、主成分提取、子空间投影、Mahalanobis距离计算等，同时确立了关键参数子空间维数k的选取依据，即k应足够大，使得子空间数据信息量达到原始数据信息量的99%以上，同时避免k过大而造成计算过程出现矩阵奇异的情况。该方法在Monte Carlo模拟实验和核电厂WBC实测γ能谱识别测试中表现出良好的准确性和可靠性。

参考文献

[1]	Canberra Industries. Genie 2000 Customization Tools Manual[R]. Meriden, CT:Canberra Industries, 2006:300-310.
[2]	Advanced Measurement Technology. GammaVision V6 Users Manual[R]. Oak Ridge, TN:Advanced Measurement Technology, 2006:180-185.
[3]	Boardman D, Reinhard M, Flynn A. Principal component analysis of gamma-ray spectra for radiation portal monitors[J]. IEEE Transactions on Nuclear Science, 2012, 59(1): 154–161. DOI:10.1109/TNS.2011.2179313
[4]	Runkle R, Tardiff F, Anderson K, et al. Analysis of spectroscopic radiation portal monitor data using principal components analysis[J]. IEEE Transactions on Nuclear Science, 2006, 53(3): 1418–1425. DOI:10.1109/TNS.2006.874883
[5]	Dragovic S, Onjia A. Classification of soil samples according to their geographic origin using gamma-ray spectrometry and principal component analysis[J]. Journal of Environmental Radioactivity, 2006, 89: 150–158. DOI:10.1016/j.jenvrad.2006.05.002
[6]	Heggemanna T, Welpa G, Amelung W, et al. Proximal gamma-ray spectrometry for site-independent in situ prediction of soil texture on ten heterogeneous fields in Germany using support vector machines[J]. Soil & Tillage Research, 2017, 168: 99–109.
[7]	Varley A, Tyler A, Smith L, et al. Mapping the spatial distribution and activity of ²²⁶Ra at legacy sites through machine learning interpretation of gamma-ray spectrometry data[J]. Science of the Total Environment, 2016, 545: 654–661.
[8]	Wang C L, Funk L, Riedel R, et al. Improved neutron-gamma discrimination for a ³He neutron detector using subspace learning method[J]. Nuclear Instruments and Methods in Physics Research A, 2017, 853: 27–35. DOI:10.1016/j.nima.2017.02.022
[9]	Hata H, Yokoyama K, Ishimori Y, et al. Application of support vector machine to rapid classification of uranium waste drums using low-resolution γ-ray spectra[J]. Applied Radiation and Isotopes, 2015, 104: 143–146. DOI:10.1016/j.apradiso.2015.06.030
[10]	杨杨. 基于模糊聚类的地层岩性测井识别研究与配套软件研制[D]. 成都: 西南石油大学, 2014. YANG Yang. Research and Development of Formation Lithology Logging Identification and Matching Software Based on Fuzzy Clustering[D]. Chengdu:Southwest Petroleum University, 2014. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10615-1014409667.htm
[11]	周志华. 机器学习[M]. 北京: 清华大学出版社, 2016.ZHOU Zhihua. Machine Learning[M]. Beijing: Tsinghua University Press, 2016. (in Chinese)
[12]	Han J, Kamber M, Pei J. Data Mining:Concepts Techniques[M]. Singapore:Elsevier (Singapore) Pte Ltd, 2012.