CN103400403A - 一种pet浓度与衰减系数的同时重建方法 - Google Patents

Abstract

本发明公开了一种PET浓度与衰减系数的同时重建方法，包括：（1）采集符合计数；（2）构建PET状态空间方程；（3）利用UKF算法估算出PET浓度分布和衰减系数分布。本发明利用鲁棒的无味卡尔曼滤波框架对衰减系数和放射性浓度进行同时估计，不需要线性化Jacobian矩阵，同时使用采样的方法来保证对均值和方差的估计；相对基于EKF的重建方法，本发明能够更好的处理系统噪声、测量数据噪声和建模的误差。

Description

一种PET浓度与衰减系数的同时重建方法

技术领域

本发明属于PET成像技术领域，具体涉及一种PET浓度与衰减系数的同时重建方法。

背景技术

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

定量分析的基础是像素计数，要求体内某部位的像素计数值能够精确反应出该处所摄取的放射性药物浓度。PET数据采集之后要进行各种各样的校正，而其中对定量准确性影响最大的就是衰减校正。当前PET的衰减校正多采用透射扫描的方式或者利用PET/CT中CT的扫描数据来进行但使用CT的扫描过程极易引入伪影，影响成像质量，目前最准确的衰减校正方法还是透射扫描。

透射扫描最大的问题在于耗时长，并且病人的状态可能发生变化，而且会增加病人的射线被曝量，所以出现了一些关于衰减系数和放射性浓度同时估计的研究，但是他们更多的集中在对SPECT（单光子发射计算机断层成像技术）采用不同的衰减方法的研究上。在以前的工作中也提到过使用状态空间方法扩展卡尔曼解法（EKF）来同时估计衰减系数和放射性浓度分布的体系，但是该体系还存在着几个问题，如Jacobian（雅克比）矩阵的导出、非线性问题的线性化等，这些都会造成估计过程的误差。同时，由于扩展了状态向量的维度，增加了状态变量和模型参数，对求解结果的标准差也会产生潜在的影响。

发明内容

针对现有技术所存在的上述技术问题，本发明提供了一种PET浓度与衰减系数的同时重建方法，相对EKF能够更好的消除重建过程中的误差。

一种PET浓度与衰减系数的同时重建方法，包括如下步骤：

（1）利用探测器对注入有放射性物质的生物组织进行探测，采集得到PET的符合计数；

（2）根据PET成像原理，建立含有衰减因素的PET状态空间方程；

（3）利用无色卡尔曼滤波算法求解所述的PET状态空间方程，得到以下迭代方程组；进而根据符合计数通过以下迭代方程组估计出PET浓度分布信息和衰减系数分布信息；

x_k＝x_k-1+K_k(y-g_k)

μ_k＝μ_k-1+H_k(y-z_k)

其中：x_k和x_k-1分别为第k次和第k-1次迭代后的PET浓度分布向量，μ_k和μ_k-1分别为第k次和第k-1次迭代后的衰减系数分布向量，K_k和H_k均为第k次迭代时的增益矩阵，g_k和z_k均为第k次迭代时的中间向量，y为符合计数。

所述的PET状态空间方程的表达式如下：

x_t+1＝x_t+v

μ_t+1＝μ_t+s

y＝GDx+e

其中：x_t和x_t+1分别为第t状态和第t+1状态下的PET浓度分布向量，μ_t和μ_t+1分别为第t状态和第t+1状态下的衰减系数分布向量，G为衰减矩阵，D为系统矩阵，v和s均为过程噪声向量；e为测量噪声向量，y为符合计数，x为PET浓度分布向量。

所述的增益矩阵K_k根据以下算式求解：

K_k＝F_kE_k ^-1

$E_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{w}_j(Y_{k,j}-g_k){(Y_{k,j}-g_k)}^T+R_0$

$F_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{w}_j(X_{k,j}-x_{k-1}){(Y_{k,j}-g_k)}^T$

其中：X_k,j和Y_k,j分别为第k次迭代时的sigma矩阵X_k和Y_k中的第j列向量，R₀为测量噪声的协方差矩阵，w_j为对应第j列向量的权重因子，n为PET浓度分布向量的维度。

所述的增益矩阵H_k根据以下算式求解：

