核脉冲反演型闪烁能谱测量技术

肖无云， 谢建明， 李京伦， 张羽中， 陈晔

(国民核生化灾害防护国家重点实验室， 北京 102205)

近半个多世纪，核脉冲信号处理技术主要是围绕解决脉冲幅度分辨率与脉冲通过率的矛盾问题展开研究。Radeka[1]根据前置放大器输出信号与噪声特性，从理论上推导出了一种最佳滤波器。根据该理论设计出了基于电阻-电容网络的核脉冲信号准高斯滤波方法，在模拟多道脉冲幅度分析器和核脉冲谱仪中得到了广泛应用。原则上讲，该理论的出发点是获得最佳能量分辨率，并没有考虑核脉冲通过率。由于准高斯脉冲较宽，尤其是下降沿有长的拖尾，不利于在高计数率下使用。Radeka[2]提出了梯形滤波成形思想。梯形脉冲作为准最佳时域响应，其前后沿陡峭，有利于提高脉冲通过率。但是采用模拟电路很难实现梯形滤波成形，这限制了该方法的推广使用。Jordanov[3]提出了一种简便的递归算法，通过数字信号处理技术设计了梯形滤波器。此后，梯形滤波器成为数字化多道与谱仪产品的核心模块。该技术推动了模拟谱仪数字化。研究表明，基于梯形滤波方法的数字化谱仪综合性能全面超越传统模拟谱仪[4]。需要指出的是，为了实现核辐射探测器输出电荷信号的完全积分，梯形滤波器的平顶宽度应大于探测器的电荷收集时间[5]。另外，由于探测器电荷收集时间存在涨落，要求梯形平顶宽度设置应有足够的裕度[6]。掺铊碘化钠[NaI(Tl)]是能谱测量中最常使用的无机闪烁体，常温下其发光衰减时间长达230 ns。为获得满意的能量分辨率，碘化钠谱仪一般使用不低于700 ns的脉冲积分时间，即使采用梯形滤波其脉宽也相当可观。在高计数率下，由此带来的脉冲堆积问题不可忽视。

近年来，微电子技术的快速发展为采用复杂数学计算来实时处理核脉冲信号提供了条件。针对以上问题，有学者提出了一些独特的核脉冲处理方法，如基于最大似然估计原理的核脉冲解析方法[7]，基于双模状态空间模型的卡尔曼平滑器[8]等，都具有很好的提升脉冲通过率的效果。但是这些方法计算量大，在现场可编程门阵列(field programmable gate array, FPGA)中实现较困难。为此，本文采用基于卷积模型[9]的核脉冲反演方法[10]，为高计数率下闪烁谱仪的核脉冲信号处理提供了简约的解决方案。本文主要探讨核脉冲反演方法在NaI(Tl)闪烁谱仪中的应用。

1 闪烁脉冲反演方法

1.1 闪烁脉冲能量分布

在谱仪核脉冲信号处理中，影响能量分辨率性能的主要因素有2个：1)单次粒子入射事件输出信号收集的完备性(无损失)与纯粹性(无堆积)；2)测量系统的信噪比。对能量分辨性能一般的无机闪烁探测器，光电收集效果是很重要的影响因素。为了平衡信号收集的完备性与纯粹性，需要深入探讨闪烁脉冲能量分布的时间特性及其处理方法。采用Φ75 mm×75 mm NaI(Tl)闪烁体配合光电倍增管以及电流型前放实测了137Cs γ射线点源产生的1万个核脉冲，图1(a)中给出了经归一化和对齐后做平均得到的标准核脉冲波形。该平均脉冲下降沿近似符合指数衰减规律，衰减时间常数约为230 ns。图1(b)展示了平均脉冲的时域积分结果，它反映了平均脉冲相应能量转换过程随时间的变化规律。从图上可以看到，欲收集95%的总电荷量，需要约1.3 μs的积分时间，这不利于高计数率应用。

