获取PET晶体能量峰值及设定能量鉴频器的方法

摘要

本发明涉及获取PET晶体能量峰值的方法及设定能量鉴频器的方法，所述获取PET晶体能量峰值的方法包括：建立PET晶体的能量分布模型并采集PET原始数据；基于所述PET原始数据统计PET晶体能量谱，以获取所述PET晶体的能量分布模型的初次拟合参数；基于所述能量分布模型和初次拟合参数，采用LM迭代算法计算所述能量分布模型的拟合参数，以获取所述PET晶体的能量分布；根据所获取的PET晶体的能量分布完成PET晶体能量峰值的识别。本发明技术方案能够寻求到PET晶体的最佳能量峰值。

说明书

技术领域

本发明涉及PET成像技术领域，特别涉及一种获取PET晶体能量峰值的方法及一种设定能量鉴频器的方法。

背景技术

PET(Positron Emission Computed Tomography，正电子发射计算机断层扫描)系统是当今核医学领域发展最快一项成像技术，在临床检测中得到广泛的应用。

在PET系统中，闪烁晶体探测器广泛用于检测gamma射线，射线通过康普顿散射或光电吸收在闪烁晶体材料中激发而发出可见光，可见光在光电倍增管(PMT)的阴极产生电流，通过光电倍增管的光电放大输出。

光电倍增管产生的脉冲信号同时被一个微分鉴别器接收处理，微分鉴别器包括低能量鉴别器(LLD)及高能量鉴别器(ULD)，结合上述两种鉴别器可用来鉴别过低和过高的脉冲信号。

低能量鉴别器需要设定其LLD阈值来区分散射光子和非散射光子，如果设置过高的LLD阈值会将一些真符合事件排斥在外，设置过低的LLD阈值，则会将过多的散射事件混入检测影响系统准确性。因此，设定恰当的低能量鉴别器LLD阈值对于PET系统的性能来说非常重要。低能量鉴别器LLD阈值与PET系统中闪烁晶体光子能量分布情况有关，尤其是晶体能量峰值。

根据已知研究，每个闪烁晶体在发生康普顿散射或光电吸收效应时，其上光子能量分布呈现一定规律：非散射光子能量集中在峰值511kev，而能量420kev以下主要是散射光子。但实际PET系统的运行中，非散射光子能量谱有一定宽度，导致PET系统的对应闪烁晶体发出光子能量的数据处理通道数目众多，一个数据处理通道对应一种光子能量，比如，在一种通用的PET系统中其闪烁晶体发出光子能量分布范围为0kev～511kev，其对应具有512个数据处理通道。PET系统的多个数据处理通道相对于闪烁晶体发出的光子，对通过相应数据处理通道的光子进行计数，计数最大的数据处理通道编号，也就是该闪烁晶体的能量峰值。

基于上述内容可知，在PET系统中，对低能量鉴别器设定恰当的LLD阈值，意味着需要基于PET晶体的光子能量寻求最佳能量峰值，现有技术无法提供一种能够获取晶体最佳能量峰值的方法。

发明内容

本发明技术方案解决的技术问题为，如何寻求PET晶体的最佳能量峰值。

为了解决上述技术问题，本发明技术方案提供了一种获取PET晶体能量峰值的方法，包括：

建立PET晶体的能量分布模型并采集PET原始数据；

基于所述PET原始数据统计PET晶体能量谱，以获取所述PET晶体的能量分布模型的初次拟合参数；

基于所述能量分布模型和初次拟合参数，采用LM迭代算法计算所述能量分布模型的拟合参数，以获取所述PET晶体的能量分布；

根据所获取的PET晶体的能量分布完成PET晶体能量峰值的识别。

可选的，所述PET晶体的能量分布模型为混合型高斯曲线模型。

可选的，基于如下算式建立所述能量分布模型：

$>f\left(x\right)={\mathrm{Ae}}^{\frac{-{(x-\mu )}^2}{2{\sigma }^2}}+{\mathrm{ax}}^3+{\mathrm{\beta x}}^2+\gamma ;$>