H_k＝V_kT_k ^-1

$V_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{w}_j(U_{k,j}-{\mathrm{\μ}}_{k-1}){(Z_{k,j}-z_k)}^T$

$T_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{w}_j(Z_{k,j}-z_k){(Z_{k,j}-z_k)}^T+R_0$

其中：U_k,j和Z_k,j分别为第k次迭代时的sigma矩阵U_k和Z_k中的第j列向量，R₀为测量噪声的协方差矩阵，w_j为对应第j列向量的权重因子，n为PET浓度分布向量的维度。

所述的中间向量g_k根据以下算式求解：

$g_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{u}_j{Y}_{k,j}$

其中：Y_k,j为第k次迭代时的sigma矩阵Y_k中的第j列向量，u_j为对应第j列向量的权重因子，n为PET浓度分布向量的维度。

所述的中间向量z_k根据以下算式求解：

$z_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{u}_j{Z}_{k,j}$

其中：Z_k,j为第k次迭代时的sigma矩阵Z_k中的第j列向量，u_j为对应第j列向量的权重因子，n为PET浓度分布向量的维度。

所述的权重因子u_j由以下算式求得：

$u_j=\frac{\mathrm{\λ}}{n+\mathrm{\λ}},j=1$

$u_j=\frac{\mathrm{\λ}}{2(n+\mathrm{\λ})},j=2,...,2n+1$

其中：λ=n(α²-1)，α为给定的计算参数。

所述的权重因子w_j由以下算式求得：

$w_j=\frac{\mathrm{\λ}}{n+\mathrm{\λ}}+(1-{\mathrm{\α}}^2+\mathrm{\β}),j=1$

$w_j=\frac{\mathrm{\λ}}{2(n+\mathrm{\λ})},j=2,...,2n+1$

其中：λ=n(α²-1)，α和β均为给定的计算参数。

所述的sigma矩阵X_k中的第j列向量X_k,j由以下关系式确定：

X_k,j＝x_k-1,j＝1

$X_{k,j}=x_{k-1}+\sqrt{n+\mathrm{\λ}}{M}_{k,j},j=2,...,n+1$

$X_{k,j}=x_{k-1}-\sqrt{n+\mathrm{\λ}}{M}_{k,j},j=n+2,...,2n+1$

其中：λ=n(α²-1)，α为给定的计算参数，M_k,j为第k次迭代时的中间矩阵M_k中的第j列向量。

所述的sigma矩阵U_k中的第j列向量U_k,j由以下关系式确定：

U_k,j＝μ_k-1,j＝1

$U_{k,j}={\mathrm{\μ}}_{k-1}+\sqrt{n+\mathrm{\λ}}{N}_{k,j},j=2,...,n+1$

$U_{k,j}={\mathrm{\μ}}_{k-1}-\sqrt{n+\mathrm{\λ}}{N}_{k,j},j=n+2,...,2n+1$

其中：λ=n(α²-1)，α为给定的计算参数，N_k,j为第k次迭代时的中间矩阵N_k中的第j列向量。

所述的中间矩阵M_k由以下关系式确定：

$M_{k,i,j}=\sqrt{A_{k,i,j}}$

A_k＝P_k-1+Q_v

P_k-1＝A_k-1-K_k-1E_k-1(K_k-1)^T

其中：M_k,i,j和A_k,i,j分别为中间矩阵M_k和中间矩阵A_k中第i行第j列的元素值，P_k-1为PET浓度分布向量x_k-1的协方差矩阵，Q_v为系统噪声的协方差矩阵。

所述的中间矩阵N_k由以下关系式确定：

$N_{k,i,j}=\sqrt{B_{k,i,j}}$

B_k＝S_k-1+Q_s

S_k-1＝B_k-1-H_k-1T_k-1(H_k-1)^T

其中：N_k,i,j和B_k,i,j分别为中间矩阵N_k和中间矩阵B_k中第i行第j列的元素值，S_k-1为衰减系数分布向量μ_k-1的协方差矩阵，Q_s为系统噪声的协方差矩阵。

所述的sigma矩阵Y_k中的第j列向量Y_k,j由以下关系式确定：

Y_k,j＝G_kDX_k,j