图1 实测碘化钠探测器平均脉冲及其能量分布Fig.1 Measured average pulse and its energy distribution from a NaI(Tl) detector

1.2 基于卷积模型的闪烁脉冲反演方法

核脉冲反演方法主要分为直接反演法和迭代法2大类。迭代法通常需要多步迭代运算，很难在FPGA中实现在线实时处理。直接反演法也有很多种，一般需要建立尽可能准确的核脉冲模型。为了实现高性能的核脉冲反演，深入分析了闪烁体与光电转换器件组合输出核脉冲信号的物理过程，提出了一种闪烁脉冲卷积模型[9]。通过比较发现，卷积模型对实测脉冲的拟合效果显著优于既有的单指数、双指数和、双指数差模型，而且这3种模型也只是卷积模型在特定条件下的简化版。

卷积模型清楚地表达了闪烁脉冲的物理意义。其中的指数卷积项体现了闪烁脉冲形成过程中物理上的衰减特性，而高斯卷积项反映了核脉冲形成过程中的随机性。卷积模型为闪烁脉冲反演提供了较清晰的解析思路。考虑到高斯函数较难反演，而实测NaI(Tl)脉冲中高斯卷积核函数的标准偏差σ的拟合值仅为9.2 ns，在本项研究中舍弃对高斯卷积核的反演，而只反演易处理的指数卷积项。反演算法具体参见文献[10]。

图2(a)给出了对碘化钠探测器平均脉冲的反演结果，其主要成分是经反演残余的高斯型窄脉冲，其脉冲宽度约为反演前的1/24。但反演后的脉冲信噪比明显下降。而且在原脉冲持续时间范围内可以观察到较明显的震荡。这反映出卷积模型与实测平均脉冲也并不完美吻合，实际脉冲中电荷收集过程并不是平滑的，并且还有其他物理因素没有在卷积模型中得到体现，如：粒子与闪烁体相互作用过程中的多次散射事件，可见光在晶体中的传输与界面反射等。图2(b)给出了反演脉冲及其积分能量分布，可见核脉冲所含有用信息已集中在约50 ns的极窄时间间隔内。因此，若将核脉冲反演技术应用于闪烁谱仪的核脉冲成形处理，可显著提升相邻脉冲的时间分辨能力。

2 反演闪烁谱仪信号处理固件设计

2.1 数字信号处理流程

根据前述核脉冲反演方法，设计了如图3所示的数字核脉冲信号处理流程。经12 bit、200 MHz模数转换器(ADC)量化得到数字核脉冲信号序列，进入FPGA后，首先经过核脉冲反演处理，成形为尖锐窄脉冲。后续处理分为2路，一路经适度低通滤波形成快成形信号，再经脉冲甄别后进行堆积判弃，产生用于幅度提取的控制信号，并计算测量过程中的死时间；另一路经过窄梯形成形和基线恢复后，通过门限积分提取脉冲幅度，最后生成幅度分布谱。采用门限积分替代梯形平顶采样，是因为反演后核脉冲信噪比差，门限积分有利于改善幅度分辨率。死时间信息和幅度谱数据通过数据接口发送到上位机，供显示和分析处理。

图2 平均脉冲反演结果及其能量分布Fig.2 Deconvoluted results and its energy distribution of the average pulse

图3 反演闪烁谱仪的数字核脉冲信号处理流程Fig.3 Digital nuclear pulse processing flow of the deconvolution scintillation spectrometer

与传统数字化多道或谱仪相比，上述信号处理流程最主要的差异有以下3点：1)在数字核脉冲信号成形处理前先做了一步反演处理，2)梯形成形时采用的是参数可调的窄梯形，3)幅度提取时不是基于梯形平顶采样，而是使用门限积分。这种设计改进有利于在缩短核脉冲宽度的前提下尽可能控制噪声干扰。

2.2 FPGA固件设计

