CN101627919A - 有限采样角度下基于卡尔曼滤波的pet浓度重建方法 - Google Patents

Abstract

本发明公开了一种有限采样角度下基于卡尔曼滤波的PET浓度重建方法，通过PET正电子发射断层扫描仪得到非完全的原始投影线的正弦图数据，结合房室模型理论建立状态空间体系，最后通过基于状态空间理论的卡尔曼滤波法求解，得出放射性活度分布，即实现重建图像。由于房室模型的先验引导以及卡尔曼滤波法对不完全数据具有很强的适应性，重建的PET图像质量并没有受到数据丢失的影响。通过与现有重建方法MLEM的实验比较，定性和定量分析结果都表明本发明方法具有一定的优越性。

Description

有限采样角度下基于卡尔曼滤波的PET浓度重建方法

技术领域

本发明涉及一种正电子发射断层影像重建方法，尤其是涉及一种在有限采样角度条件下的PET成像中的卡尔曼滤波图像重建方法。

背景技术

正电子发射断层成像(Positron emission tomography，PET)是核医学成像装置中最为重要的应用之一，已成为活体功能成像的重要手段。PET的基本思想是：向生物体内部注射同位素标记的化合物，同位素不稳定通过衰变发射正电子与人体内的电子发生湮灭，发射一对方向几乎相反的光子，通过在体外探测这对光子得到药物的放射性浓度分布。如今许多临床医学领域已经开始广泛使用PET图像进行肿瘤、心脏疾病、神经和精神系统疾病的诊断，并且在药物开发过程中也发挥着越来越重要的作用。

虽然PET在医学领域中已被普遍应用，但是其高昂的费用(主要包括回旋加速器和探测器的费用)成为它进一步发展的主要瓶颈。与单光子发射断层成像(single photon emission computed tomography，SPECT)技术相比PET拥有更高的分辨率，但是探测器的成本却远远大于SPECT。典型的PET探测器由许多独立的晶体模块组成一个圆桶状的多层探测环环绕被测物体，采集的角度覆盖0°～180°范围。通过180°内的完全投影得到的正弦图数据结构非常庞大，给后期的计算机处理和图像重建也带来许多负担。所以研究有限采样角度下的PET重建问题是具有实际意义的，不仅能够有效缩减探测器成本，也可以大大减少系统的后期处理时间，同时对于采用符合探测的SPECT的进一步研究也将具有重要的启发意义。

与传统的环形探测器得到的数据相比，有限采样角度下得到的正弦图是不完整的，存在着严重的数据丢失问题。对于现有的重建算法，包括滤波反投影法(Filtered back projection，FBP)、期望最大极大似然法(Maximum likelihood expectation maximization，MLEM)、最小二乘法(Least square，LS)以及最大后验法(Maximum a posteriori，MAP)往往不能得到最优的结果。在信息缺失的情况下，这些原本在完全采样角度下被采用的方法得到的重建图像将存在严重的缺陷和伪影。虽然在计算机断层成像(Computed tomography，CT)中有少量对有限采样角度下重建问题的研究，但是目前在PET成像领域对这个问题还没有针对的解决方法。

在对PET重建问题的研究中，已经成功地把基于状态空间理论的Kalman滤波法运用到完全采样角度条件下的PET重建中。由于该方法使用了房室模型作为先验引导，其计算效果与传统的解析法或者迭代法相比具有一定的优越性。Kalman滤波是基于最小均方误差为估计的最佳准则的递推估计的算法，其基本思想是：利用前一时刻的估计值和现时刻的观测值来更新对状态变量的估计，循环迭代直至求出最优估计。Kalman滤波法尤其适合于实时处理和计算机运算，同时对不完全数据也具有很好的适应性。当在不完全角度采样条件下，此方法也将具有良好的表现。

发明内容

本发明提出一种有限采样角度下基于卡尔曼滤波的PET浓度重建方法，在投影数据不完全的情况下能很好地重建PET浓度。

一种有限采样角度下基于卡尔曼滤波的PET浓度重建方法，包括以下步骤：

(1)在有限采样角度下PET正电子发射断层扫描仪采集生物体中的光子信号数据，经过符合判断和采集系统处理，形成原始投影线，校正后得到正弦图数据Y；

