CN103559699A - 一种基于投影估计的多能谱ct图像重建方法 - Google Patents

Abstract

本发明公开一种基于投影估计的多能谱CT图像重建方法，用以重建被测物体的多种基材料密度图像。该方法包括：利用传统单能CT重建算法由所采集的多色投影数据分别直接重建被测物体的能量图像；计算每一个能量图像沿着所有其它能谱的射线方向的线积分，估计在这些方向上的当前能谱的多色投影，求得几何参数一致的多色投影；标定多基材料分解函数，将估计出的几何参数一致的多色投影，分解为多个基材料的线积分；分别由多基材料线积分重建出所对应的多基材料密度图像。本发明方法简单实用，适用于多色投影几何参数不一致的情况下的多能谱CT图像重建。与现有技术相比，仅利用所测投影数据估计出几何参数一致的多色投影，便可重建出高质量的图像。

Description

一种基于投影估计的多能谱CT图像重建方法

技术领域

本发明涉及X射线CT成像技术领域，具体而言，涉及一种基于投影估计的多能谱CT图像重建方法。

背景技术

X射线计算机层析成像技术（简称X射线CT）是一种利用X射线与物质的相互作用原理，对物体内部信息进行成像的一种技术。X射线CT技术不仅是一种常规的医学诊疗检测手段，也是一项重要的工业支撑技术，广泛应用于多种领域的关键部件的无损检测。目前常用的X射线管发不出理想的单能量X射线束，发出的是服从一定谱分布的多色X射线束，扫描采集的数据为多色投影数据，因此传统CT也称为单能谱CT。单能谱CT仅使用一个X射线能谱对被测物体进行扫描，所重建的不同物质的CT图像可能具有相同或相近的CT值，难以区分不同的物质。此外传统CT算法中忽略了X射线的多色性，导致在重建图像中产生硬化伪影。多能谱CT扫描被测物体时，使用多个X射线能谱，能够测得比传统单能谱CT更多的被测物体信息。利用这些信息可以重建出被测物体等效原子序数和电子密度图像，或是重建出特定基材料的密度图像，基于这些信息对物质进行区分。由于多能谱CT具有更好的物质区分能力，因此具有广泛的应用前景。

在高分辨率成像情况下，多能谱CT图像问题是一个高维数的非线性反问题。理想的多能谱CT的数据采集要求为“同源、同时、同向”。所谓同源，即保证每种能谱的X射线采集时的质量完全一致；所谓同时，即保证被测物体在每种能谱的X射线采集时的位移完全一致；所谓同向，即保证每种能谱的X射线采集路径完全一致。在满足理想条件情况下，可以在投影空间将多色投影数据直接分解为各个基材料的单色投影，然后用传统CT重建算法如滤波反投影算法等即可重建出高质量的被测物体的基材料图像。然而目前常用的多能谱CT硬件***均难以满足理想数据采集条件，尤其难以满足“同向”条件，导致所采集的各种能谱的多色投影的几何参数不一致。在这种情况下，多能谱成像方程具有很高的病态性。现有的求解方法大致可以分为三类：迭代法、投影域分解法和图像域分解法。迭代法利用数值方法或者优化方法构造迭代结构，通过对图像重建结果逐步修正求出多基材料密度图像等信息。该方法虽然重建图像质量较高，但需要预先精确标定出X射线源的多个能谱以及在不同能量下的基材料的衰减系数等，且计算量大，收敛速度慢，因此难以在实际多能谱CT成像***中应用。投影域分解成像方法先在投影域分解出基材料密度图像的投影值，然后利用传统CT方法进行重建出基材料的密度图像。该方法步骤简单、成像质量高，但仅仅适用于各种能谱的多色投影几何参数一致的情况，因此不适用于当前的多能谱CT***。图像域分解法先用传统CT重建方法分别计算出物体的高能和低能衰减系数图像，然后对重建图像在图像域进行线性组合重建出基材料的密度图像。图像域分解法是一种近似成像方法，在分解成像时忽略了多能谱投影的乘积项，所重建的基材料密度图像有严重伪影。鉴于目前常用的双能谱CT成像***以及螺旋锥束扫描等多种CT扫描模式采集不到几何参数一致的多色投影，图像域分解成像方法仍然是最常用的多能谱CT成像方法，因此需要开发一种新的适用于多色投影几何参数不一致情况下的高质量多能谱CT图像重建方法。

发明内容

本发明提供一种基于投影估计的多能谱CT图像重建方法，用以提高在多能谱多色投影的几何参数不一致情况下的重建图像质量。

为达到上述目的，本发明提供了一种基于投影估计的多能谱CT图像重建方法，包括以下步骤：

