背景技术
在大部分的高能粒子探测领域,以及计算机断层成像(Computed Tomography,以下简称CT)、正电子发射断层成像(Positron Emission Tomography,以下简称PET)、单光子发射计算机断层成像(Single Photon Emission Computed Tomography,以下简称SPECT)等医疗影像领域,数据获得***所采集、处理的闪烁脉冲信号均是由高能粒子(如γ射线、X射线等)经闪烁晶体转换成可见光,再经光电转换器件转换得到可以观测的电信号。
当多个高能粒子在短时间内被同一个探测器捕获,多个闪烁脉冲的信号会叠加在一起,这种信号被称为堆积信号(pileup),如图1所示。使用传统方法测量高能粒子时,pileup会对某些参数的测量造成较大误差。
首先,pileup会使后一次事件的脉冲电荷被计入到前一次脉冲中,造成能量测量的误差,从而会造成能量谱和能量分辨率的恶化。随着计数率的增加,这种恶化会变得越来越严重,如图2所示。
其次,在部分高能粒子测量领域,如PET、SPECT等,脉冲信号的位置信息是基于测量得到的四个corner 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处理***。
典型的闪烁脉冲波形如图5所示,该波形由快速上升的上升沿和慢速下降的下降沿构成。上升沿的上升速度由闪烁晶体和光电转换器件共同决定;下降沿的衰减速度,由闪烁晶体的特性决定。
在不考虑噪声的情况下,单次闪烁脉冲可考虑成由线性上升的上升沿和指数下降的下降沿构成的理想信号模型,理想的闪烁脉冲波形如图6所示,其波形模型如式(1)所示,
其中,LineK为上升沿直线的斜率,LineK>0,LineB是上升沿截距,可为任意数值,和上升沿开始时间是线性关系;ExpK为衰减时间常数,ExpK<0,ExpB可为任意数值,与下降沿开始时间是线性关系,tp为脉冲峰值时间。因此,一个理想闪烁脉冲由四个参数LineK、LineB、ExpK、ExpB来描述,闪烁脉冲信号的开始时间、峰值时间、峰值幅值、能量、余辉常数等信息均可由这四个参数计算获得,计算公式如下,
(a)脉冲开始时间t0
(b)峰值时间tp,可解方程(3)获得近似解,
LineK×t+LineB=exp(ExpK×t+ExpB) (3)
(c)峰值幅值Vp
Vp=LineK×tp+LineB (4)
(d)能量E
(e)余辉常数τ
具体实施方式
下面结合附图和实例对本发明的技术方案作进一步的详细说明。
如图7所示,本发明方法包括以下步骤:
(1)载入采集的数字脉冲波形,从中甄别出pileup波形,具体步骤如下:
(1.1)载入采集的数字脉冲波形。该波形可以是事先采集好的,或者实时采集获得的。
(1.2)数字脉冲波形滤波。根据数字脉冲波形特征和噪声特征,可以选择在进行下一步处理前进行数字滤波,去除噪声对下一步处理的影响,可以采用均值滤波、中值滤波等方法;亦可不对数字脉冲波形滤波。
(1.3)数字脉冲波形求导。对步骤(1.2)处理后的数字脉冲波形求导,获得数字脉冲波形的导数序列。
离散序列的求导存在多种经典算法,如,采用后向求导算法,定义数字脉冲波形DigitalPulse(1:n),其中,n为数字脉冲波形长度,n为整数;定义求导步长为step,step≥l且不大于脉冲时间常数/采样频率,其中,step为整数,则其导数序列为:
其中,DeriPulse(i)表示数字脉冲波形DigitalPulse(i)对应的导数。
(1.4)根据数字脉冲波形的导数序列判断数字脉冲波形所包含的脉冲上升沿,并统计该数字脉冲波形中所包含的脉冲上升沿个数和相邻脉冲上升沿之间的时间间隔。
定义阈值Threshold,其中,Threshold>0,Threshold依据该探测器数字脉冲波形幅值和噪声水平经验设定,原则上该参数数值尽量取小,但要大于探测器噪声的导数。导数序列中大于Threshold的一段连续区域即对应脉冲的上升沿所在区域。统计导数序列中所有大于Threshold的连续区域的段数,该段数即为该数字脉冲波形所包含的脉冲上升沿个数,并计算各相邻脉冲上升沿之间的时间间隔。
(1.5)判断数字脉冲波形是否为pileup波形。若该数字脉冲波形中只存在1个脉冲上升沿,则该数字脉冲波形为非pileup波形;否则,判断导数序列中是否存在两个相邻脉冲上升沿之间的时间间隔超过某一固定时间长度T,若存在,则该数字脉冲波形为pileup波形,否则,不是pileup波形,其中,T>0,T具体数值根据探测器衰减时间常数设定,一般为探测器衰减时间常数的整数倍,如3倍等。
(2)pileup波形分割。对步骤(1)甄别出的pileup波形进行处理,根据脉冲上升沿和下降沿对pileup波形进行分割,获取pileup波形中的各个单次事件脉冲波形。具体过程为:
(2.1)在pileup波形的导数序列中,确定与脉冲上升沿区域相对应的子序列。
(a)对数字脉冲波形DigitalPulse(1:n)求导,获得其导数序列DeriPulse(1:n)。可以采用步骤(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(si:si+Δsi),其中,i=1,2,...,N,i为上升沿序号,si、si+Δsi分别为第i个上升沿在导数序列中的起始位置和终止位置,i、si、Δsi均为整数,该子序列的长度为Δsi+1;
(2.2)确定数字脉冲波形中每个脉冲上升沿的位置。根据步骤(2.1)获得的N个脉冲上升沿区域对应的子序列DeriPulse(sj:sj+Δsj),j=1,2,...,N,获得数字脉冲波形中对应的上升沿区域对应的子序列DigitalPulse(sj:sj+Δsj),对数字脉冲波形中每个子序列DigitalPulse(sj:sj+Δsj)求取其极小值对应的位置minj和极大值对应的位置maxj,两者之间的波形即为一个波形的上升沿,记DigitalPulse(minj:maxj)为该数字脉冲波形中第j个脉冲上升沿,其中,sj≤minj<maxj≤sj+Δsj,minj、maxj均为整数;
(2.3)确定数字脉冲波形中每个脉冲下降沿的位置。步骤(2.2)确定了每个脉冲的上升沿,两个脉冲上升沿之间的部分即为前一个脉冲的下降沿的波形,对于第j个脉冲,1≤j<N,其下降沿为DigitalPulse(maxj:minj+1),对第N个脉冲,从前至后在DigitalPulse(maxN:n)找到第一个数值小于0的采样点对应的位置为结束点ending,若DigitalPulse(maxN:n)中无小于0的采样点,则ending等于n+1,则第N个脉冲的下降沿为DigitalPulse(maxN:ending-1),其中,maxN<ending,ending为整数。
(3)对分割出的各个单次事件脉冲波形依次进行重建。每个单次脉冲波形可以认为是由一个快速上升的上升沿和一个指数衰减的下降沿构成。步骤(2)已经分割出每个单次事件脉冲波形的上升沿和下降沿,需要对每个单次事件脉冲波形的上升沿和下降沿分别进行重建,获得其模型参数。由于当次到达的事件脉冲波形的下降沿会叠加到下一个到达的事件脉冲波形中,在对下一个到达的事件脉冲波形进行重建时,需要去除当次的事件脉冲波形的影响。对每个单次事件脉冲波形依次进行如下具体重建步骤:
(3.1)上升沿重建。脉冲上升沿以线性公式(8)表示,
y(x)=LineK×x+LineB (8)
其中,LineK为上升沿直线的斜率,LineK>0,LineB是上升沿截距,可为任意数值,和上升沿开始时间是线性关系,x为数字脉冲序列采样位置,y为在x位置上的数字脉冲幅值。
在脉冲上升沿抽取不少于P个采样点,其中,P为整数,P≥2且小于等于上升沿长度,对这P个点按式(8)进行线性拟合获取上升沿重建参数LineK和LineB。
拟合时,可以采用最小二乘法进行拟合。(具体参见“梁家惠.用最小二乘法进行直线拟合的讨论.工科物理,pp.11-15,1995年3月”)
(3.2)下降沿重建。脉冲下降沿以指数公式(9)表示,
y(x)=exp(ExpK×x+ExpB) (9)
其中,ExpK为衰减时间常数,ExpK<0,ExpB可为任意数值,与下降沿开始时间是线性关系。
在脉冲下降沿抽取不少于Q个采样点,其中,Q为整数,Q≥2且小于等于下降沿长度,对这Q个点进行按式(9)指数拟合获取下降沿重建参数ExpK和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)进行线性拟合获取上升沿重建参数LineK和LineB,进行重建:
y(x)=LineK×x+LineB (8)
其中,LineK为上升沿直线的斜率,LineK>0,LineB是上升沿截距,可为任意值,x为数字脉冲序列采样位置,y(x)为在x位置上的数字脉冲幅值;
下降沿重建模块(310)中获得下降重建参数的过程为:
在脉冲下降沿抽取不少于Q个采样点,其中,Q为整数,Q≥2且小于等于下降沿长度,对采样点按式(9)进行指数拟合获取下降沿重建参数ExpK和ExpB,进行重建:
y(x)=exp(ExpK×x+ExpB) (9)
其中,ExpK为衰减时间常数,ExpK<0,ExpB可为任意值。
实例:
本发明提出的一种pileup波形处理方法及其***涉及到若干参数,这些参数需针对具体处理数据进行调节以达到良好的性能。此处列出本应用实例处理数据的参数:
步骤(1.1)中导入数字脉冲为使用LYSO晶体和Hamamtsu R9800 PMT采集得到闪烁脉冲,典型波形为图5所示。采样率为10GSps,能量为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,上升沿最小区域长度M=20。
图9为利用本发明提出的方法从pileup波形中恢复各个单次事件脉冲的结果。图10为本发明提出的方法对位置谱的对比结果图,其中,10(a)为未利用本发明获得的位置谱,10(b)为利用本发明获得的位置谱。
本发明的方法和***可应用于高能粒子探测数据获得***中,光子计数(Photon Counting)CT***中,正电子断层扫描(Positron Emission Tomography,PET)***中等。
本发明不仅局限于上述具体实施方式,本领域一般技术人员根据本发明公开的内容,可以采用其它多种具体实施方式实施本发明,因此,凡是采用本发明的设计结构和思路,做一些简单的变化或更改的设计,都落入本发明保护的范围。