基于改进序贯滤波的静态PET图像重建方法

摘要

一种基于改进序贯滤波的静态PET图像重建方法，包括如下步骤：1)建立PET系统的状态空间模型；2)基于改进序贯滤波的浓度重建过程如下：2.1)首先设定放射性浓度分布的初始值和初始估计协方差P0(0)＝P(0)；2.2)获取正弦图数据y(t)；2.3)将正弦图数据y(t)分成r块；2.4)利用方程(3)计算增益矩阵Ki(t)；2.5)利用量测值yi(t)及增益Ki(t)，根据状态更新方程(2)计算出空间浓度估计值并根据(4)推出相应的预估误差协方差阵Pi(t)；2.6)若i＜r,i＝i+1，跳转至步骤2.4)，否则，P(t)＝Pr(t)；2.7)若t＜N,t＝t+1，跳转至步骤2.2)，否则，算法结束，得到最终的浓度重构结果。本发明在保证重构效果的同时，降低计算成本，提高重构速度。

说明书

技术领域

本发明涉及正电子发射断层成像技术领域，具体涉及一种基于改进序贯滤波的静态PET图像重建方法。

背景技术

正电子发射断层成像(PositronEmissionTomography，PET)技术是一种通过在分子水平上对代谢变化进行追踪、监控，从而达到对疾病的早期检测、预防目的的医学成像技术。PET技术作为功能性成像技术，能反映病人体内新陈代谢的情况，因而能更早地探测出病灶。因此，PET技术被广泛应用于肿瘤、心脏疾病、神经和精神系统疾病的诊断及实验生物成像、药物筛选与开发，已显示出巨大的应用前景。

在进行PET扫描时首先使用加速器产生正电子发射核素，然后将由放射性同位核素标记的药物注入人体内，通过血液循环，这些物质在人体内各组织器官中将形成一定的分布。由于放射性同位核素的半衰期较短，且极其不稳定，将很快发生衰变。在衰变过程中，放射性同位核素会产生正电子并与附近的自由电子发生湮灭反应，产生一对方向相反、能量相等的伽玛光子。这时，可以通过体外的探测器利用符合技术对光子进行探测，再经由噪声矫正，得到投影数据。通过一些重构方法利用投影数据进行反演求解，可重建出人体的放射性物质的空间浓度分布。

目前，PET图像重建方法大致可分为两类：解析法和迭代统计法。前一类主要是滤波反投影法，计算速度快，代价小，但不能很好抑制噪声，以致重建图像质量不高。因此，出现了以极大似然法为代表的迭代统计法。由于迭代法基于统计学模型，对不完全数据适应性好，已逐渐成为PET重建算法研究的焦点。相比于解析法，迭代法获得的重建图像更加清晰。但对PET系统成像过程的准确建模是影响迭代法重建结果的关键。状态空间的引入，有力地推动了迭代法的发展。状态空间法能够将PET成像的过程模型化(给出PET成像过程的数学表达)，对PET扫描过程的统计特性与生物体的生理、结构特性进行相应的刻画，因而该模型联合一些已有的估计算法往往能够提供更好的重建效果。目前已有的状态空间求解方法主要基于H₂滤波(即卡尔曼滤波)、H_∞滤波等。然而，由于需要重建图像的清晰度往往取决于体素的个数，而后者又直接关系到滤波重构算法的计算量。因此，大多数情况下，基于滤波算法的PET图像重构方法都面临非常高维度的矩阵求逆运算，这会给基于一般滤波算法的PET图像重构带来极大的运算负担。

发明内容

为了克服现有的PET图像重建方法的运算成本较大、重构速度较慢的不足，本发明提供了一种基于改进序贯滤波的静态PET图像重建方法，该发明提出将可观测的高维正弦数据分成多个低维观测数据。当系统基于这些低维观测进行滤波时，原本要求逆的高维矩阵被一些与低维观测维数一致的低维矩阵替代了，这会极大地降低计算成本，提高重构速度，提供与基于卡尔曼滤波的PET重建方法一样的重建效果。

为了解决上述技术问题提供了如下技术方案：

一种基于改进序贯滤波的静态PET图像重建方法，包括如下步骤：

1)建立PET系统的状态空间模型：

$>\left(\begin{array}{c}x\left(t+1\right)=x\left(t\right)+v\left(t\right)\\ y\left(t\right)=Dx\left(t\right)+e\left(t\right)\end{array}\right)---\left(1\right)$>

其中，t表示时间；y(t)为观测值，就是经过噪声矫正后得到的正弦图数据；D表示人体内放射性浓度与PET扫描之间的投影关系的投影矩阵，由PET装置的固有特性决定；x(t)为放射性浓度分布，即需要重建的对象；v(t)是过程噪声；e(t)为数据采集并经符合矫正后残留的量测噪声；

2)根据下列方程得到基于改进序贯滤波的重建图像：

$>{\hat{x}}_i\left(t\right)={\hat{x}}_{i-1}\left(t\right)+K_i\left(t\right)\left(y_i\right(t)-D{\hat{x}}_{i-1}(t\left)\right),{\hat{x}}_0\left(t\right)=\hat{x}\left(0\right)---\left(2\right)$>

$>K_i\left(t\right)=P_i\left(t\right)D_i^T{(D_i{P}_i\left(t\right)D_i^T+R_i)}^{-1}---\left(3\right)$>

$>P_i\left(t\right)=P_{i-1}\left(t\right)+Q_i-K_i\left(t\right)\left(D_i{P}_i\left(t\right)D_i^T+R_i\right)K_i^T\left(t\right),P_0\left(t\right)=P\left(t\right)---\left(4\right)$>

