CN109636869B - 基于非局部全变分和低秩约束的动态pet图像重建方法 - Google Patents

Abstract

本发明公开了一种基于非局部全变分和低秩约束的动态PET图像重建方法，该方法利用PET图像的分段平滑特性和时空相关性，同时引入低秩约束和非局部全变分约束，实现了动态PET图像重建，能够移除阶梯效应并保持精细细节，有利于改善早期病灶检测。本发明通过低秩和稀疏这一个约束来优化解的时空相关性，也即是通过一个联合模型来重构出图像的背景成分和细节成分；同时在动态PET图像重建中，为充分考虑图像数据的结构光滑特性，本发明引入了非局部全变分约束，利用图像冗余特性来恢复更多的图像细节并移除阶梯效应。相比于现有技术，本发明可以给出更准确的重建图像以改善病灶检测，并且对噪声具有更好的鲁棒性。

Description

基于非局部全变分和低秩约束的动态PET图像重建方法

技术领域

本发明属于PET成像技术领域，具体涉及一种基于非局部全变分和低秩约束的动态PET图像重建方法。

背景技术

动态正电子发射断层扫描(DPET)能够监测放射性标记示踪剂的体内时空分布，其具有改善早期检测，癌症表征和治疗反应评估的潜力。DPET使用放置在物体周围的检测器系统，以便在一系列时间帧内获得许多不同的角度视图，利用这些投影数据可以重建出放射性示踪剂浓度图，即是动态PET图像重建，其通过时间序列图像能更好的评估示踪剂摄取和代谢的动态过程，在科学研究和临床应用中具有重要应用价值。

PET图像重建是一个病态的逆问题，标准的做法是使用正则项来约束问题的解，从而使得逆问题具有适定性。除了传统的ML-EM等算法，第一类算法是空间平滑约束，其中的一种是最大后验概率(MAP)算法，设计保留区域平滑度和边缘急剧变化的惩罚项一直是PET重建研究的重点；另一种方法是利用全变分约束(TV)来提升PET图像的结构平滑特性，然而基于TV的模型假设每个图像像素总是具有边缘和梯度的扩散方向，这可能会导致阶梯效应。第二类算法是同时利用时间和空间信息来提升动态PET图像重建质量，时间约束通常被用来提升空间解的附加鲁棒性，包括时空样条模型、示踪动力学、小波变换等，然而示踪动力学方法假定所有体素均由相同的模型动力学集很好地建模，这在实践中可能并非如此。另外，小波变换的方法对于选择E样条小波参数向量和合适的小波系数仍留有可选择空间，在利用PET图像进行早期病灶检测问题方面，虽然在重建期间将从CT和MRI图像得到的补充信息合并进了正则约束，但它并没有提供关于器官和病变代谢的信息，因此缺乏足够的信息来指导功能异常区域的重建。一些研究表明，使用没有准确病变轮廓的解剖学边界不会改善病变检测或量化任务。

发明内容

鉴于上述，本发明提供了一种基于非局部全变分和低秩约束的动态PET图像重建方法，该方法利用PET图像的分段平滑特性和时空相关性，同时引入低秩约束和非局部全变分约束，实现了动态PET图像重建，能够移除阶梯效应并保持精细细节，有利于改善早期病灶检测。

一种基于非局部全变分和低秩约束的动态PET图像重建方法，包括如下步骤：

(1)利用探测器对注入有放射性药剂的生物组织进行探测，动态采集得到对应各个时刻的符合计数向量，并将这些符合计数向量组合成符合计数矩阵Y；

(2)使动态的PET图像序列组合成PET浓度分布矩阵X，根据PET成像原理建立PET测量方程；

(3)通过对所述PET测量方程引入低秩约束，得到基于低秩约束的动态PET图像重建模型M1；

(4)通过对每帧图像的低秩部分和稀疏部分进行非局部全变分约束，进一步得到基于非局部全变分约束的动态PET图像重建模型M2；

(5)将动态PET图像重建模型M1和M2相结合得到动态PET重建的目标函数如下：

s.t.L+S＝X

