CN102938154B - 一种基于粒子滤波的动态pet图像重建方法 - Google Patents

Abstract

本发明公开了一种基于粒子滤波的动态PET图像重建方法，包括：(1)采集符合计数，建立状态空间方程；(2)对体素进行粒子采样，并求取体素的先验浓度值；(3)计算粒子权重值；(4)进行重采样计算；(5)迭代求取粒子浓度真值及其粒子权重真值；(6)求取体素值，遍历重建每帧PET图像。本发明方法通过利用状态空间方程与粒子滤波，将PET的数据统计特性和生理特性联合起来用于重建PET的动态图像，有效地加强了对噪声的过滤，提高了PET图像动态重建的准确率，能够将PET动态图像更好的还原出来，使重建结果在医学领域上更具有使用价值。

Description

一种基于粒子滤波的动态PET图像重建方法

技术领域

本发明属于PET成像技术领域，具体涉及一种基于粒子滤波的动态PET图像重建方法。

背景技术

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

近几年PET在实际医学领域的应用日趋广泛，但与此同时，临床上对PET成像的要求也随之变高，越来越多的医学领域需要PET能提供更高成像分辨率和能够实时的对患者进行扫描，而与此相应的图像维数的扩展和采集数据量的急剧加大对现有的重建算法提出了挑战，并且这些需求也对计算机的计算能力和存储空间有极为苛刻的要求。

目前，PET成像的重建方法大致可分为两类：解析法和迭代统计法。前一类主要是滤波反投影法(FBP)，计算速度快，但成像分辨率低并且存在伪影严重。为此，从数据统计角度进行重建的迭代统计法被提出来，由于迭代法基于统计学模型，对不完全数据适应性好，逐渐成为PET重建算法研究关注点；如ML-EM算法，其根据最大似然估计原理，将病人体素(Voxel)浓度作为待估计值、建立似然函数、求得极大值。但是，迭代统计发的计算复杂、计算量大，耗时长。为解决重建速度问题，大规模矩阵快速运算技术被不断运用到该算法体系中，如有序子集最大似然法。而针对于算法中未考虑先验估计的统计特性，提出了后验估计中加以修正的最大后验估计法。此外，将统计模型近似为高斯模型后又引入了最小二乘法，继而发展未带有惩罚加权的最小二乘、非负超松弛迭代法等。但是这些算法都只是将患者的结构信息加入重建，并没有效的对患者的生理信息进行利用。

而基于房室模型的状态空间体系则从另一个角度为PET图像的重建提供了新的思路，通过根据实际求解问题调整测量方程与状态方程的统一表达式，以实现静态、动态重建以及先验估计。通过相关方法求解，如卡尔曼滤波、H∞滤波等可适应不同噪声特性和清晰度的要求，与传统的解析法或迭代统计法相比具有优势。

但是以上这些算法在PET的实时动态重建方面都有着一些缺陷，有的无法适应动态PET成像庞大的数据采集量，从而导致重建准确率低或耗时过长，而有的则不能将患者的先验生理信息加载到重建中，使重建的信息利用率较低。

发明内容

针对现有技术所存在的上述技术缺陷，本发明提供了一种基于粒子滤波的动态PET图像重建方法，有效提高了PET图像动态重建的质量。

一种基于粒子滤波的动态PET图像重建方法，包括如下步骤：

(1)利用探测器对注入有放射性物质的生物组织进行探测，采集到不同时刻对应的多组符合计数，每一组符合计数对应一帧PET图像；根据PET成像原理，建立PET图像的状态空间方程：

x(t)＝e^AΔtx(t-1)×u

y(t)＝Dx(t)+v

其中：D为系统矩阵且为m×n维矩阵，A为状态矩阵；y(t)为第t帧PET图像对应时刻采集到的符合计数，且为m维向量；x(t)和x(t-1)分别为第t帧PET图像和第t-1帧PET图像的体素值，且为n维向量；u为放射性物质浓度变化对探测器的影响因子，v为测量噪声；Δt为第t帧PET图像对应时刻与第t-1帧PET图像对应时刻的间隔；n为PET图像的体素个数，m为符合计数中的元素个数；