其中：D为系统矩阵，G_k为第k次迭代时的衰减矩阵。

所述的衰减矩阵G_k由以下关系式确定：

G_k＝diag[exp(-Lμ_k-1)]

其中：L为投影线与像素点的交叉长度矩阵。

所述的sigma矩阵Z_k中的第j列向量Z_k,j由以下关系式确定：

Z_k,j＝W_k,jDx_k

其中：D为系统矩阵，W_k,j为第k次迭代时的第j个对角矩阵。

所述的对角矩阵W_k,j由以下关系式确定：

W_k,j＝diag[exp(-LU_k,j)]

其中：L为投影线与像素点的交叉长度矩阵。

所述的步骤（3）中，根据符合计数通过迭代方程组进行迭代计算，则迭代收敛后的PET浓度分布向量和衰减系数分布向量即为所要估计的PET浓度分布信息和衰减系数分布信息；所述的迭代收敛条件如下：

$\frac{\underset{j=1}{\overset{n^2}{\mathrm{\Σ}}}|p_{k,j}-p_{k-1,j}|}{n^2}\≤\mathrm{\ρ}$

其中：p_k,j和p_k-1,j分别为PET浓度分布向量x_k和x_k-1的协方差矩阵P_k和P_k-1中第j个元素值，ρ为给定的收敛阈值，n为PET浓度分布向量的维度。

本发明中，所述的PET浓度分布向量x以及衰减系数分布向量μ均为n维向量，符合计数y为m维向量；过程噪声向量v和s均为n维向量，表征了浓度分布状态之间转移的统计不确定性；测量噪声向量e为m维向量，代表符合事件采集过程中由随机事件、散射事件等引入的测量误差；衰减矩阵G为m×m维的对角矩阵，表示衰减系数分布对PET数据采集过程的影响，其对角元素值可通过以下算式求得：

G＝diag[exp(-Lμ)]

μ为n维的衰减系数分布向量，L为各投影线与各像素点的交叉长度矩阵，且具有m×n维，其元素L_i,j表示像素点j发出的光子对落在投影线i确定的探测器簇内的概率，由系统的几何结构确定；R₀为测量噪声向量e的协方差矩阵，且为m×m维矩阵；Q_v和Q_s分别为过程噪声向量v和s的协方差矩阵，均为n×n维矩阵。

本发明利用鲁棒的无味卡尔曼滤波（UKF）框架对衰减系数和放射性浓度进行同时估计，不需要线性化Jacobian矩阵，同时使用采样的方法来保证对均值和方差的估计；相对基于EKF的重建方法，本发明能够更好的处理系统噪声、测量数据噪声和建模的误差。

附图说明

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

具体实施方式

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

如图1所示，一种PET浓度与衰减系数的同时重建方法，包括如下步骤：

（1）采集符合计数。

利用探测器对注入有两种放射性示踪剂的生物组织进行探测，采集得到PET的符合计数y；本实施方式中，探测器采用日本滨松公司生产的型号为SHR74000的PET扫描仪。

PET扫描仪探测人体内发出的放射性信号，经过符合和采集系统处理，形成原始符合事件，PET探测器记录的符合事件包括真符合、随机符合和散射符合；本实施方式通过探测器的延时窗口和能量窗口对随机事件和散射事件进行校正，得到正弦图数据即校正后的符合计数y，其为m维向量。

（2）构建PET状态空间方程。

PET的发射采集过程为一对光子的符合采集过程。列向量x={x_j|j=1，…，n}代表生物体内的放射性药物浓度的分布，n为图像划分的总的像素数，y={y_i|i=1，…，m}为采集到的正弦图，其中、代表以探测器旋转角为序，对应投影角度下的探测器簇的排列序数，m对应各个投影角度下探测器簇的排列总数。μ={μ_j|j=1，…，n}表示第j个像素点处的生物组织对光子的线性衰减系数。当考虑衰减因素时，PET数据采集过程的系统矩阵可以分解为两部分：衰减矩阵G和不考虑衰减时的系统矩阵D。于是其采集过程可表示为：

y＝GDx+e

对于静态重建问题，放射性浓度分布认为是不变的，考虑到过程的不确定因素，其状态转移过程可表述成：

x_t+1＝x_t+v

