CN105374060A - 一种基于结构字典约束的pet图像重建方法 - Google Patents

Abstract

本发明公开了一种基于结构字典约束的PET图像重建方法，该方法通过建立重建问题的数学模型，加入结构字典约束，基于结构字典的约束来进行PET图像重建；其中结合Poisson和字典约束对PET图像进行重建的过程中，采取EM算法进行优化求解。故本发明有效利用结构字典约束，改善了计算机在进行PET图像重建的过程中产生的结果低分辨率和噪声干扰的问题；与现有重建方法的实验比较表明，本发明能获得较好的重建效果。

Description

一种基于结构字典约束的PET图像重建方法

技术领域

本发明属于PET成像技术领域，具体涉及一种基于结构字典约束的PET图像重建方法。

背景技术

正电子发射断层成像(Positronemissiontomography，PET)是一种基于核物理学、分子生物学的医学影像技术，它能够从分子水平上观察细胞的代谢活动，为早期疾病的检测和预防提供了有效依据。在进行PET测量前，首先让被测者吸入或注射一种放射性药剂，该药剂中含有与被测位置相对应的正电子核素，它们通常是由回旋加速器产生的。经过一小段时间，放射性核素到达了相应的区域，同时被目标组织吸收，此时即可进行扫描。当放射性核素衰变时，会向外发射出正电子，正电子移动很短的一段距离后，会和组织中的电子相遇并且湮灭。在湮灭时，会放出向相反方向移动的一对光子，它们的能量均为511keV。所以，我们通过探测器探测被测者体内放射出的光子对来确定湮灭事件发生的位置。通常，如果两个位置相对的光子在设定的时间窗(如10ns)之内被探测到，我们就认为这两个光子是在同一次湮灭中产生的，属于符合计数(atruecoincidence),这个正电子放射事件就会被记录下来。所有正电子放射事件的集合可以近似等于放射性核素浓度分布的线积分，记录的事件数量越多，则近似度越高。

PET图像具有高灵敏度、高特异性的优点，但由于放射性核素受人体组织的影响会发生严重的衰减，且校正衰减的方法复杂、成本高，因此由测量得到的数据重建后的图像分辨率较低，图像略模糊。传统上，放射性浓度分布重建往往采用统计迭代方法，由于迭代法基于统计学模型，对不完全数据适应性好，逐渐成为PET重建算法研究关注点，其中包括著名的MLEM(最大似然期望最大化)、MAP(最大后验)和SAGE(惩罚似然)算法。如何获得更精确、清晰的重建图像是大家研究的热点。

发明内容

针对现有技术所存在的上述技术问题，本发明提供了一种基于结构字典约束的PET图像重建方法，能够获得高质量的PET重建图像。

一种基于结构字典约束的PET图像重建方法，包括如下步骤：

(1)利用探测器对注入有放射性药剂的生物组织进行探测，采集得到对应探测器各晶体块的符合计数向量，进而构建PET的符合计数矩阵y；

(2)根据PET成像原理，建立PET的测量方程如下；通过对所述的测量方程引入Poisson噪声约束，得到PET的Poisson模型L(x)；

y＝Gx+r+v

其中：G为系统矩阵，x为PET浓度分布矩阵，r和v分别为关于反射符合事件和散射符合事件的测量噪声矩阵；

(3)通过对所述的Poisson模型L(x)引入结构字典约束，得到基于结构字典约束的PET图像重建模型如下：

$\underset{x,\mathrm{\α}}{\min }\{\mathrm{\λ}L\left(x\right)+R_{sparse}(x,\mathrm{\α})\}$

其中：λ为权重系数，R_sparse(x，α)为关于PET浓度分布矩阵x和稀疏系数矩阵α的惩罚项；

(4)对上述PET图像重建模型进行优化求解得到PET浓度分布矩阵x，进而重建获得PET图像。

所述Poisson模型L(x)的表达式如下：

$L\left(x\right)=\underset{i=1}{\overset{n_i}{\Σ}}\{{\stackrel{\‾}{y}}_i-y_i\log \left({\stackrel{\‾}{y}}_i\right)\}$ s.t.y＝Gx+r+v

其中：y_i为符合计数矩阵y中对应第i个晶体块的符合计数向量，为符合计数向量y_i中所有元素的平均值，n_i为探测器的晶体块总数。

