【技术领域】

本发明属于高能粒子探测信号处理和数字信号 处理领域， 具体涉及一种数 字化堆积 （pileup) 波形处理方法及其系统， 可应用于高能粒子探测数据获得系

【背景技术】

在大部分的高能粒子探测领域， 以及计算机断层成像 （Computed Tomography， 以下简称 CT )、 正电子发射断层成像 ( Positron Emission Tomography, 以下简称 PET )、 单光子发射计算机断层成像 （Single Photon Emission Computed Tomography, 以下简称 SPECT)等医疗影像领域， 数据获得 系统所采集、 处理的闪烁脉冲信号均是由高能粒子（如 γ射线、 X射线等）经闪 烁晶体转换成可见光， 再经光电转换器件转换得到可以观测的电信号 。

当多个高能粒子在短时间内被同一个探测器捕 获， 多个闪烁脉冲的信号会 叠加在一起， 这种信号被称为堆积信号 （pileup), 如图 1 所示。 使用传统方法 测量高能粒子时， pileup会对某些参数的测量造成较大误差。

首先， pileup会使后一次事件的脉冲电荷被计入到前一 次脉冲中，造成能量 测量的误差， 从而会造成能量谱和能量分辨率的恶化。 随着计数率的增加， 这 种恶化会变得越来越严重， 如图 2所示。

其次， 在部分高能粒子测量领域， 如 PET、 SPECT等， 脉冲信号的位置信 息是基于测量得到的四个 comer pulse的能量值， 通过计算这四个能量值之间的 相对关系获得的。 因此， pileup会造成位置测量上的错误， 会对位置谱的获得产 生较大的影响， 而且会随着计数率的增加变得越来越严重， 如图 3所示。 另外， pileup也会导致符合时间分辨率的恶化， 如图 4所示。

为了解决 pileup对测量带来的不利影响， 通常做法是通过波形能量或者波 形形状， 甄别出 pileup波形， 然后去除该波形。 这样的做法会减少高能粒子的 计数量， 增加统计噪声， 降低信噪比， 尤其对计数量较小的应用， 如 PET成像、 光子计数 （Photon Counting) CT成像会造成较大的影响。

此外， 也有一些研宄工作试图恢复 pileup 中脉冲的信息， 然而这些工作都 集中在恢复脉冲的能量信息， 以及由能量信息产生的位置信息上， 未见恢复时 间、余辉常数等其它信息。 （具体参见 " W. H. Wong, H. Li, A Scintillation Detector Signal Processing Technique with Active Pileup Prevention for Extending Scintillation Count Rates, IEEE Transaction on Nuclear Science, vol. 45, pp. 838-842. 1998")。 以上所有用于恢复 pileup 的工作均基于模拟电子技术， 对模拟脉冲波 形进行处理， 未针对数字脉冲波形。

目前存在一种 pileup信息恢复方法， 该方法根据每一帧数字脉冲中脉冲上 升沿数目判定该脉冲是否为 pileup, 因此为了提高判断准确率， 每一帧数字脉冲 长度必须得到限定。 其次， 该方法对 pileup波形的判断和对波形的分割， 是基 于整帧脉冲导数最大值的， 这会导致两个主要问题： （1 ) 该方法并不能判断每 帧数字文件中是否包含闪烁脉冲波形， 因此在该算法处理之前， 需要确保每个 数字文件都为包含闪烁脉冲的数字波形； （2 ) 当闪烁脉冲幅值较小时， 该方法 非常容易受脉冲下降沿上的噪声的影响，将正 常数字闪烁脉冲波形判别为 pileup 波形， 因此该方法在恢复脉冲位置信息上有明显的缺 陷 （位置信息的获取， 是 通过电阻网络将一个正常幅值的闪烁脉冲， 分为四个大小不一的闪烁脉冲， 通 过比较它们之间幅值 （能量） 相互比例获得的， 因此这四个脉冲幅值均较小）。 由于以上缺陷的存在， 该方法不能用于建立一套实时 pileup处理系统。 该波形由快速上升的上升沿和慢速下降 的下降沿构成。 上升沿的上升速度由闪烁晶体和光电转换器件 共同决定； 下降 沿的衰减速度， 由闪烁晶体的特性决定。

