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

Abstract

本发明公开了一种基于粒子滤波的静态PET图像重建方法，包括(1)建立状态空间方程；(2)对体素进行粒子采样；(3)求取粒子权重值；(4)对粒子进行重采样；(5)求取粒子浓度真值和粒子权重真值；(6)计算每一体素的待估计浓度值。本发明通过利用状态空间结合粒子滤波，将PET的数据统计特性和生理特性联合起来重建PET图像，提高了图像的分辨率和锐度，更好地还原了真实的PET图像，同时将PET中噪声的数据模型定为泊松分布而非高斯分布，更符合PET扫描的真实情况，从而使得重建过程中对噪声的过滤和优化更加有效，得到的重建结果比ML-EM、FBP等传统重建方法更接近PET真实情况，效果更佳。

Description

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

技术领域

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

背景技术

正电子断层成像PET(Positron Emission Tomography)，是一种核医学成像技术。与CT成像不同，PET为功能成像，能得到病人体内新陈代谢的情况，从而可以比CT、MRI更早地探测出病灶。

比如葡萄糖是人体必需的物质，能够参与人体的新陈代谢，葡萄糖在体内的分布从一个侧面反映了人体新陈代谢的情况，所以，通过追踪葡萄糖在病人体内的浓度分布，就可以得知体内哪些部位新陈代谢不正常，从而得知那些部位发生了病变。普通的葡萄糖并不能被体外的任何探测器探测到，但是，通过在葡萄糖中加入同位素，能达到这个目的。同位素在随着葡萄糖扩散的过程中，不断衰变而产生正电子(e⁺)，正电子极不稳定，它与周围的负电子(e^-)湮灭，发出一对光子，这对光子能够穿透人体组织发射到体外，从而就可以被探测器探测到。通过探测这一对光子，就可以得到一条符合线(LOR)。如果将探测器排列成环状或者其他阵列，就可以接收来自体内的大部分光子，从而得到一系列符合线，这些符合线和体内正电子湮灭的位置有一定的关系，通过这些符合线，就可以重建出病人体内同位素标记的葡萄糖的浓度分布图，即PET图像。

由正负电子湮灭产生的一对光子在飞行途中，并不一定是畅通无阻的，途中可能经历散射、衰减，受这些因素影响的光子对将不再是沿着互成180°的直线前进，当一对探测器探测到这对光子的时候，其产生的符合线未必经过正负电子湮灭的地方，即出现了定位错误，在PET成像中称为散射符合；两个毫不相关的光子如果同时被探测器探测到，那探测系统将认为这是由一对正负电子湮灭产生的，将它们的连线记录为符合线，这样的记录也是假事件，在PET成像中成为随机符合。由于以上这些原因以及病人体内的干扰，使得PET采集到的数据包含了许多噪声，这些噪声直接影响重建后的PET图像质量。所以，重建算法的鲁棒性(Robust)决定了PET图像的质量，这也是为什么重建算法成为现今PET领域的研究热点。

目前PET重建的方法大致可以分为解析法、迭代法和状态空间法三种：

解析法是以滤波反投影算法(FBP)为代表的一种方法，该算法是基于傅立叶变换来实现由投影数据得到原始数据。但是由于是解析算法，并不能很好的处理PET数据中的噪声，导致滤波反投影算法得到的图像伪影很重，质量很差。

迭代法中以ML-EM和MAP(maximum a posteriori)算法较为典型。ML-EM算法根据最大似然估计原理，将病人体素(Voxel)浓度作为待估计值、建立似然函数、求得极大值。MAP算法根据原始数据和PET的统计特性建立后验函数再求其最大值，从而得到体素浓度分布。但迭代算法对模型的依赖过大，所以把基于泊松分布的算法运用在非严格按照泊松分布的数据上，得到的结果势必不理想；除了统计特性，决定PET图像的还有生理先验、结构先验等信息，而这些信息在统计迭代算法中并不能有效地被利用，这也是ML-EM和MAP的巨大缺陷。