其中，x为PET晶体的能量值，f(x)为PET晶体的能量值为x的事件数，μ为PET晶体能量分布的峰值，σ为PET晶体能量分布的标准差，A、α、β、γ为所述能量分布模型的系数。

可选的，所述获取所述能量分布模型的初次拟合参数包括：根据PET晶体能量谱的统计结果获得PET晶体能量分布的峰值。

可选的，所述采用LM迭代算法计算所述能量分布模型的拟合参数包括：达到预先设定的迭代次数或算法收敛时迭代中止，输出所述PET晶体的能量分布曲线。

可选的，所述达到预先设定的算法收敛时迭代中止包括：达到预先设定的拟合优度时迭代中止。

可选的，所述根据所获取的PET晶体的能量分布完成PET晶体能量峰值的识别包括：

根据所述PET晶体的能量分布获取所述PET晶体的能量峰值；

对所述PET晶体的能量峰值作伪彩化处理，以获取所述PET晶体的能量峰值的伪彩化图像；

基于所述PET晶体的能量峰值的伪彩化图像，调整所述PET晶体的能量峰值。

可选的，所述对所述PET晶体的能量峰值作伪彩化处理，以获取所述PET晶体的能量峰值的伪彩化图像包括：

对PET晶体进行分组处理；

基于所拟合的PET晶体的能量分布统计每组晶体的能量分布置信区间；

将所述置信区间的边界值设定为伪彩化处理的最小阈值及最大阈值；

根据所述伪彩化处理的最小阈值及最大阈值对处于所述最小阈值及最大阈值的能量值进行上色，以获得PET晶体的能量峰值的伪彩化图像。

可选的，所述每组PET晶体的能量分布置信区间为[(μ-3σ)，(μ+3σ)]，μ为每组PET晶体能量峰值，σ为每组PET晶体能量分布的标准差。

为了解决上述技术问题，本发明技术方案还提供了一种设定能量鉴频器的方法，包括：

基于如上所述的获取PET晶体能量峰值的方法识别PET晶体能量峰值；

根据所述PET晶体能量峰值设定该PET晶体所对应能量鉴频器LLD阈值。

本发明技术方案的有益效果至少包括：

本发明技术方案通过对PET系统采集到的原始数据进行统计，预先估计晶体能量峰值，以获取基于LM算法拟合晶体能量分布模型的初始值(初次拟合参数)，使LM算法的初始值贴近真值，从而提高拟合结果的真实性及准确性，有助于寻求晶体最佳能量分布和能量峰值。

本发明技术方案使用LM迭代算法计算晶体能量分布模型的拟合参数，LM算法具有高斯-牛顿法的局部收敛性和梯度下降法的全局特性，通过自适应调整阻尼因子可达到收敛特性，具有更高的迭代收敛速度；LM迭代算法在具有较为真实的初始值的前提下，其收敛速度很快且可得到全局最优解，因此可大大加快获取系统晶体能量分布的速度，并获得最佳晶体能量分布和能量峰值。

在可选方案中，本发明技术方案还对LM算法拟合时的迭代过程进行优化控制：在达到预先设定的拟合优度时便中止迭代，并输出晶体能量分布，采用LM算法拟合优度控制的LM算法迭代过程可以避免拟合时间过长，通过及时调整阻尼因子，使公式的收敛情况符合拟合优度条件，从而获取快速准确的收敛结果，获取最佳晶体能量分布和能量峰值。

在可选方案中，本发明技术方案还对拟合计算所得晶体能量峰值数据进行伪彩化处理，通过直观的颜色分布趋势可以进一步验证和手动校正对个别晶体计算有误的能量峰值，从而实现了自动和手动的结合，使得到的最终拟合结果无限逼近真值。

附图说明

图1为本发明技术方案提供的一种获取PET晶体能量峰值的方法的流程示意图；

