CN103263275A - 一种基于h无穷滤波的pet生理参数的重构方法 - Google Patents

Abstract

本发明公开了一种基于H无穷滤波的PET生理参数的重构方法，包括：（1）采集符合计数；（2）通过血输入函数求取PET血液浓度分布数据；（3）构建PET状态空间方程；（4）利用H无穷滤波算法估算出PET生理参数。本发明把状态空间理论和PET动态重建技术相结合，建立动态的PET状态空间模型，通过H无穷滤波求解，实现了房室模型生理参数的重建。本发明一方面将生理参数重建和PET图像重建置于统一的框架之中，实现了生理参数和PET图像重建的统一求解；另一方面，H无穷理论不对系统响应和测量的统计特性进行假设，不受噪声类型的影响，有效提高了重建的鲁棒性。

Description

一种基于H无穷滤波的PET生理参数的重构方法

技术领域

本发明属于PET图像技术领域，具体涉及一种基于H无穷滤波的PET生理参数的重构方法。

背景技术

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

静态正电子发射断层成像认为在符合计数采集过程中示踪物质的分布是稳定不变的，这是一种近似的假设，从而重建的放射性活度图像可以看成是示踪物质在静态时间窗内的平均分布。事实上，示踪物质进入活体开始就不断参与生理化学分解与合成过程中，其分布也是在不断改变的，而在PET扫描过程中始终维持不变的是组织代谢水平。因此在某些情况下，我们不仅需要得到示踪物质在生物体内的空间分布，更希望对生物体组织或器官的真实代谢水平进行定量分析。

房室模型分析是运用PET动态成像进行示踪动力学研究的基础，在过去几十年中，房室参数成为活体代谢代谢定量分析中一个重要准则，包括血流定量，脑部葡萄糖代谢率定量，神经受体配基结合定量等等，房室建模在医学分析与诊断领域也得到了越来越广泛的应用。

在单组织房室模型中，示踪剂血液流向组织和组织流向血液的一阶速率常数（k₁、k₂）表征了示踪剂在活体中代谢水平，是反映活体代谢功能的重要生理参数之一。把含有反映组织功能状态的参数耦合到生理模型，结合动态测量数据追踪生理代谢过程，再通过计算机数值方法处理，就能估计得到定量参数本身。

对于结合房室模型的动态PET生理参数重建，目前主要有两类算法：直接重建和间接重建。间接重建的方法需要首先重建每个测量时间段的放射性浓度图像，接着使用该图像来估计动力学参数。间接重建的主要问题在于每个时间段采集的低信噪比数据重建的放射性活度图像通常带有较大的噪声，运用带噪的图像来估计动力学参量会引入更大的系统误差。

在直接重建的研究中，Carson和Lange在标题为Alternatives to voxels forimage representation for iterative reconstruction algorithms（Physics in Medicine andBiology，1992）的文献中提出了使用期望最大化（EM）迭代算法的直接重建框架；Kamasak等在标题为Direct reconstruction of kinetic parameter images fromdynamic PET data（IEEE Transactions on Medical Imaging，2005）的文献中提出了的基于贝叶斯估值的参量化迭代坐标下降算法（Parametric Iterative CoordinateDescent，简称PICD）。这些直接重建的方法通过PET采集的正弦数据直接重建动力学参量的空间分布；相比之下，直接重建算法难度较大，但却更为直观准确。不过这些方法都是基于对统计测量和物理系统的模型假设（例如认为测量服从高斯或者泊松分布），因而好的重建结果依赖于对统计模型和系统响应准确把握。然而测量的统计特性并不是完全确定的，同时系统响应也有其不确定性，这都降低了对动力学参数重建的准确程度，成为直接重建算法的所面临的难题。

发明内容

针对现有技术所存在的上述技术问题，本发明提供了一种基于H无穷滤波的PET生理参数的重构方法，能够克服PET成像中统计不确定性的难题，提高重构的鲁棒性和准确性。