状态空间法能够将PET成像的过程模型化，使其物理过程得到数学的表达，通过状态空间，能很好地把PET扫描中的统计特性、生理特性和结构特性结合起来，从而达到改进重建质量的目的。目前现有的状态空间求解方法为卡尔曼滤波，该方法假设数据分布符合高斯分布，这与PET数据的统计特性不符，所以重建得到的PET图像效果不是很理想，这也是目前困扰着状态空间方法的问题之一。

发明内容

本发明提供了一种基于粒子滤波的静态PET图像重建方法，利用粒子滤波很好地将PET的物理过程和统计信息结合起来，解决了现有PET重建方法所存在的上述技术缺陷，提高了PET图像重建的质量。

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

(1)根据PET成像原理，建立PET图像的状态空间方程：

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

式1中：D为系统矩阵，且为m×n维矩阵；y为采集到的符合计数，且为m维向量；x为PET图像浓度分布的待估计值，且为n维向量；e为测量噪声；x_t+1为第(t+1)帧PET图像浓度分布的待估计值，x_t为第t帧PET图像浓度分布的待估计值，A为状态矩阵；

(2)根据待估计浓度值的概率分布区间，采用线性平均采样法对PET图像中的体素进行粒子采样，获得PET图像中体素的每一粒子浓度值；

所述的线性平均采样法是通过对待估计浓度值的概率分布区间进行(s-1)等分，得到s个节点，对应的每个节点值即为PET图像中任一体素的每一粒子浓度值，s为采样个数；

(3)根据采集到的符合计数，利用FBP法求得PET图像中每一体素的先验浓度值；根据所述的状态空间方程建立以下方程表达式；根据所述的先验浓度值和粒子浓度值，通过以下方程表达式求得PET图像中每一体素的每一粒子浓度值对应的粒子权重值；

$P(j,i)=\frac{D(i,j)\×\frac{x\left[j\right]}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}x\left[j\right]}}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\{D(i,j)\×\frac{x\left[j\right]}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}x\left[j\right]}\}}---\left(2\right)$

$w_{\left(k\right)}\left[j\right]=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left\{y\right[i]\×P(j,i)\}-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\{D(i,j)\×x_{\left(k\right)}[j\left]\right\}---\left(3\right)$

式2和式3中：P(j，i)为符合计数分布矩阵，且为n×m维矩阵；x[j]为PET图像中第j体素的先验浓度值，y[i]为第i符合计数，j为体素序号，i为符合计数序号；D(i，j)为系统矩阵中第i行第j列的元素值；w_(k)[j]为PET图像中第j体素的第k粒子权重值，x_(k)[j]为PET图像中第j体素的第k粒子浓度值，k为粒子序号；

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

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

(6)根据粒子浓度真值和对应的粒子权重真值，通过以下方程表达式求得PET图像中每一体素的待估计浓度值；

$x^{*}\left[j\right]=\underset{k=1}{\overset{s}{\mathrm{\Σ}}}\left\{{x_{\left(k\right)}}^{*}\right[j]\×{w_{\left(k\right)}}^{*}[j\left]\right\}---\left(4\right)$

式4中：x_(k) ^*[j]为PET图像中第j体素的第k粒子浓度真值，w_(k) ^*[j]为PET图像中第j体素的第k粒子权重真值，x^*[j]为PET图像中第j体素的待估计浓度值，s为采样个数且为实际经验值，通常取值为500。

本发明的PET图像重建方法通过利用状态空间结合粒子滤波，将PET的数据统计特性和生理特性联合起来重建PET图像，提高了图像的分辨率和锐度，更好地还原了真实的PET图像，同时将PET中噪声的数据模型定为泊松分布而非高斯分布，更符合PET扫描的真实情况，从而使得重建过程中对噪声的过滤和优化更加有效，得到的重建结果比ML-EM、FBP等传统重建方法更接近PET真实情况，效果更佳。

附图说明

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

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

图3(a)为关于Zubal-thorax-phantom数字体模采用本发明方法重建后的PET图像。

图3(b)为关于Zubal-thorax-phantom数字体模采用ML-EM方法重建后的PET图像。

图4(a)为内嵌有六个半径不等的圆柱形容器的实验器具的实物图像。

图4(b)为内嵌有六个半径不等的圆柱形容器的实验器具的截面示意图。

