CN103230282A - 一种pet浓度均值与方差的估计方法及系统 - Google Patents

Abstract

本发明公开了一种PET浓度均值与方差的估计方法，包括：（1）采集、校正，得到符合技术；（2）建立状态空间方程；（3）卡尔曼滤波求解；（4）重建PET浓度均值与方差。本发明通过对状态空间体系进行卡尔曼滤波求解，进而对PET浓度均值与方差进行估计；对于包含了不同噪声水平的放射性浓度测量值，在给出PET重建图像的同时，也得到不同噪声水平下的误差程度；本发明还公开了一种基于上述方法的PET浓度均值与方差的估计系统，与现有重建方法相比，可确定重建结果的可靠性程度，为PET重建图像的进一步分析提供了准确依据。

Description

一种PET浓度均值与方差的估计方法及系统

技术领域

本发明属于PET成像技术领域，具体涉及一种PET浓度均值与方差的估计方法及系统。

背景技术

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

目前，PET图像重建方法大致可分为两类：解析法和迭代统计法。前一类主要是滤波反投影法，计算速度快，但伪影严重。为此，以最大似然法为代表的迭代统计法被提出来，由于迭代法基于统计学模型，对不完全数据适应性好，逐渐成为PET重建算法研究关注点。状态空间体系从一个全新的角度提供了PET图像重建的新思路，通过根据实际求解问题调整测量方程与状态方程的统一表达式，以实现静态、动态重建以及先验估计。通过相关方法求解，如卡尔曼滤波、H∞滤波、粒子滤波等可适应不同噪声特性、不同运算条件，与传统的解析法或迭代统计法相比有一定的优势。

然而，PET图像重建作为一种复杂的测量，也具有其特定的误差属性。PET图像的浓度测量值总是包含了一定的噪声，因而重建结果也具有一定的不确定性，误差属性则反映重建结果的可靠程度，为PET重建图像的进一步分析提供了准确依据。同其它种类的测量一样，如果不能确定其误差程度，PET图像重建的意义将大打折扣。

目前针对PET图像重建的研究，大都集中在对重建算法的设计上，而很少设计重建结果的误差程度。部分学者也对Map迭代算法重建图像方差的进行分析，如作者J.A.Fessler在标题为Variance images for penalized-likelihood image reconstruction（IEEE.Nuclear Science Symposium,1997）的文献中以及J.Qi和R.M.Leahy在标题为A theoretical study of the contrast recovery and variance of map reconstructions from pet data（IEEE Transactions on Medical Imaging,18(4):293–305,1999）的文献中均提出了各自关于PET均值和方差的估计方法，这些方法他们主要基于信息矩阵FIM，立足于重建图像，对重建方差进行推导；并且求解方差过程与图像重建过程相分离，其应用只局限于对重建算法本身的评价，而不是用于评价PET图像重建的误差属性。然而在状态空间体系中，通过求解，可以同时对PET浓度均值与方差进行估计，对于包含了不同噪声水平的放射性浓度测量值，在给出PET重建图像的同时，也确定不同噪声水平下的误差程度。

发明内容

针对现有技术所存在的上述技术问题，本发明提出了一种PET浓度均值与方差的估计方法及系统，能够在获取PET重建图像的同时，确定其误差程度。

一种PET浓度均值与方差的估计方法，包括如下步骤：

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

（2）根据PET成像原理，建立PET的状态空间方程；

（3）利用卡尔曼滤波算法求解所述的状态空间方程得到以下迭代方程组；根据校正后的符合计数通过以下迭代方程组估计出PET浓度的均值和方差；

x_i=z_i+K_i(y-Dz_i)

P_i=H_i-K_iDH_i

其中：x_i为第i次迭代后的PET浓度分布向量，z_i为第i次迭代后的PET浓度分布预估向量，K_i为第i次迭代后的增益矩阵，D为系统矩阵，y为校正后的符合计数，P_i为第i次迭代后的误差协方差矩阵，H_i为第i次迭代后的误差协方差预估矩阵，i为大于0的自然数。

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

y=Dx+e

x_t+1=Ax_t+v

其中：D为系统矩阵；y为校正后的符合计数，且为m维向量；x为PET浓度分布向量，且为n维向量；e为测量噪声向量；x_t+1为第t+1状态下的PET浓度分布向量，x_t为第t状态下的PET浓度分布向量，A为状态转移矩阵，v为过程噪声向量，m、n和t均为大于1的自然数。

所述的步骤（3）中，根据校正后的符合计数通过迭代方程组进行迭代计算，则迭代收敛后的PET浓度分布向量即为PET浓度的均值，迭代收敛后的误差协方差矩阵中对角线元素组成的向量即为PET浓度的方差。