其中：|| ||*表示核范数，|| ||₁表示1-范数，L为低秩部分包含了图像背景中周期变化的放射性浓度，S为稀疏部分包含了具有不同代谢率的非均匀组织的放射性浓度，J_NLTV(L)和J_NLTV(S)分别为低秩部分L和稀疏部分S经非局部全变分后的结果，λ、μ、α_L和α_S均为权重系数，Ψ(Y|X)为关于X和Y的似然函数；

(6)对上述目标函数进行最优化求解后即得到PET浓度分布矩阵X，从而还原出动态的PET图像序列。

进一步地，所述PET测量方程的表达式如下：

Y＝DX+R+S

其中：D为系统矩阵，R和S分别为反映随机事件和散射事件的测量噪声矩阵。

进一步地，所述动态PET图像重建模型M1的表达式如下：

M1＝||L||_*+λ||S||₁+μΨ(Y|X)

s.t. L+S＝X

进一步地，所述动态PET图像重建模型M2的表达式如下：

M2＝α_LJ_NLTV(L)+α_SJ_NLTV(S)+μΨ(Y|X)

s.t. L+S＝X

进一步地，所述似然函数Ψ(Y|X)的表达式如下：

其中：d_ij为系统矩阵D中第i行第j列元素值，y_im为符合计数矩阵Y中第i行第m列元素值，x_jm为PET浓度分布矩阵X中第j行第m列元素值，rim为测量噪声矩阵R中第i行第m列元素值，s_im为测量噪声矩阵S中第i行第m列元素值，i、j和m均为自然数且1≤i≤I，1≤j≤J，1≤m≤M，I为符合计数向量的维度，J为PET浓度分布矩阵X的行数即PET图像的像素点个数，M为PET浓度分布矩阵X的列数即采样时间长度。

进一步地，所述步骤(6)中采用ADMM(Alternating Direction Method ofMultipliers，交替方向乘子)算法对目标函数进行最优化求解；其中，低秩约束问题利用奇异值阈值算法和软收缩算法进行迭代优化求解，似然函数问题利用EM(ExpectationMaximizationAlgorithm，期望最大化)算法进行迭代优化求解，非局部全变分约束问题利用梯度下降法进行迭代优化求解。

不同于现有将时间先验和空间先验作为两个不同约束的方法，本发明通过低秩和稀疏这一个约束来优化解的时空相关性，也即是通过一个联合模型来重构出图像的背景成分和细节成分；同时在动态PET图像重建中，为充分考虑图像数据的结构光滑特性，本发明引入了非局部全变分(NLTV)约束，利用图像冗余特性来恢复更多的图像细节并移除阶梯效应。相比于TV，NLTV整合了图像矩阵的局部和非局部相关性，本发明可以给出更准确的重建图像以改善病灶检测，并且对噪声具有更好的鲁棒性。

本发明将非局部全变分约束引入低秩矩阵恢复分析框架中，以确保PET图像中感兴趣区域(ROIs)的结构平滑性和清晰的边界；图像序列的低秩约束通过组织内的固有平均来消除噪声，同时引入非局部全变分以改善PET图像空间分辨率也是本发明的创新点；低秩和稀疏矩阵分解可以为NLTV约束提供随机噪声成分，随着随机噪声信息的确认，NLTV约束可以提供增强的干净的PET图像，并且反过来帮助低秩矩阵和稀疏矩阵的分解，从而得到一个更加符合真实情况的重建结果。结合本发明在模拟数据实验中的表现，通过与ML-EM算法、LRTV算法(基于低秩和全变分约束)的结果对比，本发明都能获得较准确的重建结果，这对于改善早期病灶检测、评估示踪剂摄取和代谢的动态过程等具有重要的实际应用价值。

附图说明

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

图2(a)为蒙特卡洛模拟Zubal胸腔数据的模板图像。

图2(b)为Hoffman脑部数据的模板图像。

图3(a)为Zubal胸腔数据第8帧的真实图像。

图3(b)为数据计数率为1×10⁷下采用ML-EM方法对蒙特卡洛模拟Zubal胸腔数据重建的第8帧PET图像结果。