图5(a)为关于内嵌有六个半径不等的圆柱形容器的实验器具采用本发明方法重建后的PET图像。

图5(b)为关于内嵌有六个半径不等圆柱形容器的实验器具采用ML-EM方法重建后的PET图像。

具体实施方式

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

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

(1)根据PET成像原理，建立PET图像的状态空间方程：

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

式1中：D为系统矩阵，且为m×n维矩阵；y为采集到的符合计数，且为m维向量；x为PET图像浓度分布的待估计值，且为n维向量；e为测量噪声；x_t+1为第(t+1)帧PET图像浓度分布的待估计值，x_t为第t帧PET图像浓度分布的待估计值，A为状态矩阵。

(2)根据待估计浓度值的概率分布区间，采用线性平均采样法对PET图像中的体素进行粒子采样，获得PET图像中体素的每一粒子浓度值；

线性平均采样法是通过对待估计浓度值的概率分布区间进行(s-1)等分，得到s个节点，对应的每个节点值即为PET图像中任一体素的每一粒子浓度值，s为采样个数。

(3)根据采集到的符合计数，利用FBP法求得PET图像中每一体素的先验浓度值；根据状态空间方程建立以下方程表达式；根据先验浓度值和粒子浓度值，通过以下方程表达式求得PET图像中每一体素的每一粒子浓度值对应的粒子权重值；

$P(j,i)=\frac{D(i,j)\×\frac{x\left[j\right]}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}x\left[j\right]}}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\{D(i,j)\×\frac{x\left[j\right]}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}x\left[j\right]}\}}---\left(2\right)$

$w_{\left(k\right)}\left[j\right]=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left\{y\right[i]\×P(j,i)\}-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\{D(i,j)\×x_{\left(k\right)}[j\left]\right\}---\left(3\right)$

式2和式3中：P(j，i)为符合计数分布矩阵，且为n×m维矩阵；x[j]为PET图像中第j体素的先验浓度值，y[i]为第i符合计数，j为体素序号，i为符合计数序号；D(i，j)为系统矩阵中第i行第j列的元素值；w_(k)[j]为PET图像中第j体素的第k粒子权重值，x_(k)[j]为PET图像中第j体素的第k粒子浓度值，k为粒子序号。

(4)根据粒子浓度值和对应的粒子权重值，通过重采样算法求得重采样后的粒子浓度值和对应的粒子权重值；

重采样算法如下(用伪代码表示)：

a.构建一个累计分布函数(CDF)：c₁＝0， $c_i=c_{i-1}+w_k^i$

b.从一个归一化的随机分布中取一个数u₁： $U[0,N_s^{-1}]$

c.forj＝1：N_s

$u_j=u_1+N_s^{-1}(j-1)$

While u_j＞c_i

*i＝i+1

end while

对样本重新赋值： $x_k^{j*}=x_k^i$

end for

重采样过程将权重大的“粒子”保留下来，权重小的“粒子”删除，从而使整个“粒子”的分布更接近真值；经过若干次迭代后，就能得到精度范围内的结果。

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

(6)根据粒子浓度真值和对应的粒子权重真值，通过以下方程表达式求得PET图像中每一体素的待估计浓度值；

$x^{*}\left[j\right]=\underset{k=1}{\overset{s}{\mathrm{\Σ}}}\left\{{x_{\left(k\right)}}^{*}\right[j]\×{w_{\left(k\right)}}^{*}[j\left]\right\}---\left(4\right)$

式4中：x_(k) ^*[j]为PET图像中第j体素的第k粒子浓度真值，w_(k) ^*[j]为PET图像中第j体素的第k粒子权重真值，x^*[j]为PET图像中第j体素的待估计浓度值，s为采样个数且为实际经验值，本实施例取值为500。

以下通过两个实验来验证本实施方式的实用性和可靠性，第一个实验为数字体模的重建，第二个实验为真实扫描数据的重建。