一种基于H无穷滤波的PET生理参数的重构方法，包括如下步骤：

（1）利用探测器对注入有放射性物质的生物组织进行探测，动态采集得到PET的n组符合计数；n为大于1的自然数；

（2）通过对生物组织进行动脉血采样，拟合出生物组织的血输入函数，并通过血输入函数计算出放射性物质在血液中与符合计数对应的n组浓度分布数据u₁～u_n；

（3）根据PET房室模型建立含有生理参数的PET状态空间方程；

（4）利用H无穷滤波方法求解所述的PET状态空间方程得到以下迭代方程；根据所述的符合计数通过以下迭代方程估计出PET生理参数；

$\left[\begin{array}{c}{x}_k\\ {\mathrm{\θ}}_k\end{array}\right]=\left[\begin{array}{c}{x}_{k-1}\\ {\mathrm{\θ}}_{k-1}\end{array}\right]+\left[\begin{array}{cc}0& A_k\\ 0& 0\end{array}\right]\left[\begin{array}{c}{x}_{k-1}\\ {\mathrm{\θ}}_{k-1}\end{array}\right]+B_k(y_k-{\mathrm{Dx}}_{k-1})$

其中：x_k和x_k-1分别为第k次和第k-1次迭代后的PET浓度分布向量，θ_k和θ_k-1分别为第k次和第k-1次迭代后的生理参数向量；A_k为第k次迭代后的系统转移矩阵，B_k为第k次迭代后的增益矩阵，D为系统矩阵，y_k为第k组符合计数，k为自然数且1≤k≤n。

所述的步骤（4）中，根据符合计数通过迭代方程进行迭代计算，则迭代n次后输出的第n次迭代后的生理参数向量θ_n即为PET生理参数。

所述的PET状态空间方程的表达式如下：

$\left[\begin{array}{c}\mathrm{dx}\left(t\right)/\mathrm{dt}\\ \mathrm{d\θ}/\mathrm{dt}\end{array}\right]=\left[\begin{array}{c}A\left(t\right)\mathrm{\θ}\\ 0\end{array}\right]+v\left(t\right)$   A(t)＝[-x(t)  u(t)]

$y\left(t\right)=H\left[\begin{array}{c}x\left(t\right)\\ \mathrm{\θ}\end{array}\right]+e\left(t\right)$   H＝[D  0]

其中：x(t)为t时刻的PET浓度分布向量，θ为PET生理参数，A(t)为t时刻的系统转移矩阵，v(t)和e(t)分别为t时刻的过程噪声向量和测量噪声向量，u(t)为t时刻放射性物质在血液中的浓度分布数据，y(t)为t时刻PET的符合计数，t为时间。

所述的第k次迭代后的系统转移矩阵A_k通过以下算式求得：

A_k＝[-x_k-1  u_k]

其中：u_k为放射性物质在血液中对应第k组符合计数y_k的浓度分布数据。

所述的第k次迭代后的增益矩阵B_k通过以下算式求得：

B_k＝Z_k ^-1H^T  H＝[D  0]

其中：Z_k为对应第k组符合计数y_k的中间状态矩阵。

所述的中间状态矩阵Z_k通过以下方法求得：

首先，利用H无穷滤波方法求解PET状态空间方程得到关于中间状态矩阵的Riccati（黎卡提）微分方程；

然后，根据所述的Riccati微分方程求解出每一时刻对应的中间状态矩阵；

最后，根据第k组符合计数y_k的采样时刻，提取对应时刻的中间状态矩阵作为中间状态矩阵Z_k。

所述的Riccati微分方程表达式如下：