图2为本发明技术方案提供的一种适于求取PET晶体能量分布的LM算法流程示意图；

图3为一组晶体能量峰值的伪彩化图像示意图；

图4为图3中一个晶体坐标所对应晶体的拟合结果示意图；

图5为本发明技术方案提供的一种设定PET系统能量鉴频器的方法的流程示意图。

具体实施方式

本发明技术方案不同于现有技术，采用LM(Levenberg-Marquardt)算法对晶体能量分布曲线进行拟合，以获取PET晶体能量分布和能量峰值，在拟合过程中，本发明技术方案首先对晶体进行能量谱统计，基于晶体能量谱估计晶体能量分布模型的初始值(初次拟合参数)，并结合拟合优度对拟合过程进行评估，当一次拟合结果满足预设拟合优度条件即可中止迭代计算，优化迭代过程。

基于上述内容，本发明技术方案可基于如下三种手段进行实现：

选用LM算法对晶体能量分布模型进行拟合，可利用梯度求最大(小)值的算法，并结合最速下降法和高斯-牛顿法，具有高斯-牛顿法的局部收敛性和梯度下降法的全局特性，其收敛性较佳且具有较高的迭代收敛速度。

对晶体进行能量谱统计，基于晶体能量谱估计晶体能量分布模型的初始值(初次拟合参数)，能够使选择的初始值接近真实值，使LM算法的收敛速度很快且可得到全局最优解，从而避免初值远离真值导致算法过早陷入局部最优解而使的最终拟合结果完全背离真值。

使用拟合优度方法对拟合过程进行评估能够避免因从系统中采集数据本身有问题而导致拟合时间长且结果背离真值。

本发明技术方案可以是上述三种手段的结合，或利用任意一种或两种手段进行实现。

为了使本发明的目的、特征和效果能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细说明。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其他不同于在此描述的方式来实施，因此本发明不受下面公开的具体实施例的限制。

如图1所示的一种获取PET晶体能量峰值的方法，包括：

步骤S100，建立PET晶体的能量分布模型并采集PET原始数据。

在PET放射的能量过程中，PET系统的每个闪烁晶体发生康普顿散射或光电吸收效应，统计PET晶体能量分布，也就是PET系统的数据处理通道对该晶体通过相应数据处理通道的光子进行计数。在统计晶体能量分布后，可知PET晶体能量分布范围为0～511Kev，对于一个晶体而言，每个能量上的计数(Counts)之间的分布符合混合型高斯曲线模型。

在上述步骤中，建立PET晶体的能量分布模型是选择适于PET系统放射时能量过程的混合型高斯曲线模型，对于不同的PET系统或者不同的PET系统参数，所选择的PET晶体的能量分布模型可能不尽相同，本实施例基于一般PET系统，给出了如下一种确定的算式一，以给出PET系统中晶体能量分布模型的一则实施例：

$>f\left(x\right)={\mathrm{Ae}}^{\frac{-{(x-\mu )}^2}{2{\sigma }^2}}+{\mathrm{ax}}^3+{\mathrm{\beta x}}^2+\gamma ;$>

其中，x为PET晶体的能量值(PET晶体接收到的γ光子的能量值)，f(x)为PET晶体的能量值为x的事件数，μ为PET晶体能量分布的峰值，σ为PET晶体能量分布的标准差，A、α、β、γ为所述能量分布模型的系数。上述算式一所示意的PET晶体的能量分布模型也是针对一个晶体的能量计算而言的。

所述PET原始数据指的是在PET系统仪器上实时采集的数据，通过实时生成的数据可分析PET当前硬件中每个探测器上晶体能量谱。从而通过能量分析检查和校验晶体性能。

继续参考图1，本实施例所述的方法还包括：

步骤S101，基于所述PET原始数据统计PET晶体能量谱，以获取所述PET晶体的能量分布模型的初次拟合参数。

