CN102509322A - 一种基于卡尔曼滤波的pet图像重建方法 - Google Patents

Abstract

本发明涉及一种基于卡尔曼滤波的PET图像重建方法，包括采集原始投影线的正弦图数据，并基于状态空间的卡尔曼滤波方法重建PET图像。其中，利用卡尔曼滤波对PET图像进行重建的过程中，卡尔曼滤波的增益矩阵在求逆过程中采用LBFGS算法。本发明有效利用有LBFGS算法，改善了基于卡尔曼滤波的PET图像重建在处理大规模矩阵中存在的对高维矩阵存储容量及内存不足的问题；与现有重建方法的实验比较表明，在重建质量并无明显下降的前提下，计算效率有明显地提升。

Description

一种基于卡尔曼滤波的PET图像重建方法

技术领域

本发明涉及一种医学影像重建技术，特别涉及一种基于卡尔曼滤波的PET图像重建方法。

背景技术

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

近几年，PET探测器数量的增加、数据采集精度的提高以及临床上对高分辨率、大面积范围扫描图像的需求对图像重建算法提出了新的挑战：图像维数的扩展，将影响计算机处理速度、甚至造成计算机内存不足，另外，必须保证充分利用采集数据量方面的优势以提高图像分辨率。

目前，PET图像重建方法大致可分为两类：解析法和迭代统计法。前一类主要是滤波反投影法，计算速度快，但伪影严重。为此，以最大似然法为代表的迭代统计法被提出来，由于迭代法基于统计学模型，对不完全数据适应性好，逐渐成为PET重建算法研究关注点。为解决重建速度问题，大规模矩阵快速运算技术被不断运用到该算法体系中，如有序子集最大似然法。而针对于算法中未考虑先验估计的统计特性，提出了后验估计中加以修正的最大后验估计法。此外，将统计模型近似为高斯模型后又引入了最小二乘法，继而发展未带有惩罚加权的最小二乘、非负超松弛迭代法等。

状态空间体系从一个全新的角度提供了PET图像重建的新思路，通过根据实际求解问题调整测量方程与状态方程的统一表达式，以实现静态、动态重建以及先验估计。通过相关方法求解，如卡尔曼滤波、H∞滤波、粒子滤波等可适应不同噪声特性、不同运算条件，与传统的解析法或迭代统计法相比有一定的优势。但其构建的状态方程和测量方程中涉及高维度矩阵存储及运算，并且其中的求逆过程容易造成内存不足，是制约状态空间算法可行性的一大关键因素。

中国专利公告号为CN101499173B公开了一种PET成像中卡尔曼滤波图像重建方法，通过PET正电子发射断层扫描仪得到原始投影线的正弦图，然后建立状态空间体系，通过基于状态空间的卡尔曼滤波法得出放射性活度分布，重建图像。该技术在使用卡尔曼滤波对PET图像进行重建的过程中，没有解决处理大规模矩阵存在的存储量大及内存不足的问题。

发明内容

本发明提供了一种基于卡尔曼滤波的PET图像重建方法，解决了计算机在进行图像重建的过程中涉及高维度矩阵的求逆产生的内存不足的问题。

一种基于卡尔曼滤波的PET图像重建方法，包括以下步骤：

(1)采集原始投影线的正弦图数据并对其进行算法校正，得到当前的放射性浓度的测量浓度值y；

(2)基于状态空间的卡尔曼滤波方法向前预估放射性浓度的空间分布值和该空间分布值的误差协方差值；

(3)基于LBFGS算法模拟出卡尔曼滤波的增益矩阵；根据当前的测量浓度值和模拟得到的增益矩阵，对步骤(2)中的空间分布值和误差协方差值进行修正；将修正后的空间分布值和误差协方差值反馈至步骤(2)；

(4)重复交替步骤(2)和步骤(3)，不断更新放射性浓度的空间分布值，获得重建的PET图像。

在上述步骤中：

进一步包括步骤：

设计一步预测模块以用于向前预估步骤(2)中所述的空间分布值和空间分布值的误差协方差值；

所述的一步预测模块基于如下方程组：

$P_k^{-}=A{P}_{k-1}{A}^T+Q$

其中，x为PET图像中放射性浓度的空间分布值，

为向前预估的空间分布值，

为上一状态中反馈的修正的空间分布值；P为x的误差协方差，

为向前预估的误差协方差值，P_k-1为上一状态中反馈的修正的误差协方差值；k为迭代次数；A为状态转移矩阵；Q为过程噪声的方差，服从正态高斯分布。

所述的一步预测模块接收步骤(3)中反馈的修正后的空间分布值和误差协方差值，输出向前预估的空间分布值和误差协方差值。

