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

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

背景技术：
正电子发射断层成像(PositronEmissionTomography，PET)技术是一种通过在分子水平上对代谢变化进行追踪、监控，从而达到对疾病的早期检测、预防目的的医学成像技术。PET技术作为功能性成像技术，能反映病人体内新陈代谢的情况，因而能更早地探测出病灶。因此，PET技术被广泛应用于肿瘤、心脏疾病、神经和精神系统疾病的诊断及实验生物成像、药物筛选与开发，已显示出巨大的应用前景。在进行PET扫描时首先使用加速器产生正电子发射核素，然后将由放射性同位核素标记的药物注入人体内，通过血液循环，这些物质在人体内各组织器官中将形成一定的分布。由于放射性同位核素的半衰期较短，且极其不稳定，将很快发生衰变。在衰变过程中，放射性同位核素会产生正电子并与附近的自由电子发生湮灭反应，产生一对方向相反、能量相等的伽玛光子。这时，可以通过体外的探测器利用符合技术对光子进行探测，再经由噪声矫正，得到投影数据。通过一些重构方法利用投影数据进行反演求解，可重建出人体的放射性物质的空间浓度分布。目前，PET图像重建方法大致可分为两类：解析法和迭代统计法。前一类主要是滤波反投影法，计算速度快，代价小，但不能很好抑制噪声，以致重建图像质量不高。因此，出现了以极大似然法为代表的迭代统计法。由于迭代法基于统计学模型，对不完全数据适应性好，已逐渐成为PET重建算法研究的焦点。相比于解析法，迭代法获得的重建图像更加清晰。但对PET系统成像过程的准确建模是影响迭代法重建结果的关键。状态空间的引入，有力地推动了迭代法的发展。状态空间法能够将PET成像的过程模型化(给出PET成像过程的数学表达)，对PET扫描过程的统计特性与生物体的生理、结构特性进行相应的刻画，因而该模型联合一些已有的估计算法往往能够提供更好的重建效果。目前已有的状态空间求解方法主要基于H2滤波(即卡尔曼滤波)、H∞滤波等。然而，H2滤波重建需假设PET系统的噪声服从高斯分布，确知噪声的协方差阵，实际操作过程却很难准确给出噪声的协方差阵；H∞滤波重建仅假定噪声扰动能量有界，相对于实际中人们总能对PET系统噪声统计特性有一定的了解，这个假定又过于保守。因此，无论单纯使用H2滤波重建还是H∞滤波重建，都难以保证重建质量。

