CN103793929A - 一种基于h2/h∞混合滤波的静态pet图像重建方法 - Google Patents

Abstract

本发明涉及一种基于H₂/H_∞混合滤波的静态PET图像重建方法，包括采集原始投影线的正弦图数据，并基于状态空间的H₂滤波、H_∞滤波结果重建PET图像。本发明为重建过程提出了一个重建准则，用于判断并选择更适合的PET噪声模型，从而能够充分利用H₂、H_∞两种滤波重建结果的优势，得到满足鲁棒性能的最小估计误差方差重建结果。由于能预先判断了PET系统的噪声特性，找到了更适合的滤波重建方法，所以能够提高重建图像的质量。

Description

一种基于H2/H∞混合滤波的静态PET图像重建方法

技术领域

本发明涉及正电子发射断层成像技术领域，具体涉及一种基于H₂/H_∞混合滤波的静态PET图像重建方法。

背景技术

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

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

目前，PET图像重建方法大致可分为两类：解析法和迭代统计法。前一类主要是滤波反投影法，计算速度快，代价小，但不能很好抑制噪声，以致重建图像质量不高。因此，出现了以极大似然法为代表的迭代统计法。由于迭代法基于统计学模型，对不完全数据适应性好，已逐渐成为PET重建算法研究的焦点。相比于解析法，迭代法获得的重建图像更加清晰。但对PET系统成像过程的准确建模是影响迭代法重建结果的关键。状态空间的引入，有力地推动了迭代法的发展。状态空间法能够将PET成像的过程模型化（给出PET成像过程的数学表达），对PET扫描过程的统计特性与生物体的生理、结构特性进行相应的刻画，因而该模型联合一些已有的估计算法往往能够提供更好的重建效果。目前已有的状态空间求解方法主要基于H₂滤波(即卡尔曼滤波)、H_∞滤波等。然而，H₂滤波重建需假设PET系统的噪声服从高斯分布，确知噪声的协方差阵，实际操作过程却很难准确给出噪声的协方差阵；H_∞滤波重建仅假定噪声扰动能量有界，相对于实际中人们总能对PET系统噪声统计特性有一定的了解，这个假定又过于保守。因此，无论单纯使用H₂滤波重建还是H_∞滤波重建，都难以保证重建质量。

发明内容

本发明要克服现有技术的上述缺点，提供一种基于H₂/H_∞混合滤波的静态PET图像重建方法。该发明提出了一个关于实际噪声性质的判别准则，根据该准则提出：当噪声统计特性相对精确时，可直接基于H₂滤波结果重建图像；否则，基于H₂、H_∞滤波结果重建图像。该方法充分利用了H₂、H_∞滤波各自的优势，在保证PET重建图像鲁棒性能的前提下，提供了最优的重建效果。

一种基于H₂/H_∞混合滤波的静态PET图像重建方法，包括如下步骤：

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

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

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

所述的状态转移矩阵A、过程噪声v(t)在静态PET成像中一般可分别取为单位矩阵和零，但考虑到没有绝对的静态，浓度分布总存在一定的变化，因此给出更为一般的模型，v(t)的均值为零，方差为Q。

所述过程噪声v(t)和量测噪声e(t)具有双重特性。这源于：实际测量的原始PET系统符合计数过程服从Poisson分布，包含随机符合、散射符合等，需经矫正方可使用，而矫正将一定程度上破坏符合计数的Poisson分布特性，最终导致符合计数的统计特性有所偏离。但直接摒弃已知的那部分统计特性、简单地把矫正后的符合计数噪声归于能量有界噪声未免过于保守。因此，可充分利用Poisson噪声和能量有界的噪声所提供的所有信息。因为Anscombe变换能将Poisson噪声转化为高斯噪声，所以系统噪声既可描述为高斯白噪声，即v(t)～N(0,Q)，e(t)～N(0,R)，又可描述为能量有界的噪声。此时，为了得到更好的重建结果，可引入判别准则，分析噪声统计特性的准确度，基于H₂和H_∞两种滤波结果进行图像重建，获得最终PET重建图像。

(2)根据下列方程得到基于H₂滤波的重建图像：

${\hat{x}}_1\left(t\right)=A{\hat{x}}_1(t-1)+K_1\left(t\right)(y\left(t\right)-\mathrm{DA}{\hat{x}}_1(t-1))---\left(2\right)$

K₁(t)=P₁(t)D^T(DP₁(t)D^T+R)^-1  (3)

