技术领域

本发明属于 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 通过以下算式求得：

^C k = ^C k-\ ^{+ X} k ~ ^B k

其中： 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 ^T H = 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 = ^C k-\ ^{+ X} k ~ ^B k

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动态重建的准确率并且大幅提高了图 像目标区域和背景的对比度。