进一步包括步骤：

设计状态更新模块以用于修正步骤(2)中的空间分布值和误差协方差值；所述的状态更新模块基于如下方程组：

$P_k=P_k^{-}-P_k^{-}{D}^T\mathrm{HD}{P}_k^{-}$

其中，H为增益矩阵的近似矩阵，即增益矩阵的模拟值；x为PET图像中放射性浓度的空间分布值，

为当前状态修正的空间分布值，

为当前状态预估的空间分布值；P为x的误差协方差，P_k为当前状态修正的误差协方差值，为当前状态预估的误差协方差值；k为迭代次数；v为测量浓度值y与空间分布值x比较的相关量；D为系统矩阵。

所述的状态更新模块接收当前状态预估的空间分布值和当前状态预估的误差协方差值，根据模拟得到的增益矩阵和当前的测量浓度值，输出修正后的空间分布值和误差协方差值。

进一步包括步骤：

设计增益矩阵模拟模块以用于模拟步骤(3)中所述的卡尔曼滤波的增益矩阵；所述的增益矩阵模拟模块实现以下步骤：

a、初始化增益矩阵的近似矩阵为一正定矩阵，构造增益矩阵的目标函数并设定使得该目标函数最小化的收敛条件，任取该目标函数上的一点作为初始点；

b、计算该目标函数的下一个搜索方向，在该方向上确定使得目标函数极小化的步长；基于该步长，获取目标函数上的下一个点，并检测该点是否满足目标函数的收敛条件；

c、若满足收敛条件，则此时的近似矩阵即为增益矩阵的模拟值；若未满足收敛条件，则用LBFGS算法中近似矩阵的更新方程，对近似矩阵进行修正，重复步骤b和步骤c，直至满足收敛条件。

在本步骤中，所述的步长可由沃尔夫准则确定；

所述的增益矩阵模拟模块还可以包括步骤：

i、预先设定增益矩阵模拟模块中的存储量；

ii、记录步骤c中对近似矩阵进行修正的次数值，并将该次数值与该存储量进行比较；

iii、若该存储量大于或等于次数值，则增益矩阵模拟模块对其存储数据不作处理；若该存储量小于次数值，则增益矩阵模拟模块对其存储数据进行部分删除，只记录最近修正的、组数等于存储量值的存储数据。

进一步地，所述的增益矩阵模拟模块基于如下方程组：

$f\left(t\right)=\frac{1}{2}{t}^T\mathrm{Ut}-v^{T}t$ ①

d_t＝-H_lg_l                                   ②

t_l+1＝t_l+α_ld_l                              ③

$H_{l+1}=(V_{l-1}^T\·\·\·V_{l-\stackrel{\‾}{j}}^T)H_0\left(V_{l-\stackrel{\‾}{j}}{\·\·\·V}_{l-1}\right)$

${+\mathrm{\ρ}}_{l-\stackrel{\‾}{j}}(V_{l-1}^T\·\·\·V_{l-\stackrel{\‾}{j}+1}^T){\mathrm{\σ}}_{l-\stackrel{\‾}{j}}{\mathrm{\σ}}_{l-\stackrel{\‾}{j}}^T\left(V_{l-\stackrel{\‾}{j}+1}{\·\·\·V}_{l-1}\right)$

${+\mathrm{\ρ}}_{l-\stackrel{\‾}{j}+1}(V_{l-1}^T\·\·\·V_{l-\stackrel{\‾}{j}+2}^T){\mathrm{\σ}}_{l-\stackrel{\‾}{j}+1}{\mathrm{\σ}}_{l-\stackrel{\‾}{j}+1}^T\left(V_{l-\stackrel{\‾}{j}+2}{\·\·\·V}_{l-1}\right)$ ④

$+\begin{array}{c}\·\\ \·\\ \·\end{array}$

$+{\mathrm{\ρ}}_{l-1}{\mathrm{\σ}}_{l-1}{\mathrm{\σ}}_{l-1}^T$

其中，U为增益矩阵的逆矩阵，H为增益矩阵的近似矩阵，即增益矩阵的模拟值；v为测量浓度值y与空间分布值x比较的相关量；方程①为构造的U的目标函数，即为增益矩阵的目标函数；g_l为对方程①进行一阶求导所得的值，即

方程②为搜索方向的计算式；d_l为搜索方向；t为中间变量；l为迭代次数；方程④为LBFGS算法中近似矩阵H的更新方程，用以获取近似矩阵H；H_l+1为第l次迭代过程中的近似矩阵，初始值H₀为一正定矩阵；ρ_l、σ_l、δ_l和V_l均为第l次迭代过程中的中间变量，满足δ_l＝t_l+1-t_l，σ_l＝g_l+1-g_l，