图3(c)为数据计数率为1×10⁷下采用LRTV方法对蒙特卡洛模拟Zubal胸腔数据重建的第8帧PET图像结果。

图3(d)为数据计数率为1×10⁷下采用本发明方法对蒙特卡洛模拟Zubal胸腔数据重建的第8帧PET图像结果。

图4(a)为图3(a)中所框出部分放大后的图像结果。

图4(b)为图3(b)中所框出部分放大后的图像结果。

图4(c)为图3(c)中所框出部分放大后的图像结果。

图4(d)为图3(d)中所框出部分放大后的图像结果。

图5(a)为数据计数率为1×10⁷下ROI2每帧Zubal胸腔数据图像结果的偏差折线图。

图5(b)为数据计数率为1×10⁷下ROI2每帧Zubal胸腔数据图像结果的方差折线图。

图6(a)为数据计数率为1×10⁷下采用ML-EM方法对蒙特卡洛模拟Zubal胸腔数据重建的第14帧PET图像结果。

图6(b)为数据计数率为1×10⁷下采用LRTV方法对蒙特卡洛模拟Zubal胸腔数据重建的第14帧PET图像结果。

图6(c)为数据计数率为1×10⁷下采用本发明方法对蒙特卡洛模拟Zubal胸腔数据重建的第14帧PET图像结果。

图7(a)为数据计数率为1×10⁶下采用ML-EM方法对蒙特卡洛模拟Zubal胸腔数据重建的第14帧PET图像结果。

图7(b)为数据计数率为1×10⁶下采用LRTV方法对蒙特卡洛模拟Zubal胸腔数据重建的第14帧PET图像结果。

图7(c)为数据计数率为1×10⁶下采用本发明方法对蒙特卡洛模拟Zubal胸腔数据重建的第14帧PET图像结果。

图8(a)为Hoffman脑部数据第11帧的真实图像。

图8(b)为数据计数率为1×10⁷下采用ML-EM方法对蒙特卡洛模拟Hoffman脑部数据重建的第11帧PET图像结果。

图8(c)为数据计数率为1×10⁷下采用LRTV方法对蒙特卡洛模拟Hoffman脑部数据重建的第11帧PET图像结果。

图8(d)为数据计数率为1×10⁷下采用本发明方法对蒙特卡洛模拟Hoffman脑部数据重建的第11帧PET图像结果。

具体实施方式

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

如图1所示，本发明基于非局部全变分和低秩约束的动态PET图像重建方法，包括如下步骤：

S1.根据动态PET扫描的模式建立测量数据矩阵Y和系统矩阵D。

将动态PET的扫描过程按照需要划分出一定数量的时间帧，每个时间帧内探测器采集得到的符合计数向量的可以按照时间的顺序构建出一个动态PET的测量数据矩阵Y；并统计每一像素点处出射的光子被各探测器接收到的概率，从而得到系统矩阵D。

S2.建立动态PET成像模型。

对于动态PET成像，它需要对示踪剂的时间信息进行一系列连续的时间帧采样。对于每一个独立的时间帧，投影数据y代表了在该时间段内被每一个探测器所捕获的符合事件之和，即y＝{y_i,i＝1,2,…,I}，其中I是探测器的总个数；相应的放射性浓度图像被记录为x＝{x_j,j＝1,2,…,J}，其中J是像素的总个数；测量数据y与未知的浓度图x之间的关系为：

由于泊松假设的独立性，我们有y的似然函数：

为便于求解，我们对该似然函数取对数，并最小化似然函数的负对数：

对于动态PET，我们能结合所有的时间帧信息来得到一个数据矩阵。对于第m帧投影数据，向量y_m代表数据矩阵Y的第m列，

y_im代表第m帧被第i个探测器所捕获的投影数据；同样的，

x_jm代表第m帧第j个像素点的放射性浓度；我们对每一帧的负对数似然函数求和可得到：

其中：y_im代表第m帧投影数据y_m的第i个探测器的值，

是第m帧

的第i个条目。

S3.低秩和稀疏约束。

对PET图像进行低秩和稀疏分解，其中L成分包括背景中放射性强度的周期性变化，S成分指那些具有不同代谢率的非均匀的组织，基于此建立以下分解模型：

