CN103295207A - 一种基于h无穷滤波的双示踪剂pet浓度的动态重建方法 - Google Patents

Abstract

本发明公开了一种基于H无穷滤波的双示踪剂PET浓度的动态重建方法，包括：（1）采集符合计数；（2）通过血输入函数求取双示踪剂血液浓度分布数据；（3）构建PET状态空间方程；（4）利用H无穷滤波算法估算出PET浓度分布；（5）重建两种示踪剂对应的PET浓度分布数据。由于光子的能量单一的物理特性，双示踪剂同时注入的实验中的单示踪剂成分一般无法直接分离；本发明通过求解独立平行房室模型引导的PET状态空间模型，实现了对同时刻注入的双示踪剂，在两种示踪剂信号成分高度混合的情况下进行同时重建。

Description

一种基于H无穷滤波的双示踪剂PET浓度的动态重建方法

技术领域

本发明属于PET图像技术领域，具体涉及一种基于H无穷滤波的双示踪剂PET浓度的动态重建方法。

背景技术

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

为了从图像中获得尽可能多的信息，在同一活体内进行双示踪剂和多示踪剂成像研究成为PET成像研究中的重要课题。传统方法需要对每种示踪剂进行独立注射、独立扫描以及独立成像，以使得各种示踪剂成分在衰变周期内完全不产生信号干扰。但是，活体的生理状态在长时间扫描过程中会发生各种变化，因此示踪剂独立操作的方法给使用动态成像进行神经药物方面的研究带来了困难；再者，长时间操作造成设备、药物等诸多成本负担的问题，也会给病人不适，所以并不适合在实际临床应用。

双注入-单扫描是处理PET双示踪剂成像中一种比较合理的技术。由于所有正电子湮灭产生符合光子的峰值能量都为511keV，因而从硬件方面难以直接分离每种示踪剂成分产生的计数光子。常用的双注入-单扫描方法正是通过对两种示踪剂进行短间隔（10-20分钟）分别注入的方法来减弱信号重叠效应的。例如Koeppe等人在标题为Evaluation of the reliability in kinetic analysis for dualtracerinjection of FDG and flumazenil PET Study（IEEE Nuclear Science SymposiumConference Record，2001）的文献中提出了采用碳11标记的双示踪剂对脑部成像做出了初步的尝试，然而该方法却难以实现在两种示踪剂信号成分高度混合情况下的分离。

发明内容

针对现有技术所存在的上述技术问题，本发明提供了一种基于H无穷滤波的双示踪剂PET浓度的动态重建方法，能够实现对同时刻注入的双示踪剂且成分高度混合的情况下的同时重建。

一种基于H无穷滤波的双示踪剂PET浓度的动态重建方法，包括如下步骤：

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

（2）通过对生物组织进行动脉血采样，拟合出两种示踪剂对应的血输入函数，并通过血输入函数计算出双示踪剂在血液中与符合计数对应的n组浓度分布向量C₁～C_n；

（3）根据PET平行房室模型建立基于双示踪剂的PET状态空间方程；

（4）利用H无穷滤波算法求解所述的PET状态空间方程得到以下迭代方程；根据所述的符合计数通过以下迭代方程估计出对应的PET浓度分布向量X_k；

X_k＝AX_k-1+H_k(Y_k-DX_k-1)+BC_k

其中：X_k和X_k-1分别为第k次和第k-1次迭代后的PET浓度分布向量；A和B均为状态转移矩阵，H_k为第k次迭代后的增益矩阵，D为系统矩阵，Y_k为第k组符合计数，C_k为对应第k组符合计数Y_k的浓度分布向量，k为自然数且1≤k≤n；

（5）根据PET浓度分布向量X_k，动态重建出两种示踪剂对应的PET浓度分布数据。

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

$\frac{\mathrm{dX}\left(t\right)}{\mathrm{dt}}=\mathrm{AX}\left(t\right)+\mathrm{BC}\left(t\right)+v$

Y(t)＝DX(t)+e

其中：X(t)为t时刻的PET浓度分布向量，A和B均为状态转移矩阵，D为系统矩阵，C(t)为t时刻双示踪剂在血液中的浓度分布向量，Y(t)为t时刻PET的符合计数，t为时间，v和e分别为过程噪声向量和测量噪声向量。

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