$>\hat{x}\left(t\right)={\hat{x}}_r\left(t\right)---\left(5\right)$>

P(t)＝P_r(t)(6)

其中，为空间浓度的滤波重建值，为空间浓度的滤波初始值，y(t)为可测得的校正后的正弦图数据，P(t)为空间浓度的预估误差协方差阵，P(0)为初始空间浓度预估误差协方差，将y(t)分为r块，y_i(t)是y(t)的第i分块，D_i,R_i,Q_i是相应于y_i(t)的矩阵D,R,Q的分块，K_i(t)为相应于y_i(t)的滤波增益矩阵，为基于{y(0),…,y(t-1),y₁(t),…,y_i(t)}的空间浓度滤波重建值，P_i(t)是相应于{y(0),…,y(t-1),y₁(t),…,y_i(t)}的滤波误差协方差阵。迭代从初始值P(0)出发，通过量测值y(t)，经过N次迭代，最终得到放射性浓度分布重建过程如下：

2.1)首先设定放射性浓度分布的初始值和初始估计协方差P₀(0)＝P(0)；

2.2)获取正弦图数据y(t)；

2.3)将正弦图数据y(t)分成r块；

2.4)利用方程(3)计算增益矩阵K_i(t)；

2.5)利用量测值y_i(t)及增益K_i(t)，根据状态更新方程(2)计算出空间浓度估计值并根据(4)推出相应的预估误差协方差阵P_i(t)；

2.6)若i＜r,i＝i+1，跳转至步骤2.4)，否则，P(t)＝P_r(t)；

2.7)若t＜N,t＝t+1，跳转至步骤2.2)，否则，算法结束。

本发明的技术构思为：基于对PET系统使用卡尔曼滤波方法成像原理的分析，发现该重构方法重构速度缓慢的问题在于高维观测导致高维矩阵求逆，通过将高维观测分块成相对低维的观测，最终只需求解低维矩阵的逆，从而有效地解决了上述问题。

从上述技术方案可以看出，本发明的有益效果主要表现在：在保证了(与直接基于卡尔曼滤波的PET图像重构)有相同的重构效果的前提下，有效避免了高维矩阵求逆，进而提高了重构速度。

附图说明

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

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的静态PET浓度重建方法进行详细说明。

如图1所示，一种基于改进序贯滤波的静态PET图像重建方法，包括如下步骤：

1)根据PET成像原理，建立状态空间体系：

$>\left(\begin{array}{c}x\left(t+1\right)=x\left(t\right)+v\left(t\right)\\ y\left(t\right)=Dx\left(t\right)+e\left(t\right)\end{array}\right)---\left(7\right)$>

其中，t表示时间；y(t)为观测值，就是经过噪声矫正后得到的正弦图数据；D表示人体内放射性浓度与PET扫描之间的投影关系的投影矩阵，由PET装置的固有特性决定，在实验中采用单响应线模型近似计算得到；x(t)为放射性浓度分布，即需要重建的对象；v(t)是过程噪声；e(t)为数据采集并经过噪声矫正后残留的量测噪声；v(t)，e(t)服从协方差矩阵分别为对角阵Q,R的正态高斯分布；

2)根据下列方程得到基于改进序贯滤波的重建图像：

$>{\hat{x}}_i\left(t\right)={\hat{x}}_{i-1}\left(t\right)+K_i\left(t\right)\left(y_i\right(t)-D{\hat{x}}_{i-1}\left(t\right)),{\hat{x}}_0\left(t\right)=\hat{x}\left(0\right)---\left(8\right)$>

$>K_i\left(t\right)=P_i\left(t\right)D_i^T{(D_i{P}_i\left(t\right)D_i^T+R_i)}^{-1}---\left(9\right)$>

$>P_i\left(t\right)=P_{i-1}\left(t\right)+Q_i-K_i\left(t\right)\left(D_i{P}_i\left(t\right)D_i^T+R_i\right)K_i^T\left(t\right),P_0\left(t\right)=P\left(t\right)---\left(10\right)$>

P(t)＝P_r(t)(12)

其中，为空间浓度的滤波重建值，为空间浓度的滤波初始值，y(t)为可测得的校正后的正弦图数据，P(t)为空间浓度的预估误差协方差阵，P(0)为初始空间浓度预估误差协方差，下标i＝1,…,r，将y(t)分成r块，y_i(t)是y(t)的第i分块，D_i,R_i,Q_i是与y_i(t)相容的矩阵D,R,Q的分块，K_i(t)为相应于y_i(t)的滤波增益矩阵，为基于直到观测y_i(t)的空间浓度滤波重建值，P_i(t)是相应于观测直到观测y_i(t)的滤波误差协方差阵。迭代从初始值P(0)出发，通过量测值y(t)，经过多次迭代，最终得到放射性浓度分布如图1所示，基于改进序贯滤波的图像重建迭代过程如下：

2.1)首先设定放射性浓度分布的初始值和初始方差P₀(0)，这两个值由经验值给出，但在没有先验知识的情况下，可将P₀(0)分别取为零矩阵及单位矩阵；

2.2)获取正弦图数据y(t)；

2.3)将正弦图数据y(t)分成r块；

2.4)利用方程(9)计算增益矩阵K_i(t)；

2.5)利用量测值y_i(t)及增益K_i(t)，根据状态更新方程(8)计算出空间浓度估计值并根据(10)推出相应的预估误差协方差阵P_i(t)；

2.6)若i＜r,i＝i+1，跳转至步骤2.4)，否则，P(t)＝P_r(t)；

2.7)若t＜N,t＝t+1，跳转至步骤2.2)，否则，算法结束，得到最终的浓度重构结果。