α_l为第l次迭代过程中使f(t)极小化的步长；j为增益矩阵模拟模块中预先设定的存储量，

为增益矩阵模拟模块中实际记录的存储数据的组数，满足

表示增益矩阵模拟模块中记录了第l次到第

次迭代过程中、共组存储数据；所述的存储数据即为ρ_l、σ_l、δ_l以及V_l等。

本发明的有益效果为：

有效利用有LBFGS算法，改善了基于卡尔曼滤波的PET图像重建在处理大规模矩阵中存在的对高维矩阵存储容量及内存不足的问题；与现有重建方法的实验比较表明，在重建质量并无明显下降的前提下，计算效率有明显地提升。

附图说明

图1为本发明实施例中的PET图像重建方法的流程示意图；

图2为本发明实施例中设计模块的结构层次示意图；

图3为本发明实施例中设计的增益矩阵模拟模块模拟增益矩阵的流程示意图；

图4为本发明实施例中设计的增益矩阵模拟模块优化内存占用量的流程示意图；

图5为本发明实验使用的数字体模型示意图；

图6为传统的卡尔曼滤波方法的重建结果和利用本发明方法的重建结果与真值图像对比的示意图；

图7为传统的尔曼滤波方法和利用本发明方法的数据结果分析图。

具体实施方式

下面结合附图详细介绍本发明的具体实施方式。

如图1所示的一种基于卡尔曼滤波的PET图像重建方法，包括步骤：

S1采集原始投影线的正弦图数据并对其进行算法校正，得到当前的放射性浓度的测量浓度值y；

S2基于状态空间的卡尔曼滤波方法向前预估放射性浓度的空间分布值和该空间分布值的误差协方差值；

S3基于LBFGS算法模拟出卡尔曼滤波的增益矩阵；根据当前的测量浓度值和模拟得到的增益矩阵，对上述步骤中的空间分布值和误差协方差值进行修正；将修正后的空间分布值和误差协方差值反馈至S2；重复交替S2和S3，不断更新放射性浓度的空间分布值。

S4获得重建的PET图像。

结合图2：

为了实现S2中的向前预估所述的空间分布值和空间分布值的误差协方差值，设计一步预测模块D1。

一步预测模块D1基于如下方程组：

$P_k^{-}=A{P}_{k-1}{A}^T+Q$

其中，x为PET图像中放射性浓度的空间分布值，

为向前预估的空间分布值，为上一状态中反馈的修正的空间分布值；P为x的误差协方差，

为向前预估的误差协方差值，P_k-1为上一状态中反馈的修正的误差协方差值；k为迭代次数；A为状态转移矩阵；Q为过程噪声的方差，服从正态高斯分布。

一步预测模块接收上一状态中反馈的修正的空间分布值

和误差协方差值P_k-1，输出向前预估的空间分布值

和误差协方差值

为了实现S3中的修正经S2产生的空间分布值和误差协方差值，设计状态更新模块D2。

状态更新模块D2基于如下方程组：

$P_k=P_k^{-}-P_k^{-}{D}^T\mathrm{HD}{P}_k^{-}$

其中，H为增益矩阵的近似矩阵，即增益矩阵的模拟值；x为PET图像中放射性浓度的空间分布值，为当前状态修正的空间分布值，为当前状态预估的空间分布值；P为x的误差协方差，P_k为当前状态修正的误差协方差值，为当前状态预估的误差协方差值；k为迭代次数；v为测量浓度值y与空间分布值x比较的相关量；D为系统矩阵；

状态更新模块D2接收当前状态预估的空间分布值

和当前状态预估的误差协方差值

并根据增益矩阵的近似矩阵H和当前的测量浓度值y，输出修正后的空间分布值

和误差协方差值P_k。

为了模拟S3中的卡尔曼滤波的增益矩阵，设计基于LBFGS算法的增益矩阵模拟模块D3。

结合图3，增益矩阵模拟模块D3实现以下步骤：

S301初始化增益矩阵的近似矩阵H为一正定矩阵，构造增益矩阵的目标函数并设定使得该目标函数最小化的收敛条件，任取该目标函数上的一点作为初始点；

S302计算该目标函数的下一个搜索方向，在该方向上确定使得目标函数极小化的步长α；基于该步长α，获取目标函数上的下一个点，并检测该点是否满足目标函数的收敛条件；

S303若满足收敛条件，则此时的近似矩阵H即为增益矩阵的模拟值，输出H；若未满足收敛条件，则用LBFGS算法中近似矩阵H的更新方程，对近似矩阵H进行修正，重复S302和步骤S303，直至满足收敛条件。