在整个数据过程中，我们假设系统衰减系数向量μ是时不变的，同时考虑过程的不确定因素，其状态转移过程可表述成：

μ_t+1＝μ_t+s

联立以上三式，得到以下含有衰减因素的PET状态空间方程：

x_t+1＝x_t+v

μ_t+1＝μ_t+s

y＝GDx+e

其中：x_t和x_t+1分别为第t状态和第t+1状态下的PET浓度分布向量，μ_t和μ_t+1分别为第t状态和第t+1状态下的衰减系数分布向量，G为衰减矩阵，D为不考虑衰减时的系统矩阵，v和s均为过程噪声向量；e为测量噪声向量，x为PET浓度分布向量。

（3）利用无色卡尔曼滤波算法估算出PET浓度分布和衰减系数分布。

利用无色卡尔曼滤波算法求解上述PET状态空间方程，得到以下迭代方程组；进而根据符合计数y通过以下迭代方程组估计出PET浓度分布信息和衰减系数分布信息；

x_k＝x_k-1+K_k(y-g_k)

μ_k＝μ_k-1+H_k(y-z_k)

其中：x_k和x_k-1分别为第k次和第k-1次迭代后的PET浓度分布向量，μ_k和μ_k-1分别为第k次和第k-1次迭代后的衰减系数分布向量，K_k和H_k均为第k次迭代时的增益矩阵，g_k和z_k均为第k次迭代时的中间向量。

增益矩阵K_k根据以下算式求解：

K_k＝F_kE_k ^-1

$E_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{w}_j(Y_{k,j}-g_k){(Y_{k,j}-g_k)}^T+R_0$

$F_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{w}_j(X_{k,j}-x_{k-1}){(Y_{k,j}-g_k)}^T$

其中：X_k,j和Y_k,j分别为第k次迭代时的sigma矩阵X_k和Y_k中的第j列向量，R₀为测量噪声的协方差矩阵，w_j为对应第j列向量的权重因子，n为PET浓度分布向量的维度；

增益矩阵H_k根据以下算式求解：

H_k＝V_kT_k ^-1

$V_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{w}_j(U_{k,j}-{\mathrm{\μ}}_{k-1}){(Z_{k,j}-z_k)}^T$

$T_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{w}_j(Z_{k,j}-z_k){(Z_{k,j}-z_k)}^T+R_0$

其中：U_k,j和Z_k,j分别为第k次迭代时的sigma矩阵U_k和Z_k中的第j列向量；

中间向量g_k和z_k根据以下算式求解：

$g_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{u}_j{Y}_{k,j}$ $z_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{u}_j{Z}_{k,j}$

其中：u_j为对应第j列向量的权重因子；

权重因子u_j和w_j由以下算式求得：

$u_j=\frac{\mathrm{\λ}}{n+\mathrm{\λ}},j=1$

$u_j=\frac{\mathrm{\λ}}{2(n+\mathrm{\λ})},j=2,...,2n+1$

$w_j=\frac{\mathrm{\λ}}{n+\mathrm{\λ}}+(1-{\mathrm{\α}}^2+\mathrm{\β}),j=1$

$w_j=\frac{\mathrm{\λ}}{2(n+\mathrm{\λ})},j=2,...,2n+1$

其中：λ=n(α²-1)，α和β均为给定的计算参数；本实施方式中，α=0.5，β=1。

sigma矩阵X_k和U_k中的第j列向量X_k,j和U_k,j由以下关系式确定：

X_k,j＝x_k-1,j＝1

$X_{k,j}=x_{k-1}+\sqrt{n+\mathrm{\λ}}{M}_{k,j},j=2,...,n+1$

$X_{k,j}=x_{k-1}-\sqrt{n+\mathrm{\λ}}{M}_{k,j},j=n+2,...,2n+1$

U_k,j＝μ_k-1,j＝1

$U_{k,j}={\mathrm{\μ}}_{k-1}+\sqrt{n+\mathrm{\λ}}{N}_{k,j},j=2,...,n+1$

$U_{k,j}={\mathrm{\μ}}_{k-1}-\sqrt{n+\mathrm{\λ}}{N}_{k,j},j=n+2,...,2n+1$

