CN104063887A - 一种基于Low Rank的动态PET图像重建方法 - Google Patents

Abstract

本发明公开了一种基于Low Rank的动态PET图像重建方法，该方法通过建立重建问题的数学模型，等价问题模型的转化，并基于Low Rank方法重建PET图像；其中，利用Low Rank模型对PET图像进行重建的过程中，等价模型问题在求解过程中采用交替最小化算法。故本发明有效利用Low-Rank算法，改善了计算机在进行动态PET图像重建的过程中产生的结果低分辨率和噪声干扰的问题；与现有重建方法的实验比较表明，本发明能获得较好的重建效果。

Description

一种基于Low Rank的动态PET图像重建方法

技术领域

本发明属于PET成像技术领域，具体涉及一种基于Low Rank(低秩)的动态PET图像重建方法。

背景技术

正电子发射断层成像(Positron emission tomography，PET)是一种基于核物理学、分子生物学的医学影像技术，它能够从分子水平上观察细胞的代谢活动，为早期疾病的检测和预防提供了有效依据。在进行PET扫描时首先需要将由放射性同位核素标记的药物注入人体内，通过血液循环系统，这些物质在人体内各组织器官中将形成一定的分布。由于放射性同位核素不稳定，将发生衰变，衰变过程中所产生的正电子与组织中的电子发生湮灭反应，产生一对方向几乎相反的伽玛光子，经由符合采集系统对这些带有放射性药物分布信息的成对光子进行处理生成投影数据。通过相应的数学方法对投影数据进行反演求解，可重建出放射性物质的空间浓度分布。

动态正电子发射断层扫描测量放射性在生物体组织的空间分布，通过连续数据采集，生理状态成像组织可以从动态PET图像序列来获得，这使得动态PET在临床诊断的强大工具和生物医学研究。活动地图进行重构从正弦图，通过PET扫描仪收集到的原始数据。随后，图像通常被分割为多个ROI(感兴趣区域)代谢参数估计，并进一步病理分析。

传统上，放射性浓度分布重建和感兴趣区域边界分割问题被视为两个顺序步骤。放射性浓度分布重建往往采用统计迭代方法，其中包括著名的MLEM(最大似然期望最大化)、MAP(最大后验)和SAGE(惩罚似然)算法，然后再应用以聚类，形变模型或图割(Graph Cut)为基础的方法来分割PET图像。然而，由于测量数据的复杂性，所得到的放射性浓度分布存在低分辨率和噪声干扰的问题，造成后续分割的困难。

发明内容

针对现有技术所存在的上述技术问题，本发明提供了一种基于Low Rank的动态PET图像重建方法，解决了计算机在进行图像重建的过程中产生的结果低分辨率和噪声干扰的问题。

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

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

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

(3)通过对所述的测量方程引入Low Rank，得到PET的Low Rank模型如下：

X＝L+S Y＝D(L+S)+E

其中：D为系统矩阵，Y为符合计数矩阵，X为PET浓度分布矩阵，E为噪声矩阵，L为低秩的序列共享的背景，S为分割出的活动的生物组织；

(4)对所述的Low Rank模型进行最小二乘正则化，得到对应的最小二乘正则化模型；

(5)对所述的最小二乘正则化模型进行求解，得到对应的增强型拉格朗日函数；

(6)对所述的增强型拉格朗日函数进行最小化求解，同时获得背景L和生物组织S，通过使背景L与生物组织S相叠加得到PET浓度分布矩阵X，并对其进行PET成像从而获得n帧连续的PET图像。

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

y＝Dx+e

其中：y为校正后的符合计数向量，x为PET浓度分布向量，e为噪声向量。

所述的PET浓度分布矩阵X由与符合计数向量对应的n组PET浓度分布向量组合而成。

所述的最小二乘正则化模型的表达式如下：