所述惩罚项R_sparse(x，α)的表达式如下：

$R_{sparse}(x,\mathrm{\α})=\underset{s=1}{\overset{n_s}{\Σ}}\left\{\right||E_{s}x-{\mathrm{D\α}}_s|{|}_2^2+\mathrm{\μ}\left|\right|{\mathrm{\α}}_s\left|{|}_0\right\}$

其中：D为结构字典(其为利用K-SVD算法从CT图像中获取的矩阵，具体参考文献：K-SVD：Analgorithmfordesigningovercompletedictionariesforsparserepresentation)，E_s为分割算子，E_sx为PET浓度分布矩阵x中的第s个n维子矩阵，μ为权重系数，n_s＝(m-n+1)²，m为PET浓度分布矩阵x的维度，n为预设的子矩阵维度(一般取7或8)，α_s为稀疏系数矩阵α中第s列稀疏系数向量，||α_s||₀表示稀疏系数向量α_s中非零元素的个数，||||₂表示L2范数。

所述的步骤(4)中对PET图像重建模型进行优化求解的具体方法为：基于EM(ExpectationMaximization，选择期望值最大)算法通过以下迭代方程求解得到PET浓度分布矩阵x：

$x_j^{k+1}=\frac{-B_j^k+\sqrt{{\left(B_j^k\right)}^2-4A_j^k{C}_j^k}}{2A_j^k}$

$B_j^k=\underset{s=1}{\overset{n_s}{\Σ}}\underset{l=1}{\overset{n_l}{\Σ}}2e_{slj}({\[E_s{x}^k\]}_l-{\[{\mathrm{D\α}}_s^k\]}_l)-A_j^k{x}_j^k+\mathrm{\λ}\underset{i}{\overset{n_i}{\Σ}}{g}_{ij}$

$A_j^k=\underset{s=1}{\overset{n_s}{\Σ}}\underset{l=1}{\overset{n_l}{\Σ}}2e_{slj}\left(\frac{e_{slj}}{{\mathrm{\ϵ}}_{slj}}\right),C_j^k=-\mathrm{\λ}\underset{i}{\overset{n_i}{\Σ}}{\hat{w}}_{ij}^k$

其中：和均为中间变量，为第k+1次迭代PET浓度分布矩阵x^k+1中的第j个元素值，为第k次迭代PET浓度分布矩阵x^k中的第j个元素值，e_slj为分割算子E_s中的第l行第j列元素值，n_j＝m²，E_sx^k为第k次迭代PET浓度分布矩阵x^k中的第s个n维子矩阵，[E_sx^k]_l为子矩阵E_sx^k中的第l个元素值，n_l＝n²，为结构字典D与第k次迭代稀疏系数向量相乘所得到矩阵中的第l个元素值，g_ij为系统矩阵G中的第i行第j列元素值，表示从第k次迭代PET浓度分布矩阵x^k中的第j个元素中发射出被探测器第i个晶体块探测到的光子数，k为迭代次数。

所述光子数的表达式如下：

${\hat{w}}_{ij}^k=\frac{g_{ij}{x}_j^k}{{\Σ}_{j=1}^{n_j}{g}_{ij}{x}_j^k+r_i+v_i}{y}_i$

其中：r_i和v_i分别为测量噪声矩阵r和v中的第i个元素值。

所述的第k次迭代稀疏系数向量采用OMP(OrthogonalMatchingPursuit，正交匹配追踪)算法对以下目标方程求解获得：

$\begin{array}{ccc}{\mathrm{\α}}_s^k=\arg \underset{{\mathrm{\α}}_s}{\min }\left|\right|{\mathrm{\α}}_s|{|}_0& s.t.& \left|\right|E_s{x}^k-{\mathrm{D\α}}_s|{|}_2^2\≤\mathrm{\ρ}\end{array}$

其中：ρ为预设的极小值。

本发明通过建立重建问题的数学模型，加入结构字典约束，基于结构字典的约束来进行PET图像重建；其中结合Poisson和字典约束对PET图像进行重建的过程中，采取EM算法进行优化求解。故本发明有效利用结构字典约束，改善了计算机在进行PET图像重建的过程中产生的结果低分辨率和噪声干扰的问题；与现有重建方法的实验比较表明，本发明能获得较好的重建效果。

附图说明

图1(a)为CT胸腔切片图像。