在不考虑噪声的情况下， 单次闪烁脉冲可考虑成由线性上升的上升沿和 指 数下降的下降沿构成的理想信号模型， 理想的闪烁脉冲波形如图 6所示， 其波 形模型如式 （1) 所示，

LineB

0 t <

LineK

LineB

V(t) = LineK xt + LineB ≤t <tp (1)

LineK

exp(ExpK xt + ExpB) t≥tp

其中， Z « 为上升沿直线的斜率， LineK >Q， Z ^fi是上升沿截距， 可为任 意数值，和上升沿开始时间是线性关系； 为衰减时间常数， ExpK<0, ExpB 可为任意数值， 与下降沿开始时间是线性关系， tp 为脉冲峰值时间。 因此， 一 个理想闪烁脉冲由四个参数 Z « 、 LineB、 ExpK、 ; β来描述， 闪烁脉冲信号 的开始时间、 峰值时间、 峰值幅值、 能量、 余辉常数等信息均可由这四个参数

(a) 脉冲开始时间 θ

LineB

tO (2) LineK

(b) 峰值时间 ίρ， 可解方程 （3) 获得近似解，

LineK xt + LineB = Qxp(ExpK χΐ-l· ExpB) (3)

(c) 峰值幅值

Vp = LineK tp + LineB (4) LineK x tp + LineB 1

E cxp(ExpK x ί + ExpB) ( 5 )

^ , LineB、 ExpK

2(tp + )

LineK

( 6)

ExpK

【发明内容】

本发明的目的在于提供一种数字化 pileup波形处理方法， 该方法能有效甄 别 pileup脉冲波形， 精确重建 pileup中单次事件脉冲波形， 恢复 pileup中单次 事件脉冲的信息， 有效地解决 pileup对能量、 时间、 位置、 余辉常数、 计数等 信息测量的影响， 能够提高测量结果的信噪比， 提高系统的灵敏度。 本发明还 提供了实现该方法的系统。

本发明提供的用于高能粒子探测数据获得系统 的数字化 pileup波形处理方 法， 其步骤包括：

( 1 )载入采集的数字脉冲波形， 根据数字脉冲波形中所包含的脉冲上升沿 和相邻脉冲上升沿时间间隔甄别出 pileup波形；

(2 ) 对甄别出的 pileup波形按脉冲上升沿和下降沿进行分割， 获得 pileup 波形中的各个单次事件脉冲波形；

(3 ) 对分割出的各个单次事件脉冲波形依次进行重 建。

本发明提供的数字化 pileup波形处理系统， 包括 pileup波形甄别模块、 脉 冲波形分割模块和脉冲波形重建模块；

pileup 波形甄别模块用于导入数字脉冲波形， 获取数字脉冲波形的导数序 列， 根据导数序列判断数字脉冲波形中所包含的脉 冲上升沿， 进而判断该数字 脉冲波形是否为 pileup波形， 并将 pileup波形传递送给脉冲波形分割模块； 脉冲波形分割模块对 pileup波形甄别模块甄别出的 pileup波形按脉冲上升 沿和下降沿进行分割， 获得 pileup波形中的各个单次事件脉冲波形； 脉冲波形重建模块用于对脉冲波形分割模块分 割出的每个单次事件脉冲进 行脉冲波形重建， 将每个单次脉冲波形重建为一个快速线性上升 的上升沿和一 个指数衰减的下降沿。

本发明提供的数字化 pileup波形处理方法首先利用数字脉冲波形的导 数判 定脉冲上升沿区域； 根据相邻脉冲上升沿之间的时间间隔判定该脉 冲是否为 pileup波形； 将判定为 pileup的数字波形按照脉冲上升沿和下降沿进行 分割； 然 后利用分割后的波形依次重建出 pileup 中的各个单次事件脉冲波形， 各单次脉 冲信息均可从重建的波形中提取。 该方法能有效甄别 pileup脉冲波形， 精确重 建 pileup中各个单次事件脉冲波形， 准确恢复 pileup中各个单次事件脉冲的能 量、 时间、 幅值、 衰减时间常数等信息， 提高系统在高活度下的计数率， 减小 数据丢失问题， 提高测量结果的信噪比， 提高系统的灵敏度。 本发明还提供了 实现该方法的系统。