所述的光子信号是利用现有医学技术，由注射进入体内的同位素标记的化合物发出。

所述的有限采样角度，是相对于现有技术中采集的角度覆盖0°～180°范围而言的。在模拟试验中，因为可以与原图进行对比，所以有限采样角度只选区了一个范围，即30°＜θ＜150°。在真实实验中，由于缺少与原图的对比，为了说明采样角度减小对图像重建效果的影响，本发明所述的有限采样角度设置在两个范围内，分别是30°＜θ＜150°范围以及45°＜θ＜135°范围，角度覆盖分别达到120°以及90°。

所述的将采集的光子信号数据转换为正弦图数据Y的过程可以利用现有技术中PET装置及处理系统来完成，区别仅在于本发明方法的采样角度是有限角度。

所述的正弦图数据Y是通过PET正电子发射断层扫描仪在有限角度下采集得到的。PET正电子发射断层扫描仪进行透射扫描和发射扫描，透射扫描得到图像的衰减校正系数，发射扫描在非完全采样角度内进行得到不完全的原始正弦图数据。使用衰减校正系数对原始正弦图数据进行校正得到正弦图数据Y。

(2)根据正弦图数据Y以及药物动力学理论建立一个由离散的观测方程(1)和连续的发展方程(2)构成的状态空间方程：

Y(t_k)＝DX(t_k)+v(t_k)......(1)

X′(t)＝AX(t)+Bμ(t)+ω(t)......(2)

其中：

t表示时间，t_k表示离散采样时间点，下标k＝1，2，3…(自然数)；

Y是正弦图数据；

D为系统矩阵，表示发射的光子被探测器(PET装置的探测器)接收到的概率，由PET装置固有特性所决定；

X为放射性浓度的空间分布，是需要重建的对象；

X`为X的导数；

A是状态转移矩阵，B是输入矩阵，两者都由药物交换的动力学常数组成，由药物动力学理论中房室模型的结构所决定。一般不同的模型对应不同的组合，并且是已知的；

μ是血输入函数，由动脉血采样得到，在此我们可以采用经验公式 $\mathrm{\μ}\left(t\right)=(A_{1}t-A_2-A_3)e^{-{\mathrm{\λ}}_{1}t}+A_2{e}^{-{\mathrm{\λ}}_{2}t}+A_3{e}^{-{\mathrm{\λ}}_{3}t},$ 其中A₁、A₂、A₃以及λ₁、λ₂、λ₃，都是已知常量；

v为测量噪声是各种噪声的集合，服从正态高斯分布v(t)～N(0，R(t))；

ω是过程噪声，服从正态高斯分布ω(t)～N(0，Q(t))。

(3)基于步骤(2)得到的状态空间方程，利用卡尔曼滤波算法计算得到放射性浓度X值。

利用卡尔曼滤波算法计算时，采用方程(3)～(7)迭代重建实现。

${\hat{x}}^{\′}=\mathrm{Ax}\left(t\right)+\mathrm{B\μ}\left(t\right)---\left(3\right)$

P′(t)＝AP(t)+P(t)A+Q(t)          (4)

K(t_k)＝P(t_k ^-)D^T[DP(t_k ^-)D+R(t_k)]^-1 (5)

$\hat{x}\left(t_k\right)=\hat{x}\left({t_k}^{-}\right)+K\left(t_k\right)[y\left(t_k\right)-D\hat{x}\left({t_k}^{-}\right)]---\left(6\right)$

P(t_k)＝[I-K(t_k)D]P(t_k ^-)           (7)

其中，x和y与步骤(2)中的X以及Y具有相同含义，但为了递推表达方便在此使用小写形式。

D^T是系统矩阵D的转置矩阵。P是放射性浓度的空间分布X的误差协方差，P`是P的导数。

P(t_k ^-)分别是X、P的先验估计值，其中t_k ^-表示先验时刻；

和P(t_k)分别是X、P的后验估计值。

Q(t)以及R(t_k)分别表示该时刻的过程噪声以及测量噪声的方差。

