改进辛普森-SNIP算法在航空γ仪器谱本底扣除中的应用

2020-06-15葛良全张庆贤熊茂淋

核技术 2020年6期

杨津李飞葛良全张庆贤熊茂淋孙坤谷懿

（成都理工大学核技术与自动化工程学院成都 610059）

航空γ能谱测量（Airborne Gamma-ray Spectrometry Survey，AGS）是指将航空γ能谱仪系统安装在飞机、飞艇或其他飞行器上，测量地空界面的γ射线能谱，从而获得地表介质中放射性铀、钍、钾等核素的含量（比活度）。该方法也可以应用在贵金属矿、稀土矿、多金属矿矿床等与放射性核素相关的其他矿产的普查上［1］；近年来被广泛应用于岩性地质填图、地下水资源搜寻、放射性环境污染监测等诸多方面［2-4］。通过航空γ仪器谱的道址-计数，可以计算当前谱线反映的放射性核素的种类和含量。因此，航空γ仪器谱的解谱工作一直是该领域研究的重点和难点。其中对全谱的本底扣除是解谱过程中最重要的步骤之一［5］。本底扣除的可信度和精确度将直接影响到后续的寻峰、峰面积计算、核素含量反演等。

理想的扣除本底方法应满足以下特征：1）原理简单、运算快速且本底扣除的设定参数少；2）计算机能进行快速自动处理且运行成本低；3）能有效扣除因康普顿效应对谱峰计数造成的影响。针对散射本底扣除这一步骤，国内外提出了大量方法，其中包括统计敏感的非线性迭代剥峰算法（Statisticssensitive Nonlinear Iterative Peak-clipping，SNIP）及其相关衍生算法［6］、傅里叶变换法［7-8］、剥峰法［9］、迭代法［10］等。已知最简便的方法是线性法，即将能量特征峰的两波谷用一条直线连接，将其视为基线，基线以下为本底部分，但该方法会将有效计数当作本底扣除，相对误差大，只适用于粗略估计。SNIP方法被认为是其中的一种较为理想的选择，因为它具有原理简单、运行速度快、成本低等优点，国外核数据计算系统中也已将它运用于全谱本底评估和扣除。张庆贤等［11］结合SNIP本底扣除方法与傅里叶变换低通滤波方法提出了傅里叶变换本底逐步逼近法。然而，SNIP算法容易受到峰宽参数的影响，采用固定峰区宽度进行本底扣除会带来较大的误差；其次，算法迭代过程中的收敛速度缓慢，导致运行时间延长。因此，本文基于具有较高代数精度的辛普森公式的思想，针对原始算法中的不足，提出了新的改进方法，并将其应用于航空γ仪器谱本底扣除中。改进后的SNIP算法在保留原始算法简洁、高效的优点的同时，减少了迭代次数，提高了运行效率，使得谱峰的本底扣除效果和准确度进一步改善。

1 本底来源分析

通常而言，本底是指除被测量的放射源之外的其他能够在能谱区间造成计数的干扰因素，如放射性污染、宇宙射线及仪器的电子学噪声等。本底扣除是任何辐射测量中不可避免的问题，如果要提高放射性核素含量和放射性活度分析的准确性，则必须有效地识别和扣除与全能峰面积无关的计数［12］。

γ能谱中的本底主要来自康普顿散射效应和γ射线在探测器灵敏体积内产生的小角度散射以及其他干扰因素在能谱的能量区间造成的无用计数［13］。由于γ射线同物质发生光电效应、康普顿散射和电子对效应等一系列相互作用，当γ射线同周边散射和探测器自身发生散射，产生能量更低的散射光子和反冲电子［14］。当散射光子被探测器探测或者散射光子逃离探测器时，输出的γ谱中会出现康普顿平台［15］。此外一些外在的环境因素，例如环境中的放射性核素、电子仪器引起的噪声、宇宙射线、电磁干扰等会在能谱中产生环境本底［16-17］。而放射源的内在因素会导致放射性本底的出现，因为在放射源中除了待测核素外，还存在大量其他的无关核素，如果这些核素也有放射性，就会产生计数干扰［18］。在对γ能谱进行分析的过程中，散射本底加大了γ能谱的复杂程度，使得一些弱全能峰难以被识别，从而降低能谱分析数据的精度［19］。在γ能谱的分析过程中，核心问题是要提高能谱仪的测量精度，其关键在于选择合适的本底扣除方法，才能有效抑制康普顿散射对谱线的影响。