S302中的步长α可由沃尔夫准则确定；结合图4，增益矩阵模拟模块D3中还包括步骤：

S311预先设定增益矩阵模拟模块D3中的存储量；

S312记录S303中对近似矩阵H进行修正的次数值，并将该次数值与该存储量进行比较；

S313若该存储量大于或等于次数值，则增益矩阵模拟模块D3对其存储数据不作处理；若该存储量小于次数值，则增益矩阵模拟模块D3对其存储数据进行部分删除，只记录最近修正的、组数等于存储量值的存储数据。

进一步地，增益矩阵模拟模块D3基于如下方程组：

$f\left(t\right)=\frac{1}{2}{t}^T\mathrm{Ut}-v^{T}t$ ①

d_l＝-H_lg_l                                ②

t_l+1＝t_l+α_ld_l                           ③

$H_{l+1}=(V_{l-1}^T\·\·\·V_{l-\stackrel{\‾}{j}}^T)H_0\left(V_{l-\stackrel{\‾}{j}}{\·\·\·V}_{l-1}\right)$

${+\mathrm{\ρ}}_{l-\stackrel{\‾}{j}}(V_{l-1}^T\·\·\·V_{l-\stackrel{\‾}{j}+1}^T){\mathrm{\σ}}_{l-\stackrel{\‾}{j}}{\mathrm{\σ}}_{l-\stackrel{\‾}{j}}^T\left(V_{l-\stackrel{\‾}{j}+1}{\·\·\·V}_{l-1}\right)$

${+\mathrm{\ρ}}_{l-\stackrel{\‾}{j}+1}(V_{l-1}^T\·\·\·V_{l-\stackrel{\‾}{j}+2}^T){\mathrm{\σ}}_{l-\stackrel{\‾}{j}+1}{\mathrm{\σ}}_{l-\stackrel{\‾}{j}+1}^T\left(V_{l-\stackrel{\‾}{j}+2}{\·\·\·V}_{l-1}\right)$ ④

$+\begin{array}{c}\·\\ \·\\ \·\end{array}$

$+{\mathrm{\ρ}}_{l-1}{\mathrm{\σ}}_{l-1}{\mathrm{\σ}}_{l-1}^T$

其中，U为增益矩阵的逆矩阵，H为增益矩阵的近似矩阵，即增益矩阵的模拟值；v为测量浓度值y与空间分布值x比较的相关量；方程①为构造的U的目标函数，即为增益矩阵的目标函数；g_l为对方程①进行一阶求导所得的值，即

方程②为搜索方向的计算式；d_l为搜索方向；t为中间变量；l为迭代次数；方程④为LBFGS算法中近似矩阵H的更新方程，用以获取近似矩阵H；H_l+1为第l次迭代过程中的近似矩阵，初始值H₀为一正定矩阵；ρ_l、σ_l、δ_l和V_l均为第l次迭代过程中的中间变量，满足δ_l＝t_l+1-t_l，σ_l＝g_l+1-g_l， α_l为第l次迭代过程中使f(t)极小化的步长；j为增益矩阵模拟模块中预先设定的存储量，

为增益矩阵模拟模块中实际记录的存储数据的组数，满足

表示增益矩阵模拟模块中记录了第l次到第次迭代过程中、共组存储数据；所述的存储数据即为ρ_l、σ_l、δ_l以及V_l等。

基于上述方法的具体实施过程如下：

正电子发射断层扫描仪探测人体内发出的放射性信号，经过符合和采集系统处理，形成原始投影线，并以正弦图的方式存放于计算机硬盘中；对原始采集到的正弦图进行各类校正，得到当前的放射性浓度的测量浓度值y。计算机以正弦图为输入，调用相关模块，计算得出PET重建图像。

本实施例包括建立相关模块：

一步预测模块D1，用于接收上一状态中反馈的修正的空间分布值和误差协方差值P_k-1，并输出向前预估的空间分布值

和误差协方差值

状态更新模块D2，用于接收当前状态预估的空间分布值

和当前状态预估的误差协方差值

并根据增益矩阵的近似矩阵H和当前的测量浓度值y，输出修正后的空间分布值

和误差协方差值P_k；

增益矩阵模拟模块D3，用于模拟得到增益矩阵的近似值，即增益矩阵的近似矩阵H，并在此过程中实现存储数据内存占用量的优化。

传统的卡尔曼滤波方法往往基于如下迭代方程以重建PET图像：

$P_k^{-}=A{P}_{k-1}{A}^T+Q$

K＝P_kD^T[DP_kD^T+R]^-1                  ⑤

