CN103279964A

CN103279964A - 一种基于prca的pet图像动态重建方法及系统

Info

Publication number: CN103279964A
Application number: CN2013101445436A
Authority: CN
Inventors: 刘华锋; 于行建
Original assignee: Zhejiang University ZJU
Current assignee: Zhejiang University ZJU
Priority date: 2013-04-23
Filing date: 2013-04-23
Publication date: 2013-09-04
Anticipated expiration: 2033-04-23
Also published as: CN103279964B; US9576379B2; US20160078646A1; WO2014172927A1

Abstract

本发明公开了一种基于PRCA的PET图像动态重建方法，包括：（1）采集、校正，得到符合计数矩阵；（2）对PET测量方程施加修正约束；（3）迭代估计PET动态浓度数据；（4）重建出各帧PET图像。本发明通过将采集到的不同帧的数据看作一个整体来进行重建，充分利用了PET数据在时间上的相关性，使所得结果更能体现出动态PET能表现出目标区域时间变化的特点；其次，本发明使用了背景和目标区域的方法，减小了背景对目标区域的干扰，加之在重建中添加了时间和空间的修正，使得重建结果的准确度更高，目标区域与背景的对比度提高，使的重建效果比传统的FBP和ML-EM算法更加优秀，更具有医学价值。

Description

一种基于PRCA的PET图像动态重建方法及系统

技术领域

本发明属于PET成像技术领域，具体涉及一种基于PRCA（鲁棒性主元分析法）的PET图像动态重建方法及系统。

背景技术

正电子发射断层成像（Positron emission tomography，PET）是一种基于核物理学、分子生物学的医学影像技术，它能够从分子层面上观察细胞的新陈代谢活动，为早期疾病尤其是肿瘤的检测和预防提供了有效依据。PET本质上是对病人体内药物的浓度分布进行成像，被注射入病人体内的放射性同位核素标记药物通过血液进入循环系统，这些物质在人体内各组织器官中将形成一定的浓度分布。由于放射性同位核素的半衰期较短，且极其不稳定，将很快发生衰变，衰变过程中所释放的正电子与附近的自由电子发生湮灭反应，产生一对方向几乎相反、能量相等，能量大小为511kev的伽玛光子对，这些光子由探测器环接收，再经由符合采集系统对这些带有放射性药物分布信息的成对光子进行处理生成投影数据（sinogram）。之后，通过相应的数学方法对投影数据进行反演求解，可重建出人体的放射性物质的空间浓度分布。

近几年PET在实际医学领域的应用日趋广泛，但与此同时，临床上对PET成像的要求也随之变高，越来越多的医学领域需要PET能提供更高成像分辨率和能够实时的对患者进行扫描，而与此相应的图像维数的扩展和采集数据量的急剧加大对现有的重建算法提出了挑战，并且这些需求也对计算机的计算能力和存储空间有极为苛刻的要求。

目前，PET成像的重建方法大致可分为两类：解析法和传统迭代统计法。前一类主要是FBP（滤波反投影法），计算速度快，但成像分辨率低并且存在伪影严重，而后一种现在最为常用的是ML-EM（最大似然最大期望）算法，这种算法虽然很大程度上提高了图像的分辨率，但依然不能很好的解决严重的伪影问题。而且，以上这两种方法传统上都是将每帧视为独立的个体，无法在重建中将PET图像的时间关联性表现出来，加之两种方法都无法将采集时的充满噪声的背景和目标区域分离，所以不能很好的应对噪声的干扰。

发明内容

针对现有技术所存在的上述技术问题，本发明提供了一种PRCA的PET图像动态重建方法及系统，能够解决计算机在进行动态图像重建的过程中准确率低的问题。

一种基于PRCA的PET图像动态重建方法，包括如下步骤：

（1）利用探测器对注入有放射性物质的生物组织进行探测，动态采集得到PET的n组符合计数向量，并对所述的符合计数向量进行校正，进而构建PET的符合计数矩阵；n为大于1的自然数；

（2）根据PET成像原理，建立PET的测量方程；

（3）通过对所述的测量方程在时间和空间上施加修正约束得到以下迭代方程；根据所述的符合计数矩阵通过以下迭代方程估计出PET浓度动态数据；

(X_{k + 1}, Z_{k + 1}) = \underset{X, Z}{\arg \min} {{| | G (X + Z) - Y - A_{k} | |}^{2} + {| | X - B_{k} + C_{k} | |}^{2} + r {| | HZ - D_{k} + E_{k} | |}^{2}}