H_k＝Z_k(I+D^TE^-1DZ_k)^-1D^TE^-1

其中：Z_k为对应第k组符合计数Y_k的中间状态矩阵，D为系统矩阵，I为单位矩阵，E为测量噪声矩阵。

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

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

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

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

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

$\frac{\mathrm{dZ}\left(t\right)}{\mathrm{dt}}=\mathrm{AZ}\left(t\right)+Z\left(t\right)A^T+\frac{Z\left(t\right)\mathrm{QZ}\left(t\right)}{{\mathrm{\γ}}^2}+V$

其中：Z(t)为t时刻的中间状态矩阵，A为状态转移矩阵，γ为预设的噪声干扰上限值，V为过程噪声矩阵，Q为给定的权重矩阵，t为时间。

所述的步骤（5）中，根据以下算式动态重建两种示踪剂对应的PET浓度分布数据：

$X_k^1={\mathrm{\Λ}}_1(X_k-X_{k-1})$

$X_k^2={\mathrm{\Λ}}_2(X_k-X_{k-1})$

Λ₁＝[WWOO]

Λ₂＝[OOWW]

其中：

和分别为对应第k组符合计数Y_k的采样时刻两种示踪剂的PET浓度分布数据，W为m×m维的单位矩阵，O为m×m维的零矩阵。

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

所述的浓度分布向量C为由两种示踪剂在血液中对应的浓度分布数据组成的2维向量；PET浓度分布向量X为4m维向量，符合计数Y为p维向量，测量噪声向量e为p维向量，代表符合事件采集过程中由随机事件、散射事件等引入的测量误差；过程噪声向量v为4m维向量，表征了状态之间转移的统计不确定性；系统矩阵D为p×4m维矩阵，表征了发射光子被探测器接收的概率，其受探测器结构、探测效率、衰减、死时间等因素的影响；过程噪声矩阵V为过程噪声向量v的协方差矩阵，且为4m×4m维矩阵；测量噪声矩阵E为测量噪声向量e的协方差矩阵，且为p×p维矩阵；权重矩阵Q为4m×4m维矩阵，其表示估计过程中的误差因素。

所述的状态转移矩阵A为m×m维的对角分块矩阵，其每个对角线元素均为a：

$a=\left[\begin{array}{cccc}-(k_2^1+k_3^1)& k_4^1& 0& 0\\ k_3^1& -k_4^1& 0& 0\\ 0& 0& -(k_2^2+k_3^2)& k_4^2\\ 0& 0& k_3^2& -k_4^2\end{array}\right]$

所述的状态转移矩阵B为m×1维的分块矩阵，其每个元素均为b：

$b=\left[\begin{array}{cc}{k}_1^1& 0\\ 0& 0\\ 0& k_1^2\\ 0& 0\end{array}\right]$

其中：

均为PET平行房室模型中的生理参数。

由于光子的能量单一的物理特性，双示踪剂同时注入的实验中的单示踪剂成分一般无法直接分离；本发明通过求解独立平行房室模型引导的PET状态空间模型，实现了对同时刻注入的双示踪剂，在两种示踪剂信号成分高度混合的情况下进行同时重建。

附图说明

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

具体实施方式

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

如图1所示，一种基于H无穷滤波的双示踪剂PET浓度的动态重建方法，包括如下步骤：

（1）采集符合计数。

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

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

（2）通过血输入函数求取双示踪剂血液浓度分布数据。

通过对生物组织进行动脉血采样，拟合出两种示踪剂对应的血输入函数，并通过血输入函数计算出双示踪剂在血液中与符合计数对应的n组浓度分布向量C₁～C_n；浓度分布向量C为由两种示踪剂在血液中对应的浓度分布数据组成的2维向量。

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

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

当生物体内同时注入多种示踪化合物，简单的房室模型无法直接模拟每种药物的动力学代谢过程。由于房室模型具有空间独立性，我们把简单房室模型扩展为平行房室模型（parallel compartment model）来分析独立又同时进行的代谢特性。平行房室模型假设每种示踪化合物具有各自的房室空间和血输入函数，其生理代谢过程是相互独立，互不影响的。

对于双示踪剂平行两组织房室模型，其一阶微分方程写成矩阵的形式：

