本发明属于谱分析领域，具体涉及一种基于GPU和奇异值方法的γ谱仪谱分析方法。
背景技术：
γ谱仪能够对辐射与物质的相互作用产生的效应进行有效地测量和记录，对测得放射性核素的相关数据和曲线进行分析。目前，市场上还没有推出能对所有能量范围内的γ能谱进行有效测量的探测器，即现存探测器都集中在某一谱段进行探测。
核素识别是通过测量得到的γ能谱来分析现场核素的个数和种类的。传统方法一般通过对经过一段时间累积的谱线使用平滑、寻峰定位等算法，然后再利用能量刻度系数求出相应的γ能量，最后在核素库表中进行匹配来确定现场的核素种类。
当所使用的探测器分辨率较低时，如碘化钠NaI晶体探测器，所得到的谱线往往很复杂，测得的噪声会极大影响核素的弱峰信息。利用上述传统核素识别方法时，往往会产生漏报的问题。此外，市场上现存的基于NaI晶体探测器的核素识别算法计算量相对较大，故多在PC上位机上实现。而为兼顾经济能力、探测效率和操作方便等因素，市场上已经推出掌上电脑，即PDA(Personal Digital Assistant)形式的γ谱仪，以适用于对辐射源进行初步的探测和分析。本发明即在PDA形式的γ谱仪上改进现有核素识别算法的效率和精度。
另一方面，随着嵌入式处理器的发展，许多公司的新一代处理器已经集成了可编程图像处理器(Programmable Graphic Process Unit,GPU)，如Freescale公司提供的i.mx6q系列的处理器。该图像处理器具有单指令流多数据流(SIMD)的并行处理特性，并支持通用的异构开发环境，如OpenCL开发和OpenGL图形处理。利用这些新的计算能力，可以在γ谱仪上实现更复杂的计算任务。
目前，用于民用0-3MeV能量范围内的便携式γ能谱核素探测的探测设备主要有NaI探测器、高纯锗(HPGe)探测器。其中，HPGe探测器的γ能谱分辨率高，线性好，是大型核素分析设备常用的探测设备，但是其价格十分昂贵，并不适用于民用便携装置。相比较而言，NaI探测器价格低廉，十分适用于民用探测设备，但其低频段通道输出信噪比低，谱线相对较复杂。
在传统的能谱分析处理算法中，要求增长时间来满足一定的统计误差需求，来降低信噪比，之后从能谱的各个峰信息中找到相关核素的特征。此种方式会带来测量时间长的缺点。而且对于NaI探测器，长时间测量对于核素识别可信度改善也并不明显。
技术实现要素：
本发明的目的在于提供一种基于GPU和奇异值方法的γ谱仪谱分析方法，利用谱信息实现对现场核素种类数目进行估计。
本发明的技术方案如下：
一种基于GPU和奇异值方法的γ谱仪谱分析方法，包括以下步骤：
步骤一：数据采集
采用便携式γ谱仪分别对本底、样本进行采集；
设γ谱仪输出数据表示为
s(n)＝[s1(n)，s2(n)，…，sk(n)]，
其中n＝1，2，…，表示时间采样。将其转化为时间阵列序列
x(n)＝[x1(n)，x2(n)，…，xk(n)]，
满足
x(1)＝s(1)，
x(n)＝s(n)-s(n-1)，n＝2，3，…；
步骤二：数据处理
采用便携式γ谱仪嵌入式处理器的GPU单元进行并行数据处理，得到奇异值σ＝[σ1，σ2，…，σn]；
步骤三：数据分析
对得到的奇异值σ＝[σ1，σ2，…，σn]从大到小排序，可得如下序列
σ′1≥σ′2≥…σ′m≥C≥σ′m+1≥…≥σ′n≥0，m≤n，
其中，C为根据仪器设备选定的某个阈值；
之后，求得排序后的相邻奇异值之比γk＝σk-1/σk，k≥2；
求取t，使其满足
γt＝max{γ2，γ3，…，γn}≥D，
其中，D为根据仪器设备选定的某个阈值，将t作为核素数目的估计值。
步骤二中，首先利用阵列信号，构造协方差矩阵
然后利用GPU和数值分析理论中常见的QR分解方法，求取RHR的特征值λi，i＝1，2，...，n，并得到R矩阵的奇异值
步骤一中，设置采样间隔为1s，5分钟后得出γ谱样图。
本发明的显著效果在于：
(1)本发明提出一种利用特征空间分解方法来提升嵌入式数字化的γ谱仪的核素识别能力的算法，对来自数据采集探头的时间阵列序列进行综合处理与模式判断，尽可能减少对现场核素的漏报。
(2)本发明算法可以完成对现场核素种类数目进行估计，为后续改善核素种类识别能力提供有力支持。
(3)本发明利用了嵌入式处理器上的GPU单元进行数据处理，使得γ谱仪能更加快速地处理现场数据，大大减少了数据处理时间，提升了数据处理能力，同时有效降低了辐射现场对测试人员的身体损害。
附图说明
图1为通过γ谱仪采集到的数据的本底噪声波形图；
图2为采集Am241+Cs137混合物的γ谱特性曲线示意图；
图3为采集Am241+Cs137+Co60混合物的γ谱特性曲线示意图。
具体实施方式
下面结合附图及具体实施例对本发明作进一步详细说明。
一种基于GPU和奇异值方法的γ谱仪谱分析方法，包括以下步骤：
步骤一：数据采集
选择样本1为Am241+Cs137混合物，样本2为Am241+Cs137+Co60混合物。
采用γ谱仪分别对本底、样本1和样本2进行采集。其中，设置采样间隔为1s，5分钟后得出如图1-图3所示的γ谱样图，采样的阵列信号存储在γ谱仪的SD卡中。
利用γ谱仪测得的数据是时间的累积值，需要对其进行处理，得到时间阵列序列。设γ谱仪输出数据表示为
s(n)＝[s1(n)，s2(n)，…，sk(n)]，
其中n＝1，2，…，表示时间采样。将其转化为时间阵列序列
x(n)＝[x1(n)，x2(n)，…，xk(n)]，
满足
x(1)＝s(1)，
x(n)＝s(n)-s(n-1)，n＝2，3，…。
步骤二：数据处理
为了加速计算和减小谱仪的功耗，本步骤采用γ谱仪嵌入式处理器的GPU单元进行并行数据处理。
首先利用阵列信号，构造协方差矩阵：
然后利用GPU和数值分析理论中常见的QR分解方法，求取RHR的特征值λi，i＝1，2，...，n，并可以得到R矩阵的奇异值：
步骤三：数据分析
这里，假设谱仪接收到的信号阵列x(n)为信号源阵列o(n)和白噪声阵列输入N(n)之和，即其模型可以表示为x(n)＝Ao(n)+N(n)。而信号源o(n)和噪声输入N(n)不相关或者相关度较低，即E[o(n)NH(n)]小于某一个阈值。这样，得到R≈Ro+RN，其中Ro＝AE[x(n)xH(n)]AH为信号源o(n)的协方差的相似变换，RN与噪声协方差矩阵相关。当||RN||＜＜||Ro||时，有其中，为信号源的协方差矩阵特征值，为噪声协方差矩阵的特征值。故如果λi或σi较小，可估计为只有噪声；反之，则估计有信号源，即被测样本有对应的核素。
综上所述，本算法对得到的奇异值σ＝[σ1，σ2，…，σn]从大到小排序，可得如下序列
σ′1≥σ′2≥…σ′m≥C≥σ′m+1≥…≥σ′n≥0，m≤n，
其中，C为根据仪器设备选定的某个阈值。得到的结果如表1所示。
表1 本专利采样三种样品排序后的前八个奇异值
之后，求得排序后的相邻奇异值之比γk＝σk-1/σk，k≥2。求取t，使其满足
γt＝max{γ2，γ3，…，γn}≥D，
其中，D为根据仪器设备选定的某个阈值，将t作为核素数目的估计值。
对于表1中的样本，可以得到表2的计算结果。从表2可以看出，如果选择阈值D为1.1，则本底中无核素，而样本1和样本2分别在γ2和γ3取最大值，并且满足阈值条件。根据本算法，样本1有两种核素，样本2有三种核素，与实际相符合。
表2 各个样品计算的前八个奇异值之比
本发明提出一种利用信号关联矩阵奇异值比方法来提升嵌入式数字化的γ谱仪的核素识别能力的算法。对来自数据采集探头的时间阵列序列进行综合处理与模式判断，以尽可能减少对现场核素的漏报。该算法可以对现场核素种类数目进行估计，为后续改善核素种类识别提供有力支持。
同时，为了改善时间阵列序列中遇到的大数据量问题，本发明采用了嵌入式处理器上的GPU单元，运用并行运算手段提升了数据处理能力并降低手持设备的功耗。