【附图说明】

图 1为典型的 pileup波形；

图 2为 pileup对能量谱和能量分辨率的影响；

图 3为 pileup对位置谱的影响；

图 4为 pileup对时间分辨率的影响；

图 5为典型的闪烁脉冲波形；

图 6为理想的闪烁脉冲模型；

图 7为本发明数字化 pileup波形处理方法的流程图；

图 8为本发明数字化 pileup波形处理系统的系统结构图；

图 9为本发明从 pileup波形中恢复各个单次事件脉冲的示意图； 图 10为本发明对位置谱的效果， 其中 10 (a) 为未使用本发明获得的位置 谱， 10 (b) 为使用本发明获得的位置谱。

【具体实施方式】

下面结合附图和实例对本发明的技术方案作进 一步的详细说明。

如图 7所示， 本发明方法包括以下步骤：

(1) 载入采集的数字脉冲波形， 从中甄别出 pileup波形， 具体步骤如下：

(1.1) 载入采集的数字脉冲波形。 该波形可以是事先采集好的， 或者实时 采集获得的。

(1.2) 数字脉冲波形滤波。 根据数字脉冲波形特征和噪声特征， 可以选择 在进行下一步处理前进行数字滤波， 去除噪声对下一步处理的影响， 可以采用 均值滤波、 中值滤波等方法； 亦可不对数字脉冲波形滤波。

(1.3) 数字脉冲波形求导。 对步骤 （1.2) 处理后的数字脉冲波形求导， 获 得数字脉冲波形的导数序列。

离散序列的求导存在多种经典算法， 如， 采用后向求导算法， 定义数字脉 冲波形 DigitalPulse (1： n)， 其中， n为数字脉冲波形长度， n为整数； 定义求 导步长为 step, step≥l且不大于脉冲时间常数 /采样频率， 其中， step为整数， 则 其导数序列为：

DigitalPulse i + step) DigitialPulse i)

DeriPulseii) step (7)

0 i = n-step + l,-",n 其中， DeriPulse (i) 表示数字脉冲波形 DigitalPulse (i) 对应的导数。

(1.4) 根据数字脉冲波形的导数序列判断数字脉冲波 形所包含的脉冲上升 沿， 并统计该数字脉冲波形中所包含的脉冲上升沿 个数和相邻脉冲上升沿之间 定义阈值 Threshold, 其中， Threshold>0， Threshold依据该探测器数字脉冲 波形幅值和噪声水平经验设定， 原则上该参数数值尽量取小， 但要大于探测器 噪声的导数。导数序列中大于 Threshold的一段连续区域即对应脉冲的上升沿所 在区域。 统计导数序列中所有大于 Threshold的连续区域的段数， 该段数即为该 数字脉冲波形所包含的脉冲上升沿个数， 并计算各相邻脉冲上升沿之间的时间 间隔。

( 1.5 ) 判断数字脉冲波形是否为 pileup波形。 若该数字脉冲波形中只存在 1个脉冲上升沿， 则该数字脉冲波形为非 pileup波形； 否则， 判断导数序列中是 否存在两个相邻脉冲上升沿之间的时间间隔超 过某一固定时间长度 T, 若存在， 则该数字脉冲波形为 pileup波形， 否则， 不是 pileup波形， 其中， T>0， Τ具体 数值根据探测器衰减时间常数设定， 一般为探测器衰减时间常数的整数倍， 如 3 倍等。

( 2 ) pileup波形分割。 对步骤（1 )甄别出的 pileup波形进行处理， 根据脉 冲上升沿和下降沿对 pileup波形进行分割， 获取 pileup波形中的各个单次事件 脉冲波形。 具体过程为：

( 2.1 ) 在 pileup波形的导数序列中， 确定与脉冲上升沿区域相对应的子序 列。

