基于弹道亏损形状特征的堆积脉冲识别方法

王 明,周建斌,*,王怀平,2,汪雪元,2,刘 易,洪 旭

(1.成都理工大学 核技术与自动化学院,四川 成都 610059;2.东华理工大学,江西 南昌 330013)

X荧光光谱分析是一种速度快、非破坏性的测量方法,广泛应用于现场矿物筛查分类、安检等场景。X荧光光谱分析仪工作时,获得稳定的光谱需要一定的测量时间,且元素含量越低,所需的测量时间越长,工程应用中一般通过提高X射线管的照射强度减少探测时间。核辐射探测系统中,脉冲信号衰减至基线需要一定的时间,在脉冲持续时间内,若其他辐射事件继续与探测器响应,则会出现两个或多个脉冲堆积的现象,形成畸变脉冲,也称为堆积脉冲。堆积脉冲会导致能谱出现假峰(和峰)、本底提高、能量分辨率衰退等问题,常用数字脉冲成形技术及反堆积技术改善能谱性能[1]。高斯成形系统[2-3]及CR-RCn准高斯成形系统[4-7]是无限冲激响应系统,一般用于提升低计数率条件下能谱的能量分辨率。数字梯形成形算法能免疫弹道亏损、降低脉冲堆积的概率及提高信噪比[8-9],广泛用于获取高计数率与高能量分辨率能谱。此外,极零相消技术也是一种降低脉冲堆积概率的方法[10-11]。

快速成形技术是识别堆积脉冲的一种方法[12-15],该方法利用反卷积算法,将缓慢衰减的脉冲转变成冲激脉冲。冲激脉冲时间短、脉冲通过率高,利用冲激脉冲的时间间隔可判断梯形脉冲是否堆积,但该方法无法甄别时间间隔低于冲激脉冲分辨时间的堆积脉冲。此外,Imperiale等[16]利用梯形成形的上升时间来甄别堆积脉冲,Zhou等[17]采用梯形脉冲的宽度来甄别堆积脉冲,然而,梯形脉冲的上升时间及宽度随着脉冲幅度增加而增加。Yu等[18]利用双指数模型得到了平顶宽度稳定的梯形成形脉冲,并利用梯形平顶宽度来判别堆积脉冲,但测量时间特征受到ADC采样精度的限制。Jordanov[19-20]提出了一种新的脉冲形状甄别技术,该技术将原始脉冲成形为逆锯齿波脉冲和梯形脉冲,计算出其幅度比值,该比值称为时不变脉冲形状标签(TIPS),通过TIPS的大小判断脉冲是否堆积。TIPS不受脉冲幅度变化及ADC分辨时间的影响,仅与弹道亏损相关,但逆锯齿波成形器是高通滤波器,信噪比低,TIPS展宽大,不利于堆积脉冲识别。三角波成型器是带通滤波器,信噪比高,TIPS展宽小,为此,本文提出梯形-三角脉冲形状甄别技术,用于提升X荧光光谱测量系统检测堆积脉冲的能力。

1 TIPS特征描述

利用数字成形技术将原始数字脉冲分别成形为梯形脉冲及三角脉冲。辐射事件响应的梯形脉冲及三角脉冲示意图如图1所示,其中,Ⅰ为非堆积,Ⅱ为两个事件堆积,a为辐射事件序列,b、c分别为事件响应的梯形脉冲、三角脉冲。对于非堆积脉冲,梯形脉冲高度及三角脉冲高度均正比于辐射事件能量,TIPS是一个常数,与辐射事件能量无关。当脉冲堆积时,梯形脉冲能免疫弹道亏损,脉冲幅度依然正比于两个辐射事件的能量。三角脉冲不能免疫弹道亏损,脉冲幅度损失,堆积脉冲的形状特征TIPS变大。堆积时间间隔越大,TIPS越大。根据TIPS的大小,即可判断脉冲是否发生堆积。

图1 梯形脉冲及三角脉冲成形Fig.1 Trapezoidal pulse and triangular pulse

(1)

其中,PKflat、PKpin分别为梯形脉冲高度及三角脉冲高度。

2 原理

2.1 TIPS特征分析

经前端电路处理后的核脉冲信号,首先快速上升至峰值,然后缓慢衰减至基线,常用双指数衰减脉冲信号x(t)进行模拟:

(2)

其中:A为比例常数,与探头和电子学电路相关;E为辐射事件的能量;u(t)为单位阶跃信号,t为时间;τ及θ为衰减常量,τ与C-R微分电路相关,θ与探头和电路特性相关,θ通常远大于τ。τ描述的缓慢衰减指数信号用于表征核脉冲的慢成份;θ描述的快速衰减指数信号用于表征核脉冲的快成份。

梯形(三角)脉冲信号为:

(t-tb)u(t-tb)+(t-tc)u(t-tc)]

(3)

其中:yT为梯形脉冲信号;bT为比例常数;ta为上升时间,tb-ta为平顶脉冲信号时间,tc=ta+tb。当t在[ta,tb]之间时,最大值为bTE。当ta=tb时,式(3)可描述为等腰三角脉冲信号,ya为三角脉冲信号,当t=ta时,最大值为baE,ba为比例常数。

