辐射测量中最重要的两个指标, 一个是能量分辨率, 另一个则是计数率。 但能量分辨率和计数率是相互矛盾的, 改善能量分辨率将降低计数率, 相反提高计数率将损失能量分辨率^{[1, 2, 3]}。 探测器的能量分辨率由多个因素决定, 比如探测信号的统计涨落、 信号处理器的噪声、 外部干扰以及温度等。 统计涨落决定了给定探测器能量分辨率的理论极限, 而其他因素的影响则可以通过适当的噪声滤除和电子技术来降低^{[4, 5, 6]}。

以往关于能量分辨率的研究大多利用谱反卷积对获取到的能谱进行后处理, 从而降低特征峰的半高宽。 这些后处理方法是通常基于将获取到的能谱建模为输入能谱和探测器响应函数这两个随机变量的函数, 而探测器的响应函数又可以建模为输入脉冲高度与输出脉冲高度的联合概率分布。 目前至少存在三种对谱线进行后处理的方法可以提高能量分辨率, 包括使用正则化、 极大似然和最大熵。 然而, 这三种后处理方法都涉及了为光谱反褶积建模探测器的能量响应函数, 因此这些方法将是一项计算量非常大的任务^[7]。

传统的X射线光谱数据处理平台采用的是多道谱的成谱方法, 其主要特征在于谱图中每一个道址上计数的增量是基于单个脉冲幅度的, 也就是说, 每一个脉冲幅度被接收, 相应的道址上的计数值就会加一。 早在2005年, 国外就有研究人员初步提出了多脉冲成谱的相关方法^[8], 但国内尚无相关成果公布。 笔者在前期研究工作^{[9, 10]}中已经提出了脉冲剔除法和脉冲修复法处理突变脉冲得到真实可靠的脉冲幅度信息, 并因此降低了统计涨落。 本文在上述研究基础上提出一种X射线光谱数据处理平台的优化设计, 采用可在线实施的多脉冲局部平均(MPLA)成谱技术, 选择一个恰当范围的平均窗口对多个脉冲幅度进行平均, 平均后的幅度值被接收, 再将相应的道址上的计数值加一, 该技术的实施使得X射线光谱数据处理平台得以优化。 本文通过理论推导揭示了MPLA概率密度变换的基本理论特征, 也通过实验验证了MPLA算法的特征, 证明了在具有正态分布PDF的频谱峰值的典型情况下, 即使仅对两个脉冲高度进行平均, 变换后峰的FWHM也会变窄。

MPLA算法与传统多道谱的根本区别在于前者使用多个脉冲高度的平均值来确定需要增加计数的道址, 而后者则是每读取一个脉冲幅度值, 就给对应道址上的计数加一。 MPLA算法的核心在于多个脉冲幅度的平均, 但这里的平均并不是针对所有的脉冲幅度。 如果对所有测量到的脉冲幅度都进行平均, 那么最终得到的谱图就会完全失去其原有的统计特征, 也会破坏不同元素特征峰之间的计数差异, 因此MPLA算法中的平均必须是针对选定窗口内一定脉冲幅度数量进行平均, 再通过平均得到的脉冲幅度值来确定增加计数的道址。

MPLA算法涉及两项可变参数, 一是平均窗口的大小, 本文以参数r表示, 另一项参数则是每一次平均的脉冲幅度数量, 本文以参数n表示。 在下文对MPLA算法原理的描述中, 我们取平均窗口r等于2, 每一次参与平均的脉冲幅度数量为2。

在MPLA算法的执行流程包含以下几个步骤, 如图1所示。 首先读取第一个脉冲幅度时定位一个平均窗口, 当r为2时窗口内包含的道址数为5, 脉冲幅度读取成功后更新当前平均窗口的幅度和脉冲个数, 每次更新后即对平均窗口内的脉冲个数进行判断, 当其小于预设的参数n时, 则继续读取下一个脉冲幅度, 当平均窗口内的脉冲幅度数量等于参数N时则对相应平均窗口内的脉冲幅度进行平均, 得出的平均数即为需要更新计数的道址, 然后再对取平均值的窗口内脉冲幅度和脉冲个数进行清零。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 MPLA算法实现流程Fig.1 MPLA algorithm implementation flow

为了详细阐述MPLA算法的实现过程, 下文将以图文结合的方式呈现一个MPLA实现示例。 在当前示例中, MPLA算法的平均窗口r设置为2, 每一次参与平均的脉冲幅度数量设置为2。 在图2所示的MPLA算法步骤一中, 读取到的第一个脉冲幅度值为660, 平均窗口的中心则为660, 而窗口大小为2则代表窗口的范围被限定在658~622之间。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 MPLA算法实现示例步骤一Fig.2 Step 1 of MPLA algorithm implementation example