$\frac{\mathrm{dZ}\left(t\right)}{\mathrm{dt}}=-Z\left(t\right)\left[\begin{array}{cc}0& A\left(t\right)\\ 0& 0\end{array}\right]-\left[\begin{array}{cc}0& 0\\ {A\left(t\right)}^T& 0\end{array}\right]Z\left(t\right)+\left[\begin{array}{cc}{D}^{T}D& 0\\ 0& -{\mathrm{\γ}}^{2}Q\end{array}\right]-Z\left(t\right)\left[\begin{array}{cc}I& 0\\ 0& 0\end{array}\right]Z\left(t\right)$

A(t)＝[-x(t)  u(t)]

其中：Z(t)为t时刻的中间状态矩阵，A(t)为t时刻的系统转移矩阵，x(t)为t时刻的PET浓度分布向量，u(t)为t时刻放射性物质在血液中的浓度分布数据，γ为预设的噪声干扰上限值，Q为给定的权重矩阵，I为单位矩阵，t为时间。

所述的血输入函数是通过动态获取动脉的血液样本并测定其中的放射性浓度，得到动脉血液中的时间放射性浓度曲线（Time Activity Curve，TAC），作为数学模型的输入函数。由于动脉血采样是损伤性的，完整地对各个时刻的动脉血进行采样并不可取，通常采用近似的拟合函数法，根据拟合函数，只需几个时刻点的动脉血样来确定函数表达式中的每个参数，便可以获取完整的血输入函数。

本发明把状态空间理论和PET动态重建技术相结合，建立动态的PET状态空间模型，通过H无穷滤波求解，实现了房室模型生理参数的重建。本发明一方面将生理参数重建和PET图像重建置于统一的框架之中，实现了生理参数和PET图像重建的统一求解；另一方面，H无穷理论不对系统响应和测量的统计特性进行假设，不受噪声类型的影响，有效提高了重建的鲁棒性，特别适合于测量数据的统计特性以及系统不确定性都无法确定的应用情况。

附图说明

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

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

如图1所示，一种基于H无穷滤波的PET生理参数的重构方法，包括如下步骤：

（1）采集符合计数。

利用探测器对注入有放射性物质的生物组织进行探测，动态采集得到PET的n组符合计数；本实施方式中，探测器采用日本滨松公司生产的型号为SHR74000的PET扫描仪。

PET扫描仪探测人体内发出的放射性信号，经过符合和采集系统处理，形成原始符合事件。PET探测器记录的符合事件包括真符合、随机符合和散射符合。通过探测器的延时窗口和能量窗口对随机事件和散射事件进行校正，而后进行衰减校正，得到正弦图数据即校正后的符合计数。

（2）通过血输入函数求取PET血液浓度分布数据。

通过对生物组织进行动脉血采样，拟合出生物组织的血输入函数，并通过血输入函数计算出放射性物质在血液中与符合计数对应的n组浓度分布数据u₁～u_n；

血输入函数是通过动态获取动脉的血液样本并测定其中的放射性浓度，得到动脉血液中的时间放射性浓度曲线，作为数学模型的输入函数。由于动脉血采样是损伤性的，完整地对各个时刻的动脉血进行采样并不可取，通常采用近似的拟合函数法，根据拟合函数，只需几个时刻点的动脉血样来确定函数表达式中的每个参数，便可以获取完整的血输入函数。

（3）构建PET状态空间方程。

在动态PET成像中，房室模型动态交换可以通过一阶微分方程来描述。在单组织房室模型中，一阶微分方程表示为：

$\frac{{\mathrm{dC}}_T\left(t\right)}{\mathrm{dt}}=k_1{C}_P\left(t\right)-k_2{C}_T\left(t\right)$

其中，C_T(t)和C_P(t)分别表示示踪化合物在组织中和血液中的浓度，k₁和k₂分别为放射性示踪剂血液流向组织和组织流向血液的一阶速率常数，且均为m维向量。

忽略血液体积比，则可得到微分方程：

dx(t)/dt＝k₁u(t)-k₂x(t)

其中，x(t)和u(t)分别表示示踪化合物在组织中和血液中的浓度。

令[k₂  k₁]^T＝θ，根据PET房室模型建立含有生理参数的PET状态空间方程如下：