(2)对PET图像中的体素进行粒子采样，获得每帧PET图像中每一体素的每一粒子浓度初值；根据采集到的符合计数，利用FBP算法求得PET图像中每一体素的先验浓度值；

(3)根据所述的状态空间方程建立以下方程表达式；根据所述的先验浓度值和粒子浓度初值，通过以下方程表达式求得当前帧PET图像中每一体素的每一粒子浓度初值对应的粒子权重值；

${\mathrm{\ω}}_j^k\left(t\right)=f_v\{z\left(t\right)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}[D(i,j)v_j^k\left(t\right)\left]\right\}$

$z\left(t\right)=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left[y_i\left(t\right)P_{j,i}\left(t\right)\right]$

$P_{j,i}\left(t\right)=\frac{D(i,j)\×\frac{u_j\left(t\right)}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{u}_j\left(t\right)}}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\{D(i,j)\×\frac{u_j\left(t\right)}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{u}_j\left(t\right)}\}}$

其中：为第t帧PET图像中第j体素的第k粒子浓度初值，为对应的粒子权重值，u_j(t)为第t帧PET图像中第j体素的先验浓度值，y_i(t)为第t帧PET图像对应时刻采集到的符合计数中的第i元素值，D(i，j)为系统矩阵中第i行第j列的元素值；f_v为噪声的概率名都函数；k、i和j均为自然数，且1≤i≤m，1≤j≤n，k为粒子序号；

(4)根据所述的粒子浓度初值及其粒子权重值，通过重采样算法求得重采样后的粒子浓度值及其粒子权重值；

(5)将重采样后的粒子浓度值作为步骤(3)中的粒子浓度初值，重复步骤(3)和步骤(4)，直至重采样后的粒子浓度值收敛，将收敛后的粒子浓度值及其粒子权重值作为粒子浓度真值及其粒子权重真值；

(6)根据粒子浓度真值及其粒子权重真值，通过以下方程表达式求得当前帧PET图像中每一体素的体素值，进而遍历每帧PET图像。

$x_j\left(t\right)=\underset{k=1}{\overset{s}{\mathrm{\Σ}}}\{{\stackrel{\‾}{v}}_j^k\left(t\right)\×{\stackrel{\‾}{\mathrm{\ω}}}_j^k\left(t\right)\}$

其中：x_j(t)为第t帧PET图像中第j体素的体素值，为第t帧PET图像中第j体素的第k粒子浓度真值，为对应的粒子权重真值，s为粒子采样个数。

所述的步骤(2)中，对PET图像中的体素进行粒子采样的方法为：获取体素值的分布区间，所述的分布区间符合非线性泊松分布；对于任一体素通过对所述的分布区间进行采样获得s个粒子采样值，进而根据公式计算出每帧PET图像中每一体素的每一粒子浓度初值；其中，为第t帧PET图像中第j体素的第k粒子浓度初值，为第t帧PET图像中第j体素的第k粒子采样值，为对应的粒子初始权重，k为自然数且1≤k≤s。

优选地：

对于第一帧PET图像，所述的粒子初始权重为初始给定值；

对于非第一帧PET图像，所述的粒子初始权重为前一帧PET图像中第j体素的第k粒子浓度真值对应的粒子权重真值。

这样可以有效地的提升重建的准确度，避免了固定初始权重在重建不同帧，尤其是重建数据扫描时间长的帧时所产生的误差，并且加强了帧与帧之间的联系。

优选地：所述的步骤(4)中，采用粗糙重采样算法求取重采样后的粒子浓度值及其粒子权重值。粗糙重采样算法对于动态PET重建尤其是长时间扫描的动态PET重建的准确度的提高有很大的帮助，它能有效地避免重采样过程中由于权重值的过于集中而造成的采样数值贫困化的这一缺点，而采样数值贫困化是动态重建过程中降低重建结果的一个重要原因。

本发明的PET图像的动态重建方法通过利用状态空间方程与粒子滤波，将PET的数据统计特性和生理特性联合起来用于重建PET的动态图像，有效地加强了对噪声的过滤，提高了PET图像动态重建的准确率，将PET动态图像更好的还原出来，并且通过有效地运用基于房室模型的状态空间法成功的将患者的先验生理信息加入了重建当中，使重建结果在医学领域上更具有使用价值，重建结果比ML-EM、FBP等传统的重建方法的效果更佳，更具有医学的使用价值。