S1：利用单能CT重建算法由所采集的各个能谱的多色投影数据分别直接重建被测物体的能量图像；

S2：计算每一个能量图像沿着所有其它能谱的射线方向的线积分，估计在这些方向上的当前能谱的多色投影，得到几何参数一致的多色投影；

S3：标定多基材料分解函数，将几何参数一致的多色投影分解为多个基材料的线积分；

S4：利用单能CT图像重建方法分别由多个基材料的线积分重建出所对应的多基材料密度图像。

进一步地，在步骤S3中，标定多基材料分解函数的步骤具体包括：

S31：用多种基材料做成一个标定模体，标定模体中各种基材料在各个方向上的线积分组合涵盖被测物体中各种基材料在各个方向上的线积分组合，利用多能谱CT以与被测物体相同的扫描模式测量标定模体，获得多能谱多色投影数据；

S32：利用多色投影估计方法计算出标定模体的几何参数一致的多能谱多色投影；

S33：在图像域内通过最小化多色投影组合的重建图像和模体图像之间的距离，计算出多基材料分解函数的系数，得到标定了的多基材料分解函数。

进一步地，步骤S31中标定模体由多种基材料组合而成，其中基材料种数小于等于能谱个数，当基材料种数小于能谱个数时，通过组合不同的能谱投影数据重建基材料密度图像。

进一步地，在多色投影估计之前，对被测物体中的主要成分预先进行线性化校正，以提高投影估计精度。

进一步地，单能CT图像重建方法为解析方法或迭代重建方法。

进一步地，在利用单能CT图像重建方法进行图像重建中，添加物理约束条件，使当前重建的密度图像满足物理约束条件，物理约束条件包括像素值非负约束和/或分段光滑约束。

进一步地，步骤S1中利用单能CT重建算法由各个能谱的多色投影数据分别直接重建被测物体图像的公式如下：

$f_{q_m}\left(x\right)=R^{-1}\left(q_m\left(L\right)\right)$

这里q_m(L)表示第m种能谱的多色投影数据，m为自然数；R^-1是Radon逆变换算子，表示单能CT图像重建方法；

表示所重建的该种能谱下的被测物体的能量图像。

进一步地，步骤S2中通过计算每一个能量图像沿着其它能谱的射线方向的线积分，估计在这些方向上的当前能谱的多色投影的公式如下：

${q^{\′}}_m\left(L\right)={\∫}_L{f}_{q_m}\left(x\right)\mathrm{dl},L\∉{\mathrm{\Ω}}_m$

这里参数m表示能谱的序号，m=1,2,…,M，M为自然数，即共有M个能谱扫描被测物体；射线方向L是其它能谱扫描被测物体时的射线方向，不属于当前第m种能谱的射线路径集合Ω_m；q′_m(L)是沿着射线方向L穿过当前第m个能量图像

的线积分，是对第m种能谱在L方向的多色投影估计。

进一步地，步骤S3中利用所估计的几何参数一致的多色投影，分解为多个基材料线积分的公式为：

$F_m\left(L\right)=C_m({q^{\′}}_1\left(L\right),{q^{\′}}_2\left(L\right),...,{q^{\′}}_m\left(L\right))$

$=\underset{k_1,k_2,...,k_M}{\mathrm{\Σ}}{c}_{m,k_1,k_2,...,k_M}{q^{\′}}_1{\left(L\right)}^{k_1}{q^{\′}}_2{\left(L\right)}^{k_2}...{q^{\′}}_M{\left(L\right)}^{k_M}$

这里F_m(L)是其中第m个基材料的线积分，m=1,2,…,M，M为自然数；C_m是待标定的多基材料分解函数，用多次多项式函数表示，参数

是多次多项式函数的系数；k₁,k₂,…,k_M分别是所对应多色投影的指数；射线路径L∈Ω₁∪Ω₂∪…∪Ω_M，即全部M个能谱的所有X射线路径的集合；如果当前射线路径L∈Ω_m，则q′_m(L)=q_m(L)，即使用已测量得到的多色投影q_m(L)作为q′_m(L)；否则使用估计出的多色投影q′_m(L)。

进一步地，步骤S4中利用单能CT图像重建方法分别由多个基材料线积分重建出所对应的多基材料密度图像的公式是：

f_m(x)=R^-1(F_m(L))

这里f_m(x)是重建出的多基材料密度图像；F_m(L)是其中第m个基材料的线积分，m=1,2,…,M，M为自然数；R^-1是Radon逆变换算子，表示单能CT图像重建方法。

