CN105212956A - 一种基于ist的次晶体级pet系统时间修正方法 - Google Patents

Abstract

本发明公开了一种基于IST的次晶体级PET系统时间修正方法，包括：(1)扫描生物组织获取single？model数据；(2)对探测器单元进行次晶体级分割；(3)对数据进行预处理；(4)将L1范数约束加到线性方程式中；(5)使用IST方法求得时间修正序列；本发明方法通过将L1范数引入到时间修正线性过程中，并将时间修正水平从探测器级提升到次晶体级，来获取更精确的时间修正序列。这种方法有效地加强了修正算法对噪声的过滤，提高了PET时间修正的准确率，提升了PET系统的空间分辨率，使PET系统能在医学领域上具有更高的使用价值。

Description

一种基于IST的次晶体级PET系统时间修正方法

技术领域

本发明属于PET成像技术领域，具体涉及一种基于IST(IterativeShrinkage-Thresholding，迭代收缩阈值)的次晶体级PET系统时间修正方法。

背景技术

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

近几年PET在实际医学领域的应用日趋广泛，但与此同时，临床上对PET成像的要求也随之变高，越来越多的医学领域需要PET能提供更高的空间分辨率。为了实现更高的空间分辨率，一种新的PET系统被称为TOF(Time-of-Flight)-PET被广泛的应用在相关的临床医学领域。TOF-PET的基本原理是通过分别记录光子对被探测到精确时间来提升空间分辨率。因此，TOF-PET对PET系统的时间分辨率有着很高的要求，但是在实际情况中，PET系统的时间分辨率常会收到探测晶体的延时、探测器部分的延时和后端电路的延时的影响，所以PET系统的时间分辨率会因此变差。所以，PET系统的时间修正对实现高分辨率PET成像是十分必要的。

目前，PET系统时间修正方法大致可分为三类：参考探测器法、特殊散射源法和线性转化法。第一类主要是用一个快速光电倍增管作为参考探测器，通过记录同一事件在参考探测器和PET系统探测器所探测到的时间，求两者之间的时间差来估计PET系统的延时从而对系统进行时间修正的；第二种方法则是使用一个特殊设计的放射性源来获取系统时间修正序列。在这个特殊设计的放射性源中，每个被探测到的事件在源内的具体位置都是已知的，由于所有伽马光子的传播速度都是光速，因此我们可以通过计算其飞行时间来获取理论上无延时的探测时间，这个理论值和实际PET系统测量值的偏差就是我们所求的系统的时间修正序列；第三种方法则是将PET时间修正问题转化为一个线性问题，通过最小二乘法来求解这一线性方程来获取时间修正序列。

但是以上这三种方法都有着各自的局限性。第一种方法需要一个额外的参考探测器来获取时间修正序列，而这个额外的探测器会增加整个系统的构建成本，并且为了获取较为准确的估计值，这种方法往往需要较长的采集时间，这都制约了其的普遍应用。而第二种方法虽然不虚额外的添加一个参考探测器，但是它却需要特别设计的放射性源来计算时间修正序列，这使其也不便于广泛的使用。第三种方法虽然对探测器和放射性源没有特殊的要求，但是由于所采集到信号存在噪声，所以这种方法的准确度不是很好，并且还受到系统矩阵大小的制约，无法应用于大型的PET系统中。

发明内容

针对现有技术所存在的上述技术问题，本发明提供了一种基于IST的次晶体级PET系统时间修正方法，能够有效提高PET系统的时间分辨率。

一种基于IST的次晶体级PET系统时间修正方法，包括如下步骤：

(1)对PET系统中的每个探测器进行次晶体级分割，每个探测器通过分割对应得到一个由n×n个晶体单元组成的晶体单元阵列，n为大于1小于m的自然数，m为探测器原晶体阵列的维度；

(2)利用次晶体级分割后的探测器对注入放射性示踪剂的生物组织进行扫描探测，得到多组LOR(lineofresponse，响应线)数据；