图1(b)为根据图1(a)的CT胸腔切片图像训练出的结构字典矩阵图形。

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

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

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

具体实施方式

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

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

本发明基于结构字典的PET图像重建方法，包括如下步骤：

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

S2.引入Poisson和结构字典来优化问题模型；

S3.初始化，设置权重系数λ，设置初始x，α值，x＝FBP(x)；α＝0；

S4.设置的初始点开始，按照EM算法计算。E步骤中，我们先用已经估计出的第k次迭代时x^k的值和正弦图y来估计隐藏变量w_ij。将估计得的代入到Ω_x(w_ij，x，α)中，我们得到即时方程M步骤中，通过使即时方程的导数等于零，我们求出了新的x^k+1；

S5.判断是否满足迭代停止条件x＜10^-3，不满足该条件则执行步骤S4，满足则迭代停止，执行步骤S6。

S6.更新α；准备当S7判断不合格后，将更新后的参数返回步骤S4进行循环；

S7.判断是否满足迭代停止条件x＜10^-4，不满足该条件则执行步骤S4，满足则迭代停止；进而得到PET浓度分布向量X实现PET成像。

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

y＝Gx+r+v

其中：G为系统矩阵，y为校正后的符合计数向量，x为PET浓度分布向量，r和v分别表示反射符合事件和散射符合事件的测量噪声矩阵。

符合计数向量的Poisson模型的表达式如下：

s.t.y＝Gx+r+v，

其中：表示y服从均值为的泊松分布，基于独立泊松假设，y的似然方程Pr(y|x)表达式如下：

$\mathrm{Pr}\left(y\right|x)=\underset{i}{\overset{n_i}{\Π}}{e}^{-{\stackrel{\‾}{y}}_i}\frac{{{\stackrel{\‾}{y}}_i}^{yi}}{y_i!}$

$\underset{x}{\min }L\left(x\right)=\underset{x}{\min }-\log \left(\mathrm{Pr}\right(y|x\left)\right)=\underset{x}{\min }\underset{i}{\overset{n_i}{\Σ}}{\stackrel{\‾}{y}}_i-y_i\log \left({\stackrel{\‾}{y}}_i\right)$

$\begin{array}{cc}s.t.& {\stackrel{\‾}{y}}_i={\Σ}_{j=1}^{n_j}{g}_{ij}{x}_j+r_i+v_i\end{array}$

其中：L(x)为重建问题的目标函数，y_i为第i个探测器的测量值，n_i表示探测器晶体块的总数。

其加入结构字典约束的表达式如下：

$\underset{x,\mathrm{\α}}{\min }\mathrm{\Ω}(x,\mathrm{\α})=\underset{x,\mathrm{\α}}{\min }\mathrm{\λ}L\left(x\right)+\underset{s}{\overset{n_s}{\Σ}}\left|\right|E_{s}x-{\mathrm{D\α}}_s|{|}_2^2+\mathrm{\μ}|\left|{\mathrm{\α}}_s\right|{|}_0$

其中：L(x)为重建问题的目标函数，α_s是稀疏系数，D是结构字典，E_s是分割矩阵，λ是权重系数，n_s是分块数，λ和μ均为权重系数。

结构字典D是根据Aharon在2006发表的文献，利用K-SVD算法从CT图像中获取的矩阵(参考文献：K-SVD:AnAlgorithmforDesigningOvercompleteDictionariesforSparseRepresentation)。图1(a)为CT胸腔切片图像，图1(b)为根据该CT胸腔切片图像训练出的结构字典矩阵图形。

我们首先解决x子问题，将目标函数改写成如下形式：

$x=\arg \underset{x}{\min }{\mathrm{\Ω}}_x(x,\mathrm{\α})=\arg \underset{x}{\min }\mathrm{\λ}\underset{i}{\overset{n_i}{\Σ}}{\stackrel{\‾}{y}}_i-y_i\log \left({\stackrel{\‾}{y}}_i\right)+\underset{s}{\overset{n_s}{\Σ}}\left|\right|E_{s}x-{\mathrm{D\α}}_s|{|}_2^2$

$\begin{array}{cc}s.t.& {\stackrel{\‾}{y}}_i={\Σ}_{j=1}^{n_j}{g}_{ij}{x}_j+r_i+v_i\end{array}$