本发明适用于多能谱多色投影数据几何参数完全不一致的情况，仅利用所测投影数据估计出几何参数一致的多色投影，便可使图像重建质量比基于几何参数不一致投影数据的传统图像域分解法的重建质量大大提高，接近于基于几何参数一致投影数据的传统投影域分解法的重建效果。另外本发明方法不需要预先测定能谱信息，抗噪性能强，对于同一类被测物体，材料分解函数系数仅需要标定一次。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1为本发明一个实施例的多能谱CT图像重建方法流程图；

图2为实验用测试模体的断层图像；

图3为实验用标定模体的图像；

图4为包含高、低两种能谱的实验用X射线能谱示意图；

图5a为在80kV管电压下采集的投影数据图像；

图5b为在140kV管电压下采集的与低能谱多色投影几何参数不一致的高能谱多色投影数据图像；

图5c为利用本发明方法所估计的与低能谱多色投影几何参数一致的高能谱多色投影数据图像；

图6a为在所测多色投影数据无噪声的情况下采用传统图像域分解法重建的骨基材料密度图像；

图6b为在所测多色投影数据无噪声的情况下采用传统图像域分解法重建的水基材料密度图像；

图6c为在所测多色投影数据无噪声的情况下采用传统图像域分解法重建的70keV光子的线性衰减系数图像；

图6d为在所测多色投影数据无噪声的情况下采用本发明方法重建的骨基材料密度图像；

图6e为在所测多色投影数据无噪声的情况下采用本发明方法重建的水基材料密度图像；

图6f为在所测多色投影数据无噪声的情况下采用本发明方法重建的70keV光子的线性衰减系数图像；

图6g为在所测多色投影数据无噪声的情况下采用传统投影域分解法重建的骨基材料密度图像；

图6h为在所测多色投影数据无噪声的情况下采用传统投影域分解法重建的水基材料密度图像；

图6i为在所测多色投影数据无噪声的情况下采用传统投影域分解法重建的70keV光子的线性衰减系数图像；

图7a为在所测多色投影数据有噪声的情况下采用传统图像域分解法重建的骨基材料密度图像；

图7b为在所测多色投影数据有噪声的情况下采用传统图像域分解法重建的水基材料密度图像；

图7c为在所测多色投影数据有噪声的情况下采用传统图像域分解法重建的70keV光子的线性衰减系数图像；

图7d为在所测多色投影数据有噪声的情况下采用本发明方法重建的骨基材料密度图像；

图7e为在所测多色投影数据有噪声的情况下采用本发明方法重建的水基材料密度图像；

图7f为在所测多色投影数据有噪声的情况下采用本发明方法重建的70keV光子的线性衰减系数图像；

图7g为在所测多色投影数据有噪声的情况下采用传统投影域分解法重建的骨基材料密度图像；

图7h为在所测多色投影数据有噪声的情况下采用传统投影域分解法重建的水基材料密度图像；

图7i为在所测多色投影数据有噪声的情况下采用传统投影域分解法重建的70keV光子的线性衰减系数图像。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

多能谱CT成像问题可表达为式(1)的数学模型：

$q_m\left(L\right)=-\ln \∫S_m\left(E\right)e^{{-\∫}_L\mathrm{\μ}(x,E)\mathrm{dl}}\mathrm{dE},L\∈{\mathrm{\Ω}}_m.---\left(1\right)$

式(1)表示用M个能谱扫描被测物体，得到M组多色投影数据。其中参数m表示能谱的序号，m=1,2,…,M，M为自然数；q_m(L)是用第m个能谱的X射线扫描被测物体时，沿着射线路径L所探测到的X射线多色投影值；μ(x,E)为点x处关于能量E的线性衰减系数；S_m(E)为第m个归一化的X射线能谱；Ω_m为第m个能谱的所有X射线路径的集合。在理想的“同向”情况下，沿着任意一条X射线路径，都能采集到所有能谱的多色投影，即多色投影的几何参数一致。在一般的情况下，不是沿着所有的射线都能采集到每一种能谱的多色投影，甚至在一些多能谱CT扫描模式下，测得的不同能谱的多色投影在几何参数完全不一致。

当用M种X射线能谱扫描被测物体时，多能谱CT获得M套被测物体的多色投影数据。利用这些投影数据，可将被测物体作为M种物质的组合，重建出被测物体中这M种物质的分解图像。为此，将方程(1)中的μ(x,E)分解为如式(2)所示的M种与物质特性相关量的线性组合：

其中f_m(x)表示第m种基材料的密度空间分布，只与空间位置有关，即基材料密度图像，简称为基图像；表示第m种基材料的质量衰减系数，只与能量有关，可以预先测定。将式(2)代入式(1)得到：

这里F_m(L)表示沿着射线L穿过f_m(x)的线积分：

F_m(L)=∫_Lf_m(x)dl    (4)