所述的第i次迭代后的PET浓度分布预估向量z_i根据以下算式求得：

z_i=Ax_i-1

其中：x_i-1为第i-1次迭代后的PET浓度分布向量，A为状态转移矩阵。

所述的第i次迭代后的增益矩阵K_i根据以下算式求得：

K_i=P_i-1D^T[DP_i-1D^T+E]^-1

其中：P_i-1为第i-1次迭代后的误差协方差矩阵，D为系统矩阵，E为测量噪声矩阵。

所述的第i次迭代后的误差协方差预估矩阵H_i根据以下算式求得：

H_i=AP_i-1A^T+V

其中：P_i-1为第i-1次迭代后的误差协方差矩阵，A为状态转移矩阵，V为过程噪声矩阵。

所述的迭代收敛条件如下：

$\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}|p_i^j-p_{i-1}^j|}{N}\text{\≤\ρ}$

其中：

为第i次迭代后的误差协方差矩阵中第j元素值，为第i-1次迭代后的误差协方差矩阵中第j元素值，N为误差协方差矩阵中的元素总个数，ρ为给定的收敛阈值，j为自然数且1≤j≤N。

所述的测量噪声向量e为m维向量，代表符合事件采集过程中由随机事件、散射事件等引入的测量误差；过程噪声向量v为n维向量，表征了状态之间转移的统计不确定性；状态转移矩阵A为n×n维单位矩阵；系统矩阵D为m×n维矩阵，表征了发射光子被探测器接收的概率，其受探测器结构、探测效率、衰减、死时间等因素的影响；过程噪声矩阵V为过程噪声向量v的协方差矩阵，且为n×n维矩阵；测量噪声矩阵E为测量噪声向量e的协方差矩阵，且为m×m维矩阵；误差协方差矩阵P和误差协方差预估矩阵H均为n×n维矩阵。

一种PET浓度均值与方差的估计系统，包括探测器和与探测器相连的计算机；

所述的探测器用于对注入有放射性物质的生物组织进行探测，采集得到符合计数；

所述的计算机内加载有以下功能模块：

数据接收模块，用于接收所述的符合计数并对其进行校正；

均值方差估计模块，用于根据校正后的符合计数通过预设的迭代方程组估计出PET浓度的均值和方差。

所述的均值方差估计模块根据以下迭代方程组估计PET浓度的均值和方差：

x_i=z_i+K_i(y-Dz_i)

P_i=H_i-K_iDH_i

z_i=Ax_i-1

K_i=P_i-1D^T[DP_i-1D^T+E]^-1

H_i=AP_i-1A^T+V

其中：x_i为第i次迭代后的PET浓度分布向量，x_i-1为第i-1次迭代后的PET浓度分布向量，z_i为第i次迭代后的PET浓度分布预估向量，K_i为第i次迭代后的增益矩阵，D为系统矩阵，A为状态转移矩阵，V为过程噪声矩阵，E为测量噪声矩阵，y为校正后的符合计数，P_i为第i次迭代后的误差协方差矩阵，P_i-1为第i-1次迭代后的误差协方差矩阵，H_i为第i次迭代后的误差协方差预估矩阵，i为大于0的自然数；迭代收敛后的PET浓度分布向量即为PET浓度的均值，迭代收敛后的误差协方差矩阵中对角线元素组成的向量即为PET浓度的方差。

本发明的有益技术效果在于：通过对状态空间体系进行卡尔曼滤波求解，进而对PET浓度均值与方差进行估计；对于包含了不同噪声水平的放射性浓度测量值，在给出PET重建图像的同时，也得到不同噪声水平下的误差程度；与现有重建方法相比，可确定重建结果的可靠性程度，为PET重建图像的进一步分析提供了准确依据。

附图说明

图1为本发明估计系统的结构示意图。

图2为本发明估计方法的步骤流程示意图。

图3为Zubal Phantom体模型的示意图

图4(a)为Zubal Phantom体模的PET浓度真值图像。

图4(b)为Zubal Phantom体模的PET浓度均值图像。

图4(c)为Zubal Phantom体模的PET浓度标准差图像。

图5(a)为本发明PET重建结果第32行像素误差范围曲线图。

图5(b)为本发明PET重建结果第32列像素误差范围曲线图。

具体实施方式

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

如图1所示，一种PET浓度均值与方差的估计系统，包括探测器和与探测器相连的计算机；系统的操作执行流程图如图2所示；其中：