$\left[\begin{array}{c}\stackrel{\·}{C_F^1}\left(t\right)\\ \stackrel{\·}{C_B^1}\left(t\right)\\ \stackrel{\·}{C_F^2}\left(t\right)\\ \stackrel{\·}{C_B^2}\left(t\right)\end{array}\right]=\left[\begin{array}{cccc}-(k_2^1+k_3^1)& k_4^1& 0& 0\\ k_3^1& -k_4^1& 0& 0\\ 0& 0& -(k_2^2+k_3^2)& k_4^2\\ 0& 0& k_3^2& -k_4^2\end{array}\right]\left[\begin{array}{c}{C}_F^1\left(t\right)\\ C_B^1\left(t\right)\\ C_F^2\left(t\right)\\ C_B^2\left(t\right)\end{array}\right]+\left[\begin{array}{cc}{k}_1^1& 0\\ 0& 0\\ 0& k_1^2\\ 0& 0\end{array}\right]\left[\begin{array}{c}{C}_p^1\left(t\right)\\ C_p^2\left(t\right)\end{array}\right]$

其中，上标代表第一种或第二种示踪剂的成份，C_p(t)、C_F(t)和C_B(t)分别表示踪化合物在血浆中的浓度、自由或未配合示踪化合物在组织中的浓度和配合示踪化合物在组织中的浓度，

表示两组织房室与血液交换的一阶速率常数即模型中的生理参数。

将上式写成积分形式，则有：

$\stackrel{\·}{x}\left(t\right)=\mathrm{ax}\left(t\right)+b{\tilde{C}}_p\left(t\right)$

其中： $x\left(t\right)={[{\∫}_0^t{C}_F^1\left(\mathrm{\τ}\right)\mathrm{d\τ},{\∫}_0^t{C}_B^1\left(\mathrm{\τ}\right)\mathrm{d\τ},{\∫}_0^t{C}_F^2\left(\mathrm{\τ}\right)\mathrm{d\τ},{\∫}_0^t{C}_B^2\left(\mathrm{\τ}\right)\mathrm{d\τ}]}^T$

${\tilde{C}}_p\left(t\right)={[{\∫}_0^t{C}_p^1\left(\mathrm{\τ}\right)\mathrm{d\τ},{\∫}_0^t{C}_p^2\left(\mathrm{\τ}\right)\mathrm{d\τ}]}^T$

$a=\left[\begin{array}{cccc}-(k_2^1+k_3^1)& k_4^1& 0& 0\\ k_3^1& -k_4^1& 0& 0\\ 0& 0& -(k_2^2+k_3^2)& k_4^2\\ 0& 0& k_3^2& -k_4^2\end{array}\right]$ $b=\left[\begin{array}{cc}{k}_1^1& 0\\ 0& 0\\ 0& k_1^2\\ 0& 0\end{array}\right]$

故PET状态空间方程可写成如下形式：

$\frac{\mathrm{dX}\left(t\right)}{\mathrm{dt}}=\mathrm{AX}\left(t\right)+\mathrm{BC}\left(t\right)+v$

Y(t)＝DX(t)+e

其中：X(t)为t时刻的PET浓度分布向量，A和B均为状态转移矩阵，D为系统矩阵，C(t)为t时刻双示踪剂在血液中的浓度分布向量，Y(t)为t时刻PET的符合计数，t为时间，v和e分别为过程噪声向量和测量噪声向量。PET浓度分布向量X为4m维向量，状态转移矩阵A为m×m维的对角分块矩阵，其每个对角线元素均为a；状态转移矩阵B为m×1维的分块矩阵，其每个元素均为b；测量噪声向量e为p维向量，代表符合事件采集过程中由随机事件、散射事件等引入的测量误差；过程噪声向量v为4m维向量，表征了状态之间转移的统计不确定性；系统矩阵D为p×4m维矩阵，表征了发射光子被探测器接收的概率，其受探测器结构、探测效率、衰减、死时间等因素的影响。

（4）利用H无穷滤波算法估算出PET浓度分布。

利用H无穷滤波算法求解上述PET状态空间方程得到以下迭代方程；根据符合计数Y_k通过以下迭代方程估计出对应的PET浓度分布向量X_k；

X_k＝AX_k-1+H_k(Y_k-DX_k-1)+BC_k