$\begin{array}{cc}\underset{L,S}{\min }\{{\left|\right|L\left|\right|}_{*}+\mathrm{\λ}{\left|\right|S\left|\right|}_1+\frac{1}{2\mathrm{\μ}}{\left|\right|R\left|\right|}_F^2\}& R=D(L+S)-Y\end{array}$

其中：||||_*为核范数，||||₁为1-范数，||||_F为F-范数，λ和μ均为权重系数。

所述的增强型拉格朗日函数的表达式如下：

$\begin{array}{c}{L}_A(L,S,\mathrm{\γ})={\left|\right|L\left|\right|}_{*}+\mathrm{\λ}{\left|\right|S\left|\right|}_1+\frac{1}{2\mathrm{\μ}}{\left|\right|R\left|\right|}_F^2-\mathrm{tr}\left({\mathrm{\γ}}^{*}(D(L+S)-Y-R)\right)+\frac{\mathrm{\β}}{2}{\left|\right|D(L+S)-Y-R\left|\right|}_F^2\\ R=D(L+S)-Y\end{array}$

其中：L_A(L,S,γ)为关于L、S和γ的增强型拉格朗日函数，γ为拉格朗日乘子，tr()表示求迹运算；||||_*为核范数，||||₁为1-范数，||||_F为F-范数，λ和μ均为权重系数，β为惩罚系数。

所述的步骤(6)中通过以下迭代方程组对增强型拉格朗日函数进行最小化求解：

$R^{k+1}=\frac{\mathrm{\μ\β}}{\mathrm{\μ\β}+1}(D(L^k+S^k)-Y-{\mathrm{\γ}}^k/\mathrm{\β})$

