CN104101894B - 一种堆积事件侦测与重建的方法及系统 - Google Patents
一种堆积事件侦测与重建的方法及系统 Download PDFInfo
- Publication number
- CN104101894B CN104101894B CN201310115535.9A CN201310115535A CN104101894B CN 104101894 B CN104101894 B CN 104101894B CN 201310115535 A CN201310115535 A CN 201310115535A CN 104101894 B CN104101894 B CN 104101894B
- Authority
- CN
- China
- Prior art keywords
- pulse
- event
- time
- energy
- accumulation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 45
- 238000009825 accumulation Methods 0.000 title claims abstract description 31
- 238000001514 detection method Methods 0.000 title claims abstract description 19
- 238000012545 processing Methods 0.000 claims abstract description 14
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 11
- 230000004044 response Effects 0.000 claims description 15
- 238000001228 spectrum Methods 0.000 claims description 13
- 230000005284 excitation Effects 0.000 claims description 11
- 239000011159 matrix material Substances 0.000 claims description 10
- 230000011218 segmentation Effects 0.000 claims description 10
- 238000012935 Averaging Methods 0.000 claims description 8
- 238000005192 partition Methods 0.000 claims description 7
- 238000012850 discrimination method Methods 0.000 claims description 6
- 239000007787 solid Substances 0.000 claims description 6
- 238000004458 analytical method Methods 0.000 claims description 5
- 230000004397 blinking Effects 0.000 claims description 5
- 230000005855 radiation Effects 0.000 claims description 4
- 239000012634 fragment Substances 0.000 claims description 3
- 238000001208 nuclear magnetic resonance pulse sequence Methods 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 16
- 239000013078 crystal Substances 0.000 description 9
- 238000012549 training Methods 0.000 description 7
- 238000004364 calculation method Methods 0.000 description 5
- 230000008878 coupling Effects 0.000 description 5
- 238000010168 coupling process Methods 0.000 description 5
- 238000005859 coupling reaction Methods 0.000 description 5
- 230000010354 integration Effects 0.000 description 4
- 238000007476 Maximum Likelihood Methods 0.000 description 3
- 238000002591 computed tomography Methods 0.000 description 3
- 230000000630 rising effect Effects 0.000 description 3
- 238000005070 sampling Methods 0.000 description 3
- 239000004809 Teflon Substances 0.000 description 2
- 229920006362 Teflon® Polymers 0.000 description 2
- 239000002390 adhesive tape Substances 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- QWUZMTJBRUASOW-UHFFFAOYSA-N cadmium tellanylidenezinc Chemical compound [Zn].[Cd].[Te] QWUZMTJBRUASOW-UHFFFAOYSA-N 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 230000009977 dual effect Effects 0.000 description 2
- 238000003306 harvesting Methods 0.000 description 2
- 230000014759 maintenance of location Effects 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 238000002600 positron emission tomography Methods 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 238000002603 single-photon emission computed tomography Methods 0.000 description 2
- 238000000638 solvent extraction Methods 0.000 description 2
- 230000015556 catabolic process Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000006866 deterioration Effects 0.000 description 1
- 238000013213 extrapolation Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000008676 import Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 229910052751 metal Inorganic materials 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 238000009206 nuclear medicine Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
Landscapes
- Measurement Of Radiation (AREA)
- Nuclear Medicine (AREA)
Abstract
一种堆积事件侦测与重建的方法,其包括:获取低计数下的非堆积符合单事件的闪烁脉冲数据库,对各路非堆积符合单事件的闪烁脉冲数据库求平均脉冲,闪烁脉冲的形状信息由平均脉冲给定;采用MLEM算法进行退卷积,闪烁脉冲的到达时间定义为退卷积序列的峰值时间;采用MAP准则分割堆积中的单事件,并通过加和,提取事件的能量值;将上述获得的形状、时间、能量信息,存储为列表数据,完成堆积事件的重建。本发明在系统运行时,可有效地检测、分割、复原出堆积中的各个单事件,有效增加了系统计数率,提高了高计数率下的能量分辨率,特别适合于离线环境下的各种堆积事件处理。
Description
技术领域
本发明涉及数字信号处理、光电信号处理和核侦测领域,尤其涉及一种堆积事件侦测与重建的方法及系统。
背景技术
在正电子寿命谱仪、穆斯堡尔谱仪等核分析领域,能谱仪、辐射计数器等核侦测领域,以及计算机断层成像(ComputedTomography,以下简称CT)、正电子发射断层成像(PositronEmissionTomography,以下简称PET)、单光子发射断层成像(SinglePhotoEmissionComputedTomography,以下简称SPECT)等医学影像领域,侦测器部分的工作机理主要分为两种:一种是通过闪烁体将高能光子转化为能量较低的可见光子或紫外光光子,再将可见光光子通过光电器件转化为电信号;另一种是将高能光子通过碲锌镉(以下简称CZT)等半导体材料直接转化为电信号。以上两种工作机理下的侦测器输出均为电信号。
当多个高能光子或粒子在较短时间内击中侦测器,侦测器输出的电脉冲信号将发生堆积,如图3所示。
首先,当两次事件的间隔较小时,堆积会使前一次事件的电荷计入当前事件的电脉冲中,造成能量测量的误差,从而造成能量谱的扭曲和能量分辨率的恶化。随着单位时间击中侦测器的光子数的增加,或电脉冲下降沿时间常数的增加,这种恶化会变得越来越严重。
其次,堆积会使前一个脉冲的尾部叠加在当前事件脉冲的开始部分,基线的抬高将使时间标记变得不准确。一方面,基于前沿甄别器(LeadingEdge Discriminator,以下简称LED)和多电压阈值(MultipleVoltageThreshold,以下简称MVT)方法的有效阈值点将比设置阈值低。另一方面,基于恒定系数鉴别器(constantfractiondiscriminator,以下简称CFD)方法的差分量将包含下降沿的斜率成分,使得实际过零点发生偏移。
再次,在部分采用位置敏感型光电侦测器(PositionSensitive Photo-multiplierTube,以下简称PSPMT)的光子侦测设备中,脉冲信号的位置是通过计算多个角脉冲的能量值的相对关系得到的。因此,堆积带来的能量不确定因素将带来位置信息的不确定。
为了减少堆积对时间信息、能量信息、位置信息的影响,通常的做法是通过电脉冲的积分值或电脉冲的形状,甄别出堆积事件,然后剔除该事件。这样的做法将减少侦测器捕获的有效光子或粒子数。对于PET、SPECT、光子计数CT来说,计数率的减少将使图像的信噪比下降。
此外,一些工作也集中于恢复堆积事件的能量信息,未有对迁移的时间标记点进行重建。例如有人通过计算截断的电脉冲片段的积分值求能量值。截断的电脉冲能量信息的计算决定于积分时间长度,这个长度的起止点准确性决定于时间标记的准确性。未对时间标记点进行重建,影响了计算的能量值的准确性。此外,以上用于堆积恢复的技术均属于模拟技术,对模拟电脉冲进行处理,并非对数字脉冲进行处理。能量计算的方法限定在下降沿为指数形状的脉冲当中。
还有一些数字化方法采用差分+波形拟合的方法对堆积脉冲进行处理。拟合的方法对于实际数字电路系统来说,开销较大。堆积通常发生在计数率较高的情形下,故待处理的脉冲个数较多。若脉冲个数较少时,堆积发生的概率较小。且上述方法的差分操作对光电倍增管中的尖峰噪声十分敏感,造成事件到达的误触发。特别是电脉冲上升沿斜率较小时,对上升沿的捕获将十分不准确。并且,该方法受限于固定的脉冲模型。
因此,针对上述技术问题,有必要提供一种改良结构的堆积事件侦测与重建的方法及系统,以克服上述缺陷。
发明内容
有鉴于此,本发明的目的在于提供一种堆积事件侦测与重建的方法及系统,该方法及系统能有效地实现堆积事件的侦测,实现堆积事件中各个单事件的分割,并对其能量信息进行复原,增加系统的有效计数率。
为实现上述目的,本发明提供如下技术方案:
一种堆积事件侦测与重建的方法,其具体包括步骤:
S1:获取低计数下的非堆积符合单事件的闪烁脉冲数据库,对各路非堆积符合单事件的闪烁脉冲数据库求平均脉冲,闪烁脉冲的形状信息由平均脉冲给定;
S2:由步骤S1给出单位脉冲响应,采用MLEM算法进行退卷积,闪烁脉冲的到达时间定义为退卷积序列的峰值时间;
S3:处理步骤S2输出的时间序列,采用MAP准则分割堆积中的单事件,并通过加和,提取事件的能量值;
S4:将步骤S1、S2、S3给出的形状信息、时间信息、能量信息存储为列表数据,完成堆积事件的重建。。
优选的,在上述堆积事件侦测与重建的方法中,所述步骤S1具体包括:
(1.1)通过降低安置射源的辐射剂量或调整侦测器的立体角,降低每个侦测器捕获的高能光子数,每个侦测器接收到的事件是一个泊松流,其平均计数率为
其中,mi和qi分别是弱源的剂量和弱源对侦测器的立体角,i是弱源的序数,n是弱源的个数;
(1.2)按照恒定系数鉴别方法或前沿甄别方法对齐脉冲;
(1.3)对对齐后的脉冲取平均。。
优选的,在上述堆积事件侦测与重建的方法中,所述步骤S1还包括:将一能量区间选择的符合脉冲对认定为单事件脉冲对,将此数字化的电脉冲预先存储后,进行离线分析,通过对齐脉冲对,获得平均脉冲信号,将平均脉冲作为系统响应,将数字脉冲加和,得到单事件数据的能谱。
优选的,在上述堆积事件侦测与重建的方法中,所述步骤S2具体包括:
(2.1)载入一个闪烁脉冲片段S0,该闪烁脉冲片段的起点是向上过阈值Vl的时间点,终点是向下过阈值Vt的时间点,
其中,Vl是,Vt均为用于设定触发的具体电压数值;
(2.2)迭代法反卷积,建立线性脉冲响应模型:
p=f*m+n, (1)
其中,p是闪烁脉冲片段,f是激励,m是平均脉冲,n是噪声,(1)式改写成线性方程组:p=Hf+n,其中,H是常对角矩阵,H=T{m},m={k-l,k-l+1,...,k-1,k0,k1,...,kr-1,kr},卷积核的长度为l+r+1,其中l是卷积核函数在零时刻左边的非零值长度,r是卷积核函数在零时刻右边的非零值长度,因而,矩阵H的元素又表示为:
矩阵的尺寸由卷积核与载入的闪烁脉冲片段决定,
MLEM迭代表示为:
其中,是第k次迭代后的脉冲输入估计序列,是第k-1次迭代后的脉冲输入估计序列,pi是带噪声的响应序列的第i个值。
优选的,在上述堆积事件侦测与重建的方法中,所述步骤S3具体包括:
(3.1)经过MLEM退卷积后的时间序列,用均值滤波器处理;
(3.2)过阈值Vu处理,Vu为用于设定触发的具体电压数值,把每个过阈值Vu的时间点,当作S1序列中每个单事件的起点,分隔开的事件持续的时间序列为S21,S22,S23,...,S2m,m为大于1的自然数,各自计算能量和脉冲的持续时间,由
决定S2i可能的单事件数,其中E0是载入时间序列S2i的积分值,d0是载入时间序列S2i的长度,1≤i≤m,是S1给出的事件对能量的概率密度函数,而是S1给出的事件对脉冲长度的概率密度函数;
(3.3)若存在S2i的大于1,1≤i≤m,则进入MAP分割步骤,决定了该脉冲最可能的事件数,采用次MAP分割,每次分割对应定义一个最优边界点,由下式决定
其中,il和ir分别是S2i左边和右边的积分值,dl和dr分别是左边和右边的脉冲长度。pduration和pintegral分别是在每次迭代下的对脉冲长度和积分值的概率密度函数。
优选的,在上述堆积事件侦测与重建的方法中,所述步骤S4具体包括:由(4)式定义的每个事件边界内的积分值,作为该事件的能量值,并把积分值作为激励的高度置于每个事件的左边界,把每个事件的激励进行正投影可以得到光滑的脉冲序列,在不需要获取重建脉冲信息的情况下,将前面获取的时间和能量存储为列表数据,即可完成侦测和重建。
一种堆积事件侦测与重建的系统,其包括:
预采集数据训练模块,用于对预先采集好的数据进行训练,获得平均脉冲、能谱和脉冲持续时间谱;
迭代法退卷积模块,用于对输入的时间序列S0进行退卷积操作,生成退卷积后的时间序列S1,时间序列S1具有较短的脉冲持续时间;
分割堆积事件模块,用于分割堆积事件中的各个单事件,并给出各个单事件的能量信息;
重建信息整合模块,用于将得到的形状信息、时间信息、能量信息存储为列表数据,完成堆积事件的重建。
优选的,在上述堆积事件侦测与重建的系统中,所述预采集数据训练模块包括:
平均脉冲模块,用于计算单事件的平均脉冲;
能量计算模块,用于计算单事件的能量值;
脉冲长度模块,用于计算每个脉冲的持续时间。
优选的,在上述堆积事件侦测与重建的系统中,所述迭代法退卷积模块包括:
数字波形输入模块,用于导入数字化的闪烁脉冲片段,数字化的闪烁脉冲片段是单事件或者是多事件,该片段的起点是向上过阈值Vl的时间点,终点是向下过阈值Vt的时间点;
MLEM迭代模块,用于用MLEM算法对序列S0进行退卷积,输出序列S1。
优选的,在上述堆积事件侦测与重建的系统中,所述分割堆积事件模块包括:
粗脉冲分割模块,用于通过过阈值对时间序列S1进行粗分割;
MAP分割模块,用于通过MAP准则对粗分割的结果进行进一步的分割;
能量获取模块,用于计算每个事件的能量值。
从上述技术方案可以看出,本发明实施例在系统运行时,可有效地检测、分割、复原出堆积中的各个单事件;有效增加了系统计数率,提高了高计数率下的能量分辨率,特别适合于离线环境下的各种堆积事件处理。
与现有技术相比,本发明的有益效果是:
(1)较好的能量保持特点。
(2)每次迭代的信号值都为正数,稳定性好。
(3)重建的计数率多。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的有关本发明的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1为本发明堆积事件侦测与重建方法的流程图;
图2为本发明堆积事件侦测与重建系统的系统结构图;
图3为典型的双事件堆叠的示意图;
图4为本发明实施实例中平均脉冲的示意图;
图5为本发明实施实例中40事件堆积的示意图;
图6为本发明实施实例的MLEM退卷积后的时间序列S1;
图7为本发明实施实例的算法还原激励的结果示意图;
图8为本发明实施实例的再次投影的结果示意图;
图9为本发明MAP脉冲分割的示意图;
图10为本发明一种典型系统的示意图;
图11为本发明另一种典型系统的示意图。
具体实施方式
本发明公开了一种堆积事件侦测与重建的方法及系统,该方法及系统能有效地实现堆积事件的侦测,实现堆积事件中各个单事件的分割,并对其能量信息进行复原,增加系统的有效计数率。
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行详细地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例,都属于本发明保护的范围。
最大似然概率期望最大法(MaxLikelihoodExpectationMaximization,全文均简称MLEM)是一种迭代的极大似然估计方法,迭代算法中交替采用期望最大法(ExpectationMaximization,EM),利用对隐藏变量的现有估计值,计算其最大似然估计值。
最大后验概率(MaxAPosterior,全文均简称MAP)算法,是以最大后验概率为准则的优化方法。MLEM和MAP算法经常出现在逆问题领域。
本发明公开的堆积事件侦测与重建的方法首先利用预先获取好的数据库获得平均脉冲。该平均脉冲将作为电信号对时间到达的响应,对一段闪烁脉冲片段进行退卷积。退卷积后的波形具有较短的脉冲持续时间。再利用MAP准则对退卷积后的波形进行分割。在其所属区间内的加和值为该事件的能量值,
如图1所示,本发明公开的堆积事件侦测与重建的方法具体包括步骤:
S1:获取低计数下的非堆积符合单事件的闪烁脉冲数据库,对各路非堆积符合单事件的闪烁脉冲数据库求平均脉冲,闪烁脉冲的形状信息由平均脉冲给定;此步骤中一般要求脉冲数据库中的样本数大于50。脉冲数越多,其统计噪声越小。堆叠事件示意图可如图3及图5所示,图3为典型的双事件堆叠的示意图,图5为本发明实施实例中40事件堆积的示意图。
S2:由步骤S1给出单位脉冲响应,采用MLEM算法进行退卷积,闪烁脉冲的到达时间定义为退卷积序列的峰值时间;
S3:处理步骤S2输出的时间序列,采用MAP准则分割堆积中的单事件,并通过加和,提取事件的能量值;
S4:将步骤S1、S2、S3给出的形状信息、时间信息、能量信息存储为列表数据,完成堆积事件的重建。
其中,所述步骤S1具体包括:
(1.1)通过降低安置射源的辐射剂量或调整侦测器的立体角,降低每个侦测器捕获的高能光子数,每个侦测器接收到的事件是一个泊松流,其平均计数率为
其中,mi和qi分别是弱源的剂量和弱源对侦测器的立体角,i是弱源的序数,n是弱源的个数;在获得脉冲数据库时,可以令n=1。并令m1足够小。
(1.2)按照恒定系数鉴别方法或前沿甄别方法对齐脉冲;通过能量窗选择的符合脉冲对,基本上可以被认为都是单事件脉冲对。将此数字化的电脉冲预先存储后,进行离线分析。
(1.3)对对齐后的脉冲取平均。获得平均脉冲信号和脉冲残差序列,将平均脉冲作为系统响应,具体如图4所示,图4为本发明实施实例中平均脉冲的示意图。
其中,所述步骤S1还包括:将一能量区间选择的符合脉冲对认定为单事件脉冲对,将此数字化的电脉冲预先存储后,进行离线分析,通过对齐脉冲对,获得平均脉冲信号,将平均脉冲作为系统响应,将数字脉冲加和,得到单事件数据的能谱。单事件能谱将作为先验知识输入给步骤S3。如未有特别注明,本发明所述的闪烁脉冲均为正值。如实际获取的脉冲是负脉冲,则默认地经过了反相操作。因此,上升沿即是电脉冲的前沿,下降沿即是电脉冲的尾部。
在样本数不足时,对平均脉冲采用某些衍生方法,使脉冲更接近脉冲形状的期望值。衍生方法可以是滤波方法,或者拟合、外推插值等。处理原始脉冲数据库得到的所有中间信息可称作衍生脉冲数据库,该衍生脉冲数据库可以是在原始脉冲数据库基础上所做的滤波、插值、拟合、外推、组合等等。
其中,所述步骤S2具体包括:
(2.1)载入一个闪烁脉冲片段S0,该闪烁脉冲片段的起点是向上过阈值Vl的时间点,终点是向下过阈值Vt的时间点,
其中,Vl,Vt均为用于设定触发的具体电压数值;
(2.2)迭代法反卷积,建立线性脉冲响应模型:
p=f*m+n, (1)
其中,p是闪烁脉冲片段,f是激励,m是平均脉冲,n是噪声,(1)式改写成线性方程组:p=Hf+n,其中,H是常对角矩阵(又称Toeplitz矩阵),H=T{m},m={k-l,k-l+1,...,k-1,k0,k1,...kr-1,kr},卷积核的长度为l+r+1,其中l是卷积核函数在零时刻左边的非零值长度,r是卷积核函数在零时刻右边的非零值长度,因而,矩阵H的元素又表示为:
矩阵的尺寸由卷积核与载入的闪烁脉冲片段决定,
MLEM迭代表示为:
其中,是第k次迭代后的脉冲输入估计序列,f是第k-1次迭代后的脉冲输入估计序列,pi是带噪声的响应序列的第i个值。
其中,所述步骤S3具体包括:
(3.1)经过MLEM退卷积后的时间序列,用均值滤波器处理;
(3.2)过阈值Vu处理,Vu为用于设定触发的具体电压数值,把每个过阈值Vu的时间点,当作S1序列中每个单事件的起点,分隔开的事件持续的时间序列为S21,S22,S23,...,S2m,m为大于1的自然数,各自计算能量和脉冲的持续时间,由
决定S2i可能的单事件数,其中E0是载入时间序列S2i的积分值,d0是载入时间序列S2i的长度,1≤i≤m,是S1给出的事件对能量的概率密度函数,而是S1给出的事件对脉冲长度的概率密度函数,具体如图6、图7及图8所示,图6为本发明实施实例的MLEM退卷积后的时间序列S1,图7为本发明实施实例的算法还原激励的结果示意图,图8为本发明实施实例的再次投影的结果示意图;
(3.3)若存在S2i的大于1,1≤i≤m,则进入MAP分割步骤,决定了该脉冲最可能的事件数,采用次MAP分割,每次分割对应定义一个最优边界点,由下式决定
其中,il和ir分别是S2i左边和右边的积分值,dl和dr分别是左边和右边的脉冲长度。pduration和pintegral分别是在每次迭代下的对脉冲长度和积分值的概率密度函数,具体如图9所示,图9为本发明MAP脉冲分割的示意图。
其中,所述步骤S4具体包括:由(4)式定义的每个事件边界内的积分值,作为该事件的能量值,并把积分值作为激励的高度置于每个事件的左边界,把每个事件的激励进行正投影可以得到光滑的脉冲序列,在不需要获取重建脉冲信息的情况下,将前面获取的时间和能量存储为列表数据,即可完成侦测和重建。
如图2所示,本发明公开的堆积事件侦测与重建的系统,其包括预采集数据训练模块100、迭代法退卷积模块200、分割堆积事件模块300及重建信息整合模块400。预采集数据训练模块100用于对预先采集好的数据进行训练,得到平均脉冲、能谱和脉冲持续时间谱。该平均脉冲将作为卷积核提供给迭代法退卷积模块200,该能谱和脉冲持续时间谱将作为事件的能量概率密度函数和持续时间概率密度函数提供给分割堆积事件模块300。
具体讲,预采集数据训练模块100用于对预先采集好的数据进行训练,获得平均脉冲、能谱和脉冲持续时间谱;迭代法退卷积模块200用于对输入的时间序列S0进行退卷积操作,生成退卷积后的时间序列S1,时间序列S1具有较短的脉冲持续时间;分割堆积事件模块300用于分割堆积事件中的各个单事件,并给出各个单事件的能量信息;重建信息整合模块400用于将得到的形状信息、时间信息、能量信息存储为列表数据,完成堆积事件的重建。
继续如图2所示,所述预采集数据训练模块100包括三个子模块,该三个子模块分别为平均脉冲模块110、能量计算模块120及脉冲长度模块130。其中,平均脉冲模块110用于计算单事件的平均脉冲;能量计算模块120用于计算单事件的能量值;脉冲长度模块130用于计算每个脉冲的持续时间。
继续如图2所示,所述迭代法退卷积模块包括两个子模块,该两个子模块分别为数字波形输入模块210及MLEM迭代模块220。其中,数字波形输入模块210用于导入数字化的闪烁脉冲片段,数字化的闪烁脉冲片段是单事件或者是多事件,该片段的起点是向上过阈值Vl的时间点,终点是向下过阈值Vt的时间点;MLEM迭代模块220用于用MLEM算法对序列S0进行退卷积,输出序列S1。
继续如图2所示,所述分割堆积事件模块包括三个子模块,该三个子模块分别为粗脉冲分割模块310、MAP分割模块320及能量获取模块330。其中,粗脉冲分割模块310用于通过过阈值对时间序列S1进行粗分割;MAP分割模块320用于通过MAP准则对粗分割的结果进行进一步的分割;能量获取模块330用于计算每个事件的能量值。
如图10及图11所示,图10为本发明一种典型系统的示意图;图11为本发明另一种典型系统的示意图。结合图10及图11,下面将结合几个具体的实施例,对本发明的堆积事件侦测与重建方法及系统进行阐述。
实例一:
本发明提出的一种堆积事件侦测与重建方法及系统。涉及到的参数、滤波器设计需要根据与获取数据的特点进行调节以达到良好的能量分辨性能和较短的脉冲持续时间。此处列出本应用实例处理数据的参数:
步骤(1)所用的实际系统为使用LYSO晶体和Hamamatsu R9800PMT。晶体尺寸为16.5mm×16.5mm×10.0mm。晶体和PMT耦合面为100面,除开耦合面外,其他面均用特氟龙胶带包裹。数据获得系统的采样率为50Ghz,带宽16Ghz。射源为511kev的正电子湮灭伽马光子。平均脉冲上升时间约为2ns,下降沿用指数拟合后,时间常数为42.5497ns。步骤(1)通过处理3000个非堆积的闪烁脉冲,获得闪烁脉冲的响应模型。这里采用的模型是平均脉冲。
步骤(2.1)采用的Vl=0.080V,Vt=0.008V.通过步骤(2)的退卷积,使堆积数据粗分离。
步骤(3.2)采用的Vu=8Vns.进一步分离堆积数据。
实例二:
此处列出本应用实例二处理数据的参数:
步骤(1)所用的实际系统为使用LYSO晶体如图11和FM300035SIPM。晶体尺寸为2.0mm×2.0mm×10.0mm。晶体和PMT耦合面为100面,除开耦合面外,其他面均用特氟龙胶带包裹。数据获得系统的采样率为50Ghz,带宽16Ghz(如图10)。射源为511kev的正电子湮灭伽马光子。平均脉冲上升时间约为5ns。采用50个脉冲进行平均操作,再进行指数拟合,拟合的模型作为响应模型。
步骤(2.1)采用的Vl=0.060V,Vt=0.006V.通过步骤(2)的退卷积,使堆积数据粗分离。
步骤(3.2)采用的Vu=6Vns.进一步分离堆积数据。
实例三:
此处列出本应用实例三处理数据的参数:
步骤(1)所用的实际系统为使用LaBr晶体和HamamatsuR9800PMT。晶体尺寸为3.5mm×3.5mm×5.0mm。晶体和PMT耦合面为100面,除开耦合面为光导外,其他面均用金属密封。数据获得系统的采样率为50Ghz,带宽16Ghz(如图10)。射源为511kev的正电子湮灭伽马光子。平均脉冲上升时间约为2ns。
步骤(2.1)采用的Vl=0.100V,Vt=0.010V.通过步骤(2)的退卷积,使堆积数据粗分离。
步骤(3.2)采用的Vu=10Vns.进一步分离堆积数据。
本发明的方法和系统可以用于高计数率背景下的核侦测、核分析、核医学仪器。
从上述技术方案可以看出,本发明实施例在系统运行时,可有效地检测、分割、复原出堆积中的各个单事件;有效增加了系统计数率,提高了高计数率下的能量分辨率,特别适合于离线环境下的各种堆积事件处理。
与现有技术相比,本发明的有益效果是:
(1)较好的能量保持特点。
(2)每次迭代的信号值都为正数,稳定性好。
(3)重建的计数率多。
对于本领域技术人员而言,显然本发明不限于上述示范性实施例的细节,而且在不背离本发明的精神或基本特征的情况下,能够以其他的具体形式实现本发明。因此,无论从哪一点来看,均应将实施例看作是示范性的,而且是非限制性的,本发明的范围由所附权利要求而不是上述说明限定,因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
此外,应当理解,虽然本说明书按照实施方式加以描述,但并非每个实施方式仅包含一个独立的技术方案,说明书的这种叙述方式仅仅是为清楚起见,本领域技术人员应当将说明书作为一个整体,各实施例中的技术方案也可以经适当组合,形成本领域技术人员可以理解的其他实施方式。
Claims (5)
1.一种堆积事件侦测与重建的方法,其特征在于:具体包括步骤:
S1:获取低计数下的非堆积符合单事件的闪烁脉冲数据库,对各路非堆积符合单事件的闪烁脉冲数据库求平均脉冲,所述求平均脉冲的步骤为:将一能量区间选择的符合脉冲对认定为单事件脉冲对,将此数字化的电脉冲预先存储后,进行离线分析,通过对齐脉冲对,获得平均脉冲信号,将平均脉冲作为系统响应,将数字脉冲加和,得到单事件数据的能谱,闪烁脉冲的形状信息由平均脉冲给定;
S2:由步骤S1给出闪烁脉冲片段,采用MLEM算法进行退卷积,闪烁脉冲的到达时间定义为退卷积序列的峰值时间;
S3:处理步骤S2输出的时间序列,采用MAP准则分割堆积中的单事件,并通过加和,提取事件的能量值;
S4:将步骤S1、S2、S3给出的形状信息、时间信息、能量信息存储为列表数据,完成堆积事件的重建。
2.根据权利要求1所述的堆积事件侦测与重建的方法,其特征在于:所述步骤S1具体包括:
(1.1)通过降低安置射源的辐射剂量或调整侦测器的立体角,降低每个侦测器捕获的高能光子数,每个侦测器接收到的事件是一个泊松流,其平均计数率为
其中,mi和qi分别是弱源的剂量和弱源对侦测器的立体角,i是弱源的序数,n是弱源的个数;
(1.2)按照恒定系数鉴别方法或前沿甄别方法对齐脉冲;
(1.3)对对齐后的脉冲取平均。
3.根据权利要求1所述的堆积事件侦测与重建的方法,其特征在于:所述步骤S2具体包括:
(2.1)载入一个闪烁脉冲片段S0,该闪烁脉冲片段的起点是向上过阈值Vl的时间点,终点是向下过阈值Vt的时间点,
其中,Vl,Vt均为用于设定触发的具体电压数值;
(2.2)迭代法反卷积,建立线性脉冲响应模型:
p=f*m+n, (1)
其中,p是闪烁脉冲片段,f是激励,m是平均脉冲,n是噪声,(1)式改写成线性方程组:p=Hf+n,其中,H是常对角矩阵,H=T{m},m={k-l,k-l+1,…,k-1,k0,k1,…,kr-1,kr},卷积核的长度为l+r+1,其中l是卷积核函数在零时刻左边的非零值长度,r是卷积核函数在零时刻右边的非零值长度,因而,矩阵H的元素又表示为:
矩阵的尺寸由卷积核与载入的闪烁脉冲片段决定,
MLEM迭代表示为:
计序列,pi是带噪声的响应序列的第i个值。
4.根据权利要求1所述的堆积事件侦测与重建的方法,其特征在于:所述步骤S3具体包括:
(3.1)经过MLEM退卷积后的时间序列,用均值滤波器处理;
(3.2)过阈值Vu处理,Vu为用于设定触发的具体电压数值,把每个过阈值Vu的时间点,当作S1序列中每个单事件的起点,分隔开的事件持续的时间序列为S21,S22,S23,...,S2m,m为大于1的自然数,各自计算能量和脉冲的持续时间,由
决定S2i可能的单事件数,其中E0是载入时间序列S2i的积分值,d0是载入时间序列S2i的长度,1≤i≤m,是S1给出的事件对能量的概率密度函数,而是S1给出的事件对脉冲长度的概率密度函数;
(3.3)若存在S2i的大于1,1≤i≤m,则进入MAP分割步骤,决定了该脉冲最可能的事件数,采用次MAP分割,每次分割对应定义一个最优边界点 由下式决定
pduration和pintegral分别是在每次迭代下的对脉冲长度和积分值的概率密度函数。
5.根据权利要求4所述的堆积事件侦测与重建的方法,其特征在于:所述步骤S4具体包括:由(4)式定义的每个事件边界内的积分值,作为该事件的能量值,并把积分值作为激励的高度置于每个事件的左边界,把每个事件的激励进行正投影可以得到光滑的脉冲序列,在不需要获取重建脉冲信息的情况下,将前面获取的时间和能量存储为列表数据,即可完成侦测和重建。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310115535.9A CN104101894B (zh) | 2013-04-07 | 2013-04-07 | 一种堆积事件侦测与重建的方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310115535.9A CN104101894B (zh) | 2013-04-07 | 2013-04-07 | 一种堆积事件侦测与重建的方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104101894A CN104101894A (zh) | 2014-10-15 |
CN104101894B true CN104101894B (zh) | 2017-04-12 |
Family
ID=51670188
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201310115535.9A Active CN104101894B (zh) | 2013-04-07 | 2013-04-07 | 一种堆积事件侦测与重建的方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104101894B (zh) |
Families Citing this family (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP3350624B1 (en) * | 2015-09-18 | 2019-11-06 | Koninklijke Philips N.V. | Correcting photon counts in a photon counting x-ray radiation detection system |
CN105212954B (zh) * | 2015-11-05 | 2018-03-16 | 苏州瑞派宁科技有限公司 | 一种脉冲堆积事件实时处理方法与系统 |
CN105487104A (zh) * | 2015-12-31 | 2016-04-13 | 中国科学院青海盐湖研究所 | 基于原料钾矿碘化钠晶体测钾仪的测钾方法 |
CN105842544B (zh) * | 2016-03-18 | 2018-09-18 | 南京瑞派宁信息科技有限公司 | 一种迭代的闪烁脉冲时间标记及其交叉验证方法 |
CN106344060B (zh) * | 2016-09-30 | 2019-06-28 | 上海联影医疗科技有限公司 | 用于pet探测器的死时间校正方法以及死时间检测方法 |
CN109669206A (zh) * | 2019-03-03 | 2019-04-23 | 南昌华亮光电有限责任公司 | 环绕式液体闪烁体智能能谱仪系统及其信号处理方法 |
CN111245409B (zh) * | 2019-12-30 | 2021-12-24 | 中国科学院高能物理研究所 | 脉冲信号处理方法及装置 |
CN112006713B (zh) * | 2020-08-28 | 2023-07-21 | 上海联影医疗科技股份有限公司 | 动态重建、扫描对象定位、计数率显示方法和pet扫描设备 |
CN112587161B (zh) * | 2020-12-09 | 2022-09-30 | 明峰医疗系统股份有限公司 | 堆积信号恢复方法、pet成像设备及计算机可读存储介质 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101124725A (zh) * | 2004-09-16 | 2008-02-13 | 南方创新国际私人有限公司 | 用于从检测器输出数据中分解单个信号的方法和设备 |
CN102073059A (zh) * | 2010-12-31 | 2011-05-25 | 华中科技大学 | 一种数字化pileup波形处理方法及系统 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US8378310B2 (en) * | 2009-02-11 | 2013-02-19 | Prismatic Sensors Ab | Image quality in photon counting-mode detector systems |
-
2013
- 2013-04-07 CN CN201310115535.9A patent/CN104101894B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101124725A (zh) * | 2004-09-16 | 2008-02-13 | 南方创新国际私人有限公司 | 用于从检测器输出数据中分解单个信号的方法和设备 |
CN102073059A (zh) * | 2010-12-31 | 2011-05-25 | 华中科技大学 | 一种数字化pileup波形处理方法及系统 |
Non-Patent Citations (4)
Title |
---|
a digital method for separation and reconstruction of pile-up events in germanium detectors;M.Nakhostin等;《Review of Scientific Instruments》;20101025;第81卷(第10期);第5页倒数第1-6行、第6-7页,图1-7 * |
A study ofthe real-time deconvolution ofdigitized waveforms with pulse pile up for digital radiation spectroscopy;Weijun Guo等;《Nuclear Instruments and Methods in Physics Research A》;20050305;第544卷;全文 * |
PILE-UP FREE PARAMETER ESTIMATION AND DIGITAL ONLINE PEAK LOCALIZATION ALGORITHMS FOR GAMMA RAY SPECTROSCOPY;JM Noras等;《IEEETEM》;20031231;全文 * |
Pulse pile-up recovery for the front-end electronics of the PANDA Electromagnetic Calorimeter;G. Tambave等;《IOP PUBLISHING FOR SISSA MEDIALAB》;20121105;全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN104101894A (zh) | 2014-10-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104101894B (zh) | 一种堆积事件侦测与重建的方法及系统 | |
US10078009B2 (en) | Pulse processing circuit with correction means | |
JP6875277B2 (ja) | 物質特性評価のための方法 | |
JP6657398B2 (ja) | パルスパイルアップイベントのリアルタイム処理の方法及びシステム | |
AU2009230876B2 (en) | Radiation imaging method with individual signal resolution | |
EP3593171B1 (en) | Increased spatial resolution for photon-counting edge-on x-ray detectors | |
EP2748594B1 (en) | Radiographic apparatus for detecting photons with offset correction | |
CN101680954A (zh) | 谱光子计数探测器 | |
JP7317586B2 (ja) | 医用画像処理装置、方法及びプログラム | |
JP6381644B2 (ja) | タイムラベル組合せ方法及びシステム | |
CN110462442A (zh) | 实现重合的光子计数检测器 | |
JP6932250B2 (ja) | 陽電子放出断層撮影における狭いエネルギー窓カウントから再構成された放出画像推定を使用した散乱補正 | |
Luo et al. | Pulse pile-up identification and reconstruction for liquid scintillator based neutron detectors | |
CN109477903A (zh) | 光谱辐射探测器中的改进的光子计数 | |
Mohammadian-Behbahani et al. | Pile-up correction algorithm based on successive integration for high count rate medical imaging and radiation spectroscopy | |
JP6675127B2 (ja) | LaBr3シンチレーション検出器及び特定イベント排除方法 | |
CN107767427B (zh) | 一种信号波形恢复方法及装置 | |
EP4200649A1 (en) | Methods and systems for coincidence detection in x-ray detectors | |
Jordanov | Pile-up rejection using pulse-shape discrimination | |
Mohammadian-Behbahani | Non-iterative pulse tail extrapolation algorithms for correcting nuclear pulse pile-up | |
Ruiz-Gonzalez et al. | Fisher information analysis of digital pulse timing | |
US8330114B2 (en) | Systems for increasing the sensitivity of gamma-ray imagers | |
Abbaszadeh et al. | Effect of energy threshold in positioning true coincidences that undergo detector scatter for a sub-mm resolution CZT-based PET system | |
Wagner | Event classification based on spectral analysis of scintillation waveforms | |
Kao et al. | Event-time determination by waveform analysis for time-of-flight positron emission tomography |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CB03 | Change of inventor or designer information | ||
CB03 | Change of inventor or designer information |
Inventor after: Xie Qingguo Inventor after: Deng Zhenzhou Inventor after: Wang Weidong Inventor after: Xiao Peng Inventor before: Xie Qingguo Inventor before: Deng Zhenzhou Inventor before: Xiao Peng |