$P_k=P_k^{-}-\mathrm{KD}{P}_k^{-}$

其中：k表示迭代次数，迭代从方程的初始值x₀、P₀出发，通过采集数据不断更新放射性浓度空间分布x，最终收敛获得重建图像。在求卡尔曼滤波的增益矩阵时涉及大规模矩阵求逆(可见式⑤中所涉及到的大规模矩阵求逆，K为传统的卡尔曼滤波方法中的增益矩阵)，往往受到电脑内存和存储容量的限制。

本发明利用LBFGS算法解决式⑤中所涉及到的大规模矩阵求逆，将卡尔曼滤波法的增益矩阵作为参数，构造增益矩阵的目标函数以求得其近似矩阵H。U为增益矩阵的逆矩阵，以此构造的线性方程见式⑥，并将式⑥转化为式①的目标函数，即将模拟增益矩阵数值的问题转化为求式①这个严格凸二次函数的最小化问题，具体方程组如下：

Ut＝v                    ⑥

$f\left(t\right)=\frac{1}{2}{t}^T\mathrm{Ut}-v^{T}t$ ①

其中：U为增益矩阵的逆矩阵，满足U＝DP_kD^T+R；t为所构造的目标函数的中间变量；v为测量浓度值y与空间分布值x比较的相关量，满足 $v=y-D{x}_k^{-}.$

对式①进行求解，建立迭代方程组：

$g_l=\▿f\left(t_l\right)=U{t}_l-v$ ⑦

d_l＝-H_lg_l                        ⑧

t_l+1＝t_l+α_ld_l                   ⑨

$H_{l+1}=(V_{l-1}^T\·\·\·V_{l-\stackrel{\‾}{j}}^T)H_0\left(V_{l-\stackrel{\‾}{j}}{\·\·\·V}_{l-1}\right)$

${+\mathrm{\ρ}}_{l-\stackrel{\‾}{j}}(V_{l-1}^T\·\·\·V_{l-\stackrel{\‾}{j}+1}^T){\mathrm{\σ}}_{l-\stackrel{\‾}{j}}{\mathrm{\σ}}_{l-\stackrel{\‾}{j}}^T\left(V_{l-\stackrel{\‾}{j}+1}{\·\·\·V}_{l-1}\right)$

${+\mathrm{\ρ}}_{l-\stackrel{\‾}{j}+1}(V_{l-1}^T\·\·\·V_{l-\stackrel{\‾}{j}+2}^T){\mathrm{\σ}}_{l-\stackrel{\‾}{j}+1}{\mathrm{\σ}}_{l-\stackrel{\‾}{j}+1}^T\left(V_{l-\stackrel{\‾}{j}+2}{\·\·\·V}_{l-1}\right)$ ⑩

$+\begin{array}{c}\·\\ \·\\ \·\end{array}$

$+{\mathrm{\ρ}}_{l-1}{\mathrm{\σ}}_{l-1}{\mathrm{\σ}}_{l-1}^T$

其中，g_l为对方程①进行一阶求导所得的值；式⑧为搜索方向的计算式；d_l为搜索方向；式⑨为搜索目标函数下一个点的计算式，t为中间变量；l为迭代次数；方程⑩为LBFGS算法中近似矩阵H的更新方程；用以获取近似矩阵H；H_l+1为第l次迭代过程中的近似矩阵，初始值H₀为一正定矩阵；ρ_l、σ_l、δ_l和V_l均为第l次迭代过程中的中间变量，满足δ_l＝t_l+1-t_l，σ_l＝g_l+1-g_l， α_l为第l次迭代过程中使目标函数f(t)极小化的步长；j为增益矩阵模拟模块中预先设定的存储量，

为增益矩阵模拟模块中实际记录的存储数据的组数，满足

表示增益矩阵模拟模块中记录了第l次到第

次迭代过程中、共

组存储数据；所述的存储数据即为ρ_l、σ_l、δ_l以及V_l等。

为模拟得到增益矩阵U，设计的增益矩阵模拟模块D3基于式⑦、式⑧、式⑨、以及式⑩，实现以下步骤：

(1)设定使目标函数①最小化的收敛条件为||g_l||≤ε；初始化增益矩阵的近似矩阵H为一正定矩阵H₀，且H₀为一单位矩阵；t₀为任意值(即任取目标函数①上的一点作为初始点)；迭代次数l＝0；设定收敛范围ε，满足ε＝10e-6，初始化选取α₀＝1作为试验步长。

(2)通过式⑧计算目标函数①的下一个搜索方向d_l，在该方向上确定使得目标函数①极小化的步长α_l；基于该步长α_l，通过式⑨获取目标函数①上的下一个点，并通过式⑦检测该点是否满足目标函数的收敛条件||g_l||≤ε；其中，步长α_l由沃尔夫准则确定。

