CN113288189B - 一种基于ADMM-Net的PET时间校正方法 - Google Patents

Abstract

本发明公开了一种基于ADMM‑Net的PET时间校正方法，包括：(1)探测器采集来自于放射源的符合事件；(2)建立时间校正线性系统模型；(3)利用数据驱动方法的相容方程求解补偿值作为标签输入到网络中；(4)在原先传统算法求解正则项L1范数的基础上，采用神经网络的方式进行计算；(5)通过网络训练求解得到晶体补偿值。本发明方法将神经网络训练引入到时间校正线性模型中，利用现代计算机的算力，区别于传统算法的迭代过程，实现参数自适应，提高了PET时间校正的准确度，使符合时间窗口宽度变窄，减少随机噪声的影响。

Description

一种基于ADMM-Net的PET时间校正方法

技术领域

本发明属于PET成像技术领域，具体涉及一种基于ADMM-Net的PET时间校正方法。

背景技术

PET(全称为Positron emission tomography)即正电子发射断层成像，该医学影像技术基于核物理学和分子生物学，从分子层面上去观察细胞的新陈代谢活动，在某一些疾病中尤其是肿瘤、癌症的检测和早期预防提供有效依据；该方法将放射性同位核素标记的药物注入病人体内，通过血液进入循环系统，使得这些同位核素在人体内各组织器官中形成一定的浓度分布。由于放射性同位核素的半衰期较短且不稳定，很快发生衰变，衰变过程中所释放的正电子与附近的自由电子发生湮灭反应，产生一对方向几乎相反、能量相等且能量大小为511kev的伽玛光子对，这些光子对被PET系统中的探测器环接收，生成记录有光子能量、探测时间、计数率和探测器编号等相关信息，这些信息以List mode的形式储存在文件当中，并最终转化为Sinogram的形式，用于最终的图像数据重建。

为了实现更精确的诊断，医学领域对PET的性能提出了更进一步的要求，随着现代技术不断提升，TOF(Time-of-Flight)-PET应运而生，TOF-PET记录探测器探测到光子的时间信息，以实现软件符合；相较于之前的硬件符合，在商用PET中可以减少额外的硬件开支，因此TOF-PET对时间分辨率有着较高的要求。在实际情况下，有以下因素影响着PET系统的时间分辨率：PMT的延时、探测晶体的延时、后端电路的延时。

为了提高成像精度，需要对PET系统进行时间校正，PET时间校正可以分为直接校正和间接校正：

直接校正主要采用参考探测器法，即利用一个快速光电倍增管作为参考探测器，通过记录同一事件在参考探测器和PET系统探测器所记录的探测时间，来估计两者之间的时间差，从而获取PET系统的准确延时并借由此数值进行校正。

间接校正主要包括以下三种方式：

①特殊放射源法：经过特殊设计的放射源，放射源位于视场当中，由于其具体位置已知，因此事件到达的两端探测器的时间差值仍需加入位置补偿，以获取正确的校正数值；这种特殊放射源中，点源和旋转源使用较多。

②更为普通的放射源：研究者在思考是否不使用某些特殊的放射源，而是一些较为一般的放射源也可进行时间校正。因此，有在点源和旋转源使用的下，改进成圆柱体源：体积较小的点源和旋转源涉及到的探测晶体单元数目较少，通常就是对边180度附近的探测晶体单元，而圆柱体源的视场更大，涉及到的探测器对数目更多，在后面物体所占体积更大的情况下，结果相对会好一点。有研究者在此基础上进一步研究，使用普通物体，亦或是病灶本身的原始数据来进行校正；对于每一个响应线而言，以时间分布作为横轴，数目作为纵轴可以画出时间差直方图，对于某些情况：比如放射源位于中心，探测器位于两边，该曲线会有一个钟形趋势图，再通过高斯拟合，来得到更为圆滑的钟形图，取其最高峰时所对应的横坐标，即为该探测单元的时间补偿数值。上面一种情况指的是特殊源时的做法，而对于普通物体以及病灶本身的原始数据，如果直接取其响应线钟形趋势图，会得到有毛刺的上下拐动的曲线，若直接对改曲线进行高斯拟合，得到的曲线的趋势可能会与原始曲线有所差异，因此在数据的预处理上，还需要一步卷积的操作，来得到更为平滑的曲线。