其中：M_k,j为第k次迭代时的中间矩阵M_k中的第j列向量，N_k,j为第k次迭代时的中间矩阵N_k中的第j列向量；

中间矩阵M_k和N_k由以下关系式确定：

$M_{k,i,j}=\sqrt{A_{k,i,j}}$

A_k＝P_k-1+Q_v

P_k-1＝A_k-1-K_k-1E_k-1(K_k-1)^T

$N_{k,i,j}=\sqrt{B_{k,i,j}}$

B_k＝S_k-1+Q_s

S_k-1＝B_k-1-H_k-1T_k-1(H_k-1)^T

其中：M_k,i,j和A_k,i,j分别为中间矩阵M_k和中间矩阵A_k中第i行第j列的元素值，P_k-1为PET浓度分布向量x_k-1的协方差矩阵，N_k,i,j和B_k,i,j分别为中间矩阵N_k和中间矩阵B_k中第i行第j列的元素值，S_k-1为衰减系数分布向量μ_k-1的协方差矩阵，Q_v和Q_s均为系统噪声的协方差矩阵；

sigma矩阵Y_k和Z_k中的第j列向量Y_k,j和Z_k,j由以下关系式确定：

Y_k,j＝G_kDX_k,j Z_k,j＝W_k,jDx_k

其中：G_k为第k次迭代时的衰减矩阵，W_k,j为第k次迭代时的第j个对角矩阵；

衰减矩阵G_k和对角矩阵W_k,j由以下关系式确定：

G_k＝diag[exp(-Lμ_k-1)]

W_k,j＝diag[exp(-LU_k,j)]

其中：L为投影线与像素点的交叉长度矩阵；

根据上述迭代方程组进行迭代计算，迭代收敛后的PET浓度分布向量和衰减系数分布向量即为所要估计的PET浓度分布信息和衰减系数分布信息；迭代收敛条件如下：

$\frac{\underset{j=1}{\overset{n^2}{\mathrm{\Σ}}}|p_{k,j}-p_{k-1,j}|}{n^2}\≤\mathrm{\ρ}$

其中：p_k,j和p_k-1,j分别为PET浓度分布向量x_k和x_k-1的协方差矩阵P_k和P_k-1中第j个元素值，ρ为给定的收敛阈值；本实施方式中，ρ=10^-3。

以下我们采用Zubal数字胸腔体模来验证本实施方式放射性浓度分布和衰减系数分布同时估计框架的准确性和鲁棒性。Zubal数字胸腔体模分辨率为128×128，两个肺区域的平均衰减系数为0.04/cm，左肺的衰减系数是不均匀的而右肺的衰减系数为均匀的。为了使合成数据更加接近真实采集数据，我们在模拟单光子发射采集过程中加入了对衰减系数和探测器探测效率的模拟，数据中含有50%的随机和散射数据。我们共设计了两组实验：一组的示踪剂平均放射性浓度为50counts/pixel（低计数率条件），另外一组为200counts/Pixel（高计数率条件）。

对于这两组合成投影数据，放射性浓度分布图像和衰减系数分布图像分别使用如下3种重建方法重建：MLAA（maximum-likelihood reconstruction ofattenuation and activity，衰减和浓度分布的极大似然重建）、EKF和UKF（本实施方式）。

对于基于MLAA的重建方法，需要将衰减系数的先验知识加入到估计过程中，我们设定了两种不同情况下的先验知识，一组是完美条件下的先验知识，先验的阈值分别设定为0.05和0.01，这两个数值大约等同于分别在组织和肺中的衰减系数值；另一组是含有噪声的先验知识，即向第一组的完美模型先验值中加入大约30%的高斯噪声，最终结果来自大约100步的双EM（最大期望）迭代。对于本实施方式UKF求解框架，状态值（即放射性浓度分布）通过UKF状态估计器来估计得到，其中第一步的迭代中使用的衰减系数分布来自于初始的猜想值，后续的步骤中使用前面UKF参数估计步骤中估计得到的衰减系数分布值。接着，这此状态估计值被用在UKF参数估计中来估计衰减系数分布。对应估计值和真值进行一组统计分析，假设N_p为重建图像的总的像索数，

为最终的重建图像，x_tr是真值图像，那么我们可以定义如下的误差表示：