此时活动平均窗口仅有一个, 并且该窗口内可用于平均的脉冲个数小于2, 开始执行MPLA算法示例的步骤二, 如图3所示, 本次读取的脉冲幅度为650, 平均窗口中心移到650, 平均窗口范围为648~652, 该范围与前一个脉冲所处的范围并无重叠, 这就产生了一个新的平均窗口, 对该窗口内的脉冲幅度和以及脉冲个数进行更新后再次判断是否有达到平均条件的窗口, 若没有继续执行步骤三。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 MPLA算法实现示例步骤二Fig.3 Step 2 of MPLA algorithm implementation example

在执行完步骤一和步骤二后可以发现, 在执行多脉冲平均的过程中, 可能会存在多个活动的平均窗口, 如果下一个读取的脉冲高度正好在包含平均和的两个通道之间, 则会随机将脉冲高度添加到两个通道之一; 如果下一个读取的脉冲幅度所产生的平均窗口覆盖了该区域内已经存在的某个脉冲幅度和, 那么则用新的脉冲幅度与原有脉冲幅度相加, 并将该平均窗口内的脉冲个数加一。

如图4所示, 在读幅度为659的脉冲之前, 所示脉冲幅度和的区域内已经有650和660两个脉冲幅度和, 每一个脉冲幅度和对应的脉冲数量都为一, 还未达到窗口内取平均的条件。 下一个读取的脉冲幅度659产生的平均窗口范围为657~661, 刚好包含了原本存在660, 因此更新后的脉冲幅度和就为1 319, 并且该窗口内脉冲个数为2, 达到了取平均的条件。 通过脉冲幅度和除以脉冲个数计算出平均脉冲幅度为660, 已经被平均的脉冲幅度和以及脉冲个数清零, 其余窗口内的计数保持不变, 同时将取平均值得到的脉冲幅度所对应的道址上的计数加一, 如图5所示。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 MPLA算法实现示例步骤三Fig.4 Step 3 of MPLA algorithm implementation example

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 MPLA算法实现示例步骤四Fig.5 Step 4 of MPLA algorithm implementation example

理论上来说, 传统的多道谱(MCA)体现的是能量上的概率密度(pdf), 每一个脉冲幅度都对最终得到的谱图的计数率有所贡献, 而MPLA将这种概率密度作为输入, 并通过在指定范围内的平均将其转换为相同能量上新的概率密度。

设r为窗口平均参数, n为平均参数的脉冲高度数。 X₀是一个通过多脉冲幅度和本地取平均值得到脉冲幅度, 且X₀的范围为[X₀-r, X₀+r]。 令概率密度函数为f(x), 累积分布函数为F(x)。 若X₁是一个随机变量描述的落在数理统计范围内的第二脉冲幅度, 则X₁与X₀有相同的概率密度函数, 但区间为[X₁-r, X₁+r]。 道址上的计数将在采集到2个数理统计范围内的脉冲幅度时进行递增, 相应的随机变量P₀=(X₀+X₁)/2。 输入概率密度函数f(x)将完全指定X₀和X₁的分布, 以及MPLA转换产生的P₀的概率密度函数。

因为f(x)是关于μ 对称的, 所以对于任意值b可以得出f(μ -b)=f(μ +b)。 同理, 如果f(x)是关于μ 对称的, 那么f_p(y)也关于μ 对称。 如果f(x)是对称的, 并且不断的趋于平均值μ , 那么P₀的方差小于X₀。 推导过程如下:

在不损失一般性的前提下, 我们假设为μ =0。 该理论简化后为式(1)所示。

EX0+X122< E[X02](1)

表达式(1)的左侧可进一步展开为式(2)所示。

EX0+X122=E14(X02+X12+2X0X1)=14(E[X02]+E[X12]+2E[X0X1])≤14(2E[X02]+2E[X0X1])=12(E[X02]+E[X0X1])(2)

不等式存在的原因是X₁是X₀的限制量, 要证明P₀的方差小于X₀仅需要证明E[X₀X₁]< E[ X02]。

E[X0X1]=∫-∞+∞f(x0)E[x0X1|X0=x0]dx0=∫-∞+∞f(x0)x0E[X1|X0=x0]dx0(3)

由于f(x)是对称的, 并且一直朝着均值递增(此处我们假定了均值为0), 我们可以得出假如x₀=μ =0那么E[X₁|X₀=x₀]=x₀, 同样的, 假如x₀≠ μ 那么|E[X₁|X₀=x₀]|< |x₀|, 在此基础上, 式(3)可进一步推导出式(4)。