$P_1\left(t\right)=A{P}_1(t-1)A^T+Q-A{K}_1\left(t\right)({\mathrm{DP}}_1\left(t\right)D^T+R)K_1^T\left(t\right)A^T---\left(4\right)$

其中，K₁(t)为t时刻的H₂滤波增益矩阵，

为t时刻的放射性浓度的滤波重建值，y(t)为t时刻可测得的校正后的正弦图数据，P₁(t)为t时刻放射性浓度的预估误差协方差阵。迭代从初始值

P₁(0)出发，通过量测值y(t)，经过t次迭代，最终得到放射性浓度分布的估计值

具体的重建过程如下：

1)首先设定放射性浓度分布的初始值和初始估计协方差

P₁(0)；

2)利用方程(3)计算增益矩阵K₁(t)；

3)利用量测值y(t)及增益K₁(t)，根据状态更新方程(2)计算出当前时刻的放射性浓度估计值

并根据(4)推导出下一时刻的预估误差协方差阵P₁(t+1)；

4)重复步骤2)，3)，直至获得最终重建结果。

(3)利用H_∞滤波，根据下列方程得到基于H_∞滤波的重建图像：

${\hat{x}}_2\left(t\right)={\hat{x}}_2^{-}\left(t\right)+K_2\left(t\right)(y\left(t\right)-D{\hat{x}}_2^{-}\left(t\right))---\left(5\right)$

${\hat{x}}_2^{-}(t+1)=A{\hat{x}}_2^{-}\left(t\right)+A{K}_2\left(t\right)(y\left(t\right)-\mathrm{DA}{\hat{x}}_2^{-}\left(t\right))---\left(6\right)$

K₂(t)=P₂(t)D^T(I+DP₂(t)D^T)^-1  (7)

$P_2(t+1)=A{P}_2\left(t\right)A^T+I-A{P}_2\left(t\right)\left[\begin{array}{cc}{D}^T& I\end{array}\right]{\stackrel{-}{R}}_2{\left(t\right)}^{-1}\left[\begin{array}{c}D\\ I\end{array}\right]P_2\left(t\right)A^T---\left(8\right)$

${\stackrel{-}{R}}_2\left(t\right)=\left[\begin{array}{cc}I& 0\\ 0& -{\mathrm{\γ}}^{2}I\end{array}\right]+\left[\begin{array}{c}D\\ I\end{array}\right]P_2\left(t\right)\left[\begin{array}{cc}{D}^T& I\end{array}\right]---\left(9\right)$

其中，K₂(t)为t时刻的H_∞滤波增益矩阵，

为t时刻放射性浓度的H_∞滤波重建值，

为t时刻放射性浓度的H_∞预估重建值，y(t)为可测得的校正后的正弦图数据，P₂(t)为t时刻空间浓度的H_∞预估误差协方差阵，γ是给定的噪声抑制参数。迭代从初始值

P₂(0)出发，通过量测值y(t)，经过t次迭代，最终得到放射性浓度分布的滤波值

具体的图像重建迭代过程如下：

1)首先设定放射性浓度分布的初始值和初始协方差P₂(0)；

2)根据方程(7)计算增益矩阵K₂(t)；

3)利用量测值y(t)与增益矩阵K₂(t)，分别根据方程(5)(6)计算当前时刻放射性浓度的空间分布滤波值

和其在下一时刻的预测值

利用方程(8)(9)计算下一时刻的预估误差协方差阵P₂(t+1)；

4)重复步骤2)，3)直至获得最终重建结果。

(4)基于H_∞滤波重建结果和H₂滤波重建结果进行最终的图像重建。根据判别准则

$J\left(t\right)+{\left[\begin{array}{c}y\left(t\right)-D{\hat{x}}_2^{-}\left(t\right)\\ {\hat{x}}_1\left(t\right)-{\hat{x}}_2^{-}\left(t\right)\end{array}\right]}^T{\stackrel{-}{R}}_2{\left(t\right)}^{-1}\left[\begin{array}{c}y\left(t\right)-D{\hat{x}}_2^{-}\left(t\right)\\ {\hat{x}}_1\left(t\right)-{\hat{x}}_2^{-}\left(t\right)\end{array}\right]\≥0---\left(10\right)$

得到最终的图像重建

$\mathrm{\θ}=\sqrt{\frac{a\left(t\right)}{b\left(t\right){|{\hat{x}}_1\left(t\right)-{\hat{x}}_2^{-}\left(t\right)|}^2}}---\left(12\right)$