其中：X_k+1为第k+1次迭代后的PET浓度分布矩阵，Z_k+1为第k+1次迭代后的PET背景浓度分布矩阵，G为系统矩阵，H为修正矩阵，Y为符合计数矩阵，r为修正系数，A_k为第k次迭代后的第一修正系数矩阵，B_k为第k次迭代后的第二修正系数矩阵，C_k为第k次迭代后的第三修正系数矩阵，D_k为第k次迭代后的第四修正系数矩阵，E_k为第k次迭代后的第五修正系数矩阵，X和Z分别为对应X_k+1和Z_k+1的随机化矩阵，k为自然数；

（4）对所述的PET浓度动态数据进行分解，得到n帧连续的PET图像。

所述的测量方程的表达式如下：

y＝Gx

其中：G为系统矩阵；y为校正后的符合计数向量；x为PET浓度分布向量。

所述的步骤（3）中，根据符合计数矩阵通过迭代方程进行迭代计算，则迭代收敛后的PET浓度分布矩阵即为PET浓度动态数据。

所述的迭代收敛条件如下：

\frac{{| | Y - G (X_{k} + Z_{k}) | |}_{F}}{{| | Y | |}_{F}} \leq ρ

其中：X_k为第k次迭代后的PET浓度分布矩阵，Z_k为第k次迭代后的PET背景浓度分布矩阵，ρ为预设的收敛阈值，||||_F表示F范数。

所述的第k次迭代后的第一修正系数矩阵A_k通过以下算式求得：

A_k＝A_k-1+G(X_k+Z_k)-Y

其中：A_k-1为第k-1次迭代后的第一修正系数矩阵，X_k为第k次迭代后的PET浓度分布矩阵，Z_k为第k次迭代后的PET背景浓度分布矩阵。

所述的第k次迭代后的第二修正系数矩阵B_k通过以下算式求得：

B_k＝U_k·diag(max(i_k-1,0))·V_k

W_k＝U_k·diag(i_k)·V_k

其中：W_k为第k次迭代后的待分解矩阵且W_k=X_k+C_k-1，X_k为第k次迭代后的PET浓度分布矩阵，C_k-1为第k-1次迭代后的第三修正系数矩阵，U_k、diag(i_k)和V_k分别为待分解矩阵W_k经奇异值分解后得到的时间参数相关矩阵、奇异值对角矩阵和中介系数矩阵，i_k为奇异值对角矩阵diag(i_k)中的对角线元素。

所述的第k次迭代后的第三修正系数矩阵C_k通过以下算式求得：

C_k＝C_k-1+X_k-B_k

其中：X_k为第k次迭代后的PET浓度分布矩阵，C_k-1为第k-1次迭代后的第三修正系数矩阵。

所述的第k次迭代后的第四修正系数矩阵D_k通过以下算式求得：

D_k＝sgn(HZ_k+E_k-1)·max(|HZ_k+E_k-1|-r,0)

其中：Z_k为第k次迭代后的PET背景浓度分布矩阵，E_k-1为第k-1次迭代后的第五修正系数矩阵。

所述的第k次迭代后的第五修正系数矩阵E_k通过以下算式求得：

E_k＝E_k-1+HZ_k-D_k

所述的修正系数r通过以下算式求得：

r = \frac{1}{\sqrt{\max (n, m)}}

其中：m为符合计数向量的维度。

本发明中，若符合计数向量的维度为m，PET浓度分布矩阵的维度为p×n，则系统矩阵G的维度为m×p，其表征了发射光子被探测器接收的概率，其受探测器结构、探测效率、衰减、死时间等因素的影响；符合计数矩阵为由n组符合计数向量组成的维度为m×n的矩阵。修正矩阵H为m×p的矩阵且为多级紧凑型小波分解中算子，其满足H^TH＝I，I为单位矩阵。

一种基于PRCA的PET图像动态重建系统，包括探测器和与探测器相连的计算机；

所述的探测器用于对注入有放射性物质的生物组织进行探测，动态采集得到PET的n组符合计数向量，n为大于1的自然数；

所述的计算机内加载有以下功能模块：

数据接收模块，用于接收所述的符合计数向量并对其进行校正，进而构建PET的符合计数矩阵；

浓度估计模块，用于根据所述的符合计数矩阵通过预设的迭代方程组估计出PET浓度动态数据；

PET成像模块，用于对所述的PET浓度动态数据进行分解，得到n帧连续的PET图像。

所述的浓度估计模块根据以下迭代方程组估计PET浓度动态数据：