图2为著名的Zubal-thorax-phantom数字体模的模型，采样角度为180°，共有128个角度和128个半径对，最后所成图像像素为128×128；针对该模型，采用本实施方式重建的PET图像与采用ML-EM方法重建的PET图像相比较(如图3)，从图3可以看出，用本实施方式重建出来的图像分辨率和图像清晰程度都比用ML-EM方法重建出来的要高得多。为了从统计学的角度分析本发明方法的优点，表1列出了在原始sinogram(符合计数向量)中分别加入10％、20％、30％、40％噪声的情况下，用本实施方式与用ML-EM方法重建得到的图像的偏差和标准差。

表1：本实施方式与ML-EM方法重建得到的图像的偏差和标准差

从表1可以看出，本实施方式在偏差和标准差两个统计方面都明显优于现有技术常用的ML-EM方法，即证明了基于粒子滤波的静态PET图像重建方法相较于目前现有的重建方法的可靠性。

在真实数据方面，通过设计了一个实验器具可验证本实施方式的实用性。本实验使用的是一个内嵌有六个半径不等的圆柱形容器的实验器具，圆柱形容器的半径分别为37mm、28mm、22mm、17mm、13mm、10mm(如图4)，每个容器里注入107.92Bq/ml的FDG，通过日本滨松公司生产的SHR-22000PET扫描仪扫描，得到的sinogram分别用本实施方式和ML-EM方法进行重建，重建结果如图5所示，从图5可以看出由本实施方式重建的图像更加清晰、均匀。

Claims (3)

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

(1)根据PET成像原理，建立PET图像的状态空间方程：

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

式1中：D为系统矩阵，且为m×n维矩阵；y为采集到的符合计数，且为m维向量；x为PET图像浓度分布的待估计值，且为n维向量；e为测量噪声；x_t+1为第(t+1)帧PET图像浓度分布的待估计值，x_t为第t帧PET图像浓度分布的待估计值，A为状态矩阵；

(2)根据待估计浓度值的概率分布区间，采用线性平均采样法对PET图像中的体素进行粒子采样，获得PET图像中体素的每一粒子浓度值；

(3)根据采集到的符合计数，利用FBP法求得PET图像中每一体素的先验浓度值；根据所述的状态空间方程建立以下方程表达式；根据所述的先验浓度值和粒子浓度值，通过以下方程表达式求得PET图像中每一体素的每一粒子浓度值对应的粒子权重值；

$P(j,i)=\frac{D(i,j)\×\frac{x\left[j\right]}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}x\left[j\right]}}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\{D(i,j)\×\frac{x\left[j\right]}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}x\left[j\right]}\}}---\left(2\right)$

$w_{\left(k\right)}\left[j\right]=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left\{y\right[i]\×P(j,i)\}-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\{D(i,j)\×x_{\left(k\right)}[j\left]\right\}---\left(3\right)$

式2和式3中：P(j，i)为符合计数分布矩阵，且为n×m维矩阵；x[j]为PET图像中第j体素的先验浓度值，y[i]为第i符合计数，j为体素序号，i为符合计数序号；D(i，j)为系统矩阵中第i行第j列的元素值；w_(k)[j]为PET图像中第j体素的第k粒子权重值，x_(k)[j]为PET图像中第j体素的第k粒子浓度值，k为粒子序号；

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

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

(6)根据粒子浓度真值和对应的粒子权重真值，通过以下方程表达式求得PET图像中每一体素的待估计浓度值；

$x^{*}\left[j\right]=\underset{k=1}{\overset{s}{\mathrm{\Σ}}}\left\{{x_{\left(k\right)}}^{*}\right[j]\×{w_{\left(k\right)}}^{*}[j\left]\right\}---\left(4\right)$

式4中：x_(k) ^*[j]为PET图像中第j体素的第k粒子浓度真值，w_(k) ^*[j]为PET图像中第j体素的第k粒子权重真值，x^*[j]为PET图像中第j体素的待估计浓度值，s为采样个数。

2.根据权利要求1所述的基于粒子滤波的静态PET图像重建方法，其特征在于：所述的线性平均采样法是通过对待估计浓度值的概率分布区间进行(s-1)等分，得到s个节点，对应的每个节点值即为PET图像中任一体素的每一粒子浓度值。

3.根据权利要求1所述的基于粒子滤波的静态PET图像重建方法，其特征在于：所述的采样个数为500。

Images