附图说明

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

图2为粒子采样值的泊松分布区间示意图。

图3为Zubal-thorax-phantom数字体模的模型示意图。

图4(a)～(c)分别为本发明方法对于Zubal-thorax-phantom数字体模的第1、5、10帧的PET重建图像。

图4(d)～(f)分别为传统ML-EM方法对于Zubal-thorax-phantom数字体模的第1、5、10帧的PET重建图像。

具体实施方式

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

如图1所示，一种基于粒子滤波的动态PET图像重建方法，包括如下步骤：

(1)采集符合计数，建立状态空间方程。

利用探测器(正电子发射断层扫描仪)对注入有放射性物质的人体组织进行探测，采集到不同时刻对应的多组符合计数，每一组符合计数对应一帧PET图像；扫描仪探测人体内发出的放射性信号，经过符合和采集系统处理，形成原始投影线，并以正弦图(sinogram)的方式存放于计算机中；通过对原始采集到的sinogram数据进行各类校正后，得到符合计数。

根据PET成像原理，建立PET图像的状态空间方程如下：

x(t)＝e^AΔtx(t-1)×u

y(t)＝Dx(t)+v

其中：D为系统矩阵且为m×n维矩阵，A为状态矩阵；y(t)为第t帧PET图像对应时刻采集到的符合计数，且为m维向量；x(t)和x(t-1)分别为第t帧PET图像和第t-1帧PET图像的体素值，且为n维向量；u为放射性物质浓度变化对探测器的影响因子，v为测量噪声；Δt为第t帧PET图像对应时刻与第t-1帧PET图像对应时刻的间隔；n为PET图像的体素个数，m为符合计数中的元素个数。

(2)对体素进行粒子采样，并求取体素的先验浓度值。

对PET图像中的体素进行粒子采样：获取体素值的分布区间(0～5.547)，分布区间符合非线性泊松分布(如图2所示)；对于任一体素通过对分布区间进行采样获得s个粒子采样值，进而根据公式计算出每帧PET图像中每一体素的每一粒子浓度初值；其中，为第t帧PET图像中第j体素的第k粒子浓度初值，为第t帧PET图像中第j体素的第k粒子采样值，为对应的粒子初始权重，k为自然数且1≤k≤s；本实施方式中，粒子采样个数s＝500。

对于第一帧PET图像，粒子初始权重为初始给定值(本实施方式取1)；

对于非第一帧PET图像，粒子初始权重为前一帧PET图像中第j体素的第k粒子浓度真值对应的粒子权重真值。

根据采集到的符合计数，利用FBP算法求得PET图像中每一体素的先验浓度值。

(3)计算粒子权重值。

根据状态空间方程建立以下方程表达式；根据先验浓度值和粒子浓度初值，通过以下方程表达式求得当前帧PET图像中每一体素的每一粒子浓度初值对应的粒子权重值；

${\mathrm{\ω}}_j^k\left(t\right)=f_v\{z\left(t\right)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}[D(i,j)v_j^k\left(t\right)\left]\right\}$

$z\left(t\right)=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left[y_i\left(t\right)P_{j,i}\left(t\right)\right]$

$P_{j,i}\left(t\right)=\frac{D(i,j)\×\frac{u_j\left(t\right)}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{u}_j\left(t\right)}}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\{D(i,j)\×\frac{u_j\left(t\right)}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{u}_j\left(t\right)}\}}$

其中：为第t帧PET图像中第j体素的第k粒子浓度初值，为对应的粒子权重值，u_j(t)为第t帧PET图像中第j体素的先验浓度值，y_i(t)为第t帧PET图像对应时刻采集到的符合计数中的第i元素值，D(i，j)为系统矩阵中第i行第j列的元素值；f_v为噪声的概率名都函数；k、i和j均为自然数，且1≤i≤m，1≤j≤n，k为粒子序号。

(4)进行重采样计算。

根据粒子浓度初值及其粒子权重值，通过粗糙重采样算法求得重采样后的粒子浓度值及其粒子权重值；