( a)对数字脉冲波形 DigitalPulse ( 1： n)求导， 获得其导数序列 DeriPulse ( 1： n) o 可以采用步骤 （1.3 ) 的方法获得该数字脉冲波形的导数序列。

(b)设定阈值 deri— threshold,上升沿区域最小长度 M，其中， deri— threshold 根据该探测器数字脉冲波形幅值及其噪声水平 设定， 原则上该阈值数值尽量取 小， 但应大于 0。 M根据脉冲典型的上升时间和脉冲采样频率确� �， 一般地， M 不小于 0.5ns/采样频率、 不大于 10ns/采样频率。 计算导数序列 DeriPulse ( 1： n) 中所有大于 deri— threshold的连续的子序列的长度， 统计这些连续的子序列中长 度大于 M的子序列的个数 N， N为整数， 长度大于 M的子序列即对应脉冲上升 沿区域， 即 N为该导数序列中所含有的脉冲上升沿的总个� �。 记录导数序列中 所有的脉冲上升沿区域对应的子序列为 DeriPulse ( s _i: Si+Asi),其中， i=l， 2， ...， N， i为上升沿序号， _Sl 、 _§1 +八 ₈₁ 分别为第 i个上升沿在导数序列中的起始位置和 终止位置， i、 、 A _Si 均为整数， 该子序列的长度为 Asi+l ;

(2.2 )确定数字脉冲波形中每个脉冲上升沿的位置� � 根据步骤 （2.1 ) 获得 的 N个脉冲上升沿区域对应的子序列 DeriPulse ( ¾： +△ )， j=l， 2， ...， N， 获得数字脉冲波形中对应的上升沿区域对应的 子序列 DigitalPulse ( ¾： +△ )， 对数字脉冲波形中每个子序列 DigitalPulse ( ¾： ¾+Δ¾) 求取其极小值对应的位 置 minj和极大值对应的位置 m _aXj ， 两者之间的波形即为一个波形的上升沿， 记 DigitalPulse (minj: maxj)为该数字脉冲波形中第 j个脉冲上升沿， 其中， ≤ηώ¾ <maxj≤Sj+Asj， minj、 maxj均为整数；

(2.3 ) 确定数字脉冲波形中每个脉冲下降沿的位置。 步骤 （2.2 ) 确定了每 个脉冲的上升沿， 两个脉冲上升沿之间的部分即为前一个脉冲的 下降沿的波形， 对于第 j个脉冲， l≤j <N， 其下降沿为 DigitalPulse (maxj : min _j+1 ), 对第 N个 脉冲， 从前至后在 DigitalPulse (max _N : n)找到第一个数值小于 0的采样点对应 的位置为结束点 ending, 若 DigitalPulse (max _N : n) 中无小于 0的采样点， 则 ending等于 n+l， 则第 N个脉冲的下降沿为 DigitalPulse (max _N : ending-1 ), 其 中， max _N <ending, ending为整数。

( 3 )对分割出的各个单次事件脉冲波形依次进行� �建。 每个单次脉冲波形 可以认为是由一个快速上升的上升沿和一个指 数衰减的下降沿构成。 步骤 （2 ) 已经分割出每个单次事件脉冲波形的上升沿和 下降沿， 需要对每个单次事件脉 冲波形的上升沿和下降沿分别进行重建， 获得其模型参数。 由于当次到达的事 件脉冲波形的下降沿会叠加到下一个到达的事 件脉冲波形中， 在对下一个到达 的事件脉冲波形进行重建时， 需要去除当次的事件脉冲波形的影响。 对每个单 次事件脉冲波形依次进行如下具体重建步骤：

(3.1 ) 上升沿重建。 脉冲上升沿以线性公式 （8) 表示，

y(x) = LineK x x + LineB ( 8) 其中， Z « 为上升沿直线的斜率， LhieK X) , Z ^fi是上升沿截距， 可为任 意数值， 和上升沿开始时间是线性关系， X为数字脉冲序列采样位置， y为在 X 位置上的数字脉冲幅值。