假设X₀为方程(3)的初始值，P₀为方程(4)的初始值。下标k表示迭代的次数，当k＝1时，迭代从初始X₀、P₀出发，首先通过方程(3)和(4)计算得到浓度X的先验估计

以及误差协方差P的先验估计P(t_k ^-)，然后根据方程(5)计算出增益矩阵K，最后通过观测值y(t_k)来修正先验估计

和P(t_k ^-)从而得到后验估计值

和P(t_k)。

具体的迭代重建包括以下步骤：

1)首先设定放射性浓度的初始值和初始协方差X₀、P₀。初始值可以任意设定，但原则上不要选择0值。迭代从下标k＝1开始进行；

2)利用方程(3)和方程(4)推算出放射性浓度的估计值

和误差协方差的估计值P(t_k ^-)；

3)根据方程(5)计算出Kalman增益矩阵K；

4)利用离散采集时间点测得的正弦图数据y(t_k)，根据状态更新方程(6)、(7)来更新先验值和P(t_k ^-)，得到这一时刻的后验估计值和P(t_k)；

5)一次迭代完成以后，k值递增加1。重复交替步骤2)-4)，直到X值收敛，获得放射性浓度的最优估计值。

本发明的优点是：

对于有限采样角度条件下的PET重建问题，基于状态空间体系的Kalman滤波法由于使用了房室模型作为发展方程来约束浓度重建过程，使得重建结果具有生理学意义；并且Kalman法本身对不完全数据具有很好的适应能力，所以重建的PET图像质量并没有受到数据丢失的影响，结果比较理想；通过与现有重建方法的实验比较，重建结果的偏差和方差都说明了这种新方法的优越性。

附图说明

图1是模拟实验使用的胸腔体模示意图；

图2是MLEM法和利用本发明方法对图1的重建结果示意图；

图3是MLEM法和利用本发明方法重建的时间浓度曲线与真实曲线的对比效果图；

图4是真实实验使用的体模示意图；

图5是MLEM法和利用本发明方法对图4的重建结果示意图。

具体实施方式

正电子发射断层扫描仪探测人体内发出的光子信号，经过符合判断和采集系统处理，形成投影线，并以正弦图方式存放于计算机硬盘中。计算机以正弦图为输入，调用图像重建模块，计算得出人体内药物放射性浓度分布图像。

本发明方法的具体实施步骤如下：

1、在应用PET正电子发射断层扫描仪时，进行透射扫描和发射扫描。透射扫描得到图像的衰减校正系数。

在模拟试验中，因为可以与原图进行对比，所以有限采样角度只选区了一个范围，即30°＜θ＜150°。在真实实验中，由于缺少与原图的对比，为了说明采样角度减小对图像重建效果的影响，本发明所述的有限采样角度设置在两个范围内，分别是30°＜θ＜150°范围以及45°＜θ＜135°范围，角度覆盖分别达到120°以及90°。发射扫描选取有限个采样角度和径向单元进行，得到不完全的原始正弦图数据，其中系统矩阵D由具体使用的仪器性能所决定，为已知量。

2、对原始采集到的正弦图数据进行各类校正得到Y，从而建立一个离散的测量方程：

Y(t_k)＝DX(t_k)+v(t_k)

其中，Y是正弦图数据，D是系统矩阵，v为测量噪声是各种噪声的集合，t_k表示离散的采样时间点，下标k＝1，2，3...。

3、建立药物发展方程：根据药物动力学理论选择相应的房室模型描述药物在人体内的时间和空间分布过程，建立一个一阶微分方程，并且加入过程噪声ω，可得到一般的表达：

X′(t)＝AX(t)+Bμ(t)+ω(t)