重采样时是修改较小的权重保留大的权重，不过这里需要注意，当相邻两帧的采样时间差距较大时，直接按权重阈值赋予所有较小的权重，可能会导致采样权重过度集中，本实施方式选择修改替换部分较小的权重(取半数)；同时，为了防止重采样集中于大权重，导致采样种类的贫困化在各粒子中加入一定的随机噪声：

$g_j^{k+}\left(t\right)=g_j^{k-}\left(t\right)+\mathrm{\Δg}$

Δg～(0，KMS^-1/n)

其中，和分别表示重采样前后的粒子采样值，Δg是期望为0的随机变量；K代表每个粒子噪声的程度；S为粒子个数；n为状态空间维数；M为分布区间中的最大值。

(5)迭代求取粒子浓度真值及其粒子权重真值。

将重采样后的粒子浓度值作为步骤(3)中的粒子浓度初值，重复步骤(3)和步骤(4)，直至重采样后的粒子浓度值收敛，将收敛后的粒子浓度值及其粒子权重值作为粒子浓度真值及其粒子权重真值。收敛条件如下：

$\left|\frac{v_j^{k+}\left(t\right)-v_j^{k-}\left(t\right)}{v_j^{k-}\left(t\right)}\right|\≤\min (\frac{1}{s},0.002)$

其中，和分别表示前后两次迭代的粒子浓度值。

(6)求取体素值，遍历重建每帧PET图像。

根据粒子浓度真值及其粒子权重真值，通过以下方程表达式求得当前帧PET图像中每一体素的体素值，进而遍历每帧PET图像。

$x_j\left(t\right)=\underset{k=1}{\overset{s}{\mathrm{\Σ}}}\{{\stackrel{\‾}{v}}_j^k\left(t\right)\×{\stackrel{\‾}{\mathrm{\ω}}}_j^k\left(t\right)\}$

其中：x_j(t)为第t帧PET图像中第j体素的体素值，为第t帧PET图像中第j体素的第k粒子浓度真值，为对应的粒子权重真值，s为粒子采样个数。

以下我们通过实验来验证本实施方式的实用性和可靠性，图3为著名的Zubal-thorax-phantom数字体模模型，该模型包含3个组织区域(ROIa～ROIc)和一个背景区域(Background)。

我们采用的是180度内48个采样角，每个采样角下射束为32条。总扫描时间为60分钟，这60分钟被分成了18个时间段，它们被分为4个30秒、4个1分钟和4个5分钟。其中，血液中放射性浓度C_P(t)可以表示为：

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

其中：A₁＝851.1225μC_i/mL/min，A₂＝20.8113μC_i/mL，A₃＝21.8798μC_i/mL，α＝-4.133859min^-1，β＝-0.01043449min^-1，γ＝-0.1190996min^-1；所使用的放射性示踪剂为FDG，其放射性衰减常数值为0.0063min^-1。重建图像的大小为32×32，原始图像中不同区域的动力学参数如表1所示：

表1

Region	k1(min-1)	k2(min-1)	k3(min-1)	k4(min-1)
					ROIa	0.55951	2.75288	0.44793	0.01101
ROIb	0.37811	1.04746	0.13483	0.00857
					ROIc	0.78364	1.15641	0.11200	0.02706

为了能够验证本实施方式能够很好的提升动态PET成像的准确度，我们将粒子滤波的动态重建结果和用ML-EM的重建结果做了比较，具体情况如图4；其中，图4(a)～(c)分别代表了本实施方式第1、5、10帧的重建图像，而图4(d)～(f)分别表示的是ML-EM相应的第1、5、10帧的重建图像。

同时为了能够对实验结果做出定量的比较粒子滤波重建结果和ML-EM重建结果的偏差(bias)和方差(variance)，具体的公式如下：

$\mathrm{bias}=\frac{1}{n}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(X_j-\stackrel{\‾}{X})$

$\mathrm{variance}=\frac{1}{n-1}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(X_j-\stackrel{\‾}{X})}^2$

其中：X_j代表的是粒子滤波的重建结果和ML-EM的重建结果；而则代表的是实际真值。

本实施方式重建结果和ML-EM重建结果的量化比较结果如表2所示：

表2

通过以上的实验结果我们可以看出，粒子滤波的重建结果无论是在实际的图像还是偏差和方差这两个方面的表现都要优于ML-EM的重建结果，由此可以看出本实施方式有效地提高了PET动态重建的准确率。