(X_{k + 1}, Z_{k + 1}) = \underset{X, Z}{\arg \min} {{| | G (X + Z) - Y - A_{k} | |}^{2} + {| | X - B_{k} + C_{k} | |}^{2} + r {| | HZ - D_{k} + E_{k} | |}^{2}}

A_k＝A_k-1+G(X_k+Z_k)-Y

B_k＝U_k·diag(max(i_k-1,0))·V_k

W_k＝U_k·diag(i_k)·V_k

C_k＝C_k-1+X_k-B_k

D_k＝sgn(HZ_k+E_k-1)·max(|HZ_k+E_k-1|-r,0)

E_k＝E_k-1+HZ_k-D_k

r = \frac{1}{\sqrt{\max (n, m)}}

其中：其中：X_k+1为第k+1次迭代后的PET浓度分布矩阵，Z_k+1为第k+1次迭代后的PET背景浓度分布矩阵，G为系统矩阵，H为修正矩阵，Y为符合计数矩阵，r为修正系数，A_k为第k次迭代后的第一修正系数矩阵，B_k为第k次迭代后的第二修正系数矩阵，C_k为第k次迭代后的第三修正系数矩阵，D_k为第k次迭代后的第四修正系数矩阵，E_k为第k次迭代后的第五修正系数矩阵，X和Z分别为对应X_k+1和Z_k+1的随机化矩阵，A_k-1为第k-1次迭代后的第一修正系数矩阵，X_k为第k次迭代后的PET浓度分布矩阵，Z_k为第k次迭代后的PET背景浓度分布矩阵，W_k为第k次迭代后的待分解矩阵且W_k=X_k+C_k-1，C_k-1为第k-1次迭代后的第三修正系数矩阵，U_k、diag(i_k)和V_k分别为待分解矩阵W_k经奇异值分解后得到的时间参数相关矩阵、奇异值对角矩阵和中介系数矩阵，i_k为奇异值对角矩阵diag(i_k)中的对角线元素，E_k-1为第k-1次迭代后的第五修正系数矩阵，m为符合计数向量的维度，k为自然数；迭代收敛后的PET浓度分布矩阵即为PET浓度动态数据。

本发明的有益技术效果在于：通过将采集到的不同帧的数据看作一个整体来进行重建，充分利用了PET数据在时间上的相关性，使所得结果更能体现出动态PET能表现出目标区域时间变化的特点；其次，本发明使用了背景和目标区域的方法，减小了背景对目标区域的干扰，加之在重建中添加了时间和空间的修正，使得重建结果的准确度更高，目标区域与背景的对比度提高，使的重建效果比传统的FBP和ML-EM算法更加优秀，更具有医学价值。

附图说明

图1为本发明估计系统的结构示意图。

图2为本发明估计方法的步骤流程示意图。

图3(a)为本发明使用蒙特卡罗模拟在Level1等级的模拟图像。

图3(b)为本发明使用蒙特卡罗模拟在Level2等级的模拟图像。

图4(a)为Level1等级模拟图像的原始浓度结果示意图。

图4(b)为Level1等级的模拟图像采用ML-EM重建的浓度结果示意图。

图4(c)为Level1等级的模拟图像采用本发明方法重建的浓度结果示意图。

图5(a)为Level2等级模拟图像的原始浓度结果示意图。

图5(b)为Level2等级的模拟图像采用ML-EM重建的浓度结果示意图。

图5(c)为Level2等级的模拟图像采用本发明方法重建的浓度结果示意图。

图6为本发明使用的Phantom体模示意图。

图7(a)为Phantom体模的原始浓度结果示意图。

图7(b)为Phantom体模采用ML-EM重建的浓度结果示意图。

图7(c)为Phantom体模采用本发明方法重建的浓度结果示意图。

图8为Phantom体模重建结果中第40行像素值的对比示意图。

具体实施方式

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

如图1所示，一种基于PRCA的PET图像动态重建系统，包括探测器和与探测器相连的计算机；系统的操作执行流程图如图2所示；其中：

探测器用于对注入有放射性物质的生物组织进行探测，动态采集得到PET的n组符合计数向量；本实施方式中，探测器采用日本滨松公司生产的型号为SHR74000的PET扫描仪。

计算机内加载有数据接收模块、浓度估计模块和PET成像模块；其中：

数据接收模块用于接收符合计数向量并对其进行校正；正电子发射断层扫描仪探测人体内发出的放射性信号，经过符合和采集系统处理，形成原始符合事件。PET探测器记录的符合事件包括真符合、随机符合和散射符合。通过探测器的延时窗口和能量窗口对随机事件和散射事件进行校正，而后进行衰减校正，得到正弦图数据即校正后的符合计数向量y。