③线性方程：利用特殊放射源所得到的数据，获取符合事件的相关信息，对于某一个符合事件而言，其探测到该事件的两个探测单元我们分别用I和J命名，最为基础的需要记录到信息为I和J的位置信息以及该次符合事件的探测单元I和J的时间差，然后将位置信息记录到系统测量矩阵当中，时间差信息记录到时间差向量当中，建立线性方程，进行求解。在早期的线性方程的建立中，系统矩阵以响应线的方式居多，这也是因为之前的特殊源在求晶体时间补偿值时，采用的就是这样的迭代计算，后来的研究者发现系统矩阵的建立还可以以符合事件数量进行呈现；两种方法均可用于计算，适用于不同的情况下。

线性方程的求解这种方式，早期以最小二乘作为最终目标，基于对方法的鲁棒性角度出发，许多研究者加入了L1正则项，全变分TV正则项，低秩Low Rank约束核范数作为正则项，甚至求解的时候损失函数以1范数而非2范数进行求解，各种方法适用不同情况。

发明内容

鉴于上述，本发明提供了一种基于ADMM-Net的PET时间校正方法，能够有效提高PET系统的时间分辨率。

一种基于ADMM-Net的PET时间校正方法，包括如下步骤：

(1)采集放射源为点源和圆柱体源两种情况下所得到的符合事件，并求出PET探测系统的时间分辨率的FWHM(full width at half maxima，半峰全宽)图；

(2)建立PET的时间校正系统模型；

(3)基于上述时间校正系统模型中加入1范数，得到用于PET时间校正的目标函数如下：

其中：|| ||₂表示2范数，|| ||₁表示1范数，A为系统测量矩阵，x为晶体补偿值向量，b为探测单元之间所测得的时间差向量，λ为权重系数；

(4)将上述目标函数转换为约束优化问题，并采用ADMM(交替方向乘子法)结合网络训练的方式即ADMM-Net对其进行求解。

进一步地，所述步骤(1)中作为放射源的点源或圆柱体源放置于PET探测系统中心，所获得的符合事件包括即时符合事件和延时符合事件，需要将延时符合事件从即时符合事件当中去除。

进一步地，所述步骤(2)中的时间校正系统模型表达式如下：

进一步地，所述步骤(4)中将目标函数转换为约束优化问题，具体表示为：

minimize f(x)+g(z)

subject to x-z＝0

其中：

z为中间变量。

进一步地，所述步骤(4)中采用ADMM-Net对目标函数进行求解，即将ADMM计算过程中每个更新步骤以网络层的方式进行实现，总共分为X、Z、M三层，同时采用相容方程(consistency condition)求解获得网络的真值标签。

进一步地，所述网络层X、Z、M对应的前向传播算法如下：

其中：

表示以λ_l/ρ_l为比较阈值的软阈值函数，λ为权重系数，ρ为大于0的惩罚参数，η为拉格朗日乘子的更新率，β和z为中间变量，I为单位矩阵，上标n表示迭代次数，下标l表示变量中的元素值序号，T表示转置，L为设定的向量维度(一般取8)。

进一步地，根据链式法则，所述网络层X、Z、M对应的反向传播算法如下，其计算操作顺序与前向传播相反；

其中：E为损失函数，β和z为中间变量，上标n表示迭代次数，下标l表示变量中的元素值序号。

进一步地，所述损失函数E的表达式如下：

其中：|Γ|为训练集的样本数量，

表示基于网络参数Θ下且以样本x作为输入所得到的网络输出结果，训练集中的样本数据x由传统方法计算得到，x^gt为样本数据x所对应的真值标签。

进一步地，采用相容方程求解网络真值标签的实现方式如下：

首先以径向角度去观察系统中的探测器，定义一圆周的探测器总数目为2N-1，然后根据步骤(1)获得的符合事件得到两个变量

和

其中Δ_t为采集的TOF时间差，M₀和M₁为数据增量，下标j和k分别表示当前探测器编号和对应符合响应线的探测器编号，i_t表示探测器j与探测器k的时间差刻度值，

表示探测器j与探测器k之间探测响应线的数目；