基于本实施例步骤S100的相关内容可知，PET原始数据可适用于分析每个晶体能量谱，步骤S101对于晶体能量谱的统计手段主要涉及对晶体的定位及晶体能量谱中能量值的计数，本实施例给出了一种具体的统计方式，可以参考如下内容：

设本实施例基于PET系统硬件组成基本单位以及采集到的PET原始数据中，包括1152×96个PET晶体。

在统计过程中可以使用分块处理的方式：即对1152×96个晶体进行分块处理，大致分为36块，而每一块晶体个数为64×48，采用多线程技术，统计每个晶体中能量在0～511Kev的能量谱。

使用上述分块处理的方式可依每块晶体的数据基本处理单元对每块晶体中的晶体能量进行并行统计，可适用于对数据量非常大且数据结构复杂的计算数据，提高了运行效率。

按照上述PET设计的数据结构，为了定位每个晶体的具体坐标位置，需要预先得到36份晶体块的网格数据，所述网格数据把一个个的晶体框在晶格中。当统计得到的原始数据时，原始数据中携带了如下信息：

晶体位置，即晶体在上述晶格中的相对坐标值；

晶体的能量值，即0～511Kev能量的对应计数(各能量的对应计数是随机情况下产生的)。

将上述PET原始数据按照一定规则解析后，可将对应参数解析成晶体坐标值及晶体的能量值，依据解析得到的坐标值和晶格匹配后即可对该坐标位置的晶体进行能量谱统计，方便后续步骤的原始数据获取。

简言之，统计PET晶体能量分布，也就是PET系统的数据处理通道对该晶体通过相应数据处理通道的光子进行计数。本实施例的上述具体统计方式仅是给出了一种较佳的实施例，而非对于上述统计过程的限定。

当然，在步骤S101的实现过程中，基于所统计的PET晶体能量谱可以获得PET晶体的能量分布模型的初次拟合参数，参考上述算式一，所述拟合参数主要指算式一中所代表的每个晶体模型的参数，即A、μ、σ、α、β、γ。

下面给出上述拟合参数初始值的获得方式：

首先，根据PET晶体能量谱的统计结果获取PET晶体能量分布的标准差σ，σ满足：其中，i代表能量分布范围0～511kev，x_i代表晶体在能量分布范围0～511kev下的能量值(能量计数)，μ则代表晶体能量分布的峰值。

其次，根据PET晶体能量谱的统计结果获得PET晶体能量分布的峰值μ，峰值μ是通过统计所有晶体在0～511Kev范围内能量最大计数的能量值所得到的，然而对于统计所得到的PET晶体的能量谱而言，在获取晶体能量分布峰值μ的过程中，一般所获得的仅是一段区间内的局部最大值：

按照统计学和曲线分布，通常正常晶体实际能量最大值，也即上述峰值μ会分布在[μ-3σμ+3σ]区间。此时，直接认定所述局部最大值为峰值μ的初始值。

但是，由于晶体物理原因，性能有差异，某些晶体实际数据会出现一定偏差，在拟合寻峰时，则会出现如下极端情况：

当所统计的晶体能量实际数据分布在能量谱中呈现双峰，则峰值μ的初始化值会出现两个局部最大值，这时选择局部最大值最接近区间[μ-3σ，μ+3σ]的局部最大值作为峰值μ的初始值。

若所统计的晶体能量实际数据分布在能量谱中仅出现一个峰值，但该峰值出现在区间[0，μ-3σ]或区间[μ+3σ，511]，则也按照数据中所统计的局部最大值，赋值于峰值μ的初始值，但最终得到的结果可能与真值具有一定差距，得到的结果需要手动调节。

再次，系数A、α、β、γ可以依据如下方式赋予拟合初始值：

在一种PET系统中，可认为参数α、β、γ的初始值为0，且设x＝μ，从而求得A所估计的最大值，其中A可视为晶体能量谱中峰值所对应的能量谱幅度系数。