对于两个辐射事件堆积的情况,假设首先与探测器响应的辐射事件能量为El,间隔时间为ΔT时,能量为Er的辐射事件与探测器响应,且ΔT

ya(t) =ElV(t) +ErV(t-ΔT)

(4)

其中:tw为三角脉冲的宽度;V为三角脉冲信号;ya为重叠的三角脉冲信号。

设等腰三角波的宽度为tw,即tw=tc=2ta=2tb,堆积三角波的不可导位置t=tw/2或t=ΔT+tw/2的极值分别为ba(El+Er-2ΔTEr/tw)、ba(El+Er-2ΔTEl/tw)。当ba(El+Er-2ΔTEr/tw)≥ba(El+Er-2ΔTEl/tw),即El≥Er时,极大值点为ba(El+Er-2ΔTEr/tw);当ba(El+Er-2ΔTEr/tw)ΔT,堆积的梯形脉冲幅度为bT(El+Er)。

根据式(1),非堆积脉冲的TIPS为:

(5)

其中,TIPSa为非堆积脉冲的梯形-三角波形状特征,简称为三角波形状特征。

两个事件堆积时的TIPS为:

(6)

图2 TIPS-El/Er曲线Fig.2 TIPS-El/Er curve

2.2 甄别器评估方法

脉冲形状特征甄别器是一种二分类算法,即对非堆积脉冲和堆积脉冲进行分类。常用精确率、召回率以及F1得分来度量分类器的性能。精确率与召回率是矛盾的,精确率越高,召回率就越低。F1得分是精确率与召回率的调和平均数,它同时兼顾了分类模型的精确率与召回率,计算公式分别为:

(7)

(8)

(9)

其中:Pacc为精确率;Pcal为召回率;TP为正确的正例;FN为错误的反例;FP为错误的正例。本文将堆积脉冲设为正例,非堆积脉冲设为反例。精确率表示甄别系统正确识别堆积脉冲的比例;召回率表征甄别系统识别到的堆积脉冲与全部堆积脉冲的比值。

3 数字成形算法实现

对式(2)进行ADC采样量化,获得数字脉冲信号,采样时间为Ts。对数字脉冲进行z变换得z变换表达式为:

(10)

梯形(三角脉冲)传递函数HT(z)[17]为:

(11)

其中:na=ta/Ts;nb=tb/Ts;nc=tc/Ts。当na=nb时可得到三角波成形传递函数Ha(z)。

在硬件上实现数字成形算法以整数运算为基础,不存在溢出的情况下进行加法、减法以及乘法运算不会产生计算误差,除法运算容易产生截断误差。采用级联结构来设计数字成形算法可降低系统的复杂度。对于线性时不变系统,交换级联的顺序不影响最终的结果输出。差分项放在前级,累加项放在后级,能够防止溢出。

为计算TIPS,需根据式(11)将原始脉冲成形为三角脉冲和梯形脉冲。通过反卷积算法将原始脉冲x成形为冲激脉冲v,如图3a所示。为了避免浮点运算,将d1和d2改为有理分式形式。设d1≈nu1/de1、d2≈nu2/de2。其中nu1、de1、nu2、de2均为整数。级联项1-d1z-1及1-d2z-1分别改为(de1-nu1z-1)/de1以及(de2-nu2z-1)/de2。设三角脉冲的上升时间为na,缩放系数设为k1=D1/2M1,D1、M1均为整数。三角成形系统的硬件算法如图3b所示,ya为三角脉冲信号。设梯形脉冲的上升时间为na1,平顶时间为nb1-na1,缩放系数设为k2=D2/2M2,D2、M2均为整数。梯形成形系统的硬件算法如图3c所示,yT为梯形脉冲信号。

图3 算法实现逻辑图Fig.3 Logic diagram of algorithm implementation

算法部署于xc7z020CLG400-1芯片中,参数取值分别为:de1=100,nu1=95,de2=100,nu2=61,na=4,na1=8,nb1=12,D1=1,D2=800,M1=16,M2=17。算法资源消耗为:19个LUT,7个LUTRAM及359个FF。

4 实验

利用基于FAST-SDD探测器的CIT-3000MD X射线荧光光谱测量系统进行实验。FAST-SDD探测器(XR-100SDD)由Amptek生产。XR-100SDD内部封装了热电冷却固态探测器和前置放大器。XR-100SDD引出的独立电源线和信号线接入新先达(M&C Co. Ltd)的数字脉冲处理系统。光源使用科颐维公司设计的Ag靶X光管(KYW2000A),额定工作电压与额定工作电流分别为50 kV、1 mA。

梯形成形脉冲的上升时间为400 ns,平顶宽度为200 ns,脉冲持续时间为1 μs。三角脉冲的宽度为350 ns,冲激脉冲的时间分辨率约为150 ns(20 MHz 采样率条件下)[14]。采用铅黄铜样品进行实验,测量时间为100 s,光管的工作电压为33 kV,工作电流分别为3.9 μA、1 mA,冲激脉冲通过率分别为6.5 kcps/s、1.22 Mcps/s。