于是，多能谱成像问题变为一个由M套投影数q_i(L)求解M个基材料密度图像f_m(x)的问题，即需要首先求解公式(3)的逆变换公式,然后由F_m(L)利用传统单能CT重建算法重建出f_m(x)。在高分辨率成像情况下，多能谱CT图像重建问题是一个高维数的非线性反问题，由公式(3)难以解析计算得到其逆变换。可用多元多项式函数拟合其逆变换，然后通过测量标准模体对多元多项式函数进行标定，直接得到利用多色投影重建基材料图像的逆变换。

当不同能谱的多色投影的几何参数一致时，可用如下多元多项式函数拟合逆变换：

$\begin{array}{c}{F}_m\left(L\right)=C_m(q_1\left(L\right),q_2\left(L\right),...,q_M\left(L\right))\\ =\underset{k_1,k_2,...,k_M}{\mathrm{\Σ}}{c}_{m,k_1,k_2,...,k_M}{q}_1{\left(L\right)}^{k_1}{q}_2{\left(L\right)}^{k_2}...q_M{\left(L\right)}^{k_M}\end{array}---\left(5\right)$

这里函数C_m是将多能谱多色投影数据分解为基材料线积分F_m(L)的分解函数，参数

表示多元多次多项式函数的系数。公式(5)充分利用了多项式函数的性质，实验表明具有很高的拟合精度。然而由于在该公式中有不同能谱的多色投影的乘积项，仅适用于几何参数一致的多色投影。当不同能谱的多色投影在几何参数不一致时，沿着同一条X射线路径L不能采集到不同能谱的多色投影，因此不能包含不同能谱的多色投影的乘积项，这时公式(5)变为：

$\begin{array}{c}{F}_m\left(L\right)=C_m(q_1\left(L\right),q_2\left(L\right),...,q_M\left(L\right))\\ =\underset{k_1}{\mathrm{\Σ}}{c}_{m,k_1,0,...,0}{q}_1{\left(L\right)}^{k_1}+\underset{k_2}{\mathrm{\Σ}}{c}_{m,0,k_2,...,0}{q}_2{\left(L\right)}^{k_2}+...+\underset{k_M}{\mathrm{\Σ}}{c}_{m,\mathrm{0,0},...,k_M}{q}_m{\left(L\right)}^{k_M}\end{array}---\left(6\right)$

由公式(6)可见，该公式中忽略了所有的多色投影的乘积项，适用于几何参数不一致的多色投影扫描模式，但对公式(3)逆变换的拟合精度大大降低，因此该方法求解误差大，重建图像中有严重的伪影。

为此，本发明提出了一种基于投影估计的多能谱CT图像重建方法。该方法用于在多色投影几何参数不一致情况下的高质量多能谱CT图像重建。该方法的关键是先用测量到的几何参数不一致的多色投影数据集估计出几何参数一致的多色投影，然后重建图像，具体步骤如下：

1)利用传统单能CT重建算法由各个能谱的多色投影数据分别直接重建所对应的被测物体的图像：

$f_{q_m}\left(x\right)=R^{-1}\left(q_m\left(L\right)\right)---\left(7\right)$

这里q_m(L)表示第m种能谱的多色投影数据；R^-1是Radon逆变换算子，表示传统单能CT图像重建方法；

表示所重建的该种能谱下的被测物体的图像，这种由多色投影直接重建的图像称为能量图像。

2)计算每一个能量图像沿着所有其它能谱的射线方向的线积分，估计在这些方向上的当前能谱的多色投影，求得几何参数一致的多色投影。

${q^{\′}}_m\left(L\right)={\∫}_L{f}_{q_m}\left(x\right)\mathrm{dl},L\∉{\mathrm{\Ω}}_m---\left(8\right)$

这里射线方向L是其它能谱扫描被测物体时的射线方向；q′_m(L)是沿着射线方向L穿过当前第m个能量图像

的线积分，是对第m种能谱在L方向的多色投影估计。当所采集的高低能谱投影足够多时，利用公式(7)能重建出高质量的能量图像

根据Radon变换的性质，由公式(8)计算得到的线积分是对高能量能谱多色投影的一个很精确的估计，由此可以获得所有射线方向上各种能谱彼此几何参数一致的多色投影，因此本发明方法在理论上是合理的。在利用基材料标定方法进行多能谱图像重建时，需要采集大量角度的多色投影，这为估算出高精度的几何参数一致多色投影的创造了条件，因此本发明方法同时具有很强的实用性。尽管如此，由于q_m(L)是多色投影，不是理想的单能量X射线扫描下得到的单色投影，因此利用传统单能CT重建算法所重建的图像