其中，X为放射性浓度状态变量，为需要重建的对象，X`为X的导数；A是状态转移矩阵，B是输入矩阵，两者都由药物交换的动力学常数组成；μ是血输入函数，由动脉血采样得到，在此我们采用经验公式 $\mathrm{\μ}\left(t\right)=(A_{1}t-A_2-A_3)e^{-{\mathrm{\λ}}_{1}t}+A_2{e}^{-{\mathrm{\λ}}_{2}t}+A_3{e}^{-{\mathrm{\λ}}_{3}t},$ 其中A₁、A₂、A₃以及λ₁、λ₂、λ₃，都是已知常量。

4、构建状态空间体系：把离散的测量方程和连续的发展方程联合起来，建立一个混合的状态空间体系：

Y(t)＝DX(t)+v(t)

X(t)＝AX(t)+Bμ(t)+ω(t)

其中，v服从正态高斯分布v(t)～N(0，R(t))，ω服从正态高斯分布ω(t)～N(0，Q(t))。

5、利用卡尔曼滤波进行迭代，根据下列方程重建图像：

${\hat{x}}^{\′}\left(t\right)=\mathrm{Ax}\left(t\right)+\mathrm{B\μ}\left(t\right)$                              ①

P′(t)＝AP(t)+P(t)A+Q(t)           ②

K(t_k)＝P(t_k ^-)D^T[DP(t_k ^-)D+R(t_k)]^-1  ③

$\hat{x}\left(t_k\right)=\hat{x}\left({t_k}^{-}\right)+K\left(t_k\right)[y\left(t_k\right)-D\hat{x}\left({t_k}^{-}\right)]$                                              ④

P(t_k)＝[I-K(t_k)D]P(t_k ^-)            ⑤

其中，x和y是X以及Y的不同表达形式，为了递推表达方便在此使用小写形式。假设X₀为方程①的初始值，P₀为方程②的初始值。下标k表示迭代的次数，当k＝1时，迭代从初始X₀、P₀出发，通过观测值y(t_k)，不断修正放射性浓度值X，最终得出放射性活度分布，重建出图像。

采用Kalman滤波进行图像重建时，主要分成四个步骤：

1)首先设定放射性浓度的初始值和初始协方差X₀、P₀。初始值可以任意设定，但原则上不要选择0值。迭代从下标k＝1开始进行；

2)利用时间更新方程方程①和②向前推算出浓度的先验估计值P(t_k ^-)和误差协方差的先验估计值

3)根据方程③计算出Kalman增益矩阵K；

4)利用离散采集时间点测得的观测值y(t_k)，根据状态更新方程④、⑤来更新先验估计值，得到这一时刻的后验估计值

和P(t_k)；

5)一次迭代完成以后，k值递增加1。重复交替步骤2)-4)直至X值收敛，获得最优浓度估计值。

本发明技术的实验结果如下：

首先应用本发明技术进行计算机模拟实验，且和MLEM法的重建结果作比较。采用Zubal胸腔体模合成发射扫描数据，如图1所示。体模分成3个感兴趣区域，每个区域的药物分布由一组给定的动力学参数和相应的二组织房室模型来描述。通过求解房室方程得到真实的时间浓度分布曲线。模拟采样时间为60分钟，分成18个时间帧，包括30s×4，2min×4，5min×10。图像的原始分辨率为32×32像素，正弦图模拟生成过程中采用120°旋转角度下32个采样角度和34个径向采样的投影数据。为了模拟真实的采集条件，正弦图数据中加入了30％的泊松噪声，以此作为PET探测到的数据。

图2是用MLEM法和本发明的方法得到的重建结果，第一行表示在第12、14、16以及18个采样时刻的浓度真实分布情况，第二行表示MLEM的重建结果，第三行是Kalman法的结果。三个区域的平均时间浓度曲线与真实值的对比如图3所示。

进一步我们计算了两种算法的重建结果的平均偏差和方差列于表1中，MLEM法的在最高偏差达到0.0346，而本发明的偏差为0.0223；MLEM法的最高方差达到0.0099，而本发明的方差为0.0062。从图2种可以明显地看出本方法依旧能很好地还原三个感兴趣区域而MLEM法的重建结果非常模糊。图3中时间曲线的吻合情况也同样说明了本发明优于传统的MLEM法。

表1.MLEM法&Kalman法重建结果定量分析

为进一步说明本发明提出的有限角度下PET浓度重建的Kalman滤波法的有效性，我们给出临床PET系统采集得到的数据进行重建研究的结果。图4是所用的体模示意图，在六个小球体内灌注F¹⁸溶液，周围充满纯净水，一共采集25帧，分别为5×10s，5×30s，5×60s，5×120s，和5×180s。图像的分辨率为64×64像素，正弦图生成过程中采用30°＜θ＜150°和45°＜θ＜135°范围内的有效数据来模拟有限角度采样条件。图5是在上述两个采样条件下分别利用MLEM法和Kalman法对22，23，24，25帧数据的重建结果，第一行和第二行表示在30°＜θ＜150°采样条件下MLEM法和Kalman法的重建图像，第三行和第四行表示在45°＜θ＜135°采样条件下MLEM法和Kalman法的重建图像。由于所用模型比较简单，两种方法都可以重建出图像，但是可以看出，本发明在真实实验条件下的重建效果更加平滑和均匀，而用MLEM法重建出来的小球内部浓度分布则存在很明显的粗糙和不规则现象。特别是当采样角度进一步减小之后，这种差异更加明显。

Claims (4)

1、一种有限采样角度下基于卡尔曼滤波的PET浓度重建方法，其特征在于包括以下步骤：

(1)在有限角度下采集光子信号作为输出，得到投影线的正弦图数据Y；

(2)根据正弦图数据Y以及药物动力学理论建立一个由离散的观测方程(1)和连续的发展方程(2)构成的状态空间方程：

Y(t_k)＝DX(t_k)+v(t_k)......(1)

X′(t)＝AX(t)+Bμ(t)+ω(t)......(2)

其中：

t表示时间，t_k表示离散采样时间点；

D为系统矩阵；

Y是步骤(1)得到的正弦图数据；

X为放射性浓度的空间分布，为需要重建的对象；X`为X的导数；