E[X0X1]=∫-∞+∞f(x0)x0E[X1|X0=x0]dx0< ∫-∞+∞f(x0)x02dx0=E[X02](4)

上述推论表明, 不管MPLA的参数值如何, MPLA转换都不会改变分布, 并保持对称性。 如果原始分布是单个高斯峰, 则所得分布将具有相同的平均值, 维持峰值位置, 并且也将是对称的。 理想的MPLA转换降低了原始分布的半高宽, 锐化了峰形。 对于单个高斯峰, 半高宽与分布的方差有关。

MPLA算法是一种在线实时进行的脉冲处理方法, 其实施平台在硬件上包括前端的探测器、 后端的PC处理软件以及核心部分的数字信号处理器。 如前文所述, MPLA算法本身并没有使用范围的限制, 它是将传统的概率密度作为输入, 并通过在指定范围内的平均将其转换为相同能量上新的概率密度, 其本质上就是一种概率密度的转换, 并没有对探测器有任何的限定。 本文所采用的探测器是AMPTEK生产的高性能FAST-SDD探测器。 数字信号处理器包括信号转换、 多级放大、 模数转换等前端处理单元、 脉冲幅度分析单元、 MPLA单元以及脉冲成谱单元, 其结构框图如图6所示。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 数字信号处理器结构图Fig.6 Structure chart of the digital signal processor

数字信号处理器的前端处理部分不是本文讲述的重点, 因此下文中将重点介绍MPLA单元的硬件实现过程。 MPLA单元在FPGA中实现, 采用的芯片是Xilinx公司Spartan-3家族的XC3S400-4TQG144C。 如图6所示, MPLA单元是整个数字信号处理器的核心部分, 其算法的实现主要包含了算法控制、 平均窗口定位、 窗口内平均以及相应道址上的计数更新四个部分。 具体的实现流程在前文中已经有详细的描述, 此处不再赘述。 在MPLA单元处理完成后, 更新的计数则通过SPI总线传递给MCU, MCU与上位机之间通过CAN总线进行通信。 此外, MCU的另一SPI总线连接到前端处理电路中, 用于控制前端电路中的放大参数和偏移量的调整。

由图7可以看出经过MPLA算法处理后的谱线半高宽明显降低, 也就是说能量分辨率有了明显的提高, 但与此同时, 多脉冲窄窗口求平均的成谱方法也造成了计数率的降低。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 采用MPLA算法前后得到的谱图Fig.7 Spectrum obtained before and after MPLA

由于放射性衰变是随机的, 所以计数率总是在围绕一个平均值不停地摆动, 即放射源在每单位时间内发生衰变的原子数目是不同的, 时多时少, 有起有伏, 但是它比较集中地在某一范围内波动, 而这种现象就是放射性衰变的统计涨落。 由概率统计理论可知, 当统计平均值较大时, 核衰变是服从正态分布的, 这也是核衰变的重要统计特点。 在实验环节中, 以铁矿样品为测量对象, 采用MPLA算法并进行快速多脉冲成谱处理后的结果如图7所示。

图中黑色谱线是通过传统成谱方法得到的, 红色谱线是采用MPLA算法处理后得到的, 相比之下, 经MPLA算法处理后的谱线有效地降低了特征峰的半高宽。 这就验证了前文理论推导部分得出的结论, 理想的MPLA转换降低了原始分布的半高宽, 锐化了峰形。

给出了多脉冲局部平均变换的理论描述及实现原理, 证明了MPLA算法的一些基本特性, 对具有单个峰的对称分布, MPLA变换总是保持均值不变并减小FWHM。 理论推导和实验结果表明, MPLA算法具有降低谱峰的FWHM的优点。 但与此同时, MPLA的一个缺点是与传统成谱方法相比, 多脉冲平均的过程中使得计数减少。 MPLA谱中计数减少的总数受MPLA算法中平均窗口内脉冲幅度的个数控制。 有很多脉冲处理技术可以弥补这种计数损失, 例如, 可以收集大量的脉冲高度测量, 然后重复、 随机地进行采样以应用MPLA变换。 笔者将在接下来的工作中对这些脉冲处理技术进行更深入的研究, 以期为MPLA算法找到一种非常契合的脉冲处理技术来补偿计数率损失。 即便存在计数损失的缺陷, 但不可否认, MPLA算法为降低探测器噪声提供了一种实时、 高效、 通用的处理方法, 并且具有可证明的理论保证, 这对提高探测器的能量分辨率具有极大的应用意义。