本发明属于PET成像
技术领域：
，具体涉及一种基于结构字典和动力学参数字典联合稀疏约束的动态PET图像重建方法。
背景技术：
：正电子发射断层成像(PositronEmissionTomography，PET)是一种先进的核医学成像技术。PET成像技术是在生物活体内注射一种能直接或间接反映生物新陈代谢过程的放射性核素，通过PET设备中探测器环接收湮灭产生的光子对，进而计算确定正电子湮灭(发射)的位置，最后就可以由放射性元素的浓度分布反应出活体内的生理过程，从而达到诊断和分析的目的。相对于只针对静态时间窗内PET的图像重建分析，动态PET还能针对生物体的组织或器官的真实代谢水平进行定量的分析。从分子水平上观察细胞的代谢活动，为疾病的早期诊断和预防提供了有效依据；其中，尤以对癌症的发现、诊断等有着重要的应用。PET图像重建大致可以分为解析重建方法和迭代重建方法两类。解析法利用测量数据与原始图像之间的Radon变换关系，从而通过傅里叶变换重建原始图像。考虑到PET系统的径向采样模式，导致原始图像中心部分过采样，边缘部分则欠采样。因此要先对投影数据进行滤波，即增强欠采样区域，削弱过采样区域，从而在整个视场内得到相对均匀的分布。滤波反投影法(FBP)即是基于此思想的重建方法。然而由于解析法未考虑到影响PET成像质量的噪声特性，故在计数率较低以及噪声较大的情况下重建结果较差。与解析法相比，迭代法假设PET探测数据及噪声满足某种统计特性，从而能更好地描述PET系统的探测过程。如极大似然估计的期望最大法(MLEM)假设测量数据满足泊松分布，通过最大化泊松似然函数求解原始图像。然而MLEM的求解是病态，即不满足解的稳定性和唯一性。常见的解决方法是添加惩罚项，如全变分(TV)约束假设相邻像素间放射性核素浓度差异很小。利用这样的先验知识，能够有效地去除噪声和保留边界，提高成像质量。然而这种方法可能产生图像细节平滑的问题，如何权衡图像细节和信噪比的要求则较难处理。专利(一种基于结构字典约束的PET图像重建方法，申请号201510665311.4)提出了一种结构字典的方法，将PET图像分解为特定尺寸的图块，并认为这些图块可以由一个过饱和的字典稀疏表达出来。结构字典可以从组织的CT/MRI图像学习得到。由于CT/MRI图像和PET图像有着相似的边界结构，因此对于PET图像的边界区域，可以通过CT/MRI图像学习的特征元素线性组合得到；对于平滑区域，结构字典中的一个元素就足以表示它们。基于结构字典的稀疏表达方法，不仅能得到较好的图像细节如边界特性，也能有效抑制平滑区域的噪声，提高成像质量。然而这种算法一般都只针对同一帧的PET图像信息，忽略了其在时间上提供的先验知识，故在一定程度上增加了结果中的噪声。对于动态PET成像，一种简单的做法是对若干个时间点进行静态PET成像，再对时间进行插值，从而得到每个像素点的时间序列曲线。然而这种方法没有考虑到PET成像的动态响应，并且很容易受噪声的影响，因此基于动力学模型的直接动态PET成像方法成为研究的焦点。动力学模型考虑到放射性核素衰减及扩散满足微分关系，故每个房室的放射性核素浓度时间响应可以用一系列指数函数线性组合得到。将这样的先验知识加入重建算法中，从而能够在一定程度上抑制噪声，更好地估计放射性核素的浓度随时间变化关系。因此，如何充分利用动态PET成像的先验信息，获得准确的重建结果是研究动态PET成像的一个非常关键的问题。技术实现要素：基于上述，本发明提供了一种基于结构字典和动力学参数字典联合稀疏约束的动态PET图像重建方法，能够获得高质量的动态PET重建图像序列。一种基于结构字典和动力学参数字典联合稀疏约束的动态PET图像重建方法，包括如下步骤：(1)利用探测器对注入有放射性药剂的生物组织进行探测，动态采集得到nk组符合计数向量，进而组建PET的符合计数矩阵Y，nk为大于1的自然数；(2)通过使动态PET图像序列组合成待估计的PET浓度分布矩阵X，根据PET测量方程Y＝GX+R+S，建立PET的泊松分布模型L(X)：L(X)=ΣknkΣiniy‾ik-yiklogy‾ik]]>s.t.y‾ik=Σjnjgijxjk+rik+sik]]>其中：gij为系统矩阵G中第i行第j列元素值即PET图像第j个像素点发射的光子被第i对探测器所接收的概率，rik和sik分别为随机噪声矩阵R和散射噪声矩阵S中第i行第k列元素值，xjk为PET浓度分布矩阵X中第j行第k列元素值，yik为符合计数矩阵Y中第i行第k列元素值，i、j和k均为自然数且1≤i≤ni，1≤j≤nj，1≤k≤nk，ni为符合计数向量的维度，nj为PET图像的总像素个数；(3)利用生物组织的CT图像、MRI图像或高质量PET图像作为结构字典Ds的训练样本，构建PET的结构字典稀疏惩罚项Ss(X,α)：Ss(X,α)=ΣknkΣpnp||EpXk-Dsαpk||22+μs||αpk||0]]>其中：Ep为分割算子，Xk为第k帧PET图像，EpXk表示PET图像Xk中的第p个n×n维子矩阵，μs为权重系数，np＝(m-n+1)2，m为PET图像的维度，n为预设的子矩阵维度，α为待估计的结构字典稀疏编码矩阵，αpk为结构字典稀疏编码矩阵α中对应EpXk的稀疏编码，‖αpk‖0表示稀疏编码αpk中非零元素的个数，‖‖2表示2范数，p为自然数且1≤p≤np；(4)利用先验示踪动力学参数信息得到动力学参数字典Dc，进而构建PET的动力学参数字典稀疏惩罚项Sc(X,β)：Sc(X,β)=Σjnj||Xj-Dcβj||22+μc||βj||1]]>其中：Xj为PET图像第j个像素点的时间浓度序列即PET浓度分布矩阵X中的第j行横向量，μc为权重系数，β为待估计的动力学参数字典稀疏编码矩阵，βj为动力学参数字典稀疏编码矩阵β中对应Xj的稀疏编码，‖‖1表示1范数；(5)将上述泊松分布模型L(X)、结构字典稀疏惩罚项Ss(X,α)以及动力学参数字典稀疏惩罚项Sc(X,β)相加得到动态PET重建模型F(X,α,β)，并根据下式对其进行最优化求解后即得到PET浓度分布矩阵X，进而重建获得动态PET图像；minX,α,βF(X,α,β)=minX,α,β{L(X)+psSs(X,α)+pkSc(X,β)}]]>其中：ps和pk分别对应为结构字典稀疏惩罚项Ss(X,α)和动力学参数字典稀疏惩罚项Sc(X,β)的权重系数。所述的结构字典Ds通过奇异值分解算法(K-SVD)优化以下目标函数得到；minDs,αΣpnp||EpXCT-Dsαp||22+μs||αp||0]]>其中：XCT为相同生物组织的CT图像，EpXCT表示CT图像XCT中的第p个n×n维子矩阵，αp为结构字典稀疏编码矩阵α中对应EpXCT的稀疏编码，‖αp‖0表示稀疏编码αp中非零元素的个数。所述的动力学参数字典Dc的表达式如下：Dc=ψ10ψ11...ψ1nm............ψnk0ψnk1...ψnknm]]>ψk0=1tke-tks∫tkstkeCI(t)dt]]>ψkm=1tke-tks∫tkstke∫0te-θm(t-τ)CI(τ)dτdt]]>其中：CI(t)和CI(τ)分别为t时刻和τ时刻放射性药剂在血浆中的浓度值，和分别为关于第k组符合计数向量采集的开始时间和结束时间，θm对应为第m个房室组织指数函数的系数，m为自然数且1≤m≤nm，nm为大于1的自然数；的取值为在区间[θmin，θmax]内按指数级间隔进行选取，θmin和θmax分别为系数的上下限阈值，t和τ均表示时间。所述的步骤(5)对F(X,α,β)进行最优化求解过程中对于待估计量X、α和β，采用固定其中两者优化第三者的交替迭代算法进行求解；具体采用正交匹配追踪(OMP)算法优化结构字典稀疏编码矩阵α，采用基追踪降噪(BasisPursuitDenoising)算法优化动力学参数字典稀疏编码矩阵β，采用期望最大化(EM)算法优化PET浓度分布矩阵X。本发明动态PET图像重建方法将极大似然估计的期望最大法与结构字典和动力学参数字典方法相结合，既考虑每一帧图像的空间约束即每个图块可以由CT图像预训练的结构字典稀疏表达，又利用了每个像素点核素浓度随时间变化的微分模型，从而能够有效地抑制噪声，得到较好的动态PET成像结果，且与其他单独重建动态PET图像的算法相比，本发明也都能获得较好的重建结果。附图说明图1为本发明重建方法的步骤流程示意图。图2为Huffman大脑phantom及其感兴趣区域的示意图。图3(a)为本发明的第3帧重建结果示意图。图3(b)为本发明的第3帧重建误差示意图。图3(c)为本发明的第19帧重建结果示意图。图3(d)为本发明的第19帧重建误差示意图。图4(a)为ML-EM方法的第3帧重建结果示意图。图4(b)为ML-EM方法的第3帧重建误差示意图。图4(c)为ML-EM方法的第19帧重建结果示意图。图4(d)为ML-EM方法的第19帧重建误差示意图。图5(a)为结构字典单独约束方法的第3帧重建结果示意图。图5(b)为结构字典单独约束方法的第3帧重建误差示意图。图5(c)为结构字典单独约束方法的第19帧重建结果示意图。图5(d)为结构字典单独约束方法的第19帧重建误差示意图。图6(a)为动力学参数字典单独约束方法的第3帧重建结果示意图。图6(b)为动力学参数字典单独约束方法的第3帧重建误差示意图。图6(c)为动力学参数字典单独约束方法的第19帧重建结果示意图。图6(d)为动力学参数字典单独约束方法的第19帧重建误差示意图。具体实施方式为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。如图1所示，本发明基于结构字典和动力学参数字典稀疏约束的动态PET图像重建方法，包括如下步骤：S1.根据动态PET探测的原理整理测量数据矩阵Y和相应的系统矩阵G；S2.建立结构字典Ds和动力学参数字典Dc：结构字典Ds的训练方法为：a)取相同组织的CT或MRI图像XCT，将其分解为p份特定大小(如7×7)的图块EpXCT，其中Ep为图块分解矩阵。b)构建结构字典稀疏编码的目标函数：minDs,αpΣpnp||EpXCT-Dsαp||22+μs||αp||0---(1)]]>c)使用奇异值分解算法K-SVD优化目标函数，得到结构字典Ds。动力学参数字典Dc的构建方法如下：Dc=D01D11...Dnm1............D0nkD1nk...DnmnkD0k=1tke-tks∫tkstkeCI(t)dtk=1,...,nkDmk=1tke-tks∫tkstke∫0te-θm(t-τ)CI(τ)dτdtm=1,...,nm---(2)]]>其中：CI(t)和CI(τ)分别为t时刻和τ时刻放射性药剂在血浆中的浓度值，和分别为关于第k组符合计数向量采集的开始时间和结束时间，θm对应为第m个房室组织指数函数的系数，θm的取值可在区间[θmin,θmax]内按等对数间隔进行选取，θmin可选取核素衰减常数(如θmin＝λs-1)，θmax可选取一个较大的值(如θmax＝0.1s-1)。S3.初始化算法参数，设置结构字典和动力学参数字典系数惩罚项权重系数ps和pk。根据实验测试得到，ps和pk取值在0.5附近重建结果较好，并且当噪声增加的时候，算法不能准确地从字典中确定出正确的时间放射性曲线，此时需要适当增大ps和pk，增强稀疏约束，抑制噪声。S4.动力学参数字典稀疏编码β优化子问题：固定X和α，β子问题求解动力学字典稀疏编码，其优化函数为：minβΣjnj||Xj:t-Dkβj||22+μk||βj||1---(3)]]>其中：向量为像素点j在不同帧浓度的当前估计值。由于βj相互独立，故式(3)等价于也可以写作：βjt+1=argminβ||βj||1s.t.||Xj:t-Dkβj||22≤δ---(4)]]>其中：δ为动力学字典编码误差，可以用基追踪降噪(Basispursuitdenoising)求解。S5.结构字典稀疏编码α子问题：固定X和β，α子问题求解结构字典稀疏编码的优化函数可写作：minαΣknkΣpnp||EpX:kt-Dsαpk||22+μs||αpk||0---(5)]]>其中：为第k帧PET图像的当前估计值。可以看到此时αpk是相互独立的，因此式(5)等价于也可以写作：αpkt+1=argminαpk||αpk||0s.t.||EpX:kt-Dsαpk||22≤ϵ---(6)]]>其中：ε为结构字典编码误差，可以用正交匹配追踪算法(OMP)求解。S6.动态PET图像序列X子问题固定α和β，X子问题可写作：X=minXΣknkΣiniy‾ik-yiklogy‾ik+psΣknkΣpnp||EpX:k-Dsαpkt+1||22+pkΣjnj||Xj:-Dkβjt+1||22s.t.y‾ik=Σjnjgijxjk+rik+sik---(7)]]>注意到式(7)是一个经典的带有二次惩罚项的最大后验(MAP)问题，我们可以采用期望最大化(EM)来求解。根据隐藏变量公式，X子问题(7)等价于最小化下式：X=argminXΨX(C,X)=argminXΣknkΣjnjΣini(gijxjk-cijklog(gijxjk))+psΣknkΣpnp||EpX:k-Dsαpkt+1||22+pkΣjnj||Xj:-Dkβjt+1||22---(8)]]>其中：cijk为隐藏变量，表示第k帧像素点j处发射的光子对被探测器对i所接收到的事件数目。观察式(8)，如果我们能估计出隐藏变量cijk,则最小化ΨX(C,X)就可以通过导数取零来求解。因此EM算法分为两步：估计隐藏变量(E-步骤)和最小化目标函数(M-步骤)。a)E-步骤使用X的当前估计值Xt(第t次迭代结果)和观测数据Y来估计隐藏变量C其中：和分别为第k帧探测器对i的随机噪声和散射噪声的估计值。将代入式(8)中可得到优化方程b)M-步骤M-步骤则最小化式(10)来求解新的X估计值Xt+1，一个显而易见的方法是将式(10)对xjk求导取零。然而考虑到第二项对xjk求导后含有其他未知变量如xj+1,k，直接求导取零必须求解一个多元方程组，这不是件容易的事情。因此考虑使用近似的方法使用估计其导数。首先，我们先将重写为：||EpX:k-Dsαpkt+1||22=Σln([EpX:k]l-[Dsαpkt+1]l)2---(11)]]>[EpX:k]l=Σjnjηp,lj(ep,ljηp,lj(xjk-xjkt)+[EpX:kt]l)s.t.ηp,lj=ep,ljep,ljΣjnjep,lj---(12)]]>其中：[EpX:k]l为向量EpX:k第l个元素，ep,lj为矩阵Ep的第l,j个元素，为第k帧图像的当前估计值。由于相对于[EpX:k]l是凸函数，则有：([EpX:k]l-[Dsαpkt+1]l)2≤Σjnjηp,lj(ep,ljηp,lj(xjk-xjkt)+[EpX:kt]l-[Dsαpkt+1]l)2---(13)]]>将不等式(13)带入式(10)中，可以得到可分离的替代函数Φ(X；Xt)：将式(14)对xjk求导取零可得：可以看出xjk为多项式的根，其中：由于Φ(X；Xt)是一个严格的凸函数，因此为式(15)的唯一正根：xjkt+1=-Bjk+Bjk2-4AjkCjk2Ajk---(16)]]>由于我们使用估计的导数，故我们不能通过这种方法得到准确的结果，因此循环M-步骤直至X收敛。S7.判断是否满足迭代停止条件：迭代算法达到最大迭代数或满足收敛条件若不满足则执行步骤S4；若满足迭代停止条件则迭代停止，保存PET图像重建结果X。接下来我们通过对蒙特卡洛模拟的脑部phantom实验来验证本实施方式重建框架结果的准确性，图2为实验所用的脑部phantom的模板示意图，将不同的区域分为三个感兴趣的区域(regionofinterest，ROI)。实验运行环境为：8G内存，3.40GHz，64位操作系统，CPU为inteli7-3770。所模拟的PET扫描仪型号为HamamatsuSHR-22000，设定的放射性核素及药物为18F-FDG，设置sinogram为128个投影角度在每个角度下128条射束采集到的数据结果。系统矩阵G为16384×16384维实验前计算好的矩阵。本发明双字典稀疏约束的动态PET图像序列重建方法(DD)的结果和传统的ML-EM、单独的结构字典稀疏约束(DS)、单独的动力学参数字典稀疏约束(DK)的重建结果做比较。四种方法使用相同的测量数据矩阵Y和系统矩阵G以控制结果的可比性，重建第3帧和第19帧图像分别如图3～图6所示。可以直观的看到由于结构字典的作用，双字典约束方法和结构字典约束方法比ML-EM方法有更小的噪声，并且边缘部分也更好的保留下来。然而从误差图像更明显的看到，四种方法的主要误差都集中于边缘部分，而双字典约束和结构字典约束方法在非边界部分几乎不存在误差。这是由于非边界部分具有更高的稀疏度，即结构字典中的一个元素就足以表达出来，从而有效地抑制非边界部分的误差。对于边界部分，由于稀疏度较低，故可以看到边界部分存在一定的平滑。动力学参数字典约束方法只考虑到前后帧之间的动力学关联，故对于边界和非边界区域都有着相似的噪声特性。可以看出动力学方法重建结果边界特性较好，也在一定程度上抑制噪声。故我们的方法将这两者结合，不仅能提高成像图像信噪比，也较好的保持了边界的细节特性。表1和表2统计感兴趣的区域(ROI)和整体区域的相对偏差和方差数据。对于四种方法来说，由于ROI1和ROI3边界细节信息较多，故偏差和方差都较大；而ROI2处于非边界位置，其偏差和方差都小于整体值，这也与我们前面直观上的观察相符合。并且从不同帧的角度分析，由于第19帧放射性浓度低于第3帧，即第3帧信号强度大，信噪比高，因此第3帧相对偏差和方差值明显低于第19帧。双字典约束的方法由于结构字典抑制噪声的作用，故成像结果随成像系统信噪比变化不大，这也验证了我们算法的鲁棒性。表1和表2也进一步验证了本发明方法的准确性。表1表2上述对实施例的描述是为便于本
技术领域：
的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。当前第1页1 2 3