中具有一定的硬化伪影。当低能量图像中有较严重的硬化伪影时，本发明可以只选择较高能量能谱扫描时所采集的多色投影数据来估计低能谱射线方向的投影，以提高图像重建质量。这是因为利用高能谱扫描时，硬化现象不明显，所重建的图像

更接近于利用单色投影所重建的理想图像，所估计的多色投影更精确。

3)利用所估计的几何参数一致的多色投影，分解为M个基材料线积分。

$\begin{array}{c}{F}_m\left(L\right)=C_m({q^{\′}}_1\left(L\right),{q^{\′}}_2\left(L\right),...,{q^{\′}}_M\left(L\right))\\ =\underset{k_1,k_2,...,k_M}{\mathrm{\Σ}}{c}_{m,k_1,k_2,...,k_M}{q^{\′}}_1{\left(L\right)}^{k_1}{q^{\′}}_2{\left(L\right)}^{k_2}...{q^{\′}}_M{\left(L\right)}^{k_M}\end{array}---\left(9\right)$

公式(9)中，F_m(L)是其中第m个基材料的线积分；射线路径L∈Ω₁∪Ω₂∪…∪Ω_M，即全部M个能谱的所有X射线路径的集合。如果当前射线路径L∈Ω_m，则q′_m(L)=q_m(L)，即使用已测量得到的多色投影q_m(L)作为q′_m(L)；否则使用(8)式估计出的多色投影q′_m(L)。

4)利用传统单能CT图像重建方法分别由基材料线积分F_m(L)重建出所对应的多基材料密度图像f_m(x)。

f_m(x)=R^-1(F_m(L))    (10)

本发明方法中公式(9)的材料分解函数系数需要进行标定，具体的标定方法是：

1)首先用多种基材料做成一个标定模体，该模体中各种基材料在各个方向上的线积分组合应该能完全涵盖被测物体中各种基材料在各个方向上的线积分组合。利用多能谱CT以与被测物体相同的扫描模式测量标定模体，获得多能谱多色投影数据。

2)利用上述多色投影估计方法计算出标定模体的几何参数一致的多能谱多色投影；

3)利用传统单能CT重建算法重建公式(9)中各种多色投影组合的图像。利用Radon变换的线性性质，通过在图像域内通过最小化多色投影组合的重建图像和模体图像之间的距离，计算出公式(9)中的分解函数的系数得到所有标定了的分解函数C_m。

图1为本发明一个实施例的基于投影估计的多能谱CT图像重建方法流程图；如图所示，其中左侧的流程图为标定多能谱CT基材料分解函数的过程，右侧的流程图为利用基材料分解函数，进行多基材料密度图像重建的过程。该多能谱CT图像重建方法具体包括以下步骤：

S1：利用单能CT重建算法由所采集的各个能谱的多色投影数据分别直接重建被测物体的能量图像；

各个能谱的多色投影数据通过多能谱CT***采集，多能谱CT***可以是包括用两种X射线能谱扫描的CT***，即双能谱CT***，也可以是包括用两种以上X射线能谱扫描的CT***。双能谱CT***是多能谱CT成像领域内最常见的成像***，在下面的部分中本发明以双能谱CT图像重建的实施为例，进一步说明本发明方法的具体实现过程。

X射线能谱投影既包括传统X射线探测器探测到的宽X射线能谱投影，也包括光子计数探测器探测到的窄X射线能谱投影。所能适用的扫描模式包括多能谱几何参数一致的扫描以及几何参数不一致的扫描。

其中，单能CT图像重建方法可以为解析方法或迭代重建方法。

步骤S1中利用单能CT重建算法由各个能谱的多色投影数据分别直接重建被测物体图像的公式如下：

$f_{q_m}\left(x\right)=R^{-1}\left(q_m\left(L\right)\right)$

这里q_m(L)表示第m种能谱的多色投影数据，m为自然数；R^-1是Radon逆变换算子，表示单能CT图像重建方法；

表示所重建的该种能谱下的被测物体的能量图像。

S2：计算每一个能量图像沿着所有其它能谱的射线方向的线积分，估计在这些方向上的当前能谱的多色投影，得到几何参数一致的多色投影；

步骤S2中通过计算每一个能量图像沿着其它能谱的射线方向的线积分，估计在这些方向上的当前能谱的多色投影的公式如下：

${q^{\′}}_m\left(L\right)={\∫}_L{f}_{q_m}\left(x\right)\mathrm{dl},L\∉{\mathrm{\Ω}}_m$

这里参数m表示能谱的序号，m=1,2,…,M，M为自然数，即共有M个能谱扫描被测物体；射线方向L是其它能谱扫描被测物体时的射线方向，不属于当前第m种能谱的射线路径集合Ω_m；q′_m(L)是沿着射线方向L穿过当前第m个能量图像