$a\left(t\right)=J\left(t\right)+{(y\left(t\right)-D{\hat{x}}_2^{-}\left(t\right))}^T(I+{\mathrm{DP}}_2\left(t\right)D^T)(y\left(t\right)-D{\hat{x}}_2^{-}\left(t\right))---\left(13\right)$

b(t)=[γ²I-(P₂(t)^-1+D^TD)^-1]^-1  (14)

具体过程包括以下步骤：

1)将H₂滤波重建结果

代入不等式(10)的左边，验证

是否能使不等式(10)成立；

2)若

能使(10)成立，则最终重建结果

3)否则，利用式(12)(13)(14)计算权重值θ，根据θ得到混合滤波重建结果 $\hat{x}\left(t\right)=\mathrm{\θ}{\hat{x}}_1\left(t\right)+(1-\mathrm{\θ}){\hat{x}}_2\left(t\right).$

本发明的技术构思为：基于对PET系统成像原理及噪声矫正方法的分析，对噪声进行了全面、合理的刻画，充分利用了两种滤波重建方法的优势，根据所提出的判别准则，找到实际上更适合噪声的重建方法，最终保证了由于噪声特性误判而导致的PET重建质量下降问题。

从上述技术方案可以看出，本发明的有益效果主要表现在：基于状态空间模型，更客观、真实地刻画PET系统的噪声特性，并将H₂和H_∞两种滤波器的优势充分结合，有效地避免了因对噪声特性误判而导致的PET重建质量下降问题，保证了PET重建图像的可靠性。与现有重建方法相比，能够提供具有鲁棒性的最优重建结果。

附图说明

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

图2为混合滤波的步骤流程示意图。

图3为心脏模型示意图。

图4为噪声参数相对准确情况下，分别利用H₂方法，H_∞方法和本发明方法的重建结果的对比示意图。

图5为噪声参数相对不准确情况下，分别利用H₂方法，H_∞方法和本发明方法的重建结果的对比示意图。

具体实施方式

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

如图1所示，一种基于H₂/H_∞混合滤波的静态PET图像重建方法，包括如下步骤：

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

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

其中，t表示时间；y(t)为观测值，是经过噪声矫正后得到的正弦图数据；D为系统矩阵，表示人体内放射性浓度与PET扫描之间的投影关系，由PET装置的固有特性决定，在实验中采用单响应线模型近似计算得到；A是状态转移矩阵，考虑静态分布时A为单位阵；x(t)为放射性浓度分布，即需要重建的对象；v(t)是过程噪声；e(t)为数据采集并经过噪声矫正后残留的量测噪声；v(t)，e(t)服从正态高斯分布，即v(t)～N(0,Q)，e(t)～N(0,R)，已知噪声方差Q、R但它们未必准确。

(2)利用H₂滤波，根据下列方程得到基于H₂滤波的重建图像：

${\hat{x}}_1\left(t\right)=A{\hat{x}}_1(t-1)+K_1\left(t\right)(y\left(t\right)-\mathrm{DA}{\hat{x}}_1(t-1))---\left(16\right)$

K₁(t)=_P1(t)D^T(DP₁(t)D^T+R)^-1  (17)

$P_1\left(t\right)=A{P}_1(t-1)A^T+Q-A{K}_1\left(t\right)({\mathrm{DP}}_1\left(t\right)D^T+R)K_1^T\left(t\right)A^T---\left(18\right)$

其中，K₁(t)为H₂滤波增益矩阵，

为放射性浓度的滤波重建值，

为放射性浓度的滤波初始值，y(t)为可测得的校正后的正弦图数据，P₁(t)为放射性浓度的预估误差协方差阵，P₁(0)为初始放射性浓度预估误差协方差。迭代从初始值

P₁(0)出发，结合量测值y(t)，经过t次迭代，最终得到放射性浓度分布的滤波值即重建图像

如图2所示，基于H₂滤波的图像重建迭代过程如下：

1)首先设定放射性浓度分布的初始值和初始方差P₁(0)，这两个值由经验值给出，但在没有先验知识的情况下，可将

P₁(0)分别取为零矩阵及单位矩阵；

2)利用方程(17)计算增益矩阵K₁(t)；

3)利用量测值y(t)和增益矩阵K₁(t)，根据状态更新方程(16)推算出当前时刻的放射性浓度空间分布H₂滤波估计

并根据(18)计算出下一时刻的H₂预估误差协方差阵P₁(t+1)；

4)重复步骤2)，3)，4)直至获得基于H₂滤波的最终重建结果。

(3)利用H_∞滤波，根据下列方程得到基于H_∞滤波的重建图像：