当然，对于PET系统结构、物理参数或者使用性能的不同，参数α、β、γ的初始值可以有不同的设置数值，且A所估计的数值也有所不同，在其他实施例中，赋予上述拟合参数初始值的设定方式可以不受上述方式的限制，但对于参数初始值中峰值μ的寻峰方式，本实施例所公开的手段则为较优方式，其结合系统对于局部最大值选取方式，增添了基于局部最大值至峰值初始值的幅值过程，对系统寻峰的自动化及准确性都具有贡献意义。

继续参考图1，本实施例所述的方法还包括：

步骤S102，基于所述能量分布模型和初次拟合参数，采用LM迭代算法计算所述能量分布模型的拟合参数，以获取所述PET晶体的能量分布。

步骤S102描述的是：基于步骤S100所建立的能量分布模型及步骤S101所获得的有关能量分布模型参数初始值，进行LM算法计算最终能量分布模型参数的拟合过程。

所述采用LM迭代算法计算所述能量分布模型的拟合参数时，若拟合过程中，迭代次数达到预先设定的迭代次数或LM算法收敛时迭代中止，迭代中止时则输出所述PET晶体的能量分布曲线，该曲线即表征了所述能量分布模型的最终拟合结果。

使用一般意义上的LM迭代算法，利用预先设定的迭代次数及算法规定的收敛情况，则可获得晶体能量分布模型的拟合结果。

在上述方案的基础上，本实施例还提供了一种较佳的LM算法流程，其结合了拟合优度(R-square)，使用了拟合优度对一般意义上的LM算法流程进行评价，从而形成了新的LM算法。在具体实施方案中，可以设定拟合优度的拟合要求(比如设定执行区间为98％)，并同时设定迭代次数，结合上述两个指标来控制LM算法迭代的中止，即：当LM流程满足拟合优度则停止迭代，当一直未满足拟合优度条件，但达到所设定的迭代次数，也停止本地的迭代过程，停止迭代后则输出PET晶体的能量分布曲线。

本实施例涉及的结合拟合优度求取能量分布模型的LM算法过程可参考图2，包括如下步骤：

步骤S200，输入初始参数值。

基于拟合参数A、μ、σ、α、β、γ的初始参数值A0、μ0、σ0、α0、β0、γ0获取f(x_i)⁰，f(x_i)⁰为初次拟合得到的PET晶体的能量分布，x_i代表能量分布0～511Kev上接收到的能量值。

步骤S201，计算初次误差矩阵D₀。

误差矩阵D₀的方程式为：

D₀＝[F_i-f(x_i)⁰]，i＝0～511。

其中，F_i为根据PET原始数据所统计的能量谱获得晶体能量在0～511kev上的分布情况，f(x_i)⁰为根据初次拟合所得的拟合参数获取到的晶体能量在0～511kev上的分布情况。f(x_i)⁰是根据初始参数值获取的，F_i则基于晶体能量分布的统计数据获得的实际分布(参考步骤S101的相关内容)。

步骤S202，计算初次拟合误差d₀。

根据初始参数值获取到的f(x_i)⁰以及根据PET原始数据所统计的能量分布F_i，求取初次拟合误差d₀，d₀满足：

$>d_0=\underset{i=0}{\overset{511}{\Sigma }}{(F_i-f{\left(x_i\right)}^0)}^2,i=0~511.$>

步骤S203，求解雅克布矩阵。

雅克布矩阵M的求解过程可参考如下方程式，雅克布矩阵体现了一个可微方程与给出点的最优线性逼近，类似于多元函数的导数。

$>\frac{\partial f\left(x\right)}{\partial A}=e^{b{(x-\mu )}^2};$>

$>\frac{\partial f\left(x\right)}{\partial b}=A{(x-\mu )}^2{e}^{b{(x-\mu )}^2};$>

$>\frac{\partial f\left(x\right)}{\partial \mu }=-2\mathrm{Ab}(x-\mu )e^{b{(x-\mu )}^2};$>

$>\frac{\partial f\left(x\right)}{\partial a}=x^3;$>