(3)若满足收敛条件，则此时的近似矩阵H即为增益矩阵的模拟值，输出该近似矩阵H；若未满足收敛条件，则通过式⑩对近似矩阵H进行修正，重复步骤(2)和步骤(3)，直至满足收敛条件。

在式⑩中，j为预先设定增益矩阵模拟模块中的存储量，

为增益矩阵模拟模块中实际记录的存储数据的组数。实际记录的存储数据的组数

和预先设定的存储量j有如下关系的关系满足

具体如下：

增益矩阵模拟模块D3记录对近似矩阵H进行修正的次数值，若此时为第l次迭代(即第l次对近似矩阵H进行修正)，那么次数值即为l；将次数值l与存储量j进行比较，若存储量j大于或等于次数值l，则增益矩阵模拟模块D3对其存储数据(即第l次迭代之前存储的存储数据，如ρ_l、σ_l、δ_l以及V_l等)不作处理；若存储量j小于次数值l，则增益矩阵模拟模块D3对其存储数据进行部分删除，只记录最近修正的、且组数等于存储量值j的存储数据。

根据上述的解决式⑤中所涉及到的大规模矩阵求逆的LBFGS算法，改写传统的卡尔曼滤波方法基于的迭代方程为：

$P_k^{-}=A{P}_{k-1}{A}^T+Q$

$P_k=P_k^{-}-P_k^{-}{D}^T\mathrm{HD}{P}_k^{-}$

其中，

设计的一步预测模块D1基于方程组：

$P_k^{-}=A{P}_{k-1}{A}^T+Q$

其中，x为PET图像中放射性浓度的空间分布值，

为向前预估的空间分布值，

为上一状态中反馈的修正的空间分布值；P为x的误差协方差，

为向前预估的误差协方差值，P_k-1为上一状态中反馈的修正的误差协方差值；k为迭代次数；A为状态转移矩阵；Q为过程噪声的方差，服从正态高斯分布。

设计的状态更新模块D2基于如下方程组：

$P_k=P_k^{-}-P_k^{-}{D}^T\mathrm{HD}{P}_k^{-}$

其中，H为增益矩阵的近似矩阵，即增益矩阵的模拟值；x为PET图像中放射性浓度的空间分布值，

为当前状态修正的空间分布值，为当前状态预估的空间分布值；P为x的误差协方差，P_k为当前状态修正的误差协方差值，

为当前状态预估的误差协方差值；k为迭代次数；v为测量浓度值y与空间分布值x比较的相关量；D为系统矩阵。

采用基于LBFGS算法的卡尔曼滤波方法来重建PET图像时，主要包括以下几个步骤：

1)在一步预测模块D1中设定待估计空间分布值的初始值

(理论上可为任意值)，误差协方差P₀，Q，R，k＝0；

2)一步预测模块D1向前预估放射性浓度的空间分布值

和该空间分布值的误差协方差值

3)利用增益矩阵模拟模块D3模拟出当前的卡尔曼滤波的增益矩阵；状态更新模块D2根据当前的测量浓度值y和增益矩阵的近似矩阵H，对当前状态预估的空间分布值

和当前状态预估的误差协方差值

进行修正；并将修正后的空间分布值

和误差协方差值P_k反馈至一步预测模块D1。

4)不断重复步骤2)和步骤3)，直至获得最优重建结果。

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

本发明技术采用Zubal Phantom数字体模模型实验验证算法的有效性。如图五所示，该模型包含3个组织区域和一个背景区域。本实验运行环境为：3G内存，2.93GHz，32位操作系统，CPU为intel酷睿双核。

针对重建速度的验证，采用64个投影角度，每个角度下射束为64条的原始采集数据。本实验将本发明基于LBFGS算法的卡尔曼滤波方法和传统的卡尔曼滤波方法重建结果做比较，将二者的初始值保持一致以保证时间的可比性，具体参数设置如下：A＝I_n，I_n为n×n维单位矩阵，P₀＝αI_n，其中α＝1×10^-4，；Q₀＝0，Q₀为n×n维单位矩阵，R₀＝βI_m，I_m为m×m的单位矩阵，β＝0.1，空间分布值的初始值x₀＝0，n×1维向量，测量值y为m×1维向量，m＝4096；迭代次数k＝2；运用LBFGS算法求逆过程中，迭代次数l＝5，可控存储量j＝3已经能够保证完全收敛，采用沃尔夫条件进行线性搜索。

