WO2014172927A1

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

Info

Publication number: WO2014172927A1
Application number: PCT/CN2013/075593
Authority: WO
Inventors: 刘华锋; 于行健
Original assignee: 浙江大学
Priority date: 2013-04-23
Filing date: 2013-05-14
Publication date: 2014-10-30
Also published as: US9576379B2; CN103279964A; US20160078646A1; CN103279964B

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图像的时间关联性表现出来，加之两种方法都无法将采集时的充满噪声的背景和目标区域分离，所以不能 4艮好的应对噪声的干扰。发明内容

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

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

( 1 )利用探测器对注入有放射性物质的生物组织进行探测，动态采集得到

PET的 n组符合计数向量，并对所述的符合计数向量进行校正，进而构建 PET 的符合计数矩阵； n为大于 1的自然数；

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

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

(X_k+l, Z_k+l) = arg min + Ζ) - Υ - ΑΛ

+ r \\HZ - D_k + E_k 其中： 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浓度动态数据。

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

其中： X_k为第 k次迭代后的 PET浓度分布矩阵， Z_k为第 k次迭代后的 PET背景浓度分布矩阵， p为预设的收敛阈值，表示 F范数。所述的第 k次迭代后的第一修正系数矩阵 A_k通过以下算式求得：

A_k =A_k__{l +} G(X_{k +}Z_k)-Y

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

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

B_h =U_h - diag(max(i_k - 1, 0)) · V_k

W_k =U_k-diag(i_k)-V_k

其中： W_k为第 k次迭代后的待分解矩阵且

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

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

^Ck = ^Ck-\ ^{+ X}k ~^Bk

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

所述的第 k次迭代后的第四修正系数矩阵 D_k通过以下算式求得： D_k = sgn(HZ_k + E_k__x) .

+ E_k__x | - r, 0) 其中： Z_k为第 k次迭代后的 PET背景浓度分布矩阵，为第 k-1次迭代后的第五修正系数矩阵。 E_k =E_k__{l +}HZ_k-D_k

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

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

1

jmax(n,m)

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

本发明中，若符合计数向量的维度为 m, PET浓度分布矩阵的维度为 pxn, 则系统矩阵 G的维度为 mxp, 其表征了发射光子被探测器接收的概率，其受探测器结构、探测效率、衰减、死时间等因素的影响；符合计数矩阵为由 n组符合计数向量组成的维度为 mxn的矩阵。修正矩阵 H为 mxp的矩阵且为多级紧凑型小波分解中算子，其满 J H^TH = I , I为单位矩阵。一种基于 PRCA的 PET图像动态重建系统，包括探测器和与探测器相连的所述的探测器用于对注入有放射性物质的生物组织进行探测 , 动态采集得到 PET的 n组符合计数向量， n为大于 1的自然数；

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

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

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

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

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

A_k =A_k__{l +} G(X_k+Z_k)-Y

B_h =U_h- diag(max(i_k - 1, 0)) · V_k W_k = U_k - diag(i_k) - V_k

Ck = ^Ck-\ ^{+ X}k ~ ^Bk

D_k = sgn(HZ_k + E_k__x) . max( HZ_k + E_k__x ― r, 0)

E_k = E_k__{l +} HZ_k - D_t

其中：其中： 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的随机化矩阵， Aw为第 k- 1次迭代后的第一修正系数矩阵， X_k为第 k次迭代后的 PET浓度分布矩阵， Z_k为第 k次迭代后的 PET背景浓度分布矩阵， W_k为第 k次迭代后的待分解矩阵且

, 为第 k-1次迭代后的第三修正系数矩阵， U_k、 diag(i_k)和 V_k分别为待分解矩阵 W_k经奇异值分解后得到的时间参数相关矩阵、奇异值对角矩阵和中介系数矩阵， i_k为奇异值对角矩阵 diag(i_k)中的对角线元素，为第 k-1次迭代后的第五修正系数矩阵， m为符合计数向量的维度， k为自然数；迭代收敛后的 PET浓度分布矩阵即为 PET浓度动态数据。本发明的有益技术效果在于：通过将采集到的不同帧的数据看作一个整体来进行重建，充分利用了 PET数据在时间上的相关性，使所得结果更能体现出动态 PET能表现出目标区域时间变化的特点；其次，本发明使用了背景和目标区域的方法，减小了背景对目标区域的干扰，加之在重建中添加了时间和空间的修正，使得重建结果的准确度更高，目标区域与背景的对比度提高，使的重建效果比传统的 FBP和 ML-EM算法更加优秀，更具有医学价值。附图说明

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

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

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

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

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

图 4(b)为 Levell等级的模拟图像采用 ML-EM重建的浓度结果示意图。图 4(c)为 Levell等级的模拟图像采用本发明方法重建的浓度结果示意图。图 5(a)为 Level2等级模拟图像的原始浓度结果示意图。

图 5(b)为 Level2等级的模拟图像采用 ML-EM重建的浓度结果示意图。图 5(c)为 Level2等级的模拟图像采用本发明方法重建的浓度结果示意图。图 6为本发明使用的 Phantom体模示意图。