以本项目组前期设计开发的一款数字化多道脉冲幅度分析器[11]的FPGA固件源代码为基础，按照上述信号处理流程，修改其中核脉冲处理模块。采用赛灵思公司(XILINX)集成开发平台ISE，硬件描述语言为Verilog，完成核脉冲反演和门限积分，并将梯形成形器的上升沿和平顶参数调小。这样得到一台全新的反演型多道谱仪。

需要指出的是，核脉冲反演处理涉及多级浮点乘加运算。可以采用分布式查找表算法[12]实现。该方法的特点是速度快，但占用资源多，参数调整很不方便。也可以采用移位相加法[13]，该方法在浮点参数变化时需要重新安排运算结构，代码变化较大，但资源占用少，设计实现灵活方便。不过最简单、直接的方法还是使用系统自带的乘法器或乘累加器(MAC)来实现[14]，浮点参数调整和配置更为方便。门限积分器易于在FPGA中实现。只需在脉冲触发后，在设定的若干个时钟周期内，对核脉冲信号做累加即可。其结果作为核脉冲幅度，提供给后续的幅度-道址变换器。

全部信号处理模块统一采用200 MHz系统全局时钟驱动，所有操作统一采用上升沿触发模式，这样有利于避免出错。由于200 MHz的时钟频率较高，当信号处理较为复杂而FPGA芯片响应不够快时，设计可能难以正确实现。在这种情况下需要根据系统提示对出错部分的算法与逻辑进行优化设计，减少逻辑和布线的延迟，从而确保满足运行频率的要求。

在有限位宽情况下，为了避免运算过程中产生溢出错误，通常需要适当增加位宽。在本设计中，输入数字信号宽度为12 bit。全部中间处理统一采用32 bit冗余位宽。为了降低多级浮点运算过程中舍入误差的影响，在数字核脉冲信号输入级，适当左移若干位，右侧空位补零。对NaI(Tl)探测器，一般设定幅度谱的转换增益为1 024道。因此，在幅度-道址变换时，对获得的32 bit幅度值需要做位宽为10 bit的数据截取。具体截取哪一段需要根据多道谱仪的系统增益和能量测量范围要求，通过调试确定。

2.3 信号处理过程仿真

上述信号处理固件涉及到的算法和控制逻辑较为复杂，通常需要通过仿真排除设计错误，并优化控制参数的取值。常用的方法是直接使用ModelSim仿真软件或者通过ISE调用ModelSim做仿真。

图4给出了利用ModelSim对设计的Verilog代码所作的功能仿真结果。图中各信号从上到下依次为：激励信号、反演后的核脉冲信号、快甄别器输出信号、梯形成形器输出信号、基线恢复器输出信号、死时间信号、脉冲幅度值、堆积判弃控制信号。作为输入的激励信号，既可以使用用户专门制作的仿真信号，也可以使用实测核脉冲信号。这里用的是多道硬件实测信号。从图4可以看出，输入的数字核脉冲信号经反演后，去除了指数拖尾，变成了尖锐δ-形脉冲。但是，由于反演使用的是多级差分处理，使信噪比显著下降，δ-形脉冲几乎淹没在噪声中，不利于核脉冲甄别。因此需要在快成形器中施加适度的低通滤波来提高信噪比，降低噪声引起的误触发概率。主滤波器采用窄梯形成形器，在显著减小脉冲宽度的同时，也提高了信噪比，较好地兼顾了脉冲通过率和幅度分辨率。从图4中可以看到，窄梯形的宽度约为激励信号中大多数指数衰减脉冲宽度的1/5。尽管该脉冲已经很窄，但脉冲堆积概率仍然不可忽略，后续堆积判弃模块对特别靠近的堆积脉冲的判弃仍然必要。尤其是在高计数率下，这有利于降低峰堆积和尾堆积引起的谱形畸变。

图4 信号处理过程仿真结果Fig.4 Simulated result of pulse processing flow

将经仿真优化的FPGA代码编译下载到自主设计的数字化多道中，并连接NaI(Tl)探测器，得到的就是一套全新的反演闪烁谱仪。

3 能谱测量比对实验