(3)对每一组LOR数据进行预处理，以剔除每组LOR数据中时间信息有极大偏差的Singlemodel数据记录，并确定每组LOR数据的探测时延；

(4)将PET系统时间修正过程转化为线性方程，通过L1范数对该线性方程进行约束得到以下目标函数L；进而根据由各组LOR数据探测时延组成的探测时延序列ΔT，对目标函数L进行最小化求解得到PET系统的时间修正序列T_cal；

L＝||AT_cal-ΔT||²+λ||T_cal||₁

其中：A为系统矩阵，||||为L2范数，||||₁为L1范数，λ为预设的权重系数；

(5)最后根据所述的时间修正序列T_cal对预处理后的各组LOR数据Singlemodel数据记录中晶体单元的探测时间进行修正。

进一步地，所述的步骤(1)中根据计算机内存空间、系统中探测器的数量及探测器所含的晶体数对探测器进行次晶体级分割，n为除去1和m以外m的任一约数。

进一步地，所述的步骤(2)中的每组LOR数据对应一对探测到同一偶合事件且分属于不同探测器内的晶体单元J_A和J_B，每组LOR数据包含多条Singlemodel数据记录，其中每条Singlemodel数据记录对应一个探测时间差即晶体单元J_A和J_B对于同一偶合事件的探测时间之差，每条Singlemodel数据记录包含晶体单元J_A和J_B的编号、对应探测时间差精度范围内晶体单元J_A和J_B对于同一偶合事件的探测计数以及对应每次探测计数晶体单元J_A和J_B对于同一偶合事件各自的探测时间。

进一步地，所述的步骤(3)中对每组LOR数据进行预处理的具体方法为：对于任一组LOR数据，将该组LOR数据中所有Singlemodel数据记录按探测计数做成直方图，并计算所有Singlemodel数据记录探测计数的均值，进而剔除探测计数小于均值的Singlemodel数据记录。

进一步地，所述的步骤(3)中确定每组LOR数据探测时延的具体方法为：对于预处理后的任一组LOR数据，将其中所有Singlemodel数据记录对应探测时间差的最大值作为该组LOR数据的探测时延。

进一步地，所述的步骤(4)中线性方程的表达式为A·T_cal＝ΔT。

进一步地，所述的步骤(4)中采用迭代收缩阈值算法对目标函数L进行最小化求解，具体算式如下：

$T_{cal}^{k+1}=S_{\mathrm{\λ}t}(T_{cal}^k-2{\mathrm{\μA}}^T({\mathrm{AT}}_{cal}^k-\mathrm{\Δ}T))$

其中：和分别为第k次迭代和第k+1次迭代的时间修正序列，S_λt()为收缩函数，μ为收敛步长，^T表示转置，k为迭代次数。

所述收缩函数S_λt()的具体操作为：对于收缩函数S_λt()括号内向量的任一元素x，通过以下公式对元素x进行收缩运算：