图 7 为 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, 其为 mxn的矩阵，符合计数矩阵 Y中每一列对应着一组符合计数向量。

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

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

y = Gx

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

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

2 2

, Z_k+l ) = arg min G(X + Z) - 7 - A_k + X-B_{k +} C_k +r HZ-D_{k +}E_i k

A_k =A_k__{l +} G(X_k+Z_k)-Y

B_h =U_h- diag(max(i_k - 1, 0)) · V_k

W_k =U_k-diag(i_k)-V_k

D_k = sgn(HZ_k + E^) · max( HZ_k + E_k__x

E_k =E_k__{l +}HZ_k-D_k jmsLx(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的随机化矩阵， Aw为第 k-1 次迭代后的第一修正系数矩阵， X_k为第 k次迭代后的 PET浓度分布矩阵， Z_k为第 k次迭代后的 PET背景浓度分布矩阵， W_k为第 k次迭代后的待分解矩阵且 W_k=X_k+Cw , 为第 k- 1次迭代后的第三修正系数矩阵， U_k、 diag(i_k)和 V_k分别为待分解矩阵 W_k经奇异值分解后得到的时间参数相关矩阵、奇异值对角矩阵和中介系数矩阵， i_k为奇异值对角矩阵 diag(i_k)中的对角线元素，为第 k-1次迭代后的第五修正系数矩阵。

本实施方式中， PET浓度分布矩阵 X的维度为 ρχη, 该矩阵中每一列对应着一组 PET浓度分布向量 X； X。、 Z。、 A。、 B。、 C。、 D。和 E。均为初始化给定的零矩阵，修正矩阵 H为 mxp 的矩阵且为多级紧凑型小波分解中算子，其满足 Η^ΤΗ = Ι , I为单位矩阵。

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

其中： p为预设的收敛阈值， | |_f表示 F范数；本实施方式中 ρ=1(Τ³; 最大迭代次数预设为 50次。

PET成像模块用于对 PET浓度动态数据进行分解，得到 n帧连续的 PET图像；由于迭代收敛后的 PET浓度分布矩阵为 pxn的矩阵，该矩阵中每一列即对应一组 PET浓度分布向量（即一帧 PET图像数据），将每一列向量元素进行重组即可还原得到 PET图像。以下我们对本实施方式进行了两组验证实验，第一组实验采用的是蒙特卡罗模拟，第二组采用的是六个圆柱体的 phantom体模进行模拟。

在第一组试验中，我们采用的是 180度内 48个采样角，所有采样结果被分为 20帧，我们生成了两组不同噪声和计算率等级的模拟图像，如图 3所示；其中 Levell对应的是 0.53%的随机噪声和 0.03%的散射噪声，其每帧计数率约为 7xl0⁶; 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的测量方程；

( 3 )通过对所述的测量方程在时间和空间上施加^ ίι爹正约束得到以下迭代方程；根据所述的符合计数矩阵通过以下迭代方程估计出 PET浓度动态数据； d , ) = arg min ||G( + Z) - 7 - A, + X - B_k + C_k + r HZ - D_{k +} E_k χ,ζ

其中： 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图像。

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

y = Gx

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

Y - G(X_{k +} Z_k)

P

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

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

A_k =A_k__{l +} G(X_{k +}Z_k)-Y

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

^Ck = ^Ck-\ ^{+ X}k ~^Bk

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

B_h =U_h- diag(max(i_k - 1, 0)) · V_k

W_k =U_k-diag(i_k)-V_k

其中： W_k为第 k次迭代后的待分解矩阵且

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

E_k =E_k__{l +}HZ_k-D_k

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

D_k = sgn(HZ_k + E_k__x) .

+ E_k__x | - r, 0) 其中： Z_k为第 k次迭代后的 PET背景浓度分布矩阵，为第 k-1次迭代后的第五修正系数矩阵。

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

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

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

A_k =A_k__{l +} G(X_k+Z_k)-Y

B_h =U_h- diag(max(i_k - 1, 0)) · V_k

W_k =U_k-diag(i_k)-V_k

D_k = sgn(HZ_k + E^) · max( HZ_k + E_k__x

E_k =E_k__{l +}HZ_k-D_k jmsLx(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的随机化矩阵， Aw为第 k- 1次迭代后的第一修正系数矩阵， X_k为第 k次迭代后的 PET浓度分布矩阵， Z_k为第 k次迭代后的 PET背景浓度分布矩阵， W_k为第 k次迭代后的待分解矩阵且 W_k=X_k+Cw, 为第 k-1次迭代后的第三修正系数矩阵， U_k、 diag(i_k)和 V_k分别为待分解矩阵 W_k经奇异值分解后得到的时间参数相关矩阵、奇异值对角矩阵和中介系数矩阵， i_k为奇异值对角矩阵 diag(i_k)中的对角线元素，为第 k-1次迭代后的第五修正系数矩阵， m为符合计数向量的维度， k为自然数；迭代收敛后的 PET浓度分布矩阵即为 PET浓度动态数据。