$\mathrm{RMSE}={\left(\frac{1}{N_p}\mathrm{\Σ}{(\hat{x}-x_{\mathrm{tr}})}^2\right)}^{0.5}$

高低两种计数率条件下重建得到图像误差RMSE值如表1所示。使用本实施方式UKF重建的RMSE值为0.2170（低计数率条件）和0.2120（高计数率条件），尽管不是很明显，还是比EKF重建出来的结果略低，EKF重建对应的RMSE值为0.2462（低计数率条件）和0.2200（高计数率条件）。EKF重建的不理想可以认为是在EKF线性化近似时对状态转移方程和系统测量方程的一阶Taylor近似精度不够造成的。

表1

重建之后的图像和统计结果显示，如果在有精确的先验信息的条件下MLAA方法能够产生最好的重建结果，但当先验衰减系数中含有噪声时，MLAA方法将产生明显的误差。从结果中我们可以看到当先验信急含有30%的噪声影响时，重建结果将有明显的下降，而对于本实施方式UKF估计框架来说，不论在低计数率条件下还是高计数率条件下都能重建出良好质量的放射性浓度分布图像。

Claims (9)

1.一种PET浓度与衰减系数的同时重建方法，包括如下步骤：

（1）利用探测器对注入有放射性物质的生物组织进行探测，采集得到PET的符合计数；

（2）根据PET成像原理，建立含有衰减因素的PET状态空间方程；

（3）利用无色卡尔曼滤波算法求解所述的PET状态空间方程，得到以下迭代方程组；进而根据符合计数通过以下迭代方程组估计出PET浓度分布信息和衰减系数分布信息；

x_k＝x_k-1+K_k(y-g_k)

μ_k＝μ_k-1+H_k(y-z_k)

其中：x_k和x_k-1分别为第k次和第k-1次迭代后的PET浓度分布向量，μ_k和μ_k-1分别为第k次和第k-1次迭代后的衰减系数分布向量，K_k和H_k均为第k次迭代时的增益矩阵，g_k和z_k均为第k次迭代时的中间向量，y为符合计数。

2.根据权利要求1所述的同时重建方法，其特征在于：所述的PET状态空间方程的表达式如下：

x_t+1＝x_t+v

μ_t+1＝μ_t+s

y＝GDx+e

其中：x_t和x_t+1分别为第t状态和第t+1状态下的PET浓度分布向量，μ_t和μ_t+1分别为第t状态和第t+1状态下的衰减系数分布向量，G为衰减矩阵，D为系统矩阵，v和s均为过程噪声向量；e为测量噪声向量，y为符合计数，x为PET浓度分布向量。

3.根据权利要求1所述的同时重建方法，其特征在于：所述的增益矩阵K_k根据以下算式求解：

K_k＝F_kE_k ^-1

$E_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{w}_j(Y_{k,j}-g_k){(Y_{k,j}-g_k)}^T+R_0$

$F_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{w}_j(X_{k,j}-x_{k-1}){(Y_{k,j}-g_k)}^T$

其中：X_k,j和Y_k,j分别为第k次迭代时的sigma矩阵X_k和Y_k中的第j列向量，R₀为测量噪声的协方差矩阵，w_j为对应第j列向量的权重因子，n为PET浓度分布向量的维度；

所述的中间向量g_k根据以下算式求解：

$g_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{u}_j{Y}_{k,j}$

其中：u_j为对应第j列向量的权重因子；

所述的增益矩阵H_k根据以下算式求解：

H_k＝V_kT_k ^-1

$V_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{w}_j(U_{k,j}-{\mathrm{\μ}}_{k-1}){(Z_{k,j}-z_k)}^T$

$T_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{w}_j(Z_{k,j}-z_k){(Z_{k,j}-z_k)}^T+R_0$

其中：U_k,j和Z_k,j分别为第k次迭代时的sigma矩阵U_k和Z_k中的第j列向量；

所述的中间向量z_k根据以下算式求解：

$z_k=\underset{j=1}{\overset{2n+1}{\mathrm{\Σ}}}{u}_j{Z}_{k,j}.$

4.根据权利要求3所述的同时重建方法，其特征在于：所述的权重因子u_j和w_j由以下算式求得：