$x^{*}=\left\{\begin{array}{cc}x-\mathrm{\&alpha;}& x>\mathrm{\&alpha;}\\ x+\mathrm{\&alpha;}& x<\mathrm{\&alpha;}\\ 0& x=\mathrm{\&alpha;}\end{array}\right.$

其中：x^*为元素x经收缩运算后的结果，α为预设的收缩阈值。

所述收敛步长μ的表达式为μ＝β/||A^TA||，β为预设的收敛系数且为0到1之间的常数。

所述系统矩阵A的维度为M×N，M＝(N-1)N/2，N为PET系统内所有探测器晶体单元的总个数，系统矩阵A的具体表达如下：

$A=\left[\begin{array}{c}{A}_1\\ A_2\\ .\\ .\\ .\\ A_{N-1}\end{array}\right]$

其中：A₁～A_N-1均为系统矩阵A的子矩阵，对于任一子矩阵A_i，其维度为(N-i)×N，i为自然数且1≤i≤N-1；所述的子矩阵A_i中第i列向量的所有元素均为1，前i-1列向量的所有元素均为0，后N-i列向量所组成的方阵为主对角线元素均为-1的对角矩阵。

本发明PET系统时间修正方法通过利用L1范数约束与次晶体级分割方法，有效地提升了传统线性转换法对噪声的处理能力，并解决了原有方法不能应用于大型PET系统的问题，成功将PET系统时间修正提升到了次晶体级别，提高了PET系统时间修正序列估计值的准确率，使PET系统能获得更加优秀的时间分辨率，从而能更好的利用TOF信息来获得更好的空间分辨率，使PET系统能为临床诊断和药物研发提出更加精准有用的信息。

附图说明

图1为本发明PET系统时间修正方法的步骤流程示意图。

图2为本发明探测器次晶体级分割的示意图。

图3为数据采集和验证的流程示意图。

图4(a)为本发明方法与传统最小二乘法在探测器单元(普通系统，单环)下的时间修正效果对比结果示意图。

图4(b)为本发明方法与传统最小二乘法在次晶体级(普通系统，单环)下的时间修正效果对比结果示意图。

图4(c)为本发明方法与传统最小二乘法在完整PET系统下的时间修正效果对比结果示意图。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明PET系统时间修正方法进行详细说明。

如图1所示，本发明基于迭代收缩阈值法的次晶体级PET时间修正方法，包括如下步骤：

(1)利用探测器对注入放射性示踪剂(18F-FDG)的圆柱形体模(Cylinderphantom)进行扫描探测，记录探测器探测到的每一件偶合事件的探测时间和相应探测器编号，生成含有时间，探测器对和计数率信息的Singlemodel数据。

(2)根据计算机内存空间和具体的探测器数及每个探测器模块所含的晶体数，对每个探测器进行次晶体级分割；其中，对于每一个探测器模块，根据与其耦合的晶体阵列的晶体数量对其进行次晶体级分割，例如：当晶体阵列为32×32时，可选取的分割数e为32的约数即2,4,8,16，其中对于晶体阵列的行与列都进行相同分割数e的分割，如图2所示。具体的分割数取决于探测器总数、数据采集总时间和计算使用的电脑内存的大小。之后，按照分割从新记录每一条LOR上对应的新的时间，探测器对和计数率信息。

(3)对数据进行预处理，将每一条LOR上记录的事件数按探测时间做成直方图，求出直方图的均值，将计数率小于均值的数据剔除，把直方图峰值处所对应的时间记作这条LOR所对应的时间。这样变换可以剔除时间信息有极大偏差的噪声数据，获取每一条LOR上真实的时间，探测器对和计数率信息。

(4)将PET系统时间修正过程转化为线性过程，这个线性过程可表达为：

A·T_cal＝ΔT

其中：T_cal为所求的时间修正序列，其维数为系统中最小探测单元的总数；ΔT为探测单元的延时序列，其具体表达式为ΔT＝T_l-T_k，这里T_l和T_k分别代表对于某一特定事件探测单元l和k分别记录的时间信息，其维数为系统中记录的LOR的总数；A为系统矩阵，用于表示系统中每一条LOR对应的探测器单元对信息，其具体形式如下：

其中：un(·)为算法中最小的探测单元，n为系统中最小探测单元的总数，m为PET系统中记录的LOR的总数，理论上m的大小受到n大小的约束，可以表达为m＝(n-1)n/2。

为了获取精确的时间修正序列估计值，将L1范数约束加在线性方程中，则其表达式变为：

$\underset{T_{cal}}{\min }\left|\right|{\mathrm{AT}}_{cal}-\mathrm{\&Delta;}T|{|}^2+\mathrm{\&lambda;}|\left|T_{cal}\right|{|}_1$

其中：||||为L2范数，||||₁为L1范数，λ为权重系数，用于平衡式中第一项和第二项的估计精度。

(5)使用迭代收缩阈值法IST求解L1范数约束后的线性方程，得到所求时间修正序列。这里使用IST方法的目的在于避免在求解过程中涉及矩阵求逆问题，其具体求解过程如下：

$T_{cal}^{k+1}=S_{\mathrm{\&lambda;}t}(T_{cal}^k-2{\mathrm{tA}}^T({\mathrm{AT}}_{cal}^k-\mathrm{\&Delta;}T))$

其中：为第k+1次迭代后的PET时间修正序列估计值，为第k次迭代后的PET时间修正序列估计值，t为IST方法的收敛步长，其具体的表达式为t＝u/||A^TA||，这里u为0到1之间的常数。S_λt()为收缩操作，其具体表达式如下：

$S_{\mathrm{\&alpha;}}\left(x\right)=\left\{\begin{array}{cc}x-\mathrm{\&alpha;}& x>\mathrm{\&alpha;}\\ x+\mathrm{\&alpha;}& x<\mathrm{\&alpha;}\\ 0& other\end{array}\right.$

这里经过多次迭代之后，获取PET时间修正序列。迭代次数达到1000次或满足迭代停止判断是迭代结束，迭代停止判断的表达式如下：

$\frac{\left|\right|T_{cal}^{k+1}-T_{cal}^k\left|\right|}{\left|\right|T_{cal}^k\left|\right|}<{10}^{-5}$

以下我们通过实验来验证本实施方式的实用性和可靠性，图3所示了实验的基本流程。

实验中，我们使用的PET系统是日本滨松光电子有限公司的大脑PET系统(HITS-655000)，其包含5个探测器环，每个探测器环分为4的探测深度，每个探测深度由32个探测器单元组成，而每个探测器单元又和一个的32×32的晶体阵列耦合。这一系统总共包括655,360枚晶体。

我们所使用的CylinderPhantom的尺寸为直径200毫米，长度220毫米。所注射的18F-FDG的放射性浓度为7.82MBq，整体扫描时间为60分钟。

在实验中，选取的次晶体分割数为2，其具体分割方式如图2所示。这个值得选择是由于我们的电脑的内存为8G，扫描时间总长为60分钟，更大的分割数将使计算内存大于我们电脑的内存，并出现每条LOR计数率不足的现象。所以为了方便验证，我们的晶体分割数只选2。

为了检验时间修正序列的准确性，我们需要通过对比修正前后的系统时间分辨率。这里我们将放射性点源置于PET系统中心，记录其Singlemodel数据信息。之后使用之前获得的时间修正序列对数据进行修正，列出修正后的数据的时间直方图，来获取修正后的时间分辨率。本实验中所使用的放射性点源为22-Na，其直径为0.25毫米，被置于PET系统中心，总扫描时间为5分钟。

由于本发明方法主要是通过将L1范数约束引入时间修正线性方程和将时间修正过程从探测器级提升到次晶体级来提升时间修正序列估计值得准确性，通过使用IST方法来使方法适用于大数据计算。所以，为了能够充分的验证本实施方式，我们设计了三组实验来验证我们的方法，具体情况如图4所示；其中，图4(a)～(c)分别展示了本实施方式与传统最小二乘法在探测器级、小数据量(单环)的直方图结果、提升到次晶体级对时间分辨率的提升和本实施方式在次晶体级、大数据量(全系统)的直方图结果。其具体的时间分辨率如表1所示：

表1

通过以上的实验结果我们可以看出，本发明基于迭代收缩阈值法的次晶体级PET时间修正方法的修正结果的确在有效地提高了PET时间修正的准确率，提升了系统的时间分辨率，并使其可以有效的运用在大型PET系统的修正中。

Claims (10)

1.一种基于IST的次晶体级PET系统时间修正方法，包括如下步骤：

(1)对PET系统中的每个探测器进行次晶体级分割，每个探测器通过分割对应得到一个由n×n个晶体单元组成的晶体单元阵列，n为大于1小于m的自然数，m为探测器原晶体阵列的维度；

(2)利用次晶体级分割后的探测器对注入放射性示踪剂的生物组织进行扫描探测，得到多组LOR数据；

(3)对每一组LOR数据进行预处理，以剔除每组LOR数据中时间信息有极大偏差的Singlemodel数据记录，并确定每组LOR数据的探测时延；

(4)将PET系统时间修正过程转化为线性方程，通过L1范数对该线性方程进行约束得到以下目标函数L；进而根据由各组LOR数据探测时延组成的探测时延序列ΔT，对目标函数L进行最小化求解得到PET系统的时间修正序列T_cal；

L＝||AT_cal-ΔT||²+λ||T_cal||₁

其中：A为系统矩阵，||||为L2范数，||||₁为L1范数，λ为预设的权重系数；

(5)最后根据所述的时间修正序列T_cal对预处理后的各组LOR数据Singlemodel数据记录中晶体单元的探测时间进行修正。

2.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述的步骤(1)中根据计算机内存空间、系统中探测器的数量及探测器所含的晶体数对探测器进行次晶体级分割，n为除去1和m以外m的任一约数。

3.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述的步骤(2)中的每组LOR数据对应一对探测到同一偶合事件且分属于不同探测器内的晶体单元J_A和J_B，每组LOR数据包含多条Singlemodel数据记录，其中每条Singlemodel数据记录对应一个探测时间差即晶体单元J_A和J_B对于同一偶合事件的探测时间之差，每条Singlemodel数据记录包含晶体单元J_A和J_B的编号、对应探测时间差精度范围内晶体单元J_A和J_B对于同一偶合事件的探测计数以及对应每次探测计数晶体单元J_A和J_B对于同一偶合事件各自的探测时间。

4.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述的步骤(3)中对每组LOR数据进行预处理的具体方法为：对于任一组LOR数据，将该组LOR数据中所有Singlemodel数据记录按探测计数做成直方图，并计算所有Singlemodel数据记录探测计数的均值，进而剔除探测计数小于均值的Singlemodel数据记录。

5.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述的步骤(3)中确定每组LOR数据探测时延的具体方法为：对于预处理后的任一组LOR数据，将其中所有Singlemodel数据记录对应探测时间差的最大值作为该组LOR数据的探测时延。

6.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述的步骤(4)中线性方程的表达式为A·T_cal＝ΔT。

7.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述的步骤(4)中采用迭代收缩阈值算法对目标函数L进行最小化求解，具体算式如下：

$T_{cal}^{k+1}=S_{\mathrm{\&lambda;}t}(T_{cal}^k-2{\mathrm{\&mu;A}}^T({\mathrm{AT}}_{cal}^k-\mathrm{\&Delta;}T))$

其中：和分别为第k次迭代和第k+1次迭代的时间修正序列，S_λt()为收缩函数，μ为收敛步长，T表示转置，k为迭代次数。

8.根据权利要求7所述的PET系统时间修正方法，其特征在于：所述收缩函数S_λt()的具体操作为：对于收缩函数S_λt()括号内向量的任一元素x，通过以下公式对元素x进行收缩运算：

$x^{*}=\left\{\begin{array}{cc}x-\mathrm{\&alpha;}& x>\mathrm{\&alpha;}\\ x+\mathrm{\&alpha;}& x<\mathrm{\&alpha;}\\ 0& x=\mathrm{\&alpha;}\end{array}\right.$

其中：x^*为元素x经收缩运算后的结果，α为预设的收缩阈值。

9.根据权利要求7所述的PET系统时间修正方法，其特征在于：所述收敛步长μ的表达式为μ＝β/||A^TA||，β为预设的收敛系数且为0到1之间的常数。

10.根据权利要求1所述的PET系统时间修正方法，其特征在于：所述系统矩阵A的维度为M×N，M＝(N-1)N/2，N为PET系统内所有探测器晶体单元的总个数，系统矩阵A的具体表达如下：

$A=\left[\begin{array}{c}{A}_1\\ A_2\\ .\\ .\\ .\\ A_{N-1}\end{array}\right]$

其中：A₁～A_N-1均为系统矩阵A的子矩阵，对于任一子矩阵A_i，其维度为(N-i)×N，i为自然数且1≤i≤N-1；所述的子矩阵A_i中第i列向量的所有元素均为1，前i-1列向量的所有元素均为0，后N-i列向量所组成的方阵为主对角线元素均为-1的对角矩阵。