将目标函数(5)放缩为一个凸优化问题：

其中：||L||_*代表矩阵L的核范数，即L的奇异值之和。||S||₁代表S的l₁范数。

S4.非局部全变分约束。

令P＝[p₁,p₂,…,p_m,…,p_M]代表一个M帧的图像序列，其中

是第m帧图像向量。将每一列

转换至其矩阵形式

其中U×V＝J，被转换的

现在在空间

中。我们依次对每一帧图像矩阵

进行非局部全变分约束：

其中：u，v为Ω空间中的像素点，其为自然数，w(u,v)是非负的对称非局部权重函数，G_δ是标准差为δ的高斯核，h是一个滤波参数。

S5.动态PET图像的LRNLTV重建。

将非局部全变分约束分别作用于分解出的L和S，结合之前的低秩和稀疏约束，我们有以下目标函数：

s.t. L+S＝X

其中：λ，μ，α_L和α_S均为权重系数。

S6.基于增广拉格朗日乘子法的优化算法。

引入辅助变量P和Q，写出目标函数(8)的增广拉格朗日函数：

其中：Z，Z_L，Z_S为拉格朗日乘子，β、β_L、β_S是惩罚参数；将该问题分解为五个子问题进行求解，根据每个子问题的性质，将其分为三类。

6.1求解L，S子问题：

对于L子问题，令其它变量固定，仅提取出与L相关的项，并进行整理配方，有：

利用奇异值阈值法求解式(10)，则有L的迭代更新式为：

其中：A_ε(Γ)＝UB_ε(s)V^T，UsV^T是Γ的奇异值分解，其中的软收缩算法B_ε(s)为：

同样对于S子问题，化简公式后有：

利用软收缩算法求解式(12)，则有S的迭代更新式为：

6.2求解X子问题：

首先写出对于隐藏变量w的负似然函数Ψ(w|X)，其中

表示第m帧的符合线i处检测到的来自像素j的发射光子：

接下来利用EM算法来优化X。

E步：取w的条件期望，

并将其插入Ψ(w|X)；然后我们有替代函数φ(X；X^k)：

M步：计算φ(X；X^k)对于x_jm的导数，并令其等于0；经化简后发现x_jm的解是以下二次方程的根：

式(15)是凸函数，我们通过取较大的跟来更新x_jm：

s.t. a_jm＝β，

其中：[ ]_jm表示矩阵的第j个条目。

6.3求解P，Q子问题：

对于P子问题，同样让其它变量固定，仅提取出与L相关的项，并进行整理配方，有：

我们使用梯度下降法求解式(17)，首先我们写出式(7)的Euler-Lagrange：

对于固定的u，有：

我们将L矩阵的每一列

转换至其矩阵形式

Z_L的每一列也同样被转换至其矩阵形式，其中U×V＝J，并且L和Z_L的每一帧被依次处理；则对于图像的每一帧，有

的迭代更新式：

其中，

是梯度下降的步进。在所有帧都被处理之后，将矩阵

转换回至向量形式

即可重新得到L，用同样的方式恢复P和Z_L，则P子问题被解决。

对于Q子问题，由于其具有与P子问题相同的形式，故用同样的方法来解决它；其中，拉格朗日乘子照常更新。

以下是我们对蒙特卡洛模拟的Zubal胸腔和Hoffman脑部模板数据进行实验从而验证本发明系统重建结果的准确性。图2(a)和图2(b)分别为实验所用的Zubal胸腔和Hoffman脑部数据的模板示意图，将不同的区域分为三个感兴趣的区域(其中ROI4是背景区域)。实验运行环境为：8G内存，3.40GHz，64位操作系统，CPU为intel i7-3770；所模拟的PET扫描仪型号为Hamamatsu SHR-22000，设定的放射性核素及药物为¹⁸F-FDG，设置sinogram为64个投影角度在每个角度下64条射束采集到的数据结果，系统矩阵D的大小为4096×4096。在本次试验中，对1×10⁶、1×10⁷两种不同的计数率下的投影数据进行实验。