$>\frac{\partial f\left(x\right)}{\partial \beta }=x^2;$>

$>\frac{\partial f\left(x\right)}{\partial \gamma }=1;$>

其中，b为雅克布矩阵方程式的系数，x代表能量分布上接收到的连续的能量值，f(x)代表上一次拟合参数下所得的f(x)ⁿ，n为大于或等于1的自然数，对于输入的上一次拟合参数为初次拟合参数而言，则f(x)为f(x)⁰。

步骤S204，计算海瑟矩阵。

求得雅克布矩阵后依次通过雅克布矩阵求转置后相乘，得到海瑟矩阵。它是一个多元函数的二阶偏导数构成的方阵，描述了函数的局部曲率，海瑟矩阵H的方程式为：

H＝M^T×M。

步骤S205，计算混合海瑟矩阵。

为了在迭代过程加速，可引入阻尼因子并将该阻尼因子作用于海瑟矩阵，故而得到混合海瑟矩阵。阻尼因子根据每次参数的迭代逼近而做增大或减小的调整，当阻尼因子很小时，步长等于牛顿法步长，当阻尼因子很大时，步长约等于梯度下降法的步长，从而保证寻峰过程从谷底到波峰准确前进。混合海瑟矩阵H_mix的方程式为：

H_mix＝H×λ；

其中，λ为所述阻尼系数。

步骤S206，计算拟合参数步长矩阵及新的拟合参数。

计算所述拟合参数步长矩阵是为了获得能量分布模型f(x)的下一次拟合参数，也即新的拟合参数，步长矩阵S满足如下方程式：

S＝∑M^T×D×H_mix^-1；

其中，D为本次拟合参数下所得的误差矩阵，比如，若本次拟合为初次拟合，则D为D₀；H_mix^-1为混合海瑟矩阵的逆矩阵，满足H_mix^-1×H＝E，E代表主对角线元素为1，其他元素为0的矩阵。

依据步长矩阵的方程式得到的结果是有关获得下一次拟合涉及的各拟合参数的信赖步长，可以依据上述步长矩阵，得到下一次拟合的各拟合参数。

比如，若初次拟合下的步长矩阵S中，元素S[0]～S[5]分别代表下一次拟合时，拟合参数A、μ、σ、α、β、γ的信赖步长，则下一次拟合的拟合参数为：

A＝A1＝A0+S[0]；

μ＝μ1＝μ0+S[1]；

σ＝σ1＝σ0+S[2]；

α＝α1＝α0+S[3]；

β＝β1＝β0+S[4]；

γ＝γ1＝γ0+S[5]。

通过步骤S205的混合海瑟矩阵可计算出在上一次迭代的基础上，各拟合参数分别要增加的信赖步长。从而得到本次迭代时新的拟合参数。

步骤S207，基于新的拟合参数求得新的能量分布及误差矩阵。

将步骤S206所得的新的拟合参数带入能量分布模型即可得新的能量分布，新的误差矩阵D_n满足：

D_n＝[F_i-f(x_i)ⁿ]，i＝0～511；

其中，n为本次拟合次数，f(x_i)ⁿ为根据第n次拟合所得的拟合参数获取到的晶体能量在0～511kev上的分布情况。

步骤S208，计算新的拟合误差。

基于第n次能量分布f(x_i)ⁿ与基于PET原始数据所统计的能量分布(实测值)F_i，可得第n次拟合误差d_n，d₀满足：

$>d_n=\underset{i=0}{\overset{511}{\Sigma }}{(F_i-f{\left(x_i\right)}^n)}^2,i=0~511.$>

步骤S209，判断新的拟合参数是否符合拟合优度的条件。

在新的拟合结束后都计算该次拟合的拟合优度RNew，RNew满足如下方程式：

RNew＝1-(Q/P)^1/2；

$>Q=\underset{i=0}{\overset{511}{\Sigma }}{(y_i-{y^{\prime }}_i)}^2;$>