进而根据上述信息，建立一个2N×2N大小的矩阵U以及一个2N×1大小的向量Y，矩阵U中的元素值为：

其中：

为矩阵U中第j₁行第j₂列元素值，

为关于j₁和j₂的克罗内克函数，k的取值范围为以180度对边探测器作为中心向两边均匀正负方向作为对边探测器的数目，Δ_j2为探测器单元的权值(一般取1/2N)；

Y为用于求解补偿值的向量，向量Y中的元素值为：

其中：

为向量Y中第j₁行元素值，j₂＝(j₁+N+k)％(2N)，R为探测器半径，C为常数(一般取0)，k_max表示对边正方向或负方向上符合响应线的探测器数量。

最后，根据矩阵U和向量Y通过以下关系式将补偿值向量ω求出，并将其减去相应部分得到的补偿值矩阵后即作为网络的真值标签；

(U+μλ_maxI_d)·ω＝Y

其中：μ＝10-⁵，λ_max为矩阵U的最大元素值，I_d为单位矩阵。

本发明方法将神经网络训练引入到时间校正线性模型中，利用现代计算机的算力，区别于传统算法的迭代过程，实现参数自适应，提高了PET时间校正的准确度，使符合时间窗口宽度变窄，减少随机噪声的影响。

附图说明

图1(a)为本发明PET系统放射源为点源的结构示意图。

图1(b)为本发明PET系统放射源为圆柱体源的结构示意图。

图2为本发明PET时间校正方法的步骤流程示意图。

图3为相容方程求解方法的步骤流程示意图。

图4为本发明网络训练替代传统ADMM算法的步骤流程示意图。

图5为模拟补偿值的数值示意图。

图6为本发明方法与传统方法计数结果对比示意图。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

如图1(a)所示的PET探测系统中，放置于系统中心的点源向探测器发射，距离中心位置，探测器系统半径为215mm，轴向宽度为50mm，探测器晶体使用16×16的晶体阵列，晶体长度为30mm，底面边长为2.4mm，晶体之间间距为0.25mm。探测器Bank(亦称Rsector)即探测器晶体阵列，数目为32，呈圆环型包围，晶体材料采用LYSO，其中点源的半径为大小0.5mm，活度设置为1Mbq。

如图1(b)所示的情况是同一探测系统下，使用圆柱体源放置于视场中心，圆柱体源的半径20mm，长度70mm，活度设置与点源相同。

本发明基于ADMM-Net的PET时间校正方法的流程如图2所示，首先建立探测单元和探测时间差的系统矩阵模型：

其中：|| ||₂表示2范数，A为系统测量矩阵，x为晶体补偿值向量，b为探测单元之间所测得的时间差。

系统测量矩阵A的为m×n维矩阵，探测时间差向量b为m×1维向量，其中m为所采集到的符合事件的数目，晶体补偿值向量x为n×1维向量，n为系统探测单元数目，系统测量矩阵A为大型稀疏矩阵，每一行表示一次符合事件，每一行中有两个非零元素，1和-1，其他均为0，非零元素分别位于第i列和第j列，1和-1分别表示被减数和减数。向量b中所对应的矩阵A行数的数值，为该次符合事件的时间差Δt。

为了更直观的理解系统矩阵，其形式具体展示如下：

其中，u_i和u_j分别表示第i个探测单元和第j个探测单元，相应的向量b表示的是两个探测单元的时间差，在上式中，第一个事件，探测到编号为1的探测单元和编号为15的探测单元发生反应，时间差为1个单位时间差，根据不同探测器，一般以ns和ps作为时间单位。时间差的符号是由前一个探测单元减去后一个探测单元来确定下来的。同理，第二个事件表示编号为1的探测单元和编号为16的探测单元，时间差为-3。依次类推，每一个的符合事件以行向量的形式呈现，并且在向量b中有该符合事件的时间差数值。

探测单元数目n与晶体的区域划分挂钩，对于某一个特定的PET探测器(这里指的是晶体阵列耦合光学转换器件)，假设由256个晶体探测器组成的一个探测器而言，以16×16的晶体阵列为例：以参数section作为区域划分的数目，当section为1的时候，表示将整一个16×16看作一个晶体探测单元，当section为2的时候，将此16×16的晶体以两个边的中线作为划分，分成4(section数值的平方)个8×8个大小相同的晶体探测单元，那么当section为4的时候，分成16个4×4个大小相同的晶体探测单元，那么以此类推，section为4,8,…。section数值取值的最大值为最小的晶体探测单元数目，以所举例子而言是16，探测单元数目n的数值是section数值的平方。