其中：y_i为第i个探测器的测量值，α_s是稀疏系数，D是结构字典，E_s是分割矩阵，λ是权重系数，n_s是分块数，λ和μ均为权重系数。

接下来引入一个隐藏变量w_ij：

$\begin{array}{c}x=\arg \underset{x}{\min }{\mathrm{\Ω}}_x(w_{ij},x,\mathrm{\α})=\arg \underset{x}{\min }\mathrm{\λ}\underset{j=1}{\overset{n_j}{\Σ}}\underset{i=1}{\overset{n_i}{\Σ}}(g_{ij}{x}_j-w_{ij}\log (g_{ij}{x}_j\left)\right)+\\ \underset{s}{\overset{n_s}{\Σ}}\left|\right|E_{s}x-{\mathrm{D\α}}_s|{|}_2^2\end{array}$

其中：隐藏变量w_ij表示从体素j中发射出被探测器晶体块i探测到的光子数，g_ij是系统矩阵G中的第ij项。

若w_ij是可知的，则此时方程有解。之后，将EM算法分为两步进行：

步骤E：用已经估计出的第k次迭代时x^k的值和正弦图y来估计隐藏变量w_ij。将估计得的代入到Ω_x(w_ij，x，α)中，得到即时方程

步骤M：通过使即时方程的导数等于零，求出新的x。

在步骤E中，在已知目前估计值x_j和正弦图y的条件下，估计w_ij的期望式为：

${\hat{w}}_{ij}=\frac{g_{ij}{x}_j^k}{{\Σ}_j^{n_j}{g}_{ij}{x}_j^k+{\hat{r}}_i+{\hat{v}}_i}{y}_i$

在步骤M中，不能直接对这一项求导，所以用它的凸可分代理方程来代替它自己。首先将和[E_sx]_l改写成如下形式：

$\left|\right|E_{s}x-{\mathrm{D\α}}_s|{|}_2^2=\underset{l}{\overset{n_l}{\Σ}}{({\[E_{s}x\]}_l-{\[{\mathrm{D\α}}_s\]}_l)}^2$

$\begin{array}{ccc}{\[E_{s}x\]}_l=\underset{j}{\overset{n_j}{\Σ}}{\mathrm{\ϵ}}_{slj}\left(\frac{e_{slj}}{{\mathrm{\ϵ}}_{slj}}\right(x_j-x_j^k)+{\[E_s{x}^k\]}_l)& s.t.& {\mathrm{\ϵ}}_{slj}=\frac{e_{slj}}{{\Σ}_j^{n_j}{e}_{slj}}\end{array}$

其中：[E_sx]_l代表E_sx的第l项，n_l表示E_sx的矩阵中总元素个数，是当前估计出的图像x在第j个像素点的值，e_sij是矩阵E_s的第lj项，对于所有j，ε_sij均大于0。

由于([E_sx]_l-[Dα_s]_l)²是凸函数，所以：

${({\[E_{s}x\]}_l-{\[{\mathrm{D\α}}_s\]}_l)}^2\≤\underset{j}{\overset{n_j}{\Σ}}{\mathrm{\ϵ}}_{slj}{\left(\frac{e_{slj}}{{\mathrm{\ϵ}}_{slj}}\right(x_j-x_j^k)+{\[E_s{x}^k\]}_l-{\[{\mathrm{D\α}}_s\]}_l)}^2$

代入方程中，得到它的凸可分代理方程ω(x，x^k)：

$\begin{array}{c}\mathrm{\ω}(x;x^k)=\mathrm{\λ}\underset{j}{\overset{n_j}{\Σ}}\underset{i}{\overset{n_i}{\Σ}}(g_{ij}{x}_j-{\hat{w}}_{ij}\log (g_{ij}{x}_j\left)\right)\\ +\underset{s}{\overset{n_s}{\Σ}}\underset{l}{\overset{n_l}{\Σ}}\underset{j}{\overset{n_j}{\Σ}}{\mathrm{\ϵ}}_{slj}{\left(\frac{e_{slj}}{{\mathrm{\ϵ}}_{slj}}\right(x_j-x_j^k)+{\[E_s{x}^k\]}_l-{\[{\mathrm{D\α}}_s\]}_l)}^2\end{array}$

求关于x的偏导：