表1中列出了针对同一测量值，用以上两种方法分别重建出64×64维度和32×32维度图像过程中各阶段所耗费的时间。其中，T0表示对增益矩阵U的求逆过程所需时间，由于算法的差异涵盖时间的值略有不同，对本发明为计算

所需时间，而对传统卡尔曼滤波算法则计算

的值，因此前者只涉及对逆矩阵与m×1维矩阵相乘的结果值估计过程，而与高维矩阵D^T、P相乘过程的时间被包含在T3时刻；后者所包含的时间除了求逆过程外还包含与高维矩阵D^T、P相乘的时间；T1表示计算先验估计值

所需时间；T2表示计算后验估计

所需时间，包括包含高维逆矩阵运算的增益矩阵求解；T3表示计算后验估计P_k所需时间、T_ALL表示重建迭代过程总时间，各重建5次，取平均值。值得注意到是，T0秒和T3的定义对两种算法略有不同，约3秒的差异，若忽略算法运算过程的不一致性对两种算法评价结果无较大影响，从表1中可以看出，对于维数m×m恒定的矩阵求逆过程中，利用LBFGS算法求逆所需时间相对于传统的求逆方法，分别为原来的1/9.6，1/6.77；重建所需总时间降低为原来的1/1.68、1/3.1。

表1本发明技术与传统卡尔曼滤波法重建速度定量分析

针对重建图像质量的验证，采用高维度的原始采集数据128个投影角度，每个角度下射束为128条，即m＝16384，重建图像大小为128×128，即维度n＝16384；初始值设定同上。图六是真值图像、基于本发明技术方案重建的图像和传统卡尔曼滤波方法重建的图像的比对示意图，可以看直观地看出两者结果的一致性。为了进一步研究重建图像的精确度，对本发明技术方案重建的图像与传统卡尔曼滤波算法重建的图像进行了数据的分析。分别取48、64、78、96列队数据进行分析，结果如图七所示。

表2本发明技术重建结果与KF、EM、FBP重建结果数值分析

	本发明技术	KF	EM	FBP
					方差1	0.18693	0.185954	0.295927	0.317716
标准差1	0.25066	0.251748	0.383821	0.379815
					方差2	0.182603	0.180702	0.292611	0.423383
标准差2	0.312053	0.3175	0.503861	0.599736

同时，将本发明技术方案与传统的几种方法重建结果的偏差进行比较，如表2所示，KF滤波方法的重建结果优于典型的投影法——滤波反投影法(FBP)和迭代法——最大似然期望最大法(ML-EM)，而运用本发明技术所重建结果，其方差和标准差变化很小。从表2可以看出，本发明技术方案与传统的直接对增益矩阵求逆的卡尔曼滤波算法重建数据基本能够保持一致性，说明对本技术方案对增益矩阵求逆的准确性较高。

Claims (7)

1.一种基于卡尔曼滤波的PET图像重建方法，包括以下步骤：

(1)采集原始投影线的正弦图数据并对其进行算法校正，得到当前的放射性浓度的测量浓度值y；

(2)基于状态空间的卡尔曼滤波方法向前预估放射性浓度的空间分布值和该空间分布值的误差协方差值；

(3)基于LBFGS算法模拟出卡尔曼滤波的增益矩阵；根据当前的测量浓度值和模拟得到的增益矩阵，对步骤(2)中的空间分布值和误差协方差值进行修正；将修正后的空间分布值和误差协方差值反馈至步骤(2)；

(4)重复交替步骤(2)和步骤(3)，不断更新放射性浓度的空间分布值，获得重建的PET图像。

2.根据权利要求1所述的PET图像重建方法，其特征在于：包括步骤：设计一步预测模块以用于向前预估步骤(2)中所述的空间分布值和空间分布值的误差协方差值；

所述的一步预测模块基于如下方程组：

$P_k^{-}=A{P}_{k-1}{A}^T+Q$

其中，x为PET图像中放射性浓度的空间分布值，

为向前预估的空间分布值，

为上一状态中反馈的修正的空间分布值；P为x的误差协方差，

为向前预估的误差协方差值，P_k-1为上一状态中反馈的修正的误差协方差值；k为迭代次数；A为状态转移矩阵；Q为过程噪声的方差，服从正态高斯分布。

3.根据权利要求1所述的PET图像重建方法，其特征在于：包括步骤：设计状态更新模块以用于修正步骤(2)中的空间分布值和误差协方差值；所述的状态更新模块基于如下方程组：

$P_k=P_k^{-}-P_k^{-}{D}^T\mathrm{HD}{P}_k^{-}$

其中，H为增益矩阵的近似矩阵，即增益矩阵的模拟值；x为PET图像中放射性浓度的空间分布值，

为当前状态修正的空间分布值，