2 改进型辛普森-SNIP方法

2.1 SNIP算法及其改进

目前，SNIP算法被认为是最优的本底扣除算法之一。在1988年，Ryan等［20］首次提出SNIP算法扣除γ能谱本底的原理，但SNIP算法在细节描述方面，很大程度受到峰区宽度参数的影响，当选用不同的参数的时候，扣除本底的效果存在差异。随后Morháč等［21］在 1997年运用LLS运算符（taking the Square root and using the Ln operator twice）改进了在弱峰下的SNIP算法，在运用SNIP算法之前对每一道址的能谱计数进行变换。在2009年，Morháč等［22］将原始SNIP算法中的函数用高阶滤波函数替代，来缓解因将原始SNIP算法应用于不同能量区时，固定窗宽相同而产生的本底扣除误差，从而得到一个更精确的本底扣除效果。2012年，尹旺明等［23］通过借助半峰全宽（Full Width at Half Maxima，FWHM）刻度对SNIP滤波窗宽度进行自适应调整，进一步解决固定窗宽带来的影响。2019年，Li等［24］针对SNIP算法的演变给出了准确评价，提出新的算法演变思路，并指出基于SNIP的改进算法，都在避免因变换宽度而出现扣除本底过程中在峰的高能端与低能端出现的较大误差。而在SNIP光滑的本身，可以采取迭代的思想优化。

利用SNIP算法进行本底扣除，求解净峰面积主要分为4个过程。首先需要对能谱的每一道址计数变换，采用LLS算符（式（1）），通过自然对数运算压缩计数范围，计算高计数谱线，同时平方根运算可以提取弱峰信息。

式中：i表示道址；y（i）表示i道址的计数；V（i）表示经LLS运算后的变换值。

其次，利用迭代的方式，依次计算出V1（i）、V2（i）、…、Vm（i）。假设迭代到p次，则通过比较当前Vp-1（i）与Vp-1（i+p）和Vp-1（i-p）的平均值的大小，取两者中较小值作为下一次的迭代值（式（2））。

式中：p表示第p次迭代，即从1开始每次迭代运算后加1直到等于给定的m值为止（变换窗宽逐增法）；Vp（i）表示第p次迭代的LLS变换值。

迭代完毕后，再通过逆LLS变换（式（3）），求得本底谱数据。即作式（1）的逆运算：对变换值y（i）使用两次以e为底的指数运算，再进行求平方。

最后，用原始谱数据减去本底谱数据，得到γ能谱的净峰计数谱（式（4））。

本底扣除方法示意图如图1所示。

图1 SNIP本底扣除方法示意图Fig.1 The schematic diagram of SNIP background subtraction algorithm

在第二步SNIP算法迭代的过程中，实质利用的是矩形公式的思想。虽然计算简单，方法可靠，但是收敛速度太慢。本文采用辛普森公式的思想，即利用区间二等分的三个点来进行积分插值，其科特斯系数分别为1/6、4/6、1/6，由此来替代第二步的SNIP迭代计算。由于辛普森公式具有更高的代数精度，计算结果的准确度较之SNIP迭代更优秀。在迭代过程中，设定步长为单位距离，经过LLS转换的数值为对应的函数值（式（5））。

2.2 测试数据

本文采用XPG-3000A航空γ能谱勘查系统的实测数据。系统采用两箱探测器，每箱由两条40 cm×10 cm×5 cm阵列式NaI晶体探测器组成。单晶体能量分辨率在137Cs的0.661 MeV＜8.0%，能谱峰漂（208Tl的2.62 MeV）＜0.5道。能量范围最大1～10 MeV，最大线性计数率可达200 ks-1（kilo counts per second），每道最高计数264。

应用XPG-3000A系统对137Cs人工核素源进行测试响应。137Cs源活度为9 795 Bq，为了更好体现特征γ信息，采用探测器近距离测量的方式，即雷达高度为0 m。由于闪烁探测器的阵列式结构，γ射线与物质反应在探测器内部的电子沉积情况不尽相同，因此采用不同的探测点位分析，并记录谱线数据。每次实验测量时间为200 s，总测量次数为50次。

3 结果与讨论

Morháč等［21］在1997年提出单峰净峰计数保留率r和本底扣除率t两个系数，用于评价算法的本底扣除效果。

