梯形成形算法是一种广泛应用于核信号与能谱测量的数字滤波成形方法[1]，其本身有一定的噪声抑制能力，在梯形平顶时间大于探测器电荷收集时间的情况下，能够有效避免弹道亏损，且具有脉冲形状对称、脉宽窄、上升时间与平顶时间独立可调等诸多优点[2]。尤其在高计数率情况下，通过减小成形脉冲宽度，梯形成形算法能够达到一定的去除堆积、减小计数率损失的目的。
堆积效应是核信号测量中的普遍问题，也是影响幅度测量和能谱分辨率的关键问题。长期以来，中外学者对此问题进行了大量研究，目前提出的方法主要有：减小脉冲宽度以减小堆积概率[3-4]、依据参考波形外推[5-7]、极大似然估计[8-9]、能谱反卷积[10-11]等。在高计数率下，为了在现场可编程门阵列(field-programmable gate array, FPGA)中实时实现梯形成形和堆积判别等功能，减小成形脉冲宽度是一个可行且高效的方法，但此脉宽不能无限减小，需要综合考虑输入信号的上升时间、弹道亏损效应、最终对能谱的能量分辨率要求等。
X射线衍射仪是一种常用的晶体结构分析仪器。本文实验中的射线源为铜的特征X射线，主要存在Kα和Kβ两种特征，对应能量分别为8.05和8.9 keV，由于后续晶体分析实验中只用到Kα峰，故需通过单道上下阈值的方式将Kβ峰去除，防止它对分析结果产生干扰。本文实验系统结构如图 1所示。
|
| 图 1 X射线衍射实验系统结构 |
本文采用实时梯形成形算法，结合利用反卷积定位信号到达时刻的堆积脉冲识别算法，在FPGA上实现了实时的数字脉冲多道仪的大部分功能，并实际应用于X射线衍射仪上，获得良好的结果。
1 梯形成形算法在FPGA中的实现理想的梯形成形脉冲如图 2所示。
|
| 图 2 理想梯形脉冲 |
指数输入信号可表示为
| $ {V_i}\left( t \right) = A{{\rm{e}}^{-t/\tau }}u\left( t \right). $ | (1) |
其中：A为脉冲幅度，τ为指数信号的时间常数。进而可得梯形成形的系统函数[12]：
| $ H\left( z \right) = \frac{{{z^{-1}}\left( {1-d{z^{-1}}} \right)\left( {1 - {z^{ - {n_a}}}} \right)\left( {1 - {z^{ - {n_b}}}} \right)}}{{{n_a}{{\left( {1 - {z^{ - 1}}} \right)}^2}}}. $ | (2) |
| $ d = {{\rm{e}}^{-{T_{\rm{s}}}/\tau }}, $ | (3) |
| $ {n_a} = {t_a}/{T_{\rm{s}}}, {\rm{ }}{n_b} = {t_b}/{T_{\rm{s}}}, {\rm{ }}{n_c} = {t_c}/{T_{\rm{s}}}. $ | (4) |
其中Ts为系统的采样周期。
观察系统函数H(z)的表达式不难发现，此为一个无限冲激响应(infinite impulse response, IIR)数字滤波系统，具有反馈单元。为了降低系统设计的复杂度，避免量化误差与舍入误差的积累，故选择级联结构来实现。对于线性移不变系统，交换子系统的级联顺序并不会影响最终的输出，但在FPGA中具体实现时，由于存在误差与有限字长效应，系统的分解方式与级联顺序成为必须考虑的问题。
文[13]提出了一种系统的分解方式与级联顺序，并作了较为详细的分析和阐述。本文在此基础上进行了改进，利用更少的资源实现更稳定的梯形成形算法，
H(z)=H1(z)H2(z)H3(z)H4(z).
其中：
| $ \left\{ {\begin{array}{*{20}{c}} {{H_1}\left( z \right) = 1-{z^{-{n_b}}}, \;\;\;\;\;}\\ {{H_2}\left( z \right) = \frac{{1-d{z^{ - 1}}}}{{{n_a}\left( {1 - {z^{ - 1}}} \right)}}, }\\ {{H_3}\left( z \right) = \frac{{{z^{ - 1}}}}{{1 - {z^{ - 1}}}}, \;\;\;\;\;\;}\\ {{H_4}\left( z \right) = 1 - {z^{ - {n_a}}}.\;\;\;\;\;\;} \end{array}} \right. $ | (6) |
H1(z)与H4(z)相当于微分运算，H2(z)与H3(z)相当于积分运算。将微分H1放在第一级可防止溢出，微分H4放在最后一级有利于减小误差。H1与H4的区别在于存储数据的长度，前者为nb，后者为na。考虑到FPGA中滤波器系数的误差，为了避免反馈误差积累，提高运算精度，在中间运算时将系数放大了217倍，因而后端数据占用更长位宽，故将H1放在第一级以减少所需的存储空间。H3子系统在运行过程中一直在累加数据，会产生溢出，故需在无信号到来的时间段对H3的输出以及H4中存储的历史数据进行清零。
实验表明，此滤波器结构所实现的梯形成形算法，加上简单的空闲时间判定和基线恢复模块，即可准确实现长时测量。
2 堆积识别算法精确提取脉冲幅度依赖于准确判定脉冲的到达时刻。普遍的做法是利用两个通道的梯形成形，分别称为慢通道和快通道[12]。慢通道梯形脉宽较宽，用于幅度提取，以便获得较好的分辨率，但在高计数率下，为了减小堆积效应，此脉宽也不能太大。快通道梯形脉宽很窄，用于定位信号脉冲的到达时刻，取为输入信号上升时间的3~4倍即可。
当检测到相邻两个快通道信号之间间隔小于某一阈值时，即可判定信号产生了堆积。然而，快通道的梯形信号本身也有产生堆积的可能，使得最终获得的能谱产生误差。尤其在高计数率条件下，这种现象不可忽视。因而，本文研究了一种利用反卷积与快通道梯形成形相结合的脉冲定位方法，其原理如下：
归一化输入信号为
| $ {V_i}\left( t \right) = {{\rm{e}}^{-t/\tau }}u(t). $ | (7) |
将式(7) 离散化并作Z变换得到
| $ {V_i}\left( z \right) = \frac{1}{{1-d{z^{-1}}}}. $ | (8) |
理想的脉冲定位信号为δ(n)，其Z变换为 X(z)=1，故可得反卷积系统函数为
| $ H\left( z \right) = \frac{{X\left( z \right)}}{{{V_i}\left( z \right)}} = 1-d{z^{-1}}. $ | (9) |
其相应的差分方程为
| $ y\left( n \right) = x\left( n \right)-dx(n-1). $ | (10) |
可以看出，该算法可实现为一个简单的有限冲激响应(finite impulse response, FIR)滤波器，不会产生累计误差，且能够非常方便地在FPGA中实时实现。
图 3为脉冲到达时刻检测的ModelSim仿真结果，其激励信号为实采探测器输出信号，快通道成形时间为300 ns。可以看出，在输入指数脉冲与快通道梯形脉冲均发生堆积的情况下，利用反卷积算法能够准确地识别出堆积。
|
| 图 3 堆积识别算法仿真波形 |
3 实验与应用结果
本实验中所用晶体为硅粉，探测器为KETEK公司的AXAS-MH150型探测系统，含硅漂移探测器和前放，输出信号的上升时间为100 ns。数字多道部分的ADC型号为AD9255，采样频率80 MHz，FPGA为Xilinx公司的XC6SLX16。
3.1 能谱与分辨率对成形时间从600 ns到1 200 ns、计数率从约25 kcps (counts per second, 每s计数)到约270 kcps进行了能谱采集测试，数字多道设置为4 096道。在成形时间1.2 μs、180 kcps高计数率条件下，所得能谱如图 4所示。对能谱作Gauss拟合可得其均值与半高宽分别为1 429.4与35.164，则分辨率=半高宽/均值=2.46%。能谱均值对应铜的Kα射线能量8.05 keV，则等效的分辨率=8.05 keV×2.46%=198 eV。可见，本实验的能量分辨率达2.46%(198 eV)，能够有效分离Kα和Kβ特征。
|
| 图 4 8.05 keV铜Kα特征X射线能谱 |
不同输入计数率下，不同成形时间条件下所得能谱的能量分辨率曲线如图 5所示。可以看出，同一成形时间下，分辨率随计数率的升高而变差，这是由于计数率的增长使得脉冲堆积效应更加明显，对脉冲幅度的提取误差增大。在相同计数率条件下，梯形成形的上升时间越长，其对高频噪声的抑制越好；平顶时间越长，越能降低由探测器电荷收集时间所带来的弹道亏损[2]，故而能量分辨率越好。
|
| 图 5 不同计数率与成形时间下的能量分辨率 |
3.2 计数率校正
核信号通常服从Poisson分布，
| $ p\left( {x = n} \right) = \frac{{{\lambda ^n}{{\rm{e}}^{-\lambda }}}}{{n!}}. $ | (11) |
其中λ为Poisson分布参数。在时间t内无信号到来的概率为
| $ p(\Delta t \ge t) = {{\rm{e}}^{-\lambda }}. $ | (12) |
本文算法舍弃了堆积脉冲，换而言之，当且仅当脉冲的前后一定时间间隔内都没有其他脉冲的情况下，该脉冲的幅度才被记录下来，即有效记录的概率为：
| $ {p_{{\rm{record}}}} = {{\rm{e}}^{-2\lambda }}, $ | (13) |
| $ \lambda = {R_{{\rm{real}}}}{t_{{\rm{pile-up}}}}. $ | (14) |
其中：Rreal为真实计数率；tpile-up为脉冲堆积临界值，通常取为梯形的上升时间与平顶时间之和，即tb。
相应的计数率校正公式为
| $ {\rm{ }}{R_{{\rm{real}}}}.{{\rm{e}}^{-2{R_{{\rm{real}}}}{t_{{\rm{pile-up}}}}}} = {R_{{\rm{measure}}}}. $ | (15) |
实际测试中，每次增加2张铝箔遮挡X射线，能够使得计数率近似减小为原先的2/3。所用梯形成形时间为900 ns，测试结果如图 6所示。在常用计数率30~200 kcps范围内，可通过计算得到校正偏差最大为2.43%。
|
| 图 6 校正前后脉冲输入计数率与输出计数率的关系 |
4 结论
本文探讨了高计数率下的梯形成形算法，重点研究了合理的滤波器结构和脉冲到达时刻定位算法，在FPGA中实现了实时的数字脉冲多道谱仪的大部分功能。在X射线衍射仪系统上配合硅漂移探测器，在成形时间1.2 μs、180 kcps高计数率条件下，所得能量分辨率达2.46%(198 eV)。另外，本文所研究的梯形成形算法能够方便地应用于核领域不同的信号源和不同的探测器。
| [1] | Radeka V. Trapezoidal filtering of signals from large germanium detectors at high rates[J]. IEEE Transactions on Nuclear Science, 1972, 99(3): 525–539. |
| [2] | 肖无云, 魏义祥, 艾宪芸. 数字化多道脉冲幅度分析中的梯形成形算法[J]. 清华大学学报:自然科学版, 2005, 45(6): 810–812. XIAO Wuyun, WEI Yixiang, AI Xianyun. Trapezoidal shaping algorithm for digital multi-channel pulse height analysis[J]. Journal of Tsinghua University: Science and Technology, 2005, 45(6): 810–812. (in Chinese) |
| [3] | Bolic M, Drndarevic V. Digital gamma-ray spectroscopy based on FPGA technology[J]. Nuclear Instruments & Methods in Physics Research, 2002, 482(3): 761–766. |
| [4] | Bolic M, Drndarevic V. Processing architecture for high count rate spectrometry with NaI(Tl) detector [C]//IEEE Instrumentation and Measurement Technology Conference Record. Victoria, BC, Canada: IEEE Press, 2008: 274-278. |
| [5] | Cano-Ott D, Tain J L, Gadea A, et al. Pulse pileup correction of large NaI(Tl) total absorption spectra using the true pulse shape[J]. Nuclear Instruments & Methods in Physics Research, 1999, 430(2): 488–497. |
| [6] | Wong W H, Li H. A scintillation detector signal processing technique with active pileup prevention for extending scintillation count rates[J]. IEEE Transactions on Nuclear Science, 1998, 45(3): 838–842. DOI:10.1109/23.682647 |
| [7] | Haselman M D, Pasko J, Hauck S, et al. FPGA-based pulse pile-up correction with energy and timing recovery[J]. IEEE Transactions on Nuclear Science, 2012, 59(5): 1823–1830. DOI:10.1109/TNS.2012.2207403 |
| [8] | Ehrenberg J E, Ewart T E, Morris R D. Signal-processing techniques for resolving individual pulses in a multipath signal[J]. Journal of the Acoustical Society of America, 1978, 63(6): 1861–1865. DOI:10.1121/1.381926 |
| [9] | Bolic M, Drndarevic V, Gueaieb W. Pileup correction algorithms for very-high-count-rate gamma-ray spectrometry with NaI(Tl) detectors[J]. IEEE Transactions on Instrumentation & Measurement, 2010, 59(1): 122–130. |
| [10] | LIU Zhenzhou, CHEN Jinxiang. A Monte Carlo based technique for analysing gamma-ray spectra[J]. Measurement Science & Technology, 2008, 19(8): 085102. |
| [11] | Meng L J, Ramsden D. An inter-comparison of three spectral-deconvolution algorithms for gamma-ray spectroscopy[J]. IEEE Transactions on Nuclear Science, 2000, 47(4): 1329–1336. DOI:10.1109/23.872973 |
| [12] | Imperiale C, Imperiale A. On nuclear spectrometry pulses digital shaping and processing[J]. Measurement, 2001, 30(1): 49–73. DOI:10.1016/S0263-2241(00)00057-9 |
| [13] | 陈新光, 李翔宇, 孙义和. 粒子探测器读出电路数字滤波器设计[J]. 电子产品世界, 2010, 17(10): 29–31. CHEN Xinguang, LI Xiangyu, SUN Yihe. Particle detector readout circuits digital filter design[J]. Electronic Engineering & Product World, 2010, 17(10): 29–31. DOI:10.3969/j.issn.1005-5517.2010.09.006 (in Chinese) |