得到系统矩阵之后，我们建立线性模型，并且求解L1范数，得到PET时间校正的目标函数为：

其中：|| ||₂表示2范数，|| ||₁表示1范数，λ为权重系数。

采用ADMM方法对目标函数进行求解，得到以下迭代步骤：

x^k+1:＝(A^TA+ρI)-¹(A^Tb+ρ(z^k-β^k))

z^k+1:＝S_λ/ρ(x^k+1+β^k)

β^k+1:＝β^k+x^k+1-z^k+1

在z步更新的时候，S为软阈值操作(soft thresholding operator)，定义为：

在输入到网络之前，需要以某一数值作为网络的标签进行训练，我们将MichelDefrise提出的相容方程所得到的数值作为标签进行输入，因为此方法无需特定放射源，图3所示了使用相容方程求解获取网络标签的流程。

参考Michel Defrise提出的相容方程离散化求解过程，以径向角度去观察探测器，将一圆圈的探测器的总数目作为2N-1处理，对应的探测器符号α和β，以及采集的TOF时间差Δt。其中当前探测器α_j的下标j取值范围从0到2N-1，符合对应的探测器β_j,_k的下标j与α_j的下标一致，k的取值范围以180度对边探测器作为中心，向两边均匀±方向作为对边探测器的数目；Δt的取值范围则与所使用的时间精度bin数值有关，当然也是以时间差0，向两边方向取值。

α_j＝jπ/N j＝0，...，2N-1

β_j，k＝α_j+π+kπ/N-k_max≤k≤k_max

t＝i_tΔt_i i_t＝-im，...i_m

然后另外根据得到的符合事件，得到两个变量

和

用于建立之后的补偿值方程进行求解。同时每一个探测器均有其权重。在参考文献[1](Defrise M,Rezaei A,Nuyts J.Time-of-flight PET time calibrationusing data consistency[J].Physics in Medicine and Biology,2018,63(10))当中，对于普通的没有处在探测器边缘的，将其权重设为π/N(因为以径向圆形来看的话全角度2π，然后探测器的数目为2N)，对于探测器在径向上有缺失的，则将缺失部分的权值均分到它旁边的两个探测器上，既Δ_ja＝Δ_jb＝(1+w/2)π/N，其中w为缺失的探测器数目，同时将缺失的探测器的权值均设为0，因为其在修正时没有起到作用。

随后有了以上信息之后，建立大小为2N×2N的矩阵U，其元素组成为：

其中：

为克罗内克函数，当两者相等时，输出为1，否则为0；对角项

另外右边为2N×1的向量Y，其元素组成为：

其中，j₂的取值与j₁和k相关，j₂＝(j₁+N+k)％(2N)以及C为

R为探测器半径。

根据以上步骤，得到矩阵U和Y，设定补偿值向量η，但是由于U是奇异的，因此还需要加入对角矩阵来避免求解的问题，根据参考文献[1]中，加入了μλ_maxI_d到U，其中μ取值为10^-5，λ_max为矩阵U最大的特征值，随后利用高斯消元的方法来求出η：

(U+μλ_maxId)·η＝Y

将向量η求出，并减去相应部分得到补偿值矩阵，作为网络的标签进行输入。

与传统算法不同的是，利用网络训练形式的ADMM-Net将以往ADMM所需要的每个更新步骤以网络层的方式进行实现，图4所示了使用网络训练替代传统ADMM算法的步骤流程。

前向传播的算法以传统ADMM的x,z,u更新步骤进行修改，在弄好前向传播算法之后，需要反向传播算法求解相应梯度，且要求顺序与前向传播相反。

将得到的参数带入其中，得到如下形式：

有了前向传播操作就有反向传播，反向传播的操作需要与前向传播顺序相反，对于M层的反向传播：

其中

对于Z层的反向传播：

其中

对于X层的反向传播：