$\begin{array}{c}\frac{\∂\mathrm{\ω}(x;x^k)}{\∂x_j}=\mathrm{\λ}{\Σ}_i^{n_i}{g}_{ij}-\mathrm{\λ}{\Σ}_i^{n_i}{\hat{w}}_{ij}\frac{1}{x_j}+\\ {\Σ}_s^{n_s}{\Σ}_l^{n_l}2e_{slj}\left(\frac{e_{slj}}{{\mathrm{\ϵ}}_{slj}}\right(x_j-x_j^k)+{\[E_s{x}^k\]}_l-{\[{\mathrm{D\α}}_s\]}_l)=0\end{array}$

是二阶多项式的解：

$A_j{x}_j+B_j+C_j\frac{1}{x_j}=0,A_j=\underset{s}{\overset{n_s}{\Σ}}\underset{l}{\overset{n_l}{\Σ}}2e_{slj}\left(\frac{e_{slj}}{{\mathrm{\ϵ}}_{slj}}\right),C_j-\mathrm{\λ}\underset{i}{\overset{n_i}{\Σ}}{\hat{w}}_{ij}$

$B_j=\underset{s}{\overset{n_s}{\Σ}}\underset{l}{\overset{n_l}{\Σ}}2e_{slj}({\[E_s{x}^k\]}_l-{\[{\mathrm{D\α}}_s\]}_l)-A_j{x}_j^k+\mathrm{\λ}\underset{i}{\overset{n_i}{\Σ}}{g}_{ij}$

$x_j^{k+1}=\frac{-B_j+\sqrt{B_j^2-4A_j{C}_j}}{2A_j}$

从上述过程可知，x子问题是无法直接求解的，我们需要通过E步骤和M步骤的不断迭代使得收敛至一个近似解。

我们选用最佳匹配跟踪(OrthogonalMatchingPursuit，OMP)算法来求解α子问题。

$\begin{array}{ccc}\∀{\mathrm{s\α}}_s\arg \underset{{\mathrm{\α}}_s}{\min }\left|\right|{\mathrm{\α}}_s|{|}_0& s.t.& \left|\right|E_s{x}^k-{\mathrm{D\α}}_s|{|}_2^2\≤\mathrm{\ϵ}\end{array}$

其中：α_s是稀疏系数，D是结构字典，E_s是分割矩阵，||α_s||₀表示稀疏系数α_s中的非零元素的个数，∈为一个极小值。

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

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

针对重建图像质量的验证，采用高维度的原始采集数据128个投影角度，每个角度下射束为128条，即m＝16384，重建图像大小为128×128，即维度为16384；初始值设定同上。图2(a)～图2(c)分别为真值图像、基于传统ML-EM方法重建的图像和基于本实施方式重建的图像的比对示意图，可以直观地看出基于本实施方式重建的图像与ML-EM的结果相比，ML-EM算法的效果并不理想，重建图像中图像模糊、存在许多噪点且边缘十分不清晰。对比之下，我们的算法重建后的图片各区域有清晰地边界，且图像内部平滑。

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

表1

方法	偏差	方差	均方根误差
				ML-EM	0.1919	0.0543	0.2331
本实施方式	0.1153	0.0341	0.1845

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

Claims (6)

1.一种基于结构字典约束的PET图像重建方法，包括如下步骤：

(1)利用探测器对注入有放射性药剂的生物组织进行探测，采集得到对应探测器各晶体块的符合计数向量，进而构建PET的符合计数矩阵y；

(2)根据PET成像原理，建立PET的测量方程如下；通过对所述的测量方程引入Poisson噪声约束，得到PET的Poisson模型L(x)；

y＝Gx+r+v

其中：G为系统矩阵，x为PET浓度分布矩阵，r和v分别为关于反射符合事件和散射符合事件的测量噪声矩阵；

(3)通过对所述的Poisson模型L(x)引入结构字典约束，得到基于结构字典约束的PET图像重建模型如下：

$\underset{x,\mathrm{\α}}{min}\{\mathrm{\λ}L\left(x\right)+R_{sparse}(x,\mathrm{\α})\}$

其中：λ为权重系数，R_sparse(x，α)为关于PET浓度分布矩阵x和稀疏系数矩阵α的惩罚项；

(4)对上述PET图像重建模型进行优化求解得到PET浓度分布矩阵x，进而重建获得PET图像。

2.根据权利要求1所述的PET图像重建方法，其特征在于：所述Poisson模型L(x)的表达式如下：