的线积分，是对第m种能谱在L方向的多色投影估计。

进一步地，步骤S3中利用所估计的几何参数一致的多色投影，分解为多个基材料线积分的公式为：

$\begin{array}{c}{F}_m\left(L\right)=C_m({q^{\′}}_1\left(L\right),{q^{\′}}_2\left(L\right),...,{q^{\′}}_M\left(L\right))\\ =\underset{k_1,k_2,...,k_M}{\mathrm{\Σ}}{c}_{m,k_1,k_2,...,k_M}{q^{\′}}_1{\left(L\right)}^{k_1}{q^{\′}}_2{\left(L\right)}^{k_2}...{q^{\′}}_M{\left(L\right)}^{k_M}\end{array}$

这里F_m(L)是其中第m个基材料的线积分，m=1,2,…,M，M为自然数；C_m是待标定的多基材料分解函数，用多次多项式函数表示，参数

是多次多项式函数的系数；k₁,k₂,…,k_M分别是所对应多色投影的指数；射线路径L∈Ω₁∪Ω₂∪…∪Ω_M，即全部M个能谱的所有X射线路径的集合。如果当前射线路径L∈Ω_m，则q′_m(L)=q_m(L)，即使用已测量得到的多色投影q_m(L)作为q′_m(L)；否则估计出的多色投影q′_m(L)。

其中，在多色投影估计之前，还可以对被测物体中的主要成分预先进行线性化校正，以提高投影估计精度。如在医学CT中，进行“水预校正”，即求出被测物体和标定模体的等效水厚度。

投影估计时，还可以调整投影方向，生成如平行束等更易计算的投影扫描模式下的投影数据，并依此作图像重建。

S3：标定多基材料分解函数，将几何参数一致的多色投影分解为多个基材料的线积分；

其中，在步骤S3中，标定多基材料分解函数的步骤具体包括：

S31：用多种基材料做成一个标定模体，标定模体中各种基材料在各个方向上的线积分组合涵盖被测物体中各种基材料在各个方向上的线积分组合，利用多能谱CT以与被测物体相同的扫描模式测量标定模体，获得多能谱多色投影数据；

步骤S31中的标定模体由多种基材料组合而成，其中基材料种数小于等于能谱个数，当基材料种数小于能谱个数时，通过组合不同的能谱投影数据重建基材料密度图像。标定模体中的多种基材料在不同CT扫描方向的厚度组合能涵盖被测物体中这些基材料所有可能的厚度组合。

S32：利用多色投影估计方法计算出标定模体的几何参数一致的多能谱多色投影；

S33：在图像域内通过最小化多色投影组合的重建图像和模体图像之间的距离，计算出多基材料分解函数的系数，得到标定了的多基材料分解函数。

其中对于同一类被测物体，多基材料分解函数标定仅需要标定一次。

S4：利用单能CT图像重建方法分别由多个基材料的线积分重建出所对应的多基材料密度图像。

其中，在利用单能CT图像重建方法进行图像重建中，为提高图像重建质量还可以添加物理约束条件，使当前重建的密度图像满足物理约束条件，物理约束条件包括像素值非负约束和/或分段光滑约束。

步骤S4中利用单能CT图像重建方法分别由多个基材料线积分重建出所对应的多基材料密度图像的公式是：

f_m(x)=R^-1(F_m(L))

这里f_m(x)是重建出的多基材料密度图像；F_m(L)是其中第m个基材料的线积分，m=1,2,…,M，M为自然数；R^-1是Radon逆变换算子，表示单能CT图像重建方法。

本发明的上述实施例适用于多能谱多色投影数据几何参数完全不一致的情况，仅利用所测投影数据估计出几何参数一致的多色投影，便可使图像重建质量比基于几何参数不一致投影数据的传统图像域分解法的重建质量大大提高，接近于基于几何参数一致投影数据的传统投影域分解法的重建效果。另外本发明方法不需要预先测定能谱信息，抗噪性能强，对于同一类被测物体，材料分解函数系数仅需要标定一次。

