CN109893154B - 一种基于低秩约束的pet时间校正方法 - Google Patents

Abstract

本发明公开了一种基于低秩约束的PET时间校正方法，包括：(1)探测器采集来自于放射源的符合事件；(2)建立时间校正系统模型；(3)加入低秩约束；(4)使用ADMM算法求解得到晶体补偿值。本发明通过在PET时间校正线性模型的基础上加入低秩约束，提高了对大噪声的鲁棒性，同时再解决问题时采用稀疏最小二乘的方式，可以解决因为大量数据所带来的计算时间问题。同时，本发明使PET系统能够获得更好的时间分辨率，从而利用TOF信息获取更好的空间分辨率，使PET系统能够在医学成像中提供更为准确的信息，为临床诊断提供更好的帮助。

Description

一种基于低秩约束的PET时间校正方法

技术领域

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

背景技术

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

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

目前较为普遍的PET时间方法有以下几种：

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

(2)特殊放射源法：经过特殊设计的放射源，放射源位于视场当中，由于其具体位置已知，因此事件到达的两端探测器的时间差值仍需加入位置补偿，以获取正确的校正数值。例如点源，由于放置于视场中心，则到达两端探测器的距离是相同的，其位置补偿为零；旋转源(线源)相较于点源略有不同，从径向角度看，其位置与视场正中心存在一定距离，因此近端探测器中较早接收信号，远端探测器较晚接收信号，因此在计算的时候，需要加入位置补偿，来达到类似于点源的效果。其思路相较于点源大致相同，但是因为涉及到的照射面更广，即对于单个探测器而言，有着更多的对边探测器与其有关，相较点源，其参考价值更为行之有效。

(3)线性方程：在特殊放射源的基础上，获取符合事件的相关信息，我们将符合事件的两个探测器以I和J命名，将其以PET时间校正模型的方式予以呈现，然后通过求解线性方程来得到相关补偿值，用于对PET系统时间延迟的校正。当然，目前有着一些直接利用患者数据进行校正的方式，亦可写成线性方程的形式，也正因此，线性方程中系统矩阵的呈现形式分为以响应线为系统矩阵形式和以符合事件数量为系统矩阵的形式两种。

在最初的线性方程求解时，以求解方程的最小二乘解为出发点，但由于其容易受到噪声的影响，为了提高方法的鲁棒性，许多研究者分别提出了在目标函数后面加入L1范数作为正则项，TV范数正则项，或者将目标函数从原先的二范数求解改为一范数求解以减少大噪声的影响。

对于处于同一PMT上的晶体探测器而言，其PMT的延迟以及后端电路的延迟来自同一个PMT，而对于同一闪烁晶体而言，不同晶体又有着不同的特性，鉴于以上原因，引入低秩约束来对探测器进行时间校正。

发明内容

鉴于上述，本发明提供了一种基于低秩约束的PET时间校正方法，能够有效提高PET探测系统的时间分辨率。

一种基于低秩约束的PET时间校正方法，包括如下步骤：

(1)采集放射源为点源或旋转线源两种情况下PET探测系统所得到的符合事件，对于任一符合事件，记录探测器中探测到该事件的一对晶体的探测时间，进而解析出系统关于晶体对探测时差的FWHM(Full width at halfmaximum，半峰全宽)图；

(2)建立PET探测系统的时间校正模型如下：

其中：x为晶体时间补偿向量，其维度为n，n为探测器中的晶体数量，x中每个元素值即对应各晶体的时间补偿值；b为晶体对探测时差向量，其维度为m，m为符合事件的数量，b中每个元素值即对应各符合事件的晶体对探测时间差；A为m×n大小的系统矩阵，|| ||₂表示2范数；

(3)使上述时间校正模型加入低秩约束，得到PET时间校正的目标函数如下，进而对该目标函数进行优化求解得到晶体时间补偿向量x；

其中：X为晶体时间补偿矩阵，|| ||_*表示核范数，μ为设定的权重系数；

(4)根据求得的晶体时间补偿向量x，对每一晶体的探测时间进行补偿，从而得到校正后的晶体对探测时差向量b，并将其以FWHM图的形式展示。

进一步地，所述系统矩阵A中的元素值为-1、0或1，且每一行即对应一个符合事件，每一列即对应一个晶体，一行中只有两个元素值分别为-1和1，其余元素值均为0，这两个元素值对应的晶体即为探测到符合事件的晶体对，其中探测时间在先的晶体对应元素值为1，探测时间在后的晶体对应元素值为-1。

进一步地，所述晶体时间补偿矩阵X的大小为k×p，k为探测器中的晶体阵列个数，p为单个晶体阵列中的晶体数量，由于n＝k×p，故晶体时间补偿向量x中的元素值按k行p列的形式对应排布后即可得到晶体时间补偿矩阵X。

进一步地，所述步骤(3)中采用ADMM(交替方向乘子法)算法对目标函数进行优化求解得到晶体时间补偿向量x。