探测器用于对注入有放射性物质的生物组织进行探测，采集当前时刻的符合计数；本实施方式中，探测器采用日本滨松公司生产的型号为SHR74000的PET扫描仪。

计算机内加载有数据接收模块和均值方差估计模块；其中：

数据接收模块用于接收符合计数，并对其进行校正；正电子发射断层扫描仪探测人体内发出的放射性信号，经过符合和采集系统处理，形成原始符合事件。PET探测器记录的符合事件包括真符合、随机符合和散射符合。通过探测器的延时窗口和能量窗口对随机事件和散射事件进行校正，而后进行衰减校正，得到正弦图数据即校正后的符合计数y。

均值方差估计模块用于根据符合计数y通过模块内预设的迭代方程组估计出当前时刻PET浓度的均值和方差。

根据PET成像原理，可以得到PET的状态空间方程的表达式如下：

y=Dx+e

x_t+1=Ax_t+v

其中：D为系统矩阵；y为校正后的符合计数，且为m维向量；x为PET浓度分布向量，且为n维向量；e为测量噪声向量；x_t+1为第t+1状态下的PET浓度分布向量，x_t为第t状态下的PET浓度分布向量，A为状态转移矩阵，v为过程噪声向量。

利用卡尔曼滤波算法求解上述状态空间方程可得到以下迭代方程组；

x_i=z_i+K_i(y-Dz_i)

P_i=H_i-K_iDH_i

z_i=Ax_i-1

K_i=P_i-1D^T[DP_i-1D^T+E]^-1

H_i=AP_i-1A^T+V

其中：x_i为第i次迭代后的PET浓度分布向量，x_i-1为第i-1次迭代后的PET浓度分布向量，z_i为第i次迭代后的PET浓度分布预估向量，K_i为第i次迭代后的增益矩阵，D为系统矩阵，A为状态转移矩阵，V为过程噪声矩阵，E为测量噪声矩阵，P_i为第i次迭代后的误差协方差矩阵，P_i-1为第i-1次迭代后的误差协方差矩阵，H_i为第i次迭代后的误差协方差预估矩阵；PET浓度分布向量x₀和误差协方差矩阵P₀均为随机初始化生成。

根据符合计数y通过以上迭代方程组进行迭代计算，则迭代收敛后的PET浓度分布向量即为PET浓度的均值，迭代收敛后的误差协方差矩阵中对角线元素组成的向量即为PET浓度的方差。迭代收敛条件如下：

$\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}|p_i^j-p_{i-1}^j|}{N}\text{\≤\ρ}$

其中：

为第i次迭代后的误差协方差矩阵中第j元素值，

为第i-1次迭代后的误差协方差矩阵中第j元素值，N为误差协方差矩阵中的元素总个数，ρ为给定的收敛阈值；本实施方式中ρ设为10^-3。

以下我们采用Zubal Phantom体模进行蒙特卡罗模拟实验，以验证本实施方式的有效性。如图3所示，该体模包含3个组织区域和一个背景区域。测量数据的采集利用蒙特卡罗模拟方法实现，实验中所模拟的PET探测器为日本滨松公司生产的SHR74000，并将符合事件存储为64×64的正弦图数据。

实验模拟了真实环境下PET探测器对放射性符合事件的2D采集，其中包括0.02%的随机事件和8.8%的散射事件；计数总数为4,000,000。

在利用卡尔曼滤波方法进行重建时，需要给定迭代的初始参数。具体参数设置如下：

A=I，I为n×n维的单位矩阵；P₀为对角矩阵且P₀为=diag(x_fbp x_fbp ^T)，其中x_fbp为利用FBP（滤波反投影）方法重建的结果。

V=αI，其中α=1×10^-5，则V接近于零矩阵；PET浓度分布向量x₀=0，符合计数y为m维向量，n=4096，m=4096；最大迭代次数为20。

利用本实施方式重建得到PET浓度均值x和方差var；使均值x、标准差与PET浓度真值进行对比，得到图4所示的对比示意图。

取

作为重建图像的误差，得到各重建像素值的误差范围。图5为各噪声等级下误差范围曲线图（横坐标VOXEL INDEX表示像素编号，纵坐标VOXEL VALUE表示像素值），图5(a)和(b)分别取自重建图像横向第32列和纵向第32列。可以看出，真值（True Value）曲线基本落在误差范围的上下界限之内（Up bounds代表误差上限，Lower bounds代表误差下限），说明本实施方式估计出的方差能准确反映PET图像重建的误差情况。

本实施方式对PET图像的均值和方差同时进行估计，与传统的仅对图像均值进行重建的方法相比，本实施方式既提供了PET图像的重建均值，又提供了重建的误差程度，为PET重建图像的进一步分析提供了准确依据。