通过以上步骤，将传统ADMM算法计算L1范数正则项上以网络的方式进行呈现。

以下我们通过蒙特卡洛模拟的点源数据和圆柱体源数据进行试验从而验证本发明在校正方面的准确性，搭建的系统参数已在前文叙述，旋转(线)源和点源实验的采集时间为100s，其中点源实验采集共2448785个符合事件，圆柱体源实验采集2238541个符合事件。

在实验当中，根据参考文献[1]中，加入模拟补偿值，模拟补偿值由两项组成，第一项模拟补偿数值为处于区间(-25,+25)(单位:ps)均匀分布的随机数，此项为每一个晶体探测器的模拟补偿，第二项模拟补偿数值为(-30,+30)均匀分布的随机数，此项为每一个晶体阵列探测器的模拟补偿；以此上述两项的和作为模拟补偿数值，如图5所示。

以点源实验为例，原始不加入噪声的半波宽FWHM为313.09ps，在加入噪声后半波宽FWHM被影响至323.24ps；在校正噪声数据时，使用本发明方法后半波宽FWHM(316.17ps)好于传统算法的319.13ps，具体见表1，本发明与传统方法矫正后的计数结果对比如图6所示。

表1

根据以上实验，通过与传统方法比较可见，本发明基于ADMM-Net的PET系统时间校正方法的校正结果有效地提高了准确率，提升了PET系统的时间分辨率。

上述对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明，熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。

Claims (8)

1.一种基于ADMM-Net的PET时间校正方法，包括如下步骤：

(1)采集放射源为点源和圆柱体源两种情况下所得到的符合事件，并求出PET探测系统的时间分辨率的FWHM图；

(2)建立PET的时间校正系统模型；

(3)基于上述时间校正系统模型中加入1范数，得到用于PET时间校正的目标函数如下：

其中：|| ||₂表示2范数，|| ||₁表示1范数，A为系统测量矩阵，x为晶体补偿值向量，b为探测单元之间所测得的时间差向量，λ为权重系数；

(4)将上述目标函数转换为约束优化问题，并采用ADMM结合网络训练的方式即ADMM-Net对其进行求解。

2.根据权利要求1所述的PET时间校正方法，其特征在于：所述步骤(1)中作为放射源的点源或圆柱体源放置于PET探测系统中心，所获得的符合事件包括即时符合事件和延时符合事件，需要将延时符合事件从即时符合事件当中去除。

3.根据权利要求1所述的PET时间校正方法，其特征在于：所述步骤(2)中的时间校正系统模型表达式如下：

4.根据权利要求1所述的PET时间校正方法，其特征在于：所述步骤(4)中将目标函数转换为约束优化问题，具体表示为：

minimize f(x)+g(z)

subject to x-z＝0

其中：

g(z)＝λ||z||₁，z为中间变量。

5.根据权利要求1所述的PET时间校正方法，其特征在于：所述步骤(4)中采用ADMM-Net对目标函数进行求解，即将ADMM计算过程中每个更新步骤以网络层的方式进行实现，总共分为X、Z、M三层，同时采用相容方程求解获得网络的真值标签。

6.根据权利要求5所述的PET时间校正方法，其特征在于：所述网络层X、Z、M对应的前向传播算法如下：

X:

Z:

M:

其中：

表示以λ_l/ρ_l为比较阈值的软阈值函数，λ为权重系数，ρ为大于0的惩罚参数，η为拉格朗日乘子的更新率，β和z为中间变量，I为单位矩阵，上标n表示迭代次数，下标l表示变量中的元素值序号，T表示转置，L为设定的向量维度。

7.根据权利要求5所述的PET时间校正方法，其特征在于：根据链式法则，所述网络层X、Z、M对应的反向传播算法如下，其计算操作顺序与前向传播相反；

M:

Z:

X:

其中：E为损失函数，β和z为中间变量，上标n表示迭代次数，下标l表示变量中的元素值序号。

8.根据权利要求7所述的PET时间校正方法，其特征在于：所述损失函数E的表达式如下：

其中：|Γ|为训练集的样本数量，

表示基于网络参数Θ下且以样本x作为输入所得到的网络输出结果，训练集中的样本数据x由传统方法计算得到，x^gt为样本数据x所对应的真值标签。

Images