在脉冲上升沿抽取不少于 P个采样点， 其中， P为整数， P≥2且小于等于上 升沿长度，对这 P个点按式（8)进行线性拟合获取上升沿重建参 数 Z « 和 Z ^fi。

拟合时， 可以采用最小二乘法进行拟合。 （具体参见 "梁家惠. 用最小二乘 法进行直线拟合的讨论. 工科物理， pp.11-15, 1995年 3月"）

(3.2) 下降沿重建。 脉冲下降沿以指数公式 （9) 表示，

y(x) = cxp(ExpK χ χ + ExpB) (9) 其中， 为衰减时间常数， ExpK < 0， E^fi可为任意数值， 与下降沿开 始时间是线性关系。

在脉冲下降沿抽取不少于 Q个采样点， 其中， Q为整数， Q≥2且小于等于 下降沿长度， 对这 Q个点进行按式 （9) 指数拟合获取下降沿重建参数 ； ^和

ExpB。 该指数拟合可以转化为线性拟合。

(3.3 ) 去除当前重建的脉冲余辉对后续脉冲重建的影 响。 利用当前脉冲下 降沿重建参数 ExpK和 ExpB， 计算该脉冲下降沿在后续脉冲波形中的脉冲幅 值， 将之从后续每个事件脉冲中减去。 本发明提供的用于高能粒子探测数据获得系统 的数字化 pileup波形处理系 统结构图如图 8所示。该系统包括 pileup波形甄别模块 100、脉冲波形分割模块 200和脉冲波形重建模块 300。

pileup波形甄别模块 100用于导入数字脉冲波形，获取数字脉冲波形 的导数 序列， 根据导数序列判断数字脉冲波形中所包含的脉 冲上升沿， 进而判断该数 字脉冲波形是否为 pileup波形，并将 pileup波形传递送给脉冲波形分割模块 200。

该模块分为 4个子模块， 分别为波形导入模块 110， 波形滤波模块 120， 波 形求导模块 130和上升沿分析模块 140。波形导入模块 110用于导入采集的数字 脉冲波形， 并传送给波形滤波模块 120; 波形滤波模块 120用于对波形导入模块 110导入的数字脉冲波形进行滤波， 去除噪声， 并将滤波后的数字脉冲波形传送 给波形求导模块 130;波形求导模块 130用于对波形滤波模块 120滤波获得的数 字脉冲波形进行求导， 获得数字脉冲波形的导数序列， 并将该导数序列传送给 上升沿分析模块 140;上升沿分析模块 140根据波形求导模块 130获得的导数序 列判断数字脉冲波形中所包含的脉冲上升沿， 并统计脉冲上升沿个数， 根据在 某段连续时间内存在的脉冲上升沿的个数判断 数字脉冲波形是否为 pileup波形。

上升沿分析模块 （140) 中， 导数序列中大于预设的阈值 Threshold的一段 连续区域即对应脉冲的上升沿所在区域。 上升沿分析模块 （140) 中， pileup波 形的具体判断过程为： 若数字脉冲波形中只存在 1 个脉冲上升沿， 则该数字脉 冲波形为非 pileup波形； 否则， 判断是否存在任意两相邻脉冲上升沿之间的时 间间隔小于预定的固定时间长度 T, 若存在， 则该数字脉冲波形为 pileup波形。

所述脉冲波形分割模块 （200) 先确定所述甄别出的 pileup波形的上升沿和 下降沿， 并据此将所述 pileup波形分割为多个独立的单次事件脉冲， 每个单次 事件脉冲包含一个上升沿和其后的下降沿。 其中， 确定 pileup波形的上升沿和下降沿的具体过程如下：

(A)在所述 pileup波形的导数序列中， 确定与脉冲上升沿区域相对应的子 序列

计算 pileup波形的导数序列中所有大于阈值 deri— threshold的连续的子序列 的长度， 长度大于阈值 M的子序列即对应脉冲上升沿区域；

(B) 确定 pileup波形中每个脉冲上升沿的位置

对每个所述对应的子序列求取其极小值和极大 值在 pileup波形中对应的位 置， 两者之间的波形即为 pileup波形的一个上升沿；