A是状态转移矩阵，B是输入矩阵；

μ是血输入函数；

v为测量噪声；

ω是过程噪声；

(3)基于步骤(2)得到的状态空间方程，利用卡尔曼滤波算法计算得到放射性浓度X值，利用卡尔曼滤波算法计算时，采用方程(3)～(7)迭代重建实现。

${\hat{x}}^{\′}\left(t\right)=\mathrm{Ax}\left(t\right)+\mathrm{B\μ}\left(t\right)---\left(3\right)$

P′(t)＝AP(t)+P(t)A+Q(t)             (4)

K(t_k)＝P(t_k ^-)D^T[DP(t_k ^-)D+R(t_k)]^-1    (5)

$\hat{x}\left(t_k\right)=\hat{x}\left({t_k}^{-}\right)+K\left(t_k\right)[y\left(t_k\right)-D\hat{x}\left({t_k}^{-}\right)]---\left(6\right)$

P(t_k)＝[I-K(t_k)D]P(t_k ^-)              (7)

其中：

D^T是系统矩阵D的转置矩阵；

P是放射性浓度的空间分布X的误差协方差，P`是P的导数；

P(t_k ^-)分别是X、P的先验估计值，其中t_k ^-表示先验时刻；

和P(t_k)分别是X、P的后验估计值；

Q(t)以及R(t_k)分别表示该时刻的过程噪声以及测量噪声的方差。

2、根据权利要求1所述的PET浓度重建方法，其特征在于：所述的步骤(1)中所述的有限角度为采集角度θ，且30°＜θ＜150°。

3、根据权利要求1所述的PET浓度重建方法，其特征在于：所述的步骤(2)中的发展方程是根据药物动力学理论选择相应的房室模型得到。

4、根据权利要求1所述的PET浓度重建方法，其特征在于：所述的步骤(3)采用方程(3)～(7)迭代重建的步骤为：

1)首先设定放射性浓度的初始值和初始协方差X₀、P₀，迭代从下标k＝1开始进行；

2)利用方程(3)和方程(4)推算出放射性浓度的估计值和误差协方差的估计值P(t_k ^-)；

3)根据方程(5)计算出Kalman增益矩阵K；

4)利用离散采集时间点测得的正弦图数据y(t_k)，根据状态更新方程(6)、(7)来更新先验值

和P(t_k ^-)，得到这一时刻的后验估计值

和P(t_k)；

5)一次迭代完成以后，k值递增加1，重复交替步骤2)-4)，直到X值收敛，获得放射性浓度的最优估计值。

Images