H_k＝Z_k(I+D^TE^-1DZ_k)^-1DT^E-1

其中：X_k和X_k-1分别为第k次和第k-1次迭代后的PET浓度分布向量；H_k为第k次迭代后的增益矩阵，Z_k为对应Y_k的中间状态矩阵，C_k为对应Y_k的浓度分布向量，Y_k为第k组符合计数，E为测量噪声矩阵，I为单位矩阵，k为自然数且1≤k≤n；测量噪声矩阵E为测量噪声向量e的协方差矩阵，且为p×p维矩阵。

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

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

$\frac{\mathrm{dZ}\left(t\right)}{\mathrm{dt}}=\mathrm{AZ}\left(t\right)+Z\left(t\right)A^T+\frac{Z\left(t\right)\mathrm{QZ}\left(t\right)}{{\mathrm{\γ}}^2}+V$

其中：Z(t)为t时刻的中间状态矩阵，γ为预设的噪声干扰上限值，V为过程噪声矩阵，Q为给定的权重矩阵；过程噪声矩阵V为过程噪声向量v的协方差矩阵，且为4m×4m维矩阵；权重矩阵Q为4m×4m维矩阵，其表示估计过程中的误差因素；本实施方式中，γ=100。

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

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

（5）重建两种示踪剂对应的PET浓度分布数据。

根据PET浓度分布向量X_k，根据以下算式动态重建出两种示踪剂对应的PET浓度分布数据：

$X_k^1={\mathrm{\Λ}}_1(X_k-X_{k-1})$

$X_k^2={\mathrm{\Λ}}_2(X_k-X_{k-1})$

Λ₁＝[WWOO]

Λ₂＝[OOWW]

其中：

和

分别为对应第k组符合计数Y_k的采样时刻两种示踪剂的PET浓度分布数据，W为m×m维的单位矩阵，O为m×m维的零矩阵。

以下我们采用Zubal胸腔模拟实验来评估双示踪剂同时重建的精度与鲁棒性。为了减少计算量，原始图像分辨率降低为32×32像素，包含3组感兴趣区域和一组背景区域。

双示踪剂的动力模型可以分别测量两类生理状态：¹⁸F-FDG成像葡萄糖代谢水平以及¹¹C-acetate成像肿瘤增生状况。其中¹⁸F和¹¹C的半衰期分别为110.2分钟和20.4分钟；FDG和acetate动态分布都采用两房室模型分析，速率常数真值如表1所示，两种示踪剂的输入采用经验数学模型：

$C_p^{\mathrm{FDG}}\left(t\right)=(A_{1}t-A_2-A_3)e^{-{\mathrm{\λ}}_{1}t}+A_2{e}^{-{\mathrm{\λ}}_{2}t}+A_3{e}^{-{\mathrm{\λ}}_{3}t}$

$C_P^{\mathrm{acetate}}\left(t\right)=[1-0.88(1-e^{-\left(\frac{21n2}{15}\right)t})]C_p^{\mathrm{FDG}}\left(\text{t}\right)$

表1

我们模拟180度旋转角下40簇投影，每簇含48条平行投影线的PET系统，系统矩阵由Fessler教授开发的软件包线性近似计算得到。由两种示踪剂同时注入后符合信号混合模拟得到正弦数据的动态序列，包括4个0.5分钟正弦图，4个2分钟正弦图和10个5分钟正弦图。由于11C的半衰期仅有20分钟，而我们只想分析在两示踪剂信号成分高度混合情况下的分离，因此仅采用前10个时间窗的数据进行计算。

实验分别模拟不同的光子计数水平的测量数据，并加入30％泊松噪声。为了定量分析误差水平，定义重建平均百分比误差为：

其中：M=10表示总时间窗数目，算子

表示对图像所有N个像素内取算术平均。

投影模拟得到在20分钟内总计数量为10⁴、10⁵、10⁶、10⁷和10⁸等计数水平（两示踪剂放射性计数水平相同，即每种示踪剂成分的计数量各占一半）的正弦图，使用本实施方式做重建，在完美房室模型情况下，每种示踪剂重建后的聚集浓度平均百分比误差随着计数量水平增加而增大。从状态空间的观点来分析，由于采用了完美的房室参数作为模型条件，状态方程预估值与真值几乎完全一致，因此理论上可以完全不采用测量数据的修正即可实现图像重建，而计数数据本身加入了噪声成分，高计数水平使数据修正的影响增大，重建图像的精度反而降低。原始图像与10⁶计数水平下重建的图像比较,其中¹⁸F-FDG重建百分比误差为0.360%，¹¹C-acetate重建百分比误差为0.037%。