对于n组校正后的符合计数向量y，数据接收模块构建由这些符合计数向量y组成的PET符合计数矩阵Y，其为m×n的矩阵，符合计数矩阵Y中每一列对应着一组符合计数向量y。

浓度估计模块用于根据符合计数矩阵Y通过预设的迭代方程组估计出PET浓度动态数据。

根据PET成像原理，可以得到PET的测量方程的表达式如下：

y＝Gx

其中：y为校正后的符合计数向量且为m维向量；x为PET浓度分布向量且为p维向量；G为系统矩阵且为m×p维矩阵，其表征了发射光子被探测器接收的概率，其受探测器结构、探测效率、衰减、死时间等因素的影响。

通过对上述测量方程在时间和空间上施加修正约束得到以下迭代方程组：

(X_{k + 1}, Z_{k + 1}) = \underset{X, Z}{\arg \min} {{| | G (X + Z) - Y - A_{k} | |}^{2} + {| | X - B_{k} + C_{k} | |}^{2} + r {| | HZ - D_{k} + E_{k} | |}^{2}}

A_k＝A_k-1+G(X_k+Z_k)-Y

B_k＝U_k·diag(max(i_k-1,0))·V_k

W_k＝U_k·diag(i_k)·V_k

C_k＝C_k-1+X_k-B_k

D_k＝sgn(HZ_k+E_k-1)·max(|HZ_k+E_k-1|-r,0)

E_k＝E_k-1+HZ_k-D_k

r = \frac{1}{\sqrt{\max (n, m)}}

其中：X_k+1为第k+1次迭代后的PET浓度分布矩阵，Z_k+1为第k+1次迭代后的PET背景浓度分布矩阵，G为系统矩阵，H为修正矩阵，Y为符合计数矩阵，r为修正系数，A_k为第k次迭代后的第一修正系数矩阵，B_k为第k次迭代后的第二修正系数矩阵，C_k为第k次迭代后的第三修正系数矩阵，D_k为第k次迭代后的第四修正系数矩阵，E_k为第k次迭代后的第五修正系数矩阵，X和Z分别为对应X_k+1和Z_k+1的随机化矩阵，A_k-1为第k-1次迭代后的第一修正系数矩阵，X_k为第k次迭代后的PET浓度分布矩阵，Z_k为第k次迭代后的PET背景浓度分布矩阵，W_k为第k次迭代后的待分解矩阵且W_k=X_k+C_k-1，C_k-1为第k-1次迭代后的第三修正系数矩阵，U_k、diag(i_k)和V_k分别为待分解矩阵W_k经奇异值分解后得到的时间参数相关矩阵、奇异值对角矩阵和中介系数矩阵，i_k为奇异值对角矩阵diag(i_k)中的对角线元素，E_k-1为第k-1次迭代后的第五修正系数矩阵。

本实施方式中，PET浓度分布矩阵X的维度为p×n，该矩阵中每一列对应着一组PET浓度分布向量x；X₀、Z₀、A₀、B₀、C₀、D₀和E₀均为初始化给定的零矩阵，修正矩阵H为m×p的矩阵且为多级紧凑型小波分解中算子，其满足H^TH＝I，I为单位矩阵。

根据符合计数矩阵Y通过上述迭代方程组进行迭代计算，则迭代收敛后的PET浓度分布矩阵即为PET浓度动态数据；迭代收敛条件如下：

\frac{{| | Y - G (X_{k} + Z_{k}) | |}_{F}}{{| | Y | |}_{F}} \leq ρ

其中：ρ为预设的收敛阈值，||||_F表示F范数；本实施方式中ρ=10^-3；最大迭代次数预设为50次。

PET成像模块用于对PET浓度动态数据进行分解，得到n帧连续的PET图像；由于迭代收敛后的PET浓度分布矩阵为p×n的矩阵，该矩阵中每一列即对应一组PET浓度分布向量（即一帧PET图像数据），将每一列向量元素进行重组即可还原得到PET图像。

以下我们对本实施方式进行了两组验证实验，第一组实验采用的是蒙特卡罗模拟，第二组采用的是六个圆柱体的phantom体模进行模拟。