${\hat{x}}_2\left(t\right)={\hat{x}}_2^{-}\left(t\right)+K_2\left(t\right)(y\left(t\right)-D{\hat{x}}_2^{-}\left(t\right))---\left(19\right)$

${\hat{x}}_2^{-}(t+1)=A{\hat{x}}_2^{-}\left(t\right)+A{K}_2\left(t\right)(y\left(t\right)-\mathrm{DA}{\hat{x}}_2^{-}\left(t\right))---\left(20\right)$

K₂(t)=P₂(t)D^T(I+DP₂(t)D^T)^-1  (21)

$P_2(t+1)=A{P}_2\left(t\right)A^T+I-A{P}_2\left(t\right)\left[\begin{array}{cc}{D}^T& I\end{array}\right]{\stackrel{-}{R}}_2{\left(t\right)}^{-1}\left[\begin{array}{c}D\\ I\end{array}\right]P_2\left(t\right)A^T---\left(22\right)$

${\stackrel{-}{R}}_2\left(t\right)=\left[\begin{array}{cc}I& 0\\ 0& -{\mathrm{\γ}}^{2}I\end{array}\right]+\left[\begin{array}{c}D\\ I\end{array}\right]P_2\left(t\right)\left[\begin{array}{cc}{D}^T& I\end{array}\right]---\left(23\right)$

其中，K₂(t)为H_∞滤波增益矩阵，

为放射性浓度的H_∞滤波重建值，

为放射性浓度的H_∞预估重建值，y(t)为可测得的校正后的正弦图数据，P₂(t)为空间浓度的H_∞预估误差协方差阵，γ是给定的噪声抑制参数。迭代从初始值

P₂(0)出发，通过量测值y(t)，经过多次迭代，最终得到放射性浓度分布的滤波值

如流程图3所示，基于H_∞滤波的图像重建迭代过程如下：

1)首先设定放射性浓度分布的初始值和初始方差

P₂(0)以及参数γ，其中

P₂(0)由经验值给出，但在没有先验知识的情况下，可将其分别取为零矩阵及单位矩阵；

2)利用方程(21)计算增益矩阵K₂(t)；

3)利用量测值y(t)和增益矩阵K₂(t)，分别根据(19)(20)计算当前时刻放射性浓度的空间分布滤波值

和下一时刻的放射性浓度的空间分布预估值

利用方程(22)(23)计算下一时刻的H_∞预估误差协方差P₂(t+1)；

4)重复步骤2)，3)直至获得基于H_∞滤波的最终重建结果

需要注意的是，所选择的γ值应保证方程(22)有一个正定解P₂(t)、保证[γ²I-(P₂(t)^-1+D^TD)^-1正定。

(4)基于H_∞滤波结果

和H₂滤波结果

进行最终的图像重建

验证H₂滤波重建结果是否满足H_∞滤波型的约束条件(25)，获得最终混合滤波图像重建结果

$J(t+1)=J\left(t\right)+{\left[\begin{array}{c}y\left(t\right)-D{\hat{x}}_2^{-}\left(t\right)\\ {\hat{x}}_1\left(t\right)-{\hat{x}}_2^{-}\left(t\right)\end{array}\right]}^T{\stackrel{-}{R}}_2{\left(t\right)}^{-1}\left[\begin{array}{c}y\left(t\right)-D{\hat{x}}_2^{-}\left(t\right)\\ {\hat{x}}_1\left(t\right)-{\hat{x}}_2^{-}\left(t\right)\end{array}\right]\≥0,J\left(0\right)=0---\left(25\right)$

$\mathrm{\θ}=\sqrt{\frac{a\left(t\right)}{b\left(t\right){|{\hat{x}}_1\left(t\right)-{\hat{x}}_2^{-}\left(t\right)|}^2}}---\left(26\right)$

$a\left(t\right)=J\left(t\right)+{(y\left(t\right)-D{\hat{x}}_2^{-}\left(t\right))}^T(I+{\mathrm{DP}}_2\left(t\right)D^T)(y\left(t\right)-D{\hat{x}}_2^{-}\left(t\right))---\left(27\right)$

$b\left(t\right)={[{\mathrm{\γ}}^{2}I-{(P_2{\left(t\right)}^{-1}+D^{T}D)}^{-1}]}^{-1}---\left(28\right)$

如图4所示，最终的过程包括以下步骤：

1)将H₂滤波重建结果

代入（25）中的不等式左端，判断不等式是否成立；

2)若成立，则重建结果

否则，利用式(26)(27)(28)计算出权重θ，继而得到重建结果 $\hat{x}\left(t\right)=\mathrm{\θ}{\hat{x}}_1\left(t\right)+(1-\mathrm{\θ}){\hat{x}}_2\left(t\right).$