$\left[\begin{array}{c}\mathrm{dx}\left(t\right)/\mathrm{dt}\\ \mathrm{d\θ}/\mathrm{dt}\end{array}\right]=\left[\begin{array}{c}A\left(t\right)\mathrm{\θ}\\ 0\end{array}\right]+v\left(t\right)$   A(t)＝[-x(t)  u(t)]

$y\left(t\right)=H\left[\begin{array}{c}x\left(t\right)\\ \mathrm{\θ}\end{array}\right]+e\left(t\right)$   H＝[D  0]

其中：x(t)为t时刻的PET浓度分布向量，θ为PET生理参数，A(t)为t时刻的系统转移矩阵，v(t)和e(t)分别为t时刻的过程噪声向量和测量噪声向量，D为系统矩阵，u(t)为t时刻放射性物质在血液中的浓度分布数据，y(t)为t时刻PET的符合计数。

PET生理参数θ为常量，其为由k₁和k₂组成的2m维向量；PET浓度分布向量x以及放射性物质在血液中的浓度分布数据u均为m维向量，符合计数y为p维向量；系统矩阵D为p×m维矩阵，其表征了发射光子被探测器接收的概率，其受探测器结构、探测效率、衰减、死时间等因素的影响；过程噪声向量v表征了状态之间转移的统计不确定性，测量噪声向量e代表符合事件采集过程中由随机事件、散射事件等引入的测量误差；这两项噪声大小决定了Riccati微分方程中噪声干扰上限值γ和权重矩阵Q的设定，系统转移矩阵A为m×2维矩阵，其由放射性物质在组织中浓度x和血液中浓度u决定。

（4）利用H无穷滤波算法估算出PET生理参数。

利用H无穷滤波方法求解PET状态空间方程得到以下迭代方程；根据所符合计数通过以下迭代方程估计出PET生理参数；

$\left[\begin{array}{c}{x}_k\\ {\mathrm{\θ}}_k\end{array}\right]=\left[\begin{array}{c}{x}_{k-1}\\ {\mathrm{\θ}}_{k-1}\end{array}\right]+\left[\begin{array}{cc}0& A_k\\ 0& 0\end{array}\right]\left[\begin{array}{c}{x}_{k-1}\\ {\mathrm{\θ}}_{k-1}\end{array}\right]+B_k(y_k-{\mathrm{Dx}}_{k-1})$

A_k＝[-x_k-1  u_k]

B_k＝Z_k ^-1H^T  H＝[D  0]

其中：x_k和x_k-1分别为第k次和第k-1次迭代后的PET浓度分布向量，θ_k和θ_k-1分别为第k次和第k-1次迭代后的生理参数向量；A_k为第k次迭代后的系统转移矩阵，B_k为第k次迭代后的增益矩阵，D为系统矩阵，y_k为第k组符合计数，k为自然数且1≤k≤n，u_k为放射性物质在血液中对应第k组符合计数y_k的浓度分布数据，Z_k为对应第k组符合计数y_k的中间状态矩阵。

迭代n次后输出的第n次迭代后的生理参数向量θ_n即为PET生理参数，x₀、θ₀和u₀均为随机初始生成。

中间状态矩阵Z_k通过以下方法求得：

首先，利用H无穷滤波方法求解PET状态空间方程得到关于中间状态矩阵的Riccati微分方程如下：

$\frac{\mathrm{dZ}\left(t\right)}{\mathrm{dt}}=-Z\left(t\right)\left[\begin{array}{cc}0& A\left(t\right)\\ 0& 0\end{array}\right]-\left[\begin{array}{cc}0& 0\\ {A\left(t\right)}^T& 0\end{array}\right]Z\left(t\right)+\left[\begin{array}{cc}{D}^{T}D& 0\\ 0& -{\mathrm{\γ}}^{2}Q\end{array}\right]-Z\left(t\right)\left[\begin{array}{cc}I& 0\\ 0& 0\end{array}\right]Z\left(t\right)$