4.1 TIPS分布图

由于噪声以及弹道亏损等的影响,非堆积脉冲的TIPS并非固定值,而是符合一定分布,通过TIPS分布图可确定分割阈值。通过调整脉冲幅度的比例,使得TIPS谱线的均值为1 000,图4a为管流为3.9 μA时的TIPS分布,Pi为逆锯齿波形状特征TIPS分布[20]、Pa为三角波形状特征TIPS分布。低计数率时,TIPS谱线近似为非堆积脉冲的TIPS谱线,Pa展宽小,Pi展宽大。图4b为管流为1 mA时的TIPS分布。高计数率时,由于脉冲堆积的影响,TIPS分布高值区域出现长的拖尾。在高值区长尾范围内,随TIPS的增加,Pi的概率密度快速降低,Pa的概率密度缓慢衰减。

图4 TIPS分布Fig.4 TIPS distribution

表1 不同α对应的甄别阈值xαTable 1 Discrimination threshold xα corresponding to different α

4.2 测量结果及讨论

SPT谱中,铜元素Kα、Kβ特征X射线能量位于(7.6 keV,9.5 keV)能区。当脉冲通过率为6 500 s-1时,几乎没有堆积脉冲,高能区平坦。Kα峰与Kβ峰面积占总面积的95.71%。两个脉冲堆积时,一阶和峰位于(15 keV,18 keV)能区;3个脉冲堆积时,二阶和峰位于(20 keV,26 keV)能区;3个以上脉冲堆积的概率较低,且已超出ADC的测量范围,可忽略。

管流为1 mA时,梯形成形能谱如图5所示。SPT是利用传统的快成形技术[14]甄别堆积梯形脉冲获得的能谱。在快成形技术甄别堆积脉冲的基础上,利用梯形-逆锯齿波形状特征[20]或梯形-三角波形状特征进一步甄别堆积脉冲,分别获得能谱ISA-SPT[20]及TRI-SPT,其中非堆积脉冲损失率为15%。SPT中一阶和峰、二阶和峰计数率分别为26 580、2 030 s-1,ISA-SPT中一阶和峰、二阶和峰计数率分别为8 400、181 s-1,TRI-SPT中一阶和峰、二阶和峰计数率分别为6 750、116 s-1。因此,使用梯形-三角波形状特征甄别堆积脉冲,可进一步降低和峰计数率。

图5 拒绝堆积脉冲的能谱对比Fig.5 Comparison of energy spectrum rejected pile-up pulse

取Cu元素的Kα、Kβ峰射线全能峰以及其和峰计算出的精确率、召回率、F1得分列于表2。由表2可知:当α为5%时,F1得分最高;随α的增大,召回率逐渐增大,精确度逐渐降低;相同α条件下,三角波TIPS甄别法识别率高于逆锯齿波TIPS甄别法[20]。当α为15%时,三角波TIPS甄别法识别的精确率、召回率、F1得分分别为73.55%、78.75%、76.06%,分别提升了6.39%、7.66%、7.08%。

表2 不同α对应的堆积脉冲分类结果Table 2 Classification result of pile-up pulse corresponding to different α

峰总比表示全能峰的探测效率,可评价能谱的质量,峰总比表示全能峰面积与总面积的比。当电流为3.9 μA时,Cu元素Kα峰的峰总比为82.13%;当光管电流为1 mA时,Kα峰的峰总比降到62.93%。利用形状特征TIPS识别并拒绝堆积脉冲后,峰总比的测量结果如图6所示。PTi、PTa分别表示ISA-SPT、TRI-SPT能谱中Cu元素Kα峰的峰总比,峰总比随α的增加而增加。当α相同时,PTa高于PTi。当α=15%时,PTi、PTa分别为74.86%、76.60%,峰总比提升了1.74%。

图6 峰总比的测量结果Fig.6 Measurement result of peak to total ratio

5 结论

本文提出了一种基于弹道亏损形状特征的堆积脉冲甄别技术,其原理为:首先将原始脉冲成形为三角脉冲和梯形脉冲;其次计算其幅度比值,即脉冲形状特征值TIPS;最后根据TIPS判断脉冲是否发生堆积。由于TIPS与脉冲幅度及ADC采样率无关,在20 MHz采样率条件下,仍具有比较好的堆积脉冲甄别能力。通过理论计算分析了两个脉冲堆积时的TIPS特征值,堆积时间间隔越大或堆积事件的能量越接近,TIPS越大,堆积脉冲越容易甄别。

使用提出的方法识别铅黄铜样品中1.22×106s-1脉冲通过率的堆积脉冲。当非堆积脉冲的损失率为15%时,识别堆积脉冲的精确率、召回率以及F1得分分别为73.55%、78.75%以及76.06%,与文献[20]提出的TIPS甄别技术进行对比,分别提升了6.39%、7.66%以及7.08%。Cu元素的Kα峰的峰总比为76.60%,提升了1.74%。该方法减少了伪峰的计数率,在现场矿物筛查分类以及安检等快速识别微量元素场景中具有潜在优势。