本发明涉及医疗器械的信号处理领域，更具体地涉及数字pet领域的一种闪烁脉冲数字化信号的拟合方法。
背景技术：
正电子发射断层成像(positronemissiontomography，全文均简称pet)通过捕捉人体内因正电子湮灭而发出的γ光子获取以正电子核素为标记的示踪剂在人体内的分布情况，进而获取脏器功能，代谢等病理生理特征。获取γ光子所携带的能量、位置以及时间信息的准确性直接影响到系统成像的性能。闪烁探测器用于高能粒子(例如γ、x等粒子)的捕获及相应的能量沉积信息测量，由于其具有探测效率高、时间响应快和信息测量准确等特点，广泛应用于辐射探测与成像领域。闪烁探测器通常由闪烁体(例如bgo、lyso、labr等)和光电转换器(例如pmt、sipm、apd等)构成，入射的γ、x等粒子与闪烁体相互作用产生荧光，光电转换器件则将荧光转换为相应的闪烁脉冲。将闪烁脉冲数字化，利用数字信号处理技术分析数字闪烁脉冲以提取粒子的能量沉积信息，通过进一步的分析数字闪烁脉冲、优化和升级处理算法即可不断提高粒子能量沉积信息提取的精度，提高闪烁探测器的性能，进而提高pet系统的成像质量。
数字闪烁脉冲采样的一种方法是多阈值采样(multi-voltagethreshold，以下简称mvt)方法，该方法根据闪烁脉冲的特征合理的设置多个电压阈值，通过对闪烁脉冲越过电压阈值的时间的数字化最终实现闪烁脉冲的数字化采样。利用少量的阈值比较器配合时间数字转换器(time-to-digitalconverter，简称tdc)便可设计出相应的mvt采样电路进行闪烁脉冲的mvt采样，完成相应的数字闪烁脉冲的获取。mvt方法用于闪烁脉冲的数字化和粒子能量沉积信息的提取，具有性能优异、工程实现简单、硬件成本低廉等优点。现有研究结果表明，采用3～8个阈值的mvt采样电路获取的数字闪烁脉冲，配合基于先验的闪烁脉冲形状信息的数字脉冲分析方法即可实现粒子能量沉积信息的精确获取。
现有的数字脉冲分析方法之一为最小二乘算法，其包括牛顿法、梯度法等。然而，现有的大部分最小二乘算法在对pet系统获取到的数字闪烁脉冲信号进行拟合时耗时极长，加上大数量级的闪烁脉冲数目，导致实际处理过程中效率低下，从采集数据到得到医学图像的耗时过长，已经逐渐无法满足工业需求，而且这些算法具有一定的使用限定，无法在稳定性、精确性上达到令人满意的平衡。另外有一些拟合算法是运用最小二乘的思想，通过选定双指数模型中的两个参数来进行拟合，以牺牲脉冲完整性的代价来实现一定的时间效率，但是，所获取的闪烁脉冲时间、能量信息并不十分精确，最终在pet系统的时间分辨率和能量分辨率上会有一定恶化。
列文伯格-马夸尔特(levenberg-marquardt)法是最优化算法中的一种，也是使用最广泛的非线性最小二乘迭代算法，它是利用梯度求最大(小)值，介于牛顿法与梯度下降法之间的一种非线性优化方法，并且同时具有梯度法和牛顿法的优点，因此，它的应用领域非常广泛，如：经济学、管理优化、网络分析、最优设计、机械或电子设计等等。因此，为了进一步提高脉冲重建后的图像质量，如果能够将列文伯格-马夸尔特算法应用于经mvt采样获得的闪烁脉冲数字化信号的处理，将是兼顾效率、准确度、稳定性等综合因素的最佳拟合方法，能够更好的满足工业需求。
技术实现要素：
本发明的目的是提供一种闪烁脉冲数字化信号的拟合方法，从而解决现有技术的数字脉冲拟合方法不能兼顾效率、准确度和稳定性的问题。
为了解决上述技术问题，本发明的技术方案是提供一种闪烁脉冲数字化信号的拟合方法，该拟合方法包括以下步骤：
步骤s1：选定闪烁脉冲的先验模型为双指数模型v(t)＝a·ebt+c·edt，将所述双指数模型确定为目标函数；
步骤s2：选定闪烁脉冲数据库中的经典脉冲的采样点，所述采样点包括八个电压阈值序列和八个时间点序列，其中，八个所述电压阈值序列为v1、v2、v3、v4、v5、v6、v7和v8，八个所述时间点序列为t1、t2、t3、t4、t5、t6、t7和t8；
步骤s3：将所述电压阈值序列存为电压数组v＝[v1，v2，v3，v4，v5，v6，v7，v8]，对所述双指数模型的曲线进行时间上的平移，所述平移的具体步骤为：将第一个时间点t1置0，后续时间点分别与所述第一个时间点t1求差，然后将得到的数据存为时间数组t＝[0，(t2-t1)，(t3-t1)，(t4-t1)，(t5-t1)，(t6-t1)，(t7-t1)，(t8-t1)]；
步骤s4：将所述目标函数的四个初始默认参数a、b、c、d设置为1，存为初始参数数组para＝[1，1，1，1]；
步骤s5：将所述电压数组v、时间数组t、所述脉冲采样点个数n、所述初始参数数组para、所述目标函数v(t)输入到列文伯格-马夸尔特拟合函数中，进行拟合得到四个拟合参数a’、b’、c’、d’，将四个所述拟合参数存为第二参数数组paranew＝[a’，b’，c’，d’]；
步骤s6：对于闪烁脉冲数据库中的每个闪烁脉冲，除了步骤s4中用所述第二参数数组paranew代替所述初始参数数组para外，依次重复所述步骤s1-步骤s5，直至完成所有闪烁脉冲的拟合，记第k个闪烁脉冲拟合后的函数为
步骤s7：对拟合后确定拟合参数的闪烁脉冲的目标函数vk(t)积分得到第k个脉冲的能量值ek；
步骤s8：解方程vk(t)＝0得到曲线过零点tk0，得到第k个脉冲的到达时间tk＝t1+tk0，脉冲的时间信息、能量信息提取完毕。
闪烁脉冲的先验模型通过数字示波器采集得到，闪烁脉冲的曲线包括快速的上升沿和缓慢的下降沿。
闪烁脉冲数据库通过闪烁脉冲获取平台记录数字闪烁脉冲形成。
闪烁脉冲获取平台包括一对相对设置的闪烁探测器以及一个射源，射源置于一对所述闪烁探测器轴向连线的中间位置，所述闪烁探测器分别通过同轴电缆接入所述数字示波器。
闪烁探测器为si-bdm探测器，所述射源为点状18fdg射源。
数字示波器以16ghz的模拟带宽和50gsps的采样率进行数字化采样，记录采集到的闪烁脉冲。
列文伯格-马夸尔特拟合函数的原型为voidlmcurve(intn_par，double*par，intm_dat，constdouble*t，constdouble*y，double(*f)(doublet，constdouble*par)，constlm_control_struct*control，lm_status_struct*status)，其中，intn_par表示待拟合模型函数的参数个数；double*par表示参数初始值；intm_dat表示拟合采样点的个数；constdouble*t表示采样点的横坐标数组；constdouble*y表示采样点的纵坐标数组；double(*f)(doublet，constdouble*par)表示拟合依据的曲线模型；constlm_control_struct*control表示拟合算法的控制参数；lm_status_struct*status表示拟合算法的状态参数。
数据类型为包括double和float，对于性能好、配置高的服务器，数据类型设置为double，否则数据类型设置为float。
目标函数vk(t)的积分范围为0-200。
曲线过零点tk0＝[1/(dk-bk)*ln(-ak/ck)]。
本发明的闪烁脉冲数字化信号的拟合方法，利用列文伯格-马夸尔特法对基于mvt方法获取的数字闪烁脉冲信号进行非线性拟合，选用双指数模型，通过从已有的脉冲数据库中提取的双指数模型参数作为算法的初值，能大大提高算法运行时间，同时算法本身易于openmp/gpu进行加速，随着计算机的更新换代算法本身的时间性能也会逐渐提高。通过本方法对脉冲进行符合并重建后，得到的能量分辨率有显著提升，本方法不仅能够提高处理效率，而且通过本方法还能保证信息提取的准确度和稳定性，因此，本发明是在效率、准确、稳定等综合因素下的最佳拟合方法。
附图说明
图1是根据本发明的闪烁脉冲数字化信号的拟合方法的脉冲的先验模型示意图；
图2是根据本发明一个实施例的闪烁脉冲获取平台的布置示意图；
图3是根据本发明一个实施例的经典闪烁脉冲的示意图；
图4是根据本发明一个实施例的采用双指数模型的闪烁脉冲数字化信号的拟合方法应用于pet信号处理的效果图；
图5是现有技术的采用直线指数模型所获得的能量分辨率的效果图。
具体实施方式
以下结合具体实施例，对本发明做进一步说明。应理解，以下实施例仅用于说明本发明而非用于限制本发明的范围。
本发明提供的闪烁脉冲数字化信号的拟合方法，包括以下步骤：
步骤s1：选定闪烁脉冲的先验模型为双指数模型v(t)＝a·ebt+c·edt，将该双指数模型确定为目标函数；
步骤s2：选定闪烁脉冲数据库中的经典脉冲的采样点，该采样点包括八个电压阈值序列和八个时间点序列，其中，八个电压阈值序列为v1、v2、v3、v4、v5、v6、v7和v8，八个时间点序列为t1、t2、t3、t4、t5、t6、t7和t8；
步骤s3：将上述电压阈值序列存为电压数组v＝[v1，v2，v3，v4，v5，v6，v7，v8]，对双指数模型的曲线进行时间上的平移，平移的具体方法为：将第一个时间点t1置0，后续时间点分别与t1求差，将得到的数据存为时间数组t＝[0，(t2-t1)，(t3-t1)，(t4-t1)，(t5-t1)，(t6-t1)，(t7-t1)，(t8-t1)]；
步骤s4：选定目标函数的四个初始默认参数a、b、c、d为1，存为初始参数数组para＝[1，1，1，1]；
步骤s5：将电压数组v、时间数组t、脉冲采样点个数n(n＝8)、初始参数数组para、目标函数v(t)输入到拟合函数lmcurve中，根据需要设置相应的数据类型、拟合速度等参数，进行拟合，得到的输出结果为四个新的拟合参数a’、b’、c’、d’，将这四个新的拟合参数存为第二参数数组paranew＝[a’，b’，c’，d’]；
步骤s6：针对实际应用中得到的多个闪烁脉冲数据的拟合过程，除了步骤s4中用第二参数数组paranew来代替初始参数数组para外，其余的闪烁脉冲处理依次重复步骤s1-步骤s5，直至完成所有闪烁脉冲的拟合，记第k个闪烁脉冲拟合后的函数为
步骤s7：对拟合后确定拟合参数的闪烁脉冲的目标函数vk(t)积分得到第k个脉冲的能量值ek；
步骤s8：解方程vk(t)＝0得到曲线过零点tk0(为负数)，得到第k个脉冲的到达时间tk＝t1+tk0，脉冲的时间信息、能量信息提取完毕。
上述步骤s1中，脉冲的先验模型通过数字示波器采集得到，如图2所示，该先验模型的闪烁脉冲曲线包括快速的上升沿和缓慢的下降沿，该闪烁脉冲的形状信息可通过直线指数和双指数模型描述。本发明中采用了双指数模型进行描述，由于闪烁脉冲的下降沿非常快，双指数模型能有效地提高时间分辨率计算的准确性。
上述步骤s2中，闪烁脉冲数据库通过闪烁脉冲获取平台记录数字闪烁脉冲形成。图2为根据本发明一个实施例的闪烁脉冲获取平台1的布置示意图，该闪烁脉冲获取平台1包括一对相对设置的闪烁探测器10以及一个射源20，在图2的实施例中，所采用的闪烁探测器10为si-bdm探测器，射源20为点状18fdg射源，该点状18fdg射源20置于一对闪烁探测器10轴向连线的中间位置；两个闪烁探测器10分别通过同轴电缆40接入数字示波器30，比如，同轴电缆40分别接入数字示波器的通道二和通道三；数字示波器30以16ghz的模拟带宽和50gsps的采样率进行数字化采样，记录通过两个通道采集到的闪烁脉冲，从而形成闪烁脉冲数据库。
根据本发明的一个实施例，上述步骤s2中采用的经典脉冲及采样点的设置如图3所示，其中，八个电压阈值序列为v1、v2、v3、v4、v5、v6、v7和v8，八个时间点序列为t1、t2、t3、t4、t5、t6、t7和t8，图3中电压阈值的具体值的设定可以根据实际需要进行调整，在此不在赘述。
上述步骤s3中，数据的存储可通过适当的操作平台和编程语言实现，比如c++语言，在此不再赘述。
上述步骤s5中，选用的列文伯格-马夸尔特(levenberg-marquardt)拟合算法的函数原型为voidlmcurve(intn_par，double*par，intm_dat，constdouble*t，constdouble*y，double(*f)(doublet，constdouble*par)，constlm_control_struct*control，lm_status_struct*status)，该函数原型的参数中，intn_par表示待拟合模型函数的参数个数；double*par表示参数初始值；intm_dat表示拟合采样点的个数；constdouble*t表示采样点的横坐标数组；constdouble*y表示采样点的纵坐标数组；double(*f)(doublet，constdouble*par)表示拟合依据的曲线模型；constlm_control_struct*control表示拟合算法的控制参数，比如拟合速度、容错等参数，一般取默认值；lm_status_struct*status表示拟合算法的状态参数。
上述步骤s5中，根据需要设置相应的数据类型、拟合速度等参数是指根据所需拟合的闪烁脉冲的数目大小以及系统服务器的内存配置，对于性能好、配置高的服务器，可设置数据类型为double，否则设置数据类型为float；同时，该拟合算法自带的控制参数lm_control_struct*control中包含有拟合速度、步阶、拟合相对误差等参数，这些参数都可以根据当前系统服务器、硬件配置等进行调试，以实现提高效率和准确性的目的，例如，对于本发明的一个实施例，可设置数据类型为double，其余的拟合算法的控制参数设置为默认值。
上述步骤s7中，对vk(t)积分的范围0-200。
上述步骤s8中，解方程是令
图4为根据本方法对一个闪烁脉冲进行数据处理后得到的能量分辨率图谱，图5为利用直线指数模型获得的闪烁脉冲的能量分辨率的效果图，对比图4和图5可以看出，本发明采用双指数模型得到的闪烁脉冲的能量分辨率为17.9％，相比于采用直线指数模型的能量分辨率23.4％，显然本发明获得的能量分辨率具有显著提升，这对数字pet系统而言是明显的进步。
本发明基于mvt方法得到闪烁脉冲数据库，根据预设的四个电压阈值，在脉冲上升沿和下降沿都得到四个采样点，即一个脉冲用八个采样点(时间，电压序列)表示，通过将已有的脉冲数据库中提取的双指数模型参数作为拟合算法的初值，能大大提高算法运行时间，该方法本身也易于openmp/gpu进行加速，该方法做到了稳定性和精确性的平衡，在完整还原双指数模型的闪烁脉冲情况下提高了效率和算法稳定性，根据此方法最后得到的pet系统能量分辨率有显著提高。
以上所述的，仅为本发明的较佳实施例，并非用以限定本发明的范围，本发明的上述实施例还可以做出各种变化。即凡是依据本发明申请的权利要求书及说明书内容所作的简单、等效变化与修饰，皆落入本发明专利的权利要求保护范围。本发明未详尽描述的均为常规技术内容。