技术实现要素：
本发明要克服现有技术的上述缺点，提供一种基于H2/H∞混合滤波的静态PET图像重建方法。该发明提出了一个关于实际噪声性质的判别准则，根据该准则提出：当噪声统计特性相对精确时，可直接基于H2滤波结果重建图像；否则，基于H2、H∞滤波结果重建图像。该方法充分利用了H2、H∞滤波各自的优势，在保证PET重建图像鲁棒性能的前提下，提供了最优的重建效果。一种基于H2/H∞混合滤波的静态PET图像重建方法，包括如下步骤：(1)建立PET系统的状态空间模型：其中，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)，又可描述为能量有界的噪声。此时，为了得到更好的重建结果，可引入判别准则，分析噪声统计特性的准确度，基于H2和H∞两种滤波结果进行图像重建，获得最终PET重建图像。(2)根据下列方程得到基于H2滤波的重建图像：K1(t)＝P1(t)DT(DP1(t)DT+R)-1(3)其中，K1(t)为t时刻的H2滤波增益矩阵，为t时刻的放射性浓度的滤波重建值，y(t)为t时刻可测得的校正后的正弦图数据，P1(t)为t时刻放射性浓度的预估误差协方差阵。迭代从初始值P1(0)出发，通过量测值y(t)，经过t次迭代，最终得到放射性浓度分布的估计值具体的重建过程如下：1)首先设定放射性浓度分布的初始值和初始估计协方差P1(0)；2)利用方程(3)计算增益矩阵K1(t)；3)利用量测值y(t)及增益K1(t)，根据状态更新方程(2)计算出当前时刻的放射性浓度估计值并根据(4)推导出下一时刻的预估误差协方差阵P1(t+1)；4)重复步骤2)，3)，直至获得最终重建结果。(3)利用H∞滤波，根据下列方程得到基于H∞滤波的重建图像：K2(t)＝P2(t)DT(I+DP2(t)DT)-1(7)其中，K2(t)为t时刻的H∞滤波增益矩阵，为t时刻放射性浓度的H∞滤波重建值，为t时刻放射性浓度的H∞预估重建值，y(t)为可测得的校正后的正弦图数据，P2(t)为t时刻空间浓度的H∞预估误差协方差阵，γ是给定的噪声抑制参数。迭代从初始值P2(0)出发，通过量测值y(t)，经过t次迭代，最终得到放射性浓度分布的滤波值具体的图像重建迭代过程如下：1)首先设定放射性浓度分布的初始值和初始协方差P2(0)；2)根据方程(7)计算增益矩阵K2(t)；3)利用量测值y(t)与增益矩阵K2(t)，分别根据方程(5)(6)计算当前时刻放射性浓度的空间分布滤波值和其在下一时刻的预测值利用方程(8)(9)计算下一时刻的预估误差协方差阵P2(t+1)；4)重复步骤2)，3)直至获得最终重建结果。(4)基于H∞滤波重建结果和H2滤波重建结果进行最终的图像重建。根据判别准则得到最终的图像重建b(t)＝[γ2I-(P2(t)-1+DTD)-1]-1(14)具体过程包括以下步骤：1)将H2滤波重建结果代入不等式(10)的左边，验证是否能使不等式(10)成立；2)若能使(10)成立，则最终重建结果3)否则，利用式(12)(13)(14)计算权重值θ，根据θ得到混合滤波重建结果本发明的技术构思为：基于对PET系统成像原理及噪声矫正方法的分析，对噪声进行了全面、合理的刻画，充分利用了两种滤波重建方法的优势，根据所提出的判别准则，找到实际上更适合噪声的重建方法，最终保证了由于噪声特性误判而导致的PET重建质量下降问题。从上述技术方案可以看出，本发明的有益效果主要表现在：基于状态空间模型，更客观、真实地刻画PET系统的噪声特性，并将H2和H∞两种滤波器的优势充分结合，有效地避免了因对噪声特性误判而导致的PET重建质量下降问题，保证了PET重建图像的可靠性。与现有重建方法相比，能够提供具有鲁棒性的最优重建结果。附图说明图1为本发明PET图像重建方法的步骤流程示意图。图2为混合滤波的步骤流程示意图。图3为心脏模型示意图。图4为噪声参数相对准确情况下，分别利用H2方法，H∞方法和本发明方法的重建结果的对比示意图。图5为噪声参数相对不准确情况下，分别利用H2方法，H∞方法和本发明方法的重建结果的对比示意图。具体实施方式为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的静态PET浓度重建方法进行详细说明。如图1所示，一种基于H2/H∞混合滤波的静态PET图像重建方法，包括如下步骤：(1)根据PET成像原理，建立状态空间体系：其中，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)利用H2滤波，根据下列方程得到基于H2滤波的重建图像：K1(t)＝P1(t)DT(DP1(t)DT+R)-1(17)其中，K1(t)为H2滤波增益矩阵，为放射性浓度的滤波重建值，为放射性浓度的滤波初始值，y(t)为可测得的校正后的正弦图数据，P1(t)为放射性浓度的预估误差协方差阵，P1(0)为初始放射性浓度预估误差协方差。迭代从初始值P1(0)出发，结合量测值y(t)，经过t次迭代，最终得到放射性浓度分布的滤波值即重建图像如图2所示，基于H2滤波的图像重建迭代过程如下：1)首先设定放射性浓度分布的初始值和初始方差P1(0)，这两个值由经验值给出，但在没有先验知识的情况下，可将P1(0)分别取为零矩阵及单位矩阵；2)利用方程(17)计算增益矩阵K1(t)；3)利用量测值y(t)和增益矩阵K1(t)，根据状态更新方程(16)推算出当前时刻的放射性浓度空间分布H2滤波估计并根据(18)计算出下一时刻的H2预估误差协方差阵P1(t+1)；4)重复步骤2)，3)，4)直至获得基于H2滤波的最终重建结果。(3)利用H∞滤波，根据下列方程得到基于H∞滤波的重建图像：K2(t)＝P2(t)DT(I+DP2(t)DT)-1(21)其中，K2(t)为H∞滤波增益矩阵，为放射性浓度的H∞滤波重建值，为放射性浓度的H∞预估重建值，y(t)为可测得的校正后的正弦图数据，P2(t)为空间浓度的H∞预估误差协方差阵，γ是给定的噪声抑制参数。迭代从初始值P2(0)出发，通过量测值y(t)，经过多次迭代，最终得到放射性浓度分布的滤波值基于H∞滤波的图像重建迭代过程如下：1)首先设定放射性浓度分布的初始值和初始方差P2(0)以及参数γ，其中P2(0)由经验值给出，但在没有先验知识的情况下，可将其分别取为零矩阵及单位矩阵；2)利用方程(21)计算增益矩阵K2(t)；3)利用量测值y(t)和增益矩阵K2(t)，分别根据(19)(20)计算当前时刻放射性浓度的空间分布滤波值和下一时刻的放射性浓度的空间分布预估值利用方程(22)(23)计算下一时刻的H∞预估误差协方差P2(t+1)；4)重复步骤2)，3)直至获得基于H∞滤波的最终重建结果需要注意的是，所选择的γ值应保证方程(22)有一个正定解P2(t)、保证[γ2I-(P2(t)-1+DTD)-1]正定。(4)基于H2滤波结果和H∞滤波结果进行最终的图像重建验证H2滤波重建结果是否满足H∞滤波型的约束条件(25)，获得最终混合滤波图像重建结果b(t)＝[γ2I-(P2(t)-1+DTD)-1]-1(28)如图4所示，最终的过程包括以下步骤：1)将H2滤波重建结果代入(25)中的不等式左端，判断不等式是否成立；2)若成立，则重建结果否则，利用式(26)(27)(28)计算出权重θ，继而得到重建结果本发明技术采用Zubal胸腔体模模拟验证算法的有效性。心脏模型如图3所示，图像的原始分辨率为64×64像素。心脏模型正弦图是将心脏体模Zubal进行PET扫描后得到的，其分辨率也为64×64像素。针对该模型，采用本实施方式，与单纯使用H2滤波方法、H∞滤波方法的重建结果进行对比。为突出本发明方法的优点，给出两组噪声参数进行对比：一组在噪声参数相对准确的情况下重建，重建结果如图4所示；而另一组则为噪声参数相对不准确的情况，重建结果如图5所示。从两组重建结果中可以看出，不论给定的噪声统计特性是否准确，相较于单纯使用H2滤波或H∞滤波重建方法，本发明方法总能有效避免因对噪声性质误判而导致的重建质量下降问题。