$>P=\underset{i=0}{\overset{511}{\Sigma }}{y_i}^2;$>

其中，y_i＝f(x_i)⁰或f(x_i)ⁿ，y′_i＝F_i。

在求取拟合优度的方程式中，Q是基于能量0～511Kev范围上能量分布的数组的残差平方和，P是基于PET原始数据能量统计所得的0～511Kev范围上能量分布的样本数组的残差平方和。

当一次拟合的拟合误差小于该次拟合的拟合优度RNew，也即满足d_n＜RNew，则直接执行步骤S210，否则，执行重新执行步骤S203～S209。

本实施例还提供另一种终止LM迭代的方法，结合拟合优度和预设迭代次数来停止LM迭代过程。当满足拟合优度条件时停止迭代，执行步骤S210；当一直未满足拟合优度条件，但达到了预设迭代次数，也停止LM迭代过程，执行步骤S210。

步骤S210，基于本地拟合参数输出最终能量分布。

所述本地拟合参数实际上为历次拟合时拟合误差最小值所对应的拟合参数。

上述步骤S200至S201已能达到本发明技术方案的目的。其使用了拟合优度评价LM算法流程是否满足要求，若迭代时的拟合参数下的拟合误差小于拟合优度，则迭代中止，上述迭代时的拟合参数即为最终的拟合参数。

本实施例还提供另一种LM算法过程，其还基于如下思路对上述LM算法过程进行了改进：

步骤S209判断拟合参数若还不满足拟合优度的条件，则进一步进行迭代计算，该迭代计算方式与图2所示的LM算法过程有所不同：

可基于前一次迭代时的拟合参数求得的拟合方程来计算残差平方和Q_m-1，并与基于本次迭代的拟合参数所得的拟合方程计算的残差平方和Q_m做比较：

若Q_m大于Q_m-1，则进一步扩大步骤S205中的阻尼系数λ，基于新的阻尼系数λ可更新的步长矩阵，从而得到更新后的拟合参数；

若Q_m小于Q_m-1，则本次过程中缩小阻尼系数λ，且对上一次迭代得到的结果不进行更新，直到满足小于拟合优度的条件或者完成预设的迭代次数为止，终止迭代过程。得到最终拟合参数。

m为大于或等于1的自然数，其代表迭代次数。

将上述参数残差平方和Q_m、Q_m-1替换为拟合误差d_m、d_m-1来进行上述判断也是可行的。

继续参考图1，本实施例所述的方法还包括：

步骤S103，根据所获取的PET晶体的能量分布完成PET晶体能量峰值的识别。

在获取到LM算法迭代过程中最终输出的PET晶体的能量分布f(x)后，可以根据现有技术中的任意一种方式完成PET晶体能量峰值的输出，最常见的是直接输出晶体能量分布算式的参数μ的数值。

但是，基于晶体的能量分布所输出的能量峰值是针对一个晶体而言的，在对个别晶体输出的能量峰值进行识别时，还涉及PET晶体能量的自动调整或手动调整。

本实施例给出了一种较为优选的晶体能量峰值识别方式，包括如下步骤：

根据所述PET晶体的能量分布获取所述PET晶体的能量峰值；

对所述PET晶体的能量峰值作伪彩化处理，以获取所述PET晶体的能量峰值的伪彩化图像；

基于所述PET晶体的能量峰值的伪彩化图像，调整所述PET晶体的能量峰值。

上述手段基于PET晶体输出的能量峰值，将各个晶体的能量峰值转化为伪彩化图像，从而便于直观地用色彩变化规律为后续人工调节提供参照。

以下给出了一种对所述PET晶体的能量峰值作伪彩化处理的具体过程：

对PET晶体进行分组处理；

基于所拟合的PET晶体的能量分布统计每组晶体的能量分布置信区间；

将所述置信区间的边界值设定为伪彩化处理的最小阈值及最大阈值；