式中：PNC（Peak Net Counts）为扣除本底后的净峰计数；PNC0为未扣除本底的净峰计数；PTC（Peak Total Counts）为扣除本底后的总峰计数。

本文采用的是137Cs放射性核素的特征航空γ谱数据，其特征峰对应能量为0.661 MeV。在确定m值的时候，一般采用公式：

式中：W表示滤波窗的宽度，即变换宽度；m为SNIP迭代的次数。对于不同宽度W的单能峰来说，m值取在（W-1）/2附近，能获得最佳的谱峰本底扣除效果，且不会造成谱峰形状的改变。如图2所示，本文运用SNIP算法对谱线进行自适应扣除，先确定特征峰峰位道址为393道和谱峰区间为355～433道，再根据特征峰的峰区确定SNIP的迭代次数。

由图2，迭代次数会影响谱线的本底的剥离结果，当迭代次数m=7时，会导致原始谱出现“过度剥离”现象，SNIP法能削弱本底对特征峰的贡献，并减少本底干扰因素（如本底峰、屏蔽材料产生的特征X射线等）的影响，但同时过度削减了全能峰面积；当迭代次数m=12时，谱线的本底扣除会出现“欠剥离”现象，SNIP方法不仅不能有效抑制本底计数的影响，还会因本底扣除不彻底影响到净峰计数。由此，当迭代次数m=9时最为合理，能够有效抑制康普顿散射对谱线的影响。

图2 不同迭代次数下SNIP本底扣除效果Fig.2 Background subtraction effect shown on spectrum for SNIP algorithm with different iteration times

为了验证改进算法和原始的SNIP算法的本底扣除效果，本文引用了近年来常用的剥峰法［25-26］，同时作为本底扣除效果的参照。为减少剥峰法的剥离次数，提高运行效率，在剥离前对谱数据进行取对数，剥峰后，反变换可得到本底的形状。由于核素活度和元素含量的定量分析与相应特征峰的净峰面积有很直接的关系，故本文采用这三种方法对谱线进行本底扣除后得到特征峰谱线（图3），并通过对比扣除本底后得到的特征峰面积与实际净峰面积的平均相对误差来评价本底扣除效果（表1）。图3中曲线由上至下分别表示原始谱线用高斯函数拟合后的特征峰、经过原始SNIP计算得到的净计数峰、改进SNIP计算后的净计数峰以及经剥峰法处理得到的净计数峰。最后任意选取15次不同时间下测得的谱线，通过单峰净峰计数保留率r和本底扣除率t两个系数，用于进一步评价原始SNIP算法、改进SNIP算法以及剥峰法这三种算法的本底扣除效果（表2）。

表1 经本底扣除后得到的特征峰面积以及平均相对误差Table 1 The characteristic peak area and relative error table obtained after deducting from the background

通过比较三种方法得出的峰面积与实际特征峰面积的平均相对误差，由表1可知，采用原始SNIP法、改进SNIP法、剥峰法得出的峰面积的平均相对误差分别为7.60%、4.42%、-9.72%，改进方法的准确度较原始SNIP提高了41.8%。由表2可知，在本底扣除效果上，改进SNIP算法的r范围为97.79%～99.83%，原始SNIP算法的r范围为96.98%～99.45%，剥峰法的r范围为96.51%～99.18%，平均值分别为98.71%、98.29%、97.93%。改进SNIP算法的t范围为96.59%～99.47%，原始的SNIP算法的t范围为95.69%～98.78%，剥峰法的t范围为94.73%～98.76%，平均值分别为97.87%、97.15%、97.06%。通过对比，证明了在本文中改进的SNIP算法对单峰的净计数保留和本底扣除方面相对原始SNIP算法有细小提升，对特征峰的本底扣除效果最好，但两种方法的评价系数r和t均超过95%，具有足够高的精度，而剥峰法在该方面要稍次于这两种方法。

SNIP算法最大的优势是运算简单、迅速，因此需要对改进的SNIP算法的工作效率做评价。通过记录改进前后两种算法的迭代次数，其中原始的SNIP算法的最优迭代次数为12次，而基于辛普森公式迭代的改进型SNIP算法的最优迭代次数仅为9次，相对于原始的SNIP算法减少了三次，有一定的提升。进一步通过tic-toc的计时方式对两种算法的运行时间进行评价。