为当前状态预估的空间分布值；P为x的误差协方差，P_k为当前状态修正的误差协方差值，

为当前状态预估的误差协方差值；k为迭代次数；v为测量浓度值y与空间分布值x比较的相关量；D为系统矩阵。

4.根据权利要求1所述的PET图像重建方法，其特征在于：包括步骤：设计增益矩阵模拟模块以用于模拟步骤(3)中所述的卡尔曼滤波的增益矩阵；所述的增益矩阵模拟模块实现以下步骤：

a、初始化增益矩阵的近似矩阵为一正定矩阵，构造增益矩阵的目标函数并设定使得该目标函数最小化的收敛条件，任取该目标函数上的一点作为初始点；

b、计算该目标函数的下一个搜索方向，在该方向上确定使得目标函数极小化的步长；基于该步长，获取目标函数上的下一个点，并检测该点是否满足目标函数的收敛条件；

c、若满足收敛条件，则此时的近似矩阵即为增益矩阵的模拟值；若未满足收敛条件，则用LBFGS算法中近似矩阵的更新方程，对近似矩阵进行修正，重复步骤b和步骤c，直至满足收敛条件。

5.根据权利要求4所述的PET图像重建方法，其特征在于：所述的步长由沃尔夫准则确定。

6.根据权利要求4或5所述的PET图像重建方法，其特征在于：包括以下步骤：

i、预先设定增益矩阵模拟模块中的存储量；

ii、记录步骤c中对近似矩阵进行修正的次数值，并将该次数值与该存储量进行比较；

iii、若该存储量大于或等于次数值，则增益矩阵模拟模块对其存储数据不作处理；若该存储量小于次数值，则增益矩阵模拟模块对其存储数据进行部分删除，只记录最近修正的、组数等于存储量值的存储数据。

7.根据权利要求6所述的PET图像重建方法，其特征在于：所述的增益矩阵模拟模块基于如下方程组：

$f\left(t\right)=\frac{1}{2}{t}^T\mathrm{Ut}-v^{T}t$ ①

d_l＝-H_lg_l                            ②

t_l+1＝t_l+α_ld_l                       ③

$H_{l+1}=(V_{l-1}^T\·\·\·V_{l-\stackrel{\‾}{j}}^T)H_0\left(V_{l-\stackrel{\‾}{j}}{\·\·\·V}_{l-1}\right)$

${+\mathrm{\ρ}}_{l-\stackrel{\‾}{j}}(V_{l-1}^T\·\·\·V_{l-\stackrel{\‾}{j}+1}^T){\mathrm{\σ}}_{l-\stackrel{\‾}{j}}{\mathrm{\σ}}_{l-\stackrel{\‾}{j}}^T\left(V_{l-\stackrel{\‾}{j}+1}{\·\·\·V}_{l-1}\right)$

${+\mathrm{\ρ}}_{l-\stackrel{\‾}{j}+1}(V_{l-1}^T\·\·\·V_{l-\stackrel{\‾}{j}+2}^T){\mathrm{\σ}}_{l-\stackrel{\‾}{j}+1}{\mathrm{\σ}}_{l-\stackrel{\‾}{j}+1}^T\left(V_{l-\stackrel{\‾}{j}+2}{\·\·\·V}_{l-1}\right)$ ④

$+\begin{array}{c}\·\\ \·\\ \·\end{array}$

$+{\mathrm{\ρ}}_{l-1}{\mathrm{\σ}}_{l-1}{\mathrm{\σ}}_{l-1}^T$

其中，U为增益矩阵的逆矩阵，H为增益矩阵的近似矩阵，即增益矩阵的模拟值；v为测量浓度值y与空间分布值x比较的相关量；方程①为构造的U的目标函数，即为增益矩阵的目标函数；g_l为对方程①进行一阶求导所得的值，即

方程②为搜索方向的计算式；d_l为搜索方向；t为中间变量；l为迭代次数；方程④为LBFGS算法中近似矩阵H的更新方程，用以获取近似矩阵H；H_l+1为第l次迭代过程中的近似矩阵，初始值H₀为一正定矩阵；ρ_l、σ_l、δ_l和V_l均为第l次迭代过程中的中间变量，满足δ_l＝t_l+1-t_l，σ_l＝g_l+1-g_l，

α_l为第l次迭代过程中使f(t)极小化的步长；j为增益矩阵模拟模块中预先设定的存储量，

为增益矩阵模拟模块中实际记录的存储数据的组数，满足表示增益矩阵模拟模块中记录了第l次到第

次迭代过程中、共

组存储数据；所述的存储数据即为ρ_l、σ_l、δ_l以及V_l等。

Images