Claims (1)

1.一种基于粒子滤波的动态PET图像重建方法，包括如下步骤：

(1)利用探测器对注入有放射性物质的生物组织进行探测，采集到不同时刻对应的多组符合计数，每一组符合计数对应一帧PET图像；根据PET成像原理，建立PET图像的状态空间方程如下：

x(t)＝e^AΔtx(t-1)×u

y(t)＝Dx(t)+v

其中：D为系统矩阵且为m×n维矩阵，A为状态矩阵；y(t)为第t帧PET图像对应时刻采集到的符合计数，且为m维向量；x(t)和x(t-1)分别为第t帧PET图像和第t-1帧PET图像的体素值，且为n维向量；u为放射性物质浓度变化对探测器的影响因子，v为测量噪声；Δt为第t帧PET图像对应时刻与第t-1帧PET图像对应时刻的间隔；n为PET图像的体素个数，m为符合计数中的元素个数；

(2)对PET图像中的体素进行粒子采样，获得每帧PET图像中每一体素的每一粒子浓度初值，具体过程如下：

获取体素值的分布区间，所述的分布区间符合非线性泊松分布；对于任一体素通过对所述的分布区间进行采样获得s个粒子采样值，进而根据公式计算出每帧PET图像中每一体素的每一粒子浓度初值；其中，为第t帧PET图像中第j体素的第k粒子浓度初值，为第t帧PET图像中第j体素的第k粒子采样值，为对应的粒子初始权重，k为自然数且1≤k≤s；对于第一帧PET图像，所述的粒子初始权重为初始给定值；对于非第一帧PET图像，所述的粒子初始权重为前一帧PET图像中第j体素的第k粒子浓度真值对应的粒子权重真值；

根据采集到的符合计数，利用FBP算法求得PET图像中每一体素的先验浓度值；

(3)根据所述的状态空间方程建立以下方程表达式；根据所述的先验浓度值和粒子浓度初值，通过以下方程表达式求得当前帧PET图像中每一体素的每一粒子浓度初值对应的粒子权重值；

${\mathrm{\ω}}_j^k\left(t\right)=f_v\{z\left(t\right)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}[D(i,j)v_j^k\left(t\right)\left]\right\}$

$z\left(t\right)=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left[y_i\left(t\right)P_{j,i}\left(t\right)\right]$

$P_{j,i}\left(t\right)=\frac{D(i,j)\×\frac{u_j\left(t\right)}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{u}_j\left(t\right)}}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\{D(i,j)\×\frac{u_j\left(t\right)}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{u}_j\left(t\right)}\}}$

其中：为第t帧PET图像中第j体素的第k粒子浓度初值，为对应的粒子权重值，u_j(t)为第t帧PET图像中第j体素的先验浓度值，y_i(t)为第t帧PET图像对应时刻采集到的符合计数中的第i元素值，D(i,j)为系统矩阵中第i行第j列的元素值；f_v为噪声的概率密度函数；k、i和j均为自然数，且1≤i≤m，1≤j≤n，k为粒子序号，n为PET图像的体素个数，m为符合计数中的元素个数；

(4)根据所述的粒子浓度初值及其粒子权重值，通过粗糙重采样算法求得重采样后的粒子浓度值及其粒子权重值；

(5)将重采样后的粒子浓度值作为步骤(3)中的粒子浓度初值，重复步骤(3)和步骤(4)，直至重采样后的粒子浓度值收敛，将收敛后的粒子浓度值及其粒子权重值作为粒子浓度真值及其粒子权重真值；

(6)根据粒子浓度真值及其粒子权重真值，通过以下公式求得当前帧PET图像中每一体素的体素值，进而遍历每帧PET图像；

$x_j\left(t\right)=\underset{k=1}{\overset{s}{\mathrm{\Σ}}}\{{\stackrel{\‾}{v}}_j^k\left(t\right)\×{\stackrel{\‾}{\mathrm{\ω}}}_j^k\left(t\right)\}$

其中：x_j(t)为第t帧PET图像中第j体素的体素值，为第t帧PET图像中第j体素的第k粒子浓度真值，为对应的粒子权重真值，s为粒子采样个数。