进一步地，所述的步骤(1)中利用放置于PET探测系统中心的点源或距离PET探测系统中心一定位置的旋转线源来获取符合事件。

进一步地，所述的步骤(1)中获取得到的符合事件包含有即时符合事件与延时符合事件，由于延时符合事件是已知噪声信息，因此需要将延时符合事件从即时符合事件集合当中去除。

进一步地，所述的步骤(1)中当放射源为旋转线源，在解析计算FWHM图时需加入位置补偿，使得结果能以单个半波宽的形式呈现(未矫正的话波形会是双波峰)。

本发明通过在PET时间校正线性模型的基础上加入低秩约束，提高了对大噪声的鲁棒性，同时再解决问题时采用稀疏最小二乘的方式，可以解决因为大量数据所带来的计算时间问题。同时，本发明使PET系统能够获得更好的时间分辨率，从而利用TOF信息获取更好的空间分辨率，使PET系统能够在医学成像中提供更为准确的信息，为临床诊断提供更好的帮助。

附图说明

图1(a)为本发明PET探测系统放射源为点源的仿真示意图。

图1(b)为本发明PET探测系统放射源为旋转线源的仿真示意图。

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

图3为放射源为旋转线源时的位置补偿步骤流程示意图。

图4为放射源为旋转线源时添加位置补偿前后的效果示意图。

图5(a)为在随机噪声为真实数据20％时点源实验的FWHM示意图。

图5(b)为在随机噪声为真实数据20％时旋转线源实验的FWHM示意图。

图6为本发明方法与其他方法关于FWHM数值的效果对比示意图。

具体实施方式

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

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

如图1(b)所示的情况是同一探测系统下，使用旋转线源作为放射源来收集数据，旋转源距离视场中心距离为100mm，并以1deg/s的速度绕视场中心旋转，在计算旋转源的时候，选择相应半波宽能够被区分开来的对边探测器，相应地，需要加入位置补偿，其过程如图3所示，且位置补偿前后的效果如图4所示。

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

其中：|| ||₂表示2范数，A为系统矩阵(线性映射)，x为晶体补偿值，b为系统探测单元所测得的时间差。

系统矩阵(线性映射)A为m×n维矩阵，探测时间差b为m×1的向量，其中m为所采集得到的符合事件的数目，n为探测单元的数目，系统矩阵A中的每一行表示一个符合事件，每一行中在第i列为1，第j列为-1，向量b对应于该事件为Δt，表示对于此符合事件，第i个探测单元与第j个探测单元的探测时间差为Δt。

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

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

前面所讲的探测单元，其数目与晶体的区域划分有关，对于一个PET探测器，假设为32×32的晶体阵列，当划分区域section为1的时候，则将该32×32晶体看做为整体为2×2的，当section为2时，该晶体被划分为2×2，即每一个子区域块有16×16大小的晶体阵列，当section为4时，该晶体被划分为4×4，即每一个子区域块有8×8大小的晶体阵列，section为8,……32的情况下依次类推。对于PET系统，其bank数为32，晶体阵列采用16×16，因此其最小的探测单元数目为32×1×1＝32，最大的探测单元数目为32×16×16＝8192。

在得知了系统矩阵的结构之后，我们建立线性模型，并且加入低秩约束，得到PET时间校正的目标函数为：

其中：|| ||₂表示2范数，|| ||_*表示核范数，μ为权重系数，X是由x变换得到的晶体补偿矩阵。

基于上述系统模型将问题转化为约束优化问题：

minimize f(x)

其中：f为凸函数，

为非凸函数，目标函数f对应取||Ax-b||₂，约束项

取Rank(X)＝k。

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

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

其中：x步更新的时候，作为普通的最小二乘问题求解，由于系统矩阵A为大型稀疏矩阵，因此可以使用方法LSQR进行求解。

x步更新的具体形式为：

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

在z步更新的时候，是作为非凸问题的矩阵分解来进行考虑的，我们将x步更新得到的晶体补偿值向量x重新排列成晶体补偿矩阵X，

表示矩阵分解。

对于一般的矩阵分解问题而言，有以下形式：

s.t.M₁+M₂+…M_N＝X

其中，M₁,...,M_N∈R^m×n，X∈R^m×n为给定矩阵，同时，γ_i＞0为权衡参数，其目标是将矩阵X分解成各部分M_i之和，我们需要得到其中的低秩矩阵部分。

对应于我们的问题，稀疏矩阵表示的是晶体阵列中的大噪声，而低秩矩阵则代表是对应某个bank的固有延迟，由于整个晶体阵列连接在同一PMT当中，因此该晶体阵列共同拥有此PMT的固有延迟。

在我们的问题中，我们将晶体补偿矩阵X分解为小数值矩阵M₁、稀疏矩阵M₂以及低秩矩阵M₃，分别对应不同的目标项。

其中：

是F范数，该目标项是传统的最小二乘解，矩阵M₁当中的值通常较小，