$\begin{array}{cc}L\left(x\right)=\underset{i=1}{\overset{n_i}{\Σ}}\{{\stackrel{\‾}{y}}_i-y_i\log \left({\stackrel{\‾}{y}}_i\right)\}& s.t.y=Gx+r+v\end{array}$

其中：y_i为符合计数矩阵y中对应第i个晶体块的符合计数向量，为符合计数向量y_i中所有元素的平均值，n_i为探测器的晶体块总数。

3.根据权利要求2所述的PET图像重建方法，其特征在于：所述惩罚项R_sparse(x，α)的表达式如下：

$R_{sparse}(x,\mathrm{\α})=\underset{s=1}{\overset{n_s}{\mathrm{\Σ}}}\left\{\right||E_{s}x-{\mathrm{D\α}}_s|{|}_2^2+\mathrm{\μ}\left|\right|{\mathrm{\α}}_s\left|{|}_0\right\}$

其中：D为结构字典，E_s为分割算子，E_sx为PET浓度分布矩阵x中的第s个n维子矩阵，μ为权重系数，n_s＝(m-n+1)²，m为PET浓度分布矩阵x的维度，n为预设的子矩阵维度，α_s为稀疏系数矩阵α中第s列稀疏系数向量，||α_s||₀表示稀疏系数向量α_s中非零元素的个数，||||₂表示L2范数。

4.根据权利要求3所述的PET图像重建方法，其特征在于：所述的步骤(4)中对PET图像重建模型进行优化求解的具体方法为：基于EM算法通过以下迭代方程求解得到PET浓度分布矩阵x：

$x_j^{k+1}=\frac{-B_j^k+\sqrt{{\left(B_j^k\right)}^2-4A_j^k{C}_j^k}}{2A_j^k}$

$B_j^k=\underset{s=1}{\overset{n_s}{\Σ}}\underset{l=1}{\overset{n_l}{\Σ}}2e_{elj}({\[E_s{x}^k\]}_l-{\[{\mathrm{D\α}}_s^k\]}_l)-A_j^k{x}_j^k+\mathrm{\λ}\underset{i}{\overset{n_i}{\Σ}}{g}_{ij}$

$\begin{array}{cc}{A}_i^k=\underset{s=1}{\overset{n_s}{\Σ}}\underset{l=1}{\overset{n_l}{\Σ}}2e_{slj}\left(\frac{e_{slj}}{{\mathrm{\ϵ}}_{slj}}\right),& C_j^k=-\mathrm{\λ}\underset{i}{\overset{n_i}{\Σ}}{\hat{w}}_{ij}^k\end{array}$

其中：和均为中间变量，为第k+1次迭代PET浓度分布矩阵x^k+1中的第j个元素值，为第k次迭代PET浓度分布矩阵x^k中的第j个元素值，e_slj为分割算子E_s中的第l行第j列元素值，n_j＝m²，E_sx^k为第k次迭代PET浓度分布矩阵x^k中的第s个n维子矩阵，[E_sx^k]_l为子矩阵E_sx^k中的第l个元素值，n_l＝n²，为结构字典D与第k次迭代稀疏系数向量相乘所得到矩阵中的第l个元素值，g_ij为系统矩阵G中的第i行第j列元素值，表示从第k次迭代PET浓度分布矩阵x^k中的第j个元素中发射出被探测器第i个晶体块探测到的光子数，k为迭代次数。

5.根据权利要求4所述的PET图像重建方法，其特征在于：所述光子数的表达式如下：

${\hat{w}}_{ij}^k=\frac{g_{ij}{x}_j^k}{{\mathrm{\Σ}}_{j=1}^{n_j}{g}_{ij}{x}_j^k+r_i+v_i}{y}_i$

其中：r_i和v_i分别为测量噪声矩阵r和v中的第i个元素值。

6.根据权利要求4所述的PET图像重建方法，其特征在于：所述的第k次迭代稀疏系数向量采用OMP算法对以下目标方程求解获得：

$\begin{array}{cc}{\mathrm{\α}}_s^k=\arg \underset{{\mathrm{\α}}_s}{min}\left|\right|{\mathrm{\α}}_s|{|}_0& s.t.\left|\right|E_s{x}^k-{\mathrm{D\α}}_s|{|}_2^2\≤\mathrm{\ρ}\end{array}$

其中：ρ为预设的极小值。