本发明技术采用Zubal胸腔体模模拟验证算法的有效性。心脏模型如图5所示，图像的原始分辨率为64×64像素。心脏模型正弦图是将心脏体模Zubal进行PET扫描后得到的，其分辨率也为64×64像素。针对该模型，采用本实施方式，与单纯使用H₂滤波方法、H_∞滤波方法的重建结果进行对比。为突出本发明方法的优点，给出两组噪声参数进行对比：一组在噪声参数相对准确的情况下重建，重建结果如图4所示；而另一组则为噪声参数相对不准确的情况，重建结果如图5所示。从两组重建结果中可以看出，不论给定的噪声统计特性是否准确，相较于单纯使用H₂滤波或H_∞滤波重建方法，本发明方法总能有效避免因对噪声性质误判而导致的重建质量下降问题。

Claims (5)

1.一种基于H₂/H_∞混合滤波的静态PET图像重建方法，包括以下步骤： 

（1）通过数据采集及校正，得到系统的观测数据，建立相应的状态空间模型； 

（2）基于状态空间H₂滤波方法重构放射性浓度分布； 

（3）基于状态空间H_∞滤波方法重构放射性浓度分布； 

（4）提出H_∞型的判别准则，若基于H₂滤波的浓度重构结果满足H_∞型的判别准则，最终将基于H₂的滤波结果进行图像重建；否则，基于H₂、H_∞滤波的加权结果进行图像重建。 

2.根据权利要求1所述的PET图像重建方法，其特征在于：步骤（1）建立了如下的状态空间模型 

其中，t表示时间；y(t)是经过噪声矫正后得到的正弦图数据；D为系统矩阵，表示人体内放射性浓度与PET扫描之间的投影关系；A是状态转移矩阵；x(t)为放射性浓度分布，即需要重建的对象；v(t)是过程噪声；e(t)为数据采集并经噪声矫正后残留的噪声；v(t)，e(t)相互独立且服从均值为0、方差分别为Q、R的正态分布，即

考虑到实际操作中无法确知噪声的真实统计特性Q、R，同时作出合理假设，有限时间内噪声能量有界。 

3.根据权利要求1所述的PET图像重建方法，其特征在于：步骤（2）根据下列方程得到基于H₂滤波的重建结果： 

K₁(t)=P₁(t)D^T(DP₁(t)D^T+R)^-1  (3) 

其中，K₁(t)为t时刻的H₂滤波增益矩阵，为t时刻放射性浓度的滤波重建值， 

为放射性浓度的滤波初始值，y(t)为可测得的校正后的正弦图数据，P₁(t)为放射性浓度的预估误差协方差阵，P₁(0)为初始放射性浓度的预估误差协方差。迭代从初始值

P₁(0)出发，结合量测值y(t)，经过t次迭代，最终得到放射性浓度分布的估计值

。

4.根据权利要求1所述的PET图像重建方法，其特征在于：步骤（3）根据下列方程得到基于H_∞滤波的重建结果： 

K₂(t)=P₂(t)D^T(I+DP₂(t)D^T)^-1  (7) 

其中，K₂(t)为t时刻的H_∞滤波增益矩阵，

为t时刻放射性浓度的H_∞滤波重建值， 

为t时刻放射性浓度的H_∞预估重建值，

为放射性浓度的H_∞初始预估值，y(t)为可测得的校正后的正弦图数据，P₂(t)为空间浓度的H_∞预估误差协方差阵，P₂(0)为初始空间浓度的H_∞预估误差协方差阵，γ是给定的噪声抑制参数。迭代从初始值

P₂(0)出发，通过量测值y(t)，经过t次迭代，最终得到放射性浓度分布的估计值

5.根据权利要求1所述的PET图像重建方法，其特征在于：步骤（4）中考虑到实际操作中很难准确选取PET系统噪声参数、但仅考虑噪声能量有界又可能导致估计结果的保守性；因此，结合H₂、H_∞滤波重构各自的优势，提出了判别准则J(t+1)≥0  (10)其中， 

并根据判别准则按如下方式给出放射性浓度的最终重建结果： 

b(t)=[γ²I-(P₂(t)^-1+D^TD)^-1]^-1  (15) 

即先判断不等式（10）是否成立；如果成立，则得到最终的放射性浓度重建结果 

否则，根据(13)(14)(15)计算θ值，对两种滤波重建结果进行加权，得到最终的放射性浓度重建结果

Images