该目标项使得M₂为稀疏，

该目标项使得M₃为低秩。将矩阵M₁用M₁＝X-(M₂+M₃)进行表示，则矩阵分解等价于以下式子：

其求解步骤为：

其中：

为M₁,...,M_N的均值，

则为近端算子。

近端算子

对应不同目标函数有不同解法：

其中：V和W为正交向量，(M₃)＝Vdiag(σ_s(M₃))W为M₃的奇异值分解，现在我们得到了低秩矩阵M₃，将其再变换成补偿值向量，作为下一次迭代的输入。

因此，最终的整体步骤可以由以下表示：

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

迭代收敛条件如下：

其中：n为探测单元数目。

另外，在旋转线源的实验当中，由于旋转源的位置已知，对于探测得到的某条响应线LOR，对应探测到的探测器分为近端探测器和远端探测器。为了使探测得到的时间差值能够用于校正，需要计算出对应的响应线相交旋转源的圆形轨迹的弦长。

之所以在旋转线源的实验中加入补偿，是因为我们在选取对边探测器的时候，其对边探测器所探测得到的波峰能够被较好的区分开，在这种情况下的响应线是用于校正的数据，若响应线与旋转线源的运动轨迹相切亦或是不相交，则这样的List数据是不被用于校正的。

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

在点源实验当中，对于每个探测器阵列而言，加入了固有延迟；另外，在实验中加入一定数量的随机事件，作为噪声影响结果。

由图5(a)和图5(b)可见，在噪声为20％的情况下，相较于其他方法，本发明方法所得到的结果最接近于原始没有噪声的FWHM。

由表1可见，在噪声分别为0～50％的情况下，本发明方法相较于其他方法能够得到在同等噪声情况下最好的FWHM数值，并且随着噪声的增长，受影响程度较低，说明方法鲁棒性较好；图6为表1更为直观的表达形式。

表1

通过以上的实验结果我们可以看出，本发明基于低秩约束法的PET系统时间校正方法的校正结果有效地提高了准确率，提升了PET系统的时间分辨率。

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

Claims (7)

1.一种基于低秩约束的PET时间校正方法，包括如下步骤：

(1)采集放射源为点源或旋转线源两种情况下PET探测系统所得到的符合事件，对于任一符合事件，记录探测器中探测到该事件的一对晶体的探测时间，进而解析出系统关于晶体对探测时差的FWHM图；

(2)建立PET探测系统的时间校正模型如下：

其中：x为晶体时间补偿向量，其维度为n，n为探测器中的晶体数量，x中每个元素值即对应各晶体的时间补偿值；b为晶体对探测时差向量，其维度为m，m为符合事件的数量，b中每个元素值即对应各符合事件的晶体对探测时间差；A为m×n大小的系统矩阵，||||₂表示2范数；

(3)使上述时间校正模型加入低秩约束，得到PET时间校正的目标函数如下，进而对该目标函数进行优化求解得到晶体时间补偿向量x；

其中：X为晶体时间补偿矩阵，|| ||_*表示核范数，μ为设定的权重系数；

(4)根据求得的晶体时间补偿向量x，对每一晶体的探测时间进行补偿，从而得到校正后的晶体对探测时差向量b，并将其以FWHM图的形式展示。

2.根据权利要求1所述的PET时间校正方法，其特征在于：所述系统矩阵A中的元素值为-1、0或1，且每一行即对应一个符合事件，每一列即对应一个晶体，一行中只有两个元素值分别为-1和1，其余元素值均为0，这两个元素值对应的晶体即为探测到符合事件的晶体对，其中探测时间在先的晶体对应元素值为1，探测时间在后的晶体对应元素值为-1。

3.根据权利要求1所述的PET时间校正方法，其特征在于：所述晶体时间补偿矩阵X的大小为k×p，k为探测器中的晶体阵列个数，p为单个晶体阵列中的晶体数量，由于n＝k×p，故晶体时间补偿向量x中的元素值按k行p列的形式对应排布后即可得到晶体时间补偿矩阵X。

4.根据权利要求1所述的PET时间校正方法，其特征在于：所述步骤(3)中采用ADMM算法对目标函数进行优化求解得到晶体时间补偿向量x。

5.根据权利要求1所述的PET时间校正方法，其特征在于：所述的步骤(1)中利用放置于PET探测系统中心的点源或距离PET探测系统中心一定位置的旋转线源来获取符合事件。

6.根据权利要求1所述的PET时间校正方法，其特征在于：所述的步骤(1)中获取得到的符合事件包含有即时符合事件与延时符合事件，由于延时符合事件是已知噪声信息，因此需要将延时符合事件从即时符合事件集合当中去除。

7.根据权利要求1所述的PET时间校正方法，其特征在于：所述的步骤(1)中当放射源为旋转线源，在解析计算FWHM图时需加入位置补偿，使得结果能以单个半波宽的形式呈现。

Images