将本发明重建框架的PET图像结果分别与ML-EM和LRTV两种重建方法的图像结果进行比较，二者使用相同的测量数据矩阵Y和系统矩阵D以控制结果的可比性。从图3(a)～图3(d)中可以看出，本发明重建框架的图像结果在区域内的噪声明显小于另两种方法所重建的图像，在保证边缘对比的情况下功能区域内更加平滑。图4(a)～图4(d)分别是图3(a)～图3(d)所框出部分放大的图像结果，明显可以看出，本发明所重建的图像结果具有更清晰的ROI边界信息，对于早期的病灶检测具有明显的改善作用。图5(a)～图5(b)展示了Zubal胸腔数据的ROI2的每一帧的量化误差，进一步说明了本发明重建结果的准确性。图6(a)～图6(c)和图7(a)～图7(c)分别是计数率为1×10⁷和1×10⁶下三种重建方法对Zubal胸腔数据第14帧的重建结果，验证了本发明的重建结果对噪声具有鲁棒性，表1是其进一步的量化结果分析。图8(a)～图8(d)是在Hoffman脑部数据上进行的重建实验，表明了本发明的重建结果对不同的模板数据具有鲁棒性。

表1

上述对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。

Claims (6)

1.一种基于非局部全变分和低秩约束的动态PET图像重建方法，包括如下步骤：

(1)利用探测器对注入有放射性药剂的生物组织进行探测，动态采集得到对应各个时刻的符合计数向量，并将这些符合计数向量组合成符合计数矩阵Y；

(2)使动态的PET图像序列组合成PET浓度分布矩阵X，根据PET成像原理建立PET测量方程；

(3)通过对所述PET测量方程引入低秩约束，得到基于低秩约束的动态PET图像重建模型M1；

(4)通过对每帧图像的低秩部分和稀疏部分进行非局部全变分约束，进一步得到基于非局部全变分约束的动态PET图像重建模型M2；

(5)将动态PET图像重建模型M1和M2相结合得到动态PET重建的目标函数如下：

s.t.L+S＝X

其中：|| ||_*表示核范数，|| ||₁表示1-范数，L为低秩部分包含了图像背景中周期变化的放射性浓度，S为稀疏部分包含了具有不同代谢率的非均匀组织的放射性浓度，J_NLTV(L)和J_NLTV(S)分别为低秩部分L和稀疏部分S经非局部全变分后的结果，λ、μ、α_L和α_S均为权重系数，Ψ(Y|X)为关于X和Y的似然函数；

(6)对上述目标函数进行最优化求解后即得到PET浓度分布矩阵X，从而还原出动态的PET图像序列。

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

Y＝DX+R+S

其中：D为系统矩阵，R和S分别为反映随机事件和散射事件的测量噪声矩阵。

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

M1＝||L||_*+λ||S||₁+μΨ(Y|X)

s.t.L+S＝X。

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

M2＝α_LJ_NLTV(L)+α_SJ_NLTV(S)+μΨ(Y|X)

s.t.L+S＝X。

5.根据权利要求2所述的动态PET图像重建方法，其特征在于：所述似然函数Ψ(Y|X)的表达式如下：

其中：d_ij为系统矩阵D中第i行第j列元素值，y_im为符合计数矩阵Y中第i行第m列元素值，x_jm为PET浓度分布矩阵X中第j行第m列元素值，r_im为测量噪声矩阵R中第i行第m列元素值，s_im为测量噪声矩阵S中第i行第m列元素值，i、j和m均为自然数且1≤i≤I，1≤j≤J，1≤m≤M，I为符合计数向量的维度，J为PET浓度分布矩阵X的行数即PET图像的像素点个数，M为PET浓度分布矩阵X的列数即采样时间长度。

6.根据权利要求1所述的动态PET图像重建方法，其特征在于：所述步骤(6)中采用ADMM算法对目标函数进行最优化求解；其中，低秩约束问题利用奇异值阈值算法和软收缩算法进行迭代优化求解，似然函数问题利用EM算法进行迭代优化求解，非局部全变分约束问题利用梯度下降法进行迭代优化求解。

Images