( C)确定数字脉冲波形中每个脉冲下降沿的位置� �� 两个脉冲上升沿之间的 部分即为这两个脉冲中前一个脉冲的下降沿波 形。

脉冲波形重建模块 300用于对脉冲波形分割模块 200分割出的每个单次事 件脉冲进行脉冲波形重建， 将每个单次脉冲波形重建为一个快速线性上升 的上 升沿和一个指数衰减的下降沿。

该模块分为 3个模块， 分别是上升沿重建模块 310， 下降沿重建模块 320和 余辉消除模块 330。上升沿重建模块 310用于对脉冲波形分割模块 200分割出的 每个单次事件脉冲上升沿部分进行线性拟合， 获得上升沿重建参数； 下降沿重 建模块 320用于对脉冲波形分割模块 200分割出的每个单次事件脉冲下降沿部 分进行指数拟合， 获得下降沿重建参数； 余辉消除模块 330利用下降沿重建模 块 320获得的下降沿重建参数， 计算该脉冲下降沿在后续脉冲波形中的脉冲幅 值， 将该脉冲幅值从后续的每个单次事件脉冲中减 去， 消除当前脉冲波形下降 沿对随后的各个单次事件脉冲波形重建的影响 。

上升沿重建模块 （310) 中获得上升重建参数的过程为：

在脉冲上升沿抽取不少于 P个采样点， 其中， P为整数， P≥2且小于等于上 升沿长度， 对采样点按式（8)进行线性拟合获取上升沿重� ��参数 Z « 和 Z ^fi , 进行重建：

y(x) = LineK x x + LineB (8) 其中， Z « 为上升沿直线的斜率， LhieKX), Z ^fi是上升沿截距， 可为任 意值， 为数字脉冲序列采样位置， 为在 位置上的数字脉冲幅值；

下降沿重建模块 （310) 中获得下降重建参数的过程为：

在脉冲下降沿抽取不少于 Q个采样点， 其中， Q为整数， Q≥2且小于等于 下降沿长度，对采样点按式（9)进行指数拟合� ��取下降沿重建参数 ¾ ^和 E _X pB， 进行重建：

y(x) = Qxp(ExpK χ + ExpB) (9) 其中， Ej ^T为衰减时间常数， ΕχρΚ<ϋ， ψβ可为任意值。

实例：

本发明提出的一种 pileup波形处理方法及其系统涉及到若干参数， 这些参 数需针对具体处理数据进行调节以达到良好的 性能。 此处列出本应用实例处理 数据的参数：

步骤 （1.1) 中导入数字脉冲为使用 LYSO晶体和 HamamtsuR9800 PMT采 集得到闪烁脉冲， 典型波形为图 5所示。 采样率为 lOGSps, 能量为 511KeV高 能光子产生脉冲平均峰值约为 300mV，平均脉冲上升沿时间约为 2ns，探测器衰 减时间常数为 47ns。

步骤 （1.2) 采用均值滤波， 滤波窗尺寸为 100;

步骤 （1.3) 采用后向求导， 求导步长 step=50;

步骤 （1.4) 对上升沿区域进行判断， 阈值 Threshold=0.0002;

步骤 （1.5) 判断数字脉冲波形是否为 pileup波形， 设定常数 T=140ns; 步骤 (2.2 ) (b) 阈值 deri— threshold 设定为 0.0001， 上升沿最小区域长度 图 9为利用本发明提出的方法从 pileup波形中恢复各个单次事件脉冲的结 果。 图 10为本发明提出的方法对位置谱的对比结果图� �� 其中， 10(a)为未利用本 发明获得的位置谱， 10(b)为利用本发明获得的位置谱。

本发明的方法和系统可应用于高能粒子探测数 据获得系统中， 光子计数 ( Photon Counting) CT系统中，正电子断层扫描 ( Positron Emission Tomography, PET) 系统中等。

本发明不仅局限于上述具体实施方式， 本领域一般技术人员根据本发明公 开的内容， 可以采用其它多种具体实施方式实施本发明， 因此， 凡是采用本发 明的设计结构和思路， 做一些简单的变化或更改的设计， 都落入本发明保护的 范围。