根据所述伪彩化处理的最小阈值及最大阈值对处于所述最小阈值及最大阈值的能量值进行上色(也即彩色映射处理)，以获得PET晶体的能量峰值的伪彩化图像。

所述每组PET晶体的能量分布置信区间为[(μ′-3σ′)，(μ′+3σ′)]，μ′为每组PET晶体能量峰值，σ′为每组PET晶体能量分布的标准差。

在一个实例中，仍以PET系统具有1152×128个晶体为例，将上述晶体分为36组处理，每一组晶体个数为64×64，在步骤S103中，完成全部晶体的能量峰值识别后，可形成基于每组晶体的伪彩化处理，并依据每组晶体能量分布置信区间进行能量值的上色，对于一组晶体可产生一个对应该组晶体峰值的文件(Peak文件)。可对相应峰值调整工具添加自动播放模式，输入各组对应的峰值文件，以便形成直观的伪彩化图像。

在调整时，可对完成的36组晶体的峰值文件进行遍历浏览进一步确认是否有漏洞，也可点击任意一组晶体的峰值文件，以形成伪彩化图像来检查。上述手段可利用伪彩化图像的整体分布来快速预览各晶体峰值情况，不需要一个个的遍历晶体去检查其峰值情况，从而有利于峰值调整的效率，也便于手动调整个别晶体的峰值情况。

如图3所示，是一组晶体峰值文件的伪彩化图像，其显示了64×64晶体的峰值情况，图3中用灰度示意了相关色彩映射，在伪彩化处理中，其可根据置信区间的边界值对晶体峰值能量值进行上色，该色彩映射是任意的。其中，图3右侧显示了灰度示意下灰度与峰值能量值之间的关系。

结合图4，通过图3中的行项输入栏及列项输入栏，可定位图3示意的伪彩化图像中具体的一个晶体坐标(X0，Y0)，其中，X0表示行项输入的具体数值，Y0表示列项输入的具体数值，并于图4中显示关于该晶体坐标所对应晶体的拟合结果，所述拟合结果包括该晶体的能量峰值分布情况及预计该能量峰值所得LLD阈值。

从图4中可知，实线示意的是坐标(X0，Y0)所对应晶体的能量峰值分布情况的示意图，而点线则代表基于该能量峰值所得到的能量鉴频器LLD阈值，可手动拖动实线所示意的峰值曲线来手动修正峰值分布情况，此时点线所代表的LLD阈值也依据LLD阈值与峰值的对应关系表达式而被修正。

也可以基于图3整体的色彩分布是否均匀来发现可能存在拟合准确性问题的晶体，对于伪彩化图像上有色彩突变的点或区域，可通过直接点击的方式来查看对于该点或区域的晶体的拟合结果(所示拟合结果类似图4所示)，进而对拟合结果所示意的峰值分布及LLD阈值情况进行直接调节。

对拟合计算所得数据进行伪彩化处理后，能通过直观的颜色分布趋势进一步验证晶体能量峰值的输出情况，并可对个别晶体计算有误的能量峰值进行手动校正，从而实现了自动和手动结合的方法，得到最终拟合结果无限逼近真值。

基于本实施例所述获取PET晶体能量峰值的方法，如图5所示，本实施例还提供了一种设定PET系统能量鉴频器的方法，包括：

步骤S300，识别PET晶体能量峰值。

可依据步骤S100～S103的相关内容识别所述PET晶体能量峰值。

步骤S301，根据所述PET晶体能量峰值设定该PET晶体所对应能量鉴频器LLD阈值。

可依据如下算式计算能量鉴频器LLD阈值T：

T＝Pe×420/511：

其中，Pe为所述PET晶体能量峰值。

本发明虽然已以较佳实施例公开如上，但其并不是用来限定本发明，任何本领域技术人员在不脱离本发明的精神和范围内，都可以利用上述揭示的方法和技术内容对本发明技术方案做出可能的变动和修改，因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化及修饰，均属于本发明技术方案的保护范围。