在第一组试验中，我们采用的是180度内48个采样角，所有采样结果被分为20帧，我们生成了两组不同噪声和计算率等级的模拟图像，如图3所示；其中Level1对应的是0.53%的随机噪声和0.03%的散射噪声，其每帧计数率约为7×10⁶；Level2对应的是1.28%的随机噪声和22.9%的散射噪声，其每帧计数率约为2.5×10⁵。我们对比了传统ML-EM和本实施方式的重建结果，如图4和图5所示，我们主要对比的是第8帧的图像结果。为了能更好地量化结果我们还分别求出了ML-EM和本实施方式重建结果的偏差和方差，结果如表1所示：

表1

在第二组中我们对实际的phantom体模数据进行了重建，我们采用的圆柱型体模，整个体模的大小为直径200mm，深度290mm，如图6所示，其中包含有六个不同直径的圆形区域，其直径分别为：37mm，28mm，22mm，17mm，13mm，10mm。在其中我们注入放射性浓度为107.92Bq/ml的F-18示踪剂，采样时间为120分钟，其对比结果如图7所示；同理，为了能更好的量化结果，我们选取图中第40行的像素值进行比对，其结果如图8所示。

通过以上的实验结果我们可以看出，本实施方式的重建结果无论是在实际的图像还是偏差和方差这两个方面的表现都要优于ML-EM的重建结果，由此我们可以看出RPCA算法有效地提高了PET动态重建的准确率并且大幅提高了图像目标区域和背景的对比度。

Claims

1.一种基于PRCA的PET图像动态重建方法，包括如下步骤：

（2）根据PET成像原理，建立PET的测量方程；

(X_{k + 1}, Z_{k + 1}) = \underset{X, Z}{\arg \min} {{| | G (X + Z) - Y - A_{k} | |}^{2} + {| | X - B_{k} + C_{k} | |}^{2} + r {| | HZ - D_{k} + E_{k} | |}^{2}}

2.根据权利要求1所述的PET图像动态重建方法，其特征在于：所述的测量方程的表达式如下：

y＝Gx

3.根据权利要求1所述的PET图像动态重建方法，其特征在于：所述的步骤（3）中，根据符合计数矩阵通过迭代方程进行迭代计算，则迭代收敛后的PET浓度分布矩阵即为PET浓度动态数据；迭代收敛条件如下：

\frac{{| | Y - G (X_{k} + Z_{k}) | |}_{F}}{{| | Y | |}_{F}} \leq ρ

4.根据权利要求1所述的PET图像动态重建方法，其特征在于：所述的第k次迭代后的第一修正系数矩阵A_k通过以下算式求得：

A_k＝A_k-1+G(X_k+Z_k)-Y

5.根据权利要求1所述的PET图像动态重建方法，其特征在于：所述的第k次迭代后的第三修正系数矩阵C_k通过以下算式求得：

C_k＝C_k-1+X_k-B_k

6.根据权利要求1或5所述的PET图像动态重建方法，其特征在于：所述的第k次迭代后的第二修正系数矩阵B_k通过以下算式求得：

B_k＝U_k·diag(max(i_k-1,0))·V_k

W_k＝U_k·diag(i_k)·V_k

7.根据权利要求1所述的PET图像动态重建方法，其特征在于：所述的第k次迭代后的第五修正系数矩阵E_k通过以下算式求得：

E_k＝E_k-1+HZ_k-D_k

8.根据权利要求1或7所述的PET图像动态重建方法，其特征在于：所述的第k次迭代后的第四修正系数矩阵D_k通过以下算式求得：

D_k＝sgn(HZ_k+E_k-1)·max(|HZ_k+E_k-1|-r,0)

9.一种基于PRCA的PET图像动态重建系统，包括探测器和与探测器相连的计算机；所述的探测器用于对注入有放射性物质的生物组织进行探测，动态采集得到PET的n组符合计数向量，n为大于1的自然数；其特征在于：

所述的计算机内加载有以下功能模块：

10.根据权利要求9所述的PET图像动态重建系统，其特征在于：所述的浓度估计模块根据以下迭代方程组估计PET浓度动态数据：

(X_{k + 1}, Z_{k + 1}) = \underset{X, Z}{\arg \min} {{| | G (X + Z) - Y - A_{k} | |}^{2} + {| | X - B_{k} + C_{k} | |}^{2} + r {| | HZ - D_{k} + E_{k} | |}^{2}}

A_k＝A_k-1+G(X_k+Z_k)-Y

B_k＝U_k·diag(max(i_k-1,0))·V_k

W_k＝U_k·diag(i_k)·V_k

C_k＝C_k-1+X_k-B_k

D_k＝sgn(HZ_k+E_k-1)·max(|HZ_k+E_k-1|-r,0)

E_k＝E_k-1+HZ_k-D_k

r = \frac{1}{\sqrt{\max (n, m)}}