$u_j=\frac{\mathrm{\λ}}{n+\mathrm{\λ}},j=1$

$u_j=\frac{\mathrm{\λ}}{2(n+\mathrm{\λ})},j=2,...,2n+1$

$w_j=\frac{\mathrm{\λ}}{n+\mathrm{\λ}}+(1-{\mathrm{\α}}^2+\mathrm{\β}),j=1$

$w_j=\frac{\mathrm{\λ}}{2(n+\mathrm{\λ})},j=2,...,2n+1$

其中：λ=n(α²-1)，α和β均为给定的计算参数。

5.根据权利要求3所述的同时重建方法，其特征在于：所述的sigma矩阵X_k和U_k中的第j列向量X_k,j和U_k,j由以下关系式确定：

X_k,j＝x_k-1,j＝1

$X_{k,j}=x_{k-1}+\sqrt{n+\mathrm{\λ}}{M}_{k,j},j=2,...,n+1$

$X_{k,j}=x_{k-1}-\sqrt{n+\mathrm{\λ}}{M}_{k,j},j=n+2,...,2n+1$

U_k,j＝μ_k-1,j＝1

$U_{k,j}={\mathrm{\μ}}_{k-1}+\sqrt{n+\mathrm{\λ}}{N}_{k,j},j=2,...,n+1$

$U_{k,j}={\mathrm{\μ}}_{k-1}-\sqrt{n+\mathrm{\λ}}{N}_{k,j},j=n+2,...,2n+1$

其中：λ=n(α²-1)，α为给定的计算参数，M_k,j为第k次迭代时的中间矩阵M_k中的第j列向量，N_k,j为第k次迭代时的中间矩阵N_k中的第j列向量。

6.根据权利要求5所述的同时重建方法，其特征在于：所述的中间矩阵M_k和N_k由以下关系式确定：

$M_{k,i,j}=\sqrt{A_{k,i,j}}$

A_k＝P_k-1+Q_v

P_k-1＝A_k-1-K_k-1E_k-1(K_k-1)^T

$N_{k,i,j}=\sqrt{B_{k,i,j}}$

B_k＝S_k-1+Q_s

S_k-1＝B_k-1-H_k-1T_k-1(H_k-1)^T

其中：M_k,i,j和A_k,i,j分别为中间矩阵M_k和中间矩阵A_k中第i行第j列的元素值，P_k-1为PET浓度分布向量x_k-1的协方差矩阵，N_k,i,j和B_k,i,j分别为中间矩阵N_k和中间矩阵B_k中第i行第j列的元素值，S_k-1为衰减系数分布向量μ_k-1的协方差矩阵，Q_v和Q_s均为系统噪声的协方差矩阵。

7.根据权利要求3所述的同时重建方法，其特征在于：所述的sigma矩阵Y_k和Z_k中的第j列向量Y_k,j和Z_k,j由以下关系式确定：

Y_k,j＝G_kDX_k,j

Z_k,j＝W_k,jDx_k

其中：D为系统矩阵，G_k为第k次迭代时的衰减矩阵，W_k,j为第k次迭代时的第j个对角矩阵。

8.根据权利要求7所述的同时重建方法，其特征在于：所述的衰减矩阵G_k和对角矩阵W_k,j由以下关系式确定：

G_k＝diag[exp(-Lμ_k-1)]

W_k,j＝diag[exp(-LU_k,j)]

其中：L为投影线与像素点的交叉长度矩阵。

9.根据权利要求1所述的同时重建方法，其特征在于：所述的步骤（3）中，根据符合计数通过迭代方程组进行迭代计算，则迭代收敛后的PET浓度分布向量和衰减系数分布向量即为所要估计的PET浓度分布信息和衰减系数分布信息；所述的迭代收敛条件如下：

$\frac{\underset{j=1}{\overset{n^2}{\mathrm{\Σ}}}|p_{k,j}-p_{k-1,j}|}{n^2}\≤\mathrm{\ρ}$

其中：p_k,j和p_k-1,j分别为PET浓度分布向量x_k和x_k-1的协方差矩阵P_k和P_k-1中第j个元素值，ρ为给定的收敛阈值，n为PET浓度分布向量的维度。

Images