3.1 能谱测量实验

为了验证自主设计的反演闪烁谱仪在实际能谱测量应用中的有效性，采用该谱仪实测了3.7×105Bq137Cs和3.7×105Bq133Ba放射源的混合γ能谱，并在探测器、放射源和几何条件完全相同的情况下，更换另外2种多道测量能谱以用于比较。如图5所示，其中标识为“Lynx”的是美国Canberra公司生产的通用数字化多道Lynx的测量结果，其梯形滤波器上升沿设置为1.6 μs，平顶宽度为0.6 μs。标有“RX1200”的测量结果来自美国Avicenna公司的RX1200型数字化多道，它采用标准的部分积分方法直接处理核脉冲，积分时间为1 μs。标有“本文设备”的是本文设计的反演谱仪的测量结果。其中使用的窄梯形成形的上升沿为125 ns，平顶宽度为75 ns。脉冲积分时间设定为500 ns。从实测谱可直观看到，在约2 000道以下低能区，反演谱仪的计数率相对最高；而在2 000道以上高能区，反演谱仪给出的计数率最低，表明其受脉冲堆积效应的影响最轻微。

图5 实测能谱比较Fig.5 Comparison of three measured spectra

3.2 综合性能参数比较

表1列出了对上述3个实测谱的分析结果。

表1 3种实测谱的分析结果Table 1 Analysis results of three measured spectra

可见基于经典梯形滤波法的Lynx多道的脉冲通过性能最差，但能量分辨率最好；反演谱仪测得的全谱总计数率和662 keV全能峰净计数率都是最高的，峰总比也最大，而能量分辨率则介于Lynx和RX1200的测量结果之间。

3.3 能量分辨性能随成形参数的变化

设计的反演闪烁谱仪具有良好的参数可调特性，通过改变梯形成形参数和幅度提取环节的脉冲积分时间，能够获得不同的能量分辨率。图6给出了在完全相同的探测器、放射源、几何条件和模拟电路参数配置下，采用上述反演闪烁谱仪，实测的137Cs γ能谱662 keV处能量分辨率随梯形成形参数和门限积分时间参数的变化情况。所用放射源活度为3.7×104Bq。

从图6中可以看出，在不同梯形成形参数下，相对能量分辨率都随门限积分时间的增加呈现出先急剧减小后缓慢增大的趋势。即在一定的梯形成形参数下，存在一个最优门限积分时间使得能量分辨率最佳。另外，当梯形脉宽相对较小时脉冲通过率高，只是所能达到的最佳能量分辨率差于探测器的标称能量分辨率(7%)。但只要适当增加梯形脉宽就能达到探测器的标称能量分辨率。因此，反演闪烁谱仪优越的参数可调特性使之能够覆盖高脉冲通过率和高能量分辨应用的需求。在一定能量分辨率要求下，反演闪烁谱仪能够提供更高的脉冲通过率；而在一定脉冲通过率下，反演闪烁谱仪也能够获得优异的能量分辨率。可以说，基于核脉冲反演的能谱测量技术对缓解脉冲通过率和能量分辨率的矛盾具有明显的实用价值，有利于解决大体积闪烁探测器或者高辐射水平下的能谱测量难题。

图6 实测662 keV处能量分辨率随不同梯形成形参数和门限积分时间的变化Fig.6 Measured energy resolution varying with different trapezoidal shaping parameters and gated integral time at 662 keV

4 结论

1)基于卷积模型的核脉冲反演方法与算法适于在FPGA中实时在线处理闪烁脉冲，具有显著缩短脉冲宽度的实际效果。

2)核脉冲反演型闪烁谱仪和普通闪烁谱仪相比，具有更高的核脉冲通过率，在强γ辐射场中，全谱以及全能峰计数率都有成倍提高。

3)反演闪烁谱仪的梯形成形与积分时间参数可独立调节，能很好地兼顾计数率和分辨率指标，为满足不同应用对能谱测量性能的要求提供了更通用和优化的解决方案。