图2为本发明一个实施例中所用的测试模体的断层图像，该测试模体是一个标准的Forbild胸腔模型。选定水和骨组织为双基材料，水的密度为1.0g/cm³，骨组织的密度为1.92g/cm³。图3为本实施例中所用的标定模体，用于标定分解函数的系数。该模体包含一个较大的圆盘和一个较小的圆盘，其中较大圆盘的直径为120mm，较小的直径为70mm。大圆盘中填充了水，而小圆盘为骨组织。双能谱CT***发出高、低X射线能谱如图4所示。这两个能谱是X光管的管电压分别为80kV和140kV（1mm铜滤波片）下发出的。在数值实现时，能谱的采样间隔设为1keV。在双能谱CT扫描中，采用扇束扫描。***的扫描参数为：射线源到探测器的距离为1200mm，射线源到转台中心的距离为1000mm，探测器由1024个0.3mm的探测单元组成，视野的直径为256mm。射线源旋转一周，交错均匀采集1440个角度的高能和低能投影数据各一组，即相邻的高能投影和低能投影有0.125度的夹角，这样得到了几何参数不一致的高、低能多色投影数据。重建图像大小为512×512。在80kV和140kV管电压下分别采集到的测试模体的多色投影数据如图5所示，其中图5a为80kV管电压下采集的低能谱多色投影数据；图5b为140kV管电压下采集的与低能谱多色投影几何参数不一致的高能谱多色投影数据；图5c为利用本发明方法所估计的与低能谱多色投影几何参数一致的高能谱多色投影数据。用归一化均方距离（NMSB）和归一化平均绝对距离(NMAD)与该方向真实的高能谱多色投影数据作量化比较，NMSB和NMAD的值分别为0.012356和0.009745，表明本发明方法估计的高能谱多色投影数据精度很高。

图6为在所测多色投影数据无噪声的情况下，传统图像域分解法、本发明方法和传统投影域分解法这三种方法重建效果的比较，其中图6a、6b、6c为传统图像域分解法重建的图像，图6d、6e、6f为本方法重建的图像，图6g、6h、6i为传统投影域分解法重建的图像，图6a、6d、6g为水基材料密度图像，图6b、6e、6h是骨基材料密度图像，图6c、6f、6i是能量为70keV的光子的线性衰减系数图像。这里线性衰减系数图像是通过将重建的水基图像与水对指定能量的X射线光子质量衰减系数的乘积以及重建的骨基图像与骨对指定能量的X射线光子质量衰减系数的乘积求和得到的。为了凸显所重建图像的特征，图6中用了窄显示灰度窗，其中图6a、6d、6g的显示灰度窗为[0.951.15],图6b、6e、6h的显示灰度窗为[0.71.2],图6c、6f、6i的显示灰度窗为[-150HU150HU]。图6中传统图像域分解法和本发明方法所用的投影数据都是几何参数不一致的多能谱多色投影，而传统投影域分解方法所用的投影数据是几何参数一致的多能谱多色投影，作为一个理想重建结果的参照。比较图6中三种方法分别重建的图像，可以看出本发明实施例的方法可以很好地分解出双基材料密度图像，基材料密度图像中仅有一些轻微的伪影。远远优于传统图像域分解法的重建结果，特别是水基材料密度图像。和基于几何参数一致投影的传统投影域分解法相比，本发明方法的重建效果虽然与之非常接近，特别是所合成的线性衰减系数图像的质量基本相当。

图7为在所测多色投影数据有噪声的情况下，传统图像域分解法、本发明方法和传统投影域分解法这三种方法重建效果的比较。同样，图7a、7b、7c为传统图像域分解法重建的图像，图7d、7e、7f为本方法重建的图像，图7g、7h、7i为传统投影域分解法重建的图像，图7a、7d、7g为骨基材料密度图像，图7b、7e、7h是水基材料密度图像，图7c、7f、7i是能量为70keV的光子的线性衰减系数图像。在多色投影数据中所添加的噪声为对应每条射线1,000,000个光子的泊松噪声。为了凸显所重建图像的特征，图7中仍然使用和图6中相同的窄显示灰度窗。从图7中可以看到，在所测投影数据被噪声污染情况下，本发明方法仍然能够重建出高质量的基材料图像和线性衰减系数图像，重建结果仍然优于传统图像域分解法的重建结果，接近于传统投影域分解法的重建结果。但是三种方法由于受到噪声的干扰，被测模体中的一些细小的低对比度的特征都重建不出来。

本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。

本领域普通技术人员可以理解：实施例中的装置中的模块可以按照实施例描述分布于实施例的装置中，也可以进行相应变化位于不同于本实施例的一个或多个装置中。上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围。

Claims (10)

1.一种基于投影估计的多能谱CT图像重建方法，其特征在于，包括以下步骤：

S1：利用单能CT重建算法由所采集的各个能谱的多色投影数据分别直接重建被测物体的能量图像；

S2：计算每一个能量图像沿着所有其它能谱的射线方向的线积分，估计在这些方向上的当前能谱的多色投影，得到几何参数一致的多色投影；

S3：标定多基材料分解函数，将所述几何参数一致的多色投影分解为多个基材料的线积分；

S4：利用单能CT图像重建方法分别由多个基材料的线积分重建出所对应的多基材料密度图像。