A(t)＝[-x(t)  u(t)]

其中：Z(t)为t时刻的中间状态矩阵，A(t)为t时刻的系统转移矩阵，x(t)为t时刻的PET浓度分布向量，u(t)为t时刻放射性物质在血液中的浓度分布数据，γ为预设的噪声干扰上限值（γ>1），Q为给定的权重矩阵（Q为(m+2)×(m+2)维矩阵），I为单位矩阵。

然后，根据Riccati微分方程求解出每一时刻对应的中间状态矩阵；

最后，根据第k组符合计数y_k的采样时刻，提取对应时刻的中间状态矩阵作为中间状态矩阵Z_k。

以下我们采用合成图像模拟实验来验证本实施方式重建PET生理参数的可行性和鲁棒性；原始图像分辨率为128×128象素，其具有四个感兴趣区域ROI。

我们模拟180度旋转角下90簇投影，每簇含64条平行投影线的PET系统。正弦图为35分钟内30个离散时间窗数据，包括10个30秒窗间隔数据，10个60秒窗间隔数据和10个120秒窗间隔数据；血输入函数使用经验数学方程：

u(t)＝(A₁t-A₂-A₃)e^αt+A₂e^βt+A₃e^γt

其中：A₁=851.125μCi/mL/min，A₂=20.8113μCi/mL/min，A₃=21.8798μCi/mL/min，α=-4.133859/min，β=-0.04043449/min，γ=0.1190996/min。

单组织房室模型假设的速率常数真值如表1所示。本实施方式模拟PET测量数据的步骤为：每个时间窗的活度图像分别投影成正弦图中的真符合事件；然后分别加入40％的随机符合噪声和-10%的散射符合噪声；最后通过总计数事件和噪声符合时间相减获得预校正的正弦图。

表1

	ROI1	ROI2	ROI3	ROI4
					k1	0.95	0.75	0.55	0.35
k2	0.96	0.76	0.56	0.36

我们分别以本实施方式与加权最小二乘法重建作为比较，为了定量比较两种方法的误差，定义相对误差函数为：

$\mathrm{RMSE}=\sqrt{\frac{\frac{1}{N_p}\mathrm{\Σ}{(k-k_{\mathrm{tr}})}^2}{\frac{1}{N_p}\mathrm{\Σ}{(k_{\mathrm{wls}}-k_{\mathrm{tr}})}^2}}$

其中：N_p表示总像素数目，k_tr、k和k_wls分别表示真值、本实施方式重建值和加权最小二乘法重建值。通过对k₁和k₂计算的RMSE分别为0.19和0.43，这两个数值都比较小（小于1），故说明了本实施方式对生理参数的重建值较加权最小二乘法方法更接近于真实值，具有较好的准确性和鲁棒性。

Claims (7)

1.一种基于H无穷滤波的PET生理参数的重构方法，包括如下步骤：

（1）利用探测器对注入有放射性物质的生物组织进行探测，动态采集得到PET的n组符合计数，n为大于1的自然数；

（2）通过对生物组织进行动脉血采样，拟合出生物组织的血输入函数，并通过血输入函数计算出放射性物质在血液中与符合计数对应的n组浓度分布数据u₁～u_n；

（3）根据PET房室模型建立含有生理参数的PET状态空间方程；

（4）利用H无穷滤波方法求解所述的PET状态空间方程得到以下迭代方程；根据所述的符合计数通过以下迭代方程估计出PET生理参数；

$\left[\begin{array}{c}{x}_k\\ {\mathrm{\θ}}_k\end{array}\right]=\left[\begin{array}{c}{x}_{k-1}\\ {\mathrm{\θ}}_{k-1}\end{array}\right]+\left[\begin{array}{cc}0& A_k\\ 0& 0\end{array}\right]\left[\begin{array}{c}{x}_{k-1}\\ {\mathrm{\θ}}_{k-1}\end{array}\right]+B_k(y_k-{\mathrm{Dx}}_{k-1})$