L^k+1＝D_τ/β(L^k-τD(D(L^k+S^k)-b^k)

S^k+1＝S_τ/β(S^k-τD(D(L^k+S^k)-b^k)

γ^k+1＝γ^k-β(D(L^k+1+S^k+1)-Y-R^k+1)

其中：L^k和L^k+1分别为第k次和第k+1次迭代的背景，S^k和S^k+1分别为第k次和第k+1次迭代的生物组织，γ^k和γ^k+1分别为第k次和第k+1次迭代的拉格朗日乘子，D_τ/β()表示奇异值阈值算子，S_τ/β()表示软收缩算子，μ为权重系数，β为惩罚系数，b^k为第k次迭代的中间矩阵，τ为接近变量且ρ(D*D)表示矩阵D*D的谱半径。

所述的中间矩阵b^k的表达式如下：

b^k＝Y+R^k+1+γ^k/β

通过所述的迭代方程组进行迭代计算，则使达到最大迭代次数或迭代收敛后的背景与生物组织相叠加即作为估计得到的PET浓度分布矩阵X；迭代收敛条件如下：

$\frac{{\left|\right|L^{k+1}-L^k\left|\right|}_F}{{\left|\right|L^k\left|\right|}_F}<{10}^{-3}$ 且 $\frac{{\left|\right|S^{k+1}-S^k\left|\right|}_F}{{\left|\right|S^k\left|\right|}_F}<{10}^{-3}$

其中：||||_F为F-范数。

本发明通过建立重建问题的数学模型，等价问题模型的转化，并基于Low Rank方法重建PET图像；其中，利用Low Rank模型对PET图像进行重建的过程中，等价模型问题在求解过程中采用交替最小化算法。故本发明有效利用Low-Rank算法，改善了计算机在进行动态PET图像重建的过程中产生的结果低分辨率和噪声干扰的问题；与现有重建方法的实验比较表明，本发明能获得较好的重建效果。

附图说明

图1为本发明PET图像重建方法的流程示意图。

图2(a)为关于肺部体模的真值图像。

图2(b)为采用EM算法重建肺部体模的PET图像。

图2(c)为采用本发明方法重建肺部体模的PET图像。

图2(d)为采用本发明方法分割出的肺部组织的PET图像。

具体实施方式

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

如图1所示，一种基于Low Rank的动态PET图像重建方法，包括步骤：

正电子发射断层扫描仪探测人体内发出的放射性信号，经过符合和采集系统处理，形成原始投影线，并以正弦图的方式存放于计算机硬盘中；对原始采集到的sinogram，以sinogram和已知的系统矩阵D为输入项，调用相关模块。

S1.根据PET探测的原理建立重建问题的基本模型；

S2.引入Low-Rank来优化问题模型；

S3.初始化，设置初始点γ⁰，L⁰和S⁰；设置β、λ、μ，

S4.设置的初始点开始，按照交替最小化方法计算增强版本的拉格朗日方程中的噪声r^k+1和重建的图像的Low-Rank部分L^k+1以及Sparse部分S^k+1；

S5.更新γ^k+1；准备将更新后的参数返回步骤S4进行循环；

S6.判断是否满足迭代停止条件，不满足该条件则执行步骤S5，满足则迭代停止；进而使L+S＝X得到的PET浓度分布向量X实现PET成像。

为了完成PET图像的重建，PET检测过程的基本模型基于如下方程：

y＝Dx+e

其中，y为符合计数向量，这里测量所得的数据一般为sinogram，x为PET浓度分布向量，e为噪声向量，D为系统矩阵(系统矩阵表征了发射光子被探测器接收的概率，其受探测器结构、探测效率、衰减、死时间等因素的决定)。重建的过程即在采集到y后，通过已知的稀疏的D来计算x。

由于D的稀疏性，不能够直接通过求逆来完成重建，因此考虑转化问题模型。故进一步地，通过Low Rank模型来解决上述问题，模型如下：

X＝L+S Y＝D(L+S)+E

其中，X为PET浓度分布矩阵(由与符合计数向量y对应的多组PET浓度分布向量x组合而成)，L为低秩的序列共享的背景，S为分割出的活动的生物组织。重建的过程即在采集到Y(由多组符合计数向量组合成的符合计数矩阵)后，通过已知的稀疏的D来计算L和S，进一步求得X。

进一步地对Low Rank模型进行最小二乘正则化，将问题优化为：

$\begin{array}{cc}\underset{L,S}{\min }\{{\left|\right|L\left|\right|}_{*}+\mathrm{\&lambda;}{\left|\right|S\left|\right|}_1+\frac{1}{2\mathrm{\&mu;}}{\left|\right|R\left|\right|}_F^2\}& R=D(L+S)-Y\end{array}$

其中，λ为权重系数，||L||_*表示核范数，表示L的奇异值之和；||S||₁表示S的L1范数，||S||₁＝Σ_i,j|S_i,j|；μ为随噪声等级调整的权重参数；||||_F是弗罗贝尼乌斯范数。

进一步地对上式进行求解，得到对应的增强型拉格朗日函数为：

$L_A(L,S,\mathrm{\&gamma;})={\left|\right|L\left|\right|}_{*}+\mathrm{\&lambda;}{\left|\right|S\left|\right|}_1+\frac{1}{2\mathrm{\&mu;}}{\left|\right|R\left|\right|}_F^2-\mathrm{tr}\left({\mathrm{\&gamma;}}^{*}(D(L+S)-Y-R)\right)+\frac{\mathrm{\&beta;}}{2}{\left|\right|D(L+S)-Y-R\left|\right|}_F^2$

其中：L_A(L,S,γ)为关于L、S和γ的增强型拉格朗日函数，γ为拉格朗日乘子，tr() 表示求迹运算，β为惩罚系数。

通过以下迭代方程组对上述拉格朗日目标函数进行最小化求解：

$R^{k+1}=\arg \underset{R}{\min }{L}_A(L^k,S^k,R,{\mathrm{\&gamma;}}^k)$

$L^{k+1}=\arg \underset{L}{\min }{L}_A(L,S^k,R^k,{\mathrm{\&gamma;}}^k)$

$S^{k+1}=\arg \underset{S}{\min }{L}_A(L^{k+1},S,R^{k+1},{\mathrm{\&gamma;}}^k)$

γ^k+1＝γ^k-β(D(L^k+1+S^k+1)-Y-R^k+1)

其中：R^k+1的计算表达式如下：

$R^{k+1}=\frac{\mathrm{\&mu;\&beta;}}{\mathrm{\&mu;\&beta;}+1}(D(L^k+S^k)-Y-{\mathrm{\&gamma;}}^k/\mathrm{\&beta;})$

L^k+1的计算表达式如下：

L^k+1＝D_τ/β(L^k-τD(D(L^k+S^k)-b^k)

S^k+1的计算表达式如下：

S^k+1＝S_τ/β(S^k-τD(D(L^k+S^k)-b^k)

b^k＝Y+R^k+1+γ^k/β

τ为接近变量，ρ(D*D)为矩阵D*D的谱半径；D_τ/β()表示奇异值阈值算子，S_τ/β()表示软收缩算子。

通过上述迭代方程组进行迭代计算，则当达到最大迭代次数或迭代收敛后的背景与生物组织相叠加即作为估计得到的PET浓度分布矩阵X；迭代收敛条件如下：

$\frac{{\left|\right|L^{k+1}-L^k\left|\right|}_F}{{\left|\right|L^k\left|\right|}_F}<{10}^{-3}$ 且 $\frac{{\left|\right|S^{k+1}-S^k\left|\right|}_F}{{\left|\right|S^k\left|\right|}_F}<{10}^{-3}$

以下我们采用脑部数字体模模型实验来验证本实施方式的有效性，该模型包含一些高浓度区域。本实验运行环境为：8G内存，3.40GHz，64位操作系统， CPU为intel酷睿双核。

将本实施方式基于Low-Rank的PET重建方法和传统的EM(最大期望算法)方法重建结果做比较，二者使用相同的观测值Y和相同的系统矩阵D以保证结果的可比性，具体参数设置如下：Y为n×n维采集到的sinogram，令m＝n×n，D为m×m维事先计算好的系统矩阵；这里n＝64，即m＝4096。

针对重建图像质量的验证，采用高维度的原始采集数据64个投影角度，每个角度下射束为64条，即m＝4096，重建图像大小为64×64，即维度n＝4096；初始值设定同上。图2(a)～(c)是真值图像、传统EM方法重建的图像和基于本实施方式重建的图像的比对示意图，可以直观地看出基于本实施方式重建的图像与EM的结果相比，能够恢复出更多的结构，且能够在重建的同时完成图像分割，获得活动组织的图像；图2(d)是分割出来的活动的生物组织图像。

对于相同的数据，分别采用本实施方式和传统的EM方法进行比较，如表1所示；应用本实施方式重建结果在与真值的偏差，方差和均方根误差均小于传统的EM方法，说明对本发明技术方案在提高精确度和降低噪声方面的可行性。

表1

方法	偏差	方差	均方根误差
				EM	0.0599	0.0136	0.1168
本实施方法	0.0361	0.0061	0.0781

Claims (8)

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

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

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

(3)通过对所述的测量方程引入Low Rank，得到PET的Low Rank模型如下：

X＝L+S Y＝D(L+S)+E

其中：D为系统矩阵，Y为符合计数矩阵，X为PET浓度分布矩阵，E为噪声矩阵，L为低秩的序列共享的背景，S为分割出的活动的生物组织；

(4)对所述的Low Rank模型进行最小二乘正则化，得到对应的最小二乘正则化模型；

(5)对所述的最小二乘正则化模型进行求解，得到对应的增强型拉格朗日函数；

(6)对所述的增强型拉格朗日函数进行最小化求解，同时获得背景L和生物组织S，通过使背景L与生物组织S相叠加得到PET浓度分布矩阵X，并对其进行PET成像从而获得n帧连续的PET图像。

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

y＝Dx+e

其中：y为校正后的符合计数向量，x为PET浓度分布向量，e为噪声向量。

3.根据权利要求2所述的动态PET图像重建方法，其特征在于：所述的PET浓度分布矩阵X由与符合计数向量对应的n组PET浓度分布向量组合而成。

4.根据权利要求1所述的动态PET图像重建方法，其特征在于：所述的最小二乘正则化模型的表达式如下：

$\begin{array}{cc}\underset{L,S}{\min }\{{\left|\right|L\left|\right|}_{*}+\mathrm{\&lambda;}{\left|\right|S\left|\right|}_1+\frac{1}{2\mathrm{\&mu;}}{\left|\right|R\left|\right|}_F^2\}& R=D(L+S)-Y\end{array}$

其中：||||_*为核范数，||||₁为1-范数，||||_F为F-范数，λ和μ均为权重系数。

5.根据权利要求1所述的动态PET图像重建方法，其特征在于：所述的增强型拉格朗日函数的表达式如下：

$\begin{array}{c}{L}_A(L,S,\mathrm{\&gamma;})={\left|\right|L\left|\right|}_{*}+\mathrm{\&lambda;}{\left|\right|S\left|\right|}_1+\frac{1}{2\mathrm{\&mu;}}{\left|\right|R\left|\right|}_F^2-\mathrm{tr}\left({\mathrm{\&gamma;}}^{*}(D(L+S)-Y-R)\right)+\frac{\mathrm{\&beta;}}{2}{\left|\right|D(L+S)-Y-R\left|\right|}_F^2\\ R=D(L+S)-Y\end{array}$

其中：L_A(L,S,γ)为关于L、S和γ的增强型拉格朗日函数，γ为拉格朗日乘子，tr()表示求迹运算；||||_*为核范数，||||₁为1-范数，||||_F为F-范数，λ和μ均为权重系数，β为惩罚系数。

6.根据权利要求1所述的动态PET图像重建方法，其特征在于：所述的步骤(6)中通过以下迭代方程组对增强型拉格朗日函数进行最小化求解：

$R^{k+1}=\frac{\mathrm{\&mu;\&beta;}}{\mathrm{\&mu;\&beta;}+1}(D(L^k+S^k)-Y-{\mathrm{\&gamma;}}^k/\mathrm{\&beta;})$

L^k+1＝D_τ/β(L^k-τD(D(L^k+S^k)-b^k)

S^k+1＝S_τ/β(S^k-τD(D(L^k+S^k)-b^k)

γ^k+1＝γ^k-β(D(L^k+1+S^k+1)-Y-R^k+1)

其中：L^k和L^k+1分别为第k次和第k+1次迭代的背景，S^k和S^k+1分别为第k次和第k+1次迭代的生物组织，γ^k和γ^k+1分别为第k次和第k+1次迭代的拉格朗日乘子，D_τ/β()表示奇异值阈值算子，S_τ/β()表示软收缩算子，μ为权重系数，β为惩罚系数，b^k为第k次迭代的中间矩阵，τ为接近变量且ρ(D*D)表示矩阵D*D的谱半径。

7.根据权利要求6所述的动态PET图像重建方法，其特征在于：所述的中间矩阵b^k的表达式如下：

b^k＝Y+R^k+1+γ^k/β。

8.根据权利要求6所述的动态PET图像重建方法，其特征在于：通过所述的迭代方程组进行迭代计算，则使达到最大迭代次数或迭代收敛后的背景与生物组织相叠加即作为估计得到的PET浓度分布矩阵X；迭代收敛条件如下：

$\frac{{\left|\right|L^{k+1}-L^k\left|\right|}_F}{{\left|\right|L^k\left|\right|}_F}<{10}^{-3}$ 且 $\frac{{\left|\right|S^{k+1}-S^k\left|\right|}_F}{{\left|\right|S^k\left|\right|}_F}<{10}^{-3}$

其中：||||_F为F-范数。