Claims (9)

1.一种PET浓度均值与方差的估计方法，包括如下步骤：

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

（2）根据PET成像原理，建立PET的状态空间方程；

（3）利用卡尔曼滤波算法求解所述的状态空间方程得到以下迭代方程组；根据校正后的符合计数通过以下迭代方程组估计出PET浓度的均值和方差；

x_i=z_i+K_i(y-Dz_i)

P_i=H_i-K_iDH_i

其中：x_i为第i次迭代后的PET浓度分布向量，z_i为第i次迭代后的PET浓度分布预估向量，K_i为第i次迭代后的增益矩阵，D为系统矩阵，y为校正后的符合计数，P_i为第i次迭代后的误差协方差矩阵，H_i为第i次迭代后的误差协方差预估矩阵，i为大于0的自然数。

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

y=Dx+e

x_t+1=Ax_t+v

其中：D为系统矩阵；y为校正后的符合计数，且为m维向量；x为PET浓度分布向量，且为n维向量；e为测量噪声向量；x_t+1为第t+1状态下的PET浓度分布向量，x_t为第t状态下的PET浓度分布向量，A为状态转移矩阵，v为过程噪声向量，m、n和t均为大于1的自然数。

3.根据权利要求1所述的估计方法，其特征在于：所述的步骤（3）中，根据校正后的符合计数通过迭代方程组进行迭代计算，则迭代收敛后的PET浓度分布向量即为PET浓度的均值，迭代收敛后的误差协方差矩阵中对角线元素组成的向量即为PET浓度的方差。

4.根据权利要求1所述的估计方法，其特征在于：所述的第i次迭代后的PET浓度分布预估向量z_i根据以下算式求得：

z_i=Ax_i-1

其中：x_i-1为第i-1次迭代后的PET浓度分布向量，A为状态转移矩阵。

5.根据权利要求1所述的估计方法，其特征在于：所述的第i次迭代后的增益矩阵K_i根据以下算式求得：

K_i=P_i-1D^T[DP_i-1D^T+E]^-1

其中：P_i-1为第i-1次迭代后的误差协方差矩阵，D为系统矩阵，E为测量噪声矩阵。

6.根据权利要求1所述的估计方法，其特征在于：所述的第i次迭代后的误差协方差预估矩阵H_i根据以下算式求得：

H_i=AP_i-1A^T+V

其中：P_i-1为第i-1次迭代后的误差协方差矩阵，A为状态转移矩阵，V为过程噪声矩阵。

7.根据权利要求3所述的估计方法，其特征在于：所述的迭代收敛条件如下：

$\frac{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}|p_i^j-p_{i-1}^j|}{N}\text{\≤\ρ}$

其中：

为第i次迭代后的误差协方差矩阵中第j元素值，

为第i-1次迭代后的误差协方差矩阵中第j元素值，N为误差协方差矩阵中的元素总个数，ρ为给定的收敛阈值，j为自然数且1≤j≤N。

8.一种PET浓度均值与方差的估计系统，包括探测器和与探测器相连的计算机；所述的探测器用于对注入有放射性物质的生物组织进行探测，采集得到符合计数；

所述的计算机内加载有以下功能模块：

数据接收模块，用于接收所述的符合计数并对其进行校正；

均值方差估计模块，用于根据校正后的符合计数通过预设的迭代方程组估计出PET浓度的均值和方差。

9.根据权利要求8所述的估计系统，其特征在于：所述的均值方差估计模块根据以下迭代方程组估计PET浓度的均值和方差：

x_i=z_i+K_i(y-Dz_i)

P_i=H_i-K_iDH_i

z_i=Ax_i-1

K_i=P_i-1D^T[DP_i-1D^T+E]^-1

H_i=AP_i-1A^T+V

其中：x_i为第i次迭代后的PET浓度分布向量，x_i-1为第i-1次迭代后的PET浓度分布向量，z_i为第i次迭代后的PET浓度分布预估向量，K_i为第i次迭代后的增益矩阵，D为系统矩阵，A为状态转移矩阵，V为过程噪声矩阵，E为测量噪声矩阵，y为校正后的符合计数，P_i为第i次迭代后的误差协方差矩阵，P_i-1为第i-1次迭代后的误差协方差矩阵，H_i为第i次迭代后的误差协方差预估矩阵，i为大于0的自然数；迭代收敛后的PET浓度分布向量即为PET浓度的均值，迭代收敛后的误差协方差矩阵中对角线元素组成的向量即为PET浓度的方差。

Images