2.根据权利要求1所述的多能谱CT图像重建方法，其特征在于，在步骤S3中，标定多基材料分解函数的步骤具体包括：

S31：用多种基材料做成一个标定模体，标定模体中各种基材料在各个方向上的线积分组合涵盖被测物体中各种基材料在各个方向上的线积分组合，利用多能谱CT以与被测物体相同的扫描模式测量标定模体，获得多能谱多色投影数据；

S32：利用多色投影估计方法计算出标定模体的几何参数一致的多能谱多色投影；

S33：在图像域内通过最小化多色投影组合的重建图像和模体图像之间的距离，计算出多基材料分解函数的系数，得到标定了的多基材料分解函数。

3.根据权利要求2所述的多能谱CT图像重建方法，其特征在于，步骤S31中所述标定模体由多种基材料组合而成，其中基材料种数小于等于能谱个数，当基材料种数小于能谱个数时，通过组合不同的能谱投影数据重建基材料密度图像。

4.根据权利要求1所述的多能谱CT图像重建方法，其特征在于，在多色投影估计之前，对被测物体中的主要成分预先进行线性化校正，以提高投影估计精度。

5.根据权利要求1或2所述的多能谱CT图像重建方法，其特征在于，所述单能CT图像重建方法为解析方法或迭代重建方法。

6.根据权利要求1所述的多能谱CT图像重建方法，其特征在于，在利用单能CT图像重建方法进行图像重建中，添加物理约束条件，使当前重建的密度图像满足所述物理约束条件，所述物理约束条件包括像素值非负约束和/或分段光滑约束。

7.根据权利要求1所述的多能谱CT图像重建方法，其特征在于，步骤S1中利用单能CT重建算法由各个能谱的多色投影数据分别直接重建被测物体图像的公式如下：

$f_{q_m}\left(x\right)=R^{-1}\left(q_m\left(L\right)\right)$

这里q_m(L)表示第m种能谱的多色投影数据，m为自然数；R^-1是Radon逆变换算子，表示单能CT图像重建方法；

表示所重建的该种能谱下的被测物体的能量图像。

8.根据权利要求1所述的多能谱CT图像重建方法，其特征在于，步骤S2中通过计算每一个能量图像沿着所有其它能谱的射线方向的线积分，估计在这些方向上的当前能谱的多色投影的公式如下：

${q^{\′}}_m\left(L\right)={\∫}_L{f}_{q_m}\left(x\right)\mathrm{dl},L\∉{\mathrm{\Ω}}_m$

这里参数m表示能谱的序号，m=1,2,…,M，M为自然数，即共有M个能谱扫描被测物体；射线方向L是其它能谱扫描被测物体时的射线方向，不属于当前第m种能谱的射线路径集合Ω_m；q′_m(L)是沿着射线方向L穿过当前第m个能量图像

的线积分，是对第m种能谱在L方向的多色投影估计。

9.根据权利要求1所述的多能谱CT图像重建方法，其特征在于，步骤S3中利用所估计的几何参数一致的多色投影，分解为多个基材料线积分的公式为：

$F_m\left(L\right)=C_m({q^{\′}}_1\left(L\right),{q^{\′}}_2\left(L\right),...,{q^{\′}}_m\left(L\right))$

$=\underset{k_1,k_2,...,k_M}{\mathrm{\Σ}}{c}_{m,k_1,k_2,...,k_M}{q^{\′}}_1{\left(L\right)}^{k_1}{q^{\′}}_2{\left(L\right)}^{k_2}...{q^{\′}}_M{\left(L\right)}^{k_M}$

这里F_m(L)是其中第m个基材料的线积分，m=1,2,…,M，M为自然数；C_m是待标定的多基材料分解函数，用多次多项式函数表示，参数

是多次多项式函数的系数；k₁,k₂,…,k_M分别是所对应多色投影的指数；射线路径L∈Ω₁∪Ω₂∪…∪Ω_M，即全部M个能谱的所有X射线路径的集合；如果当前射线路径L∈Ω_m，则q′_m(L)=q_m(L)，即使用已测量得到的多色投影q_m(L)作为q′_m(L)；否则使用估计出的多色投影q′_m(L)。

10.根据权利要求1所述的多能谱CT图像重建方法，其特征在于，步骤S4中利用单能CT图像重建方法分别由多个基材料线积分重建出所对应的多基材料密度图像的公式是：

f_m(x)=R^-1(F_m(L))

这里f_m(x)是重建出的多基材料密度图像；F_m(L)是其中第m个基材料的线积分，m=1,2,…,M，M为自然数；R^-1是Radon逆变换算子，表示单能CT图像重建方法。

Images