其中：x_k和x_k-1分别为第k次和第k-1次迭代后的PET浓度分布向量，θ_k和θ_k-1分别为第k次和第k-1次迭代后的生理参数向量；A_k为第k次迭代后的系统转移矩阵，B_k为第k次迭代后的增益矩阵，D为系统矩阵，y_k为第k组符合计数，k为自然数且1≤k≤n。

2.根据权利要求1所述的重构方法，其特征在于：所述的步骤（4）中，根据符合计数通过迭代方程进行迭代计算，则迭代n次后输出的第n次迭代后的生理参数向量θ_n即为PET生理参数。

3.根据权利要求1所述的重构方法，其特征在于：所述的PET状态空间方程的表达式如下：

$\left[\begin{array}{c}\mathrm{dx}\left(t\right)/\mathrm{dt}\\ \mathrm{d\θ}/\mathrm{dt}\end{array}\right]=\left[\begin{array}{c}A\left(t\right)\mathrm{\θ}\\ 0\end{array}\right]+v\left(t\right)$   A(t)＝[-x(t)  u(t)]

$y\left(t\right)=H\left[\begin{array}{c}x\left(t\right)\\ \mathrm{\θ}\end{array}\right]+e\left(t\right)$   H＝[D  0]

其中：x(t)为t时刻的PET浓度分布向量，θ为PET生理参数，A(t)为t时刻的系统转移矩阵，v(t)和e(t)分别为t时刻的过程噪声向量和测量噪声向量，u(t)为t时刻放射性物质在血液中的浓度分布数据，y(t)为t时刻PET的符合计数，t为时间。

4.根据权利要求1所述的重构方法，其特征在于：所述的第k次迭代后的系统转移矩阵A_k通过以下算式求得：

A_k＝[-x_k-1  u_k]

其中：u_k为放射性物质在血液中对应第k组符合计数y_k的浓度分布数据。

5.根据权利要求1所述的重构方法，其特征在于：所述的第k次迭代后的增益矩阵B_k通过以下算式求得：

B_k＝Z_k ^-1H^T  H＝[D  0]

其中：Z_k为对应第k组符合计数y_k的中间状态矩阵。

6.根据权利要求5所述的重构方法，其特征在于：所述的中间状态矩阵Z_k通过以下方法求得：

首先，利用H无穷滤波方法求解PET状态空间方程得到关于中间状态矩阵的Riccati微分方程；

然后，根据所述的Riccati微分方程求解出每一时刻对应的中间状态矩阵；

最后，根据第k组符合计数y_k的采样时刻，提取对应时刻的中间状态矩阵作为中间状态矩阵Z_k。

7.根据权利要求6所述的重构方法，其特征在于：所述的Riccati微分方程表达式如下：

$\frac{\mathrm{dZ}\left(t\right)}{\mathrm{dt}}=-Z\left(t\right)\left[\begin{array}{cc}0& A\left(t\right)\\ 0& 0\end{array}\right]-\left[\begin{array}{cc}0& 0\\ {A\left(t\right)}^T& 0\end{array}\right]Z\left(t\right)+\left[\begin{array}{cc}{D}^{T}D& 0\\ 0& -{\mathrm{\γ}}^{2}Q\end{array}\right]-Z\left(t\right)\left[\begin{array}{cc}I& 0\\ 0& 0\end{array}\right]Z\left(t\right)$

A(t)＝[-x(t)  u(t)]

其中：Z(t)为t时刻的中间状态矩阵，A(t)为t时刻的系统转移矩阵，x(t)为t时刻的PET浓度分布向量，u(t)为t时刻放射性物质在血液中的浓度分布数据，γ为预设的噪声干扰上限值，Q为给定的权重矩阵，I为单位矩阵，t为时间。

Images