Claims (7)

1.一种基于H无穷滤波的双示踪剂PET浓度的动态重建方法，包括如下步骤：

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

（2）通过对生物组织进行动脉血采样，拟合出两种示踪剂对应的血输入函数，并通过血输入函数计算出双示踪剂在血液中与符合计数对应的n组浓度分布向量C₁～C_n；

（3）根据PET平行房室模型建立基于双示踪剂的PET状态空间方程；

（4）利用H无穷滤波算法求解所述的PET状态空间方程得到以下迭代方程；根据所述的符合计数通过以下迭代方程估计出对应的PET浓度分布向量X_k；

X_k＝AX_k-1+H_k(Y_k-DX_k-1)+BC_k

其中：X_k和X_k-1分别为第k次和第k-1次迭代后的PET浓度分布向量；A和B均为状态转移矩阵，H_k为第k次迭代后的增益矩阵，D为系统矩阵，Y_k为第k组符合计数，C_k为对应第k组符合计数Y_k的浓度分布向量，k为自然数且1≤k≤n；

（5）根据PET浓度分布向量X_k，动态重建出两种示踪剂对应的PET浓度分布数据。

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

$\frac{\mathrm{dX}\left(t\right)}{\mathrm{dt}}=\mathrm{AX}\left(t\right)+\mathrm{BC}\left(t\right)+v$

Y(t)＝DX(t)+e

其中：X(t)为t时刻的PET浓度分布向量，A和B均为状态转移矩阵，D为系统矩阵，C(t)为t时刻双示踪剂在血液中的浓度分布向量，Y(t)为t时刻PET的符合计数，t为时间，v和e分别为过程噪声向量和测量噪声向量。

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

H_k＝Z_k(I+D^TE^-1DZ_k)^-1D^TE^-1

其中：Z_k为对应第k组符合计数Y_k的中间状态矩阵，D为系统矩阵，I为单位矩阵，E为测量噪声矩阵。

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

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

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

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

5.根据权利要求4所述的动态重建方法，其特征在于：所述的Riccati微分方程表达式如下：

$\frac{\mathrm{dZ}\left(t\right)}{\mathrm{dt}}=\mathrm{AZ}\left(t\right)+Z\left(t\right)A^T+\frac{Z\left(t\right)\mathrm{QZ}\left(t\right)}{{\mathrm{\γ}}^2}+V$

其中：Z(t)为t时刻的中间状态矩阵，A为状态转移矩阵，γ为预设的噪声干扰上限值，V为过程噪声矩阵，Q为给定的权重矩阵，t为时间。

6.根据权利要求1所述的动态重建方法，其特征在于：所述的步骤（5）中，根据以下算式动态重建两种示踪剂对应的PET浓度分布数据：

$X_k^1={\mathrm{\Λ}}_1(X_k-X_{k-1})$

$X_k^2={\mathrm{\Λ}}_2(X_k-X_{k-1})$

Λ₁＝[WWOO]

Λ₂＝[OOWW]

其中：

和

分别为对应第k组符合计数Y_k的采样时刻两种示踪剂的PET浓度分布数据，W为m×m维的单位矩阵，O为m×m维的零矩阵。

7.根据权利要求1或2所述的动态重建方法，其特征在于：所述的状态转移矩阵A为m×m维的对角分块矩阵，其每个对角线元素均为a：

$a=\left[\begin{array}{cccc}-(k_2^1+k_3^1)& k_4^1& 0& 0\\ k_3^1& -k_4^1& 0& 0\\ 0& 0& -(k_2^2+k_3^2)& k_4^2\\ 0& 0& k_3^2& -k_4^2\end{array}\right]$

所述的状态转移矩阵B为m×1维的分块矩阵，其每个元素均为b：

$b=\left[\begin{array}{cc}{k}_1^1& 0\\ 0& 0\\ 0& k_1^2\\ 0& 0\end{array}\right]$

其中：

均为PET平行房室模型中的生理参数。

Images