CN103559729A - 一种双能谱ct图像迭代重建方法 - Google Patents

Abstract

本发明公开一种双能谱CT图像迭代重建方法，用以重建被测物体的基材料密度图像。该方法包括：1)根据被测物体的尺寸和材料，选定一个包含两种不同基材料的标定模体；通过双能谱CT系统测量该模体，标定表示多色投影值与基材料厚度组合之间曲线关系的双能谱多色正投影方程；2)利用双能谱CT系统对被测物体扫描，采集被测物体的双能谱多色投影数据；(3)根据所标定的多色正投影方程，构造迭代格式，通过与所实测多色投影数据作比较，并加入物理约束条件，迭代重建被测物体的双基材料密度图像。与现有方法相比，本发明方法不需要预先测定双能CT系统的高、低能谱以及基材料的衰减系数，而且同时适用于高、低能投影在几何上相匹配和不匹配的情况。

Description

一种双能谱CT图像迭代重建方法

技术领域

本发明涉及X射线CT成像技术领域，具体而言，涉及一种双能谱CT图像迭代重建方法。

背景技术

X射线计算机层析成像技术（简称X射线CT）是一种利用X射线与物质的相互作用原理，对物体内部信息进行成像的一种技术。自从1967年第一台X射线医学CT问世以来，X射线CT技术得到了迅速发展。各种CT设备为医学诊疗带来了革命性变化，成为众多疾病的常规检测手段。CT技术也成为一项重要的工业支撑技术，广泛应用于汽车、船舶、飞机、电子、石油等行业关键部件的无损检测。传统的单能谱X射线CT使用一个X射线能谱对被测物体进行扫描，不同物质的CT重建图像可能具·有相同或相似的CT值，难以区分不同的物质。双能谱CT使用两个X射线能谱对被测物体进行扫描，能够获得比传统单能谱CT更多的被测物体信息。利用这些信息可以重建出被测物体等效原子序数和电子密度图像，或是重建出特定基材料的密度图像，基于这些信息对物质进行区分。由于双能谱CT具有更好的物质区分能力，因此具有广泛的应用前景，比如：骨密度测量、PET的衰减矫正、骨髓成分分析、石油岩心分析和安检等。另外，双能谱CT还能有效去除重建图像中的硬化伪影。

双能谱CT图像重建问题是一个非线性反问题，具有非线性性、多解性与不适定性、高维数等特点，难以直接求解。由于X射线源发出的是服从一定谱分布的多色X射线，而不是单能X射线，因此双能谱CT扫描被测物体采集到的是多色投影，而不是理想的单色投影。相关重建方法的研究始于上个世纪七十年代。现有的成像方法大致可以分为三类：图像域分解法、投影域分解法和迭代法。图像域分解法通常先用传统CT重建方法分别计算出物体的高能和低能衰减系数图像，然后对重建图像在图像域进行线性组合以获得基材料的密度图像。图像域分解法是一种近似成像方法，在分解成像时忽略了高、低能谱投影的乘积项，所重建的基材料密度图像有严重伪影。投影域分解成像方法先用测得的高、低能谱投影在投影域估计出基材料密度图像的投影值，然后利用传统CT方法进行重建以获得基材料的密度图像。该方法的成像质量虽然优于图像域分解成像方法，但仅仅适用于高、低能谱几何参数一致扫描模式，即沿着每一条X射线，均可以采集到高、低能谱两种多色投影。然而当前临床医用的双能谱CT成像系统，包括双源双探测器系统和快速能量切换的单源单探测器系统等，采集的高、低能谱投影都是几何参数不一致的投影，且常用的多种双能谱CT扫描模式也采集不到几何参数一致的投影。因此目前图像域分解成像方法仍然是一种双能谱CT成像系统中最常用的方法。迭代法基于双能谱投影的非线性解析方程，利用数值方法或者优化方法构造迭代结构，通过对图像重建结果逐步修正求出基材料密度图像等信息。该方法虽然重建图像质量较高，但需要预先精确标定出X射线源的高、低能谱以及在不同能量下的基材料的衰减系数等，且计算量大，收敛速度慢，耗时长，因此难以在实际双能谱CT成像系统中应用。另外该方法忽略了X射线散射的影响，导致计算结果准确性下降。

发明内容

本发明提供一种双能谱CT图像迭代重建方法，用以解决现有技术中当所采集的高、低能谱投影几何参数不一致时，传统的图像域分解法重建的双能谱CT图像质量显著下降的问题。

为达到上述目的，本发明提供了一种双能谱CT图像迭代重建方法，包括以下步骤：

S1：根据多色投影值与基材料厚度组合之间的曲线关系获取双能谱多色正投影方程；

S2：利用双能谱CT系统扫描被测物体，采集被测物体的高、低能谱多色投影数据；

S3：根据所述双能谱多色正投影方程，构造迭代格式，利用所测得的高、低能谱多色投影数据，迭代重建被测物体的双基材料密度图像。

进一步地，其中所述步骤S1包括：

S11：根据被测物体的尺寸和材料，选定一个包含两种不同的基材料的标定模体；

S12：利用双能谱CT系统扫描标定模体，采集标定模体在各个角度下的双能谱多色投影数据；

S13：通过最小化标定模体中不同双基材料厚度组合与其双能谱多色投影值之间的距离，标定双能谱CT系统的双能谱多色正投影方程。

进一步地，其中所述步骤S3包括：

S31：根据所述双能谱多色正投影方程，构造求双基材料密度图像的迭代格式；

S32：为待重建的被测物体的双基材料密度图像赋初始值；

S33：利用所述双能谱多色正投影方程，获取当前双基材料密度图像的多色投影的估计值，并与被测物体实测的多色投影的真实值相比较，得到所重建的双基材料密度图像估计值的误差；

S34：利用计算得到的双基材料密度图像估计值的误差，更新当前重建的双基材料密度图像；

S35：加入物理约束条件，使当前重建的双基材料密度图像满足所述物理约束条件；

S36：重复步骤S33-S35，直到满足迭代结束条件，重建出双基材料密度图像。

进一步地，所述双能谱CT系统的扫描模式包括圆周扫描、螺旋扫描、有限角度扫描、稀疏角度扫描及一般扫描模式中的至少一种。

进一步地，步骤S3中所述迭代格式是每个角度修正一次或是用有序子集方法修正。

进一步地，所述物理约束条件包括像素值非负约束和/或分段光滑约束。

进一步地，步骤S12中所述双能谱多色正投影方程用多项式方程、分段多项式方程或样条曲线方程表示；在用多项式方程表示多色正投影方程的情况下：

$q_{i,L}=C_i(F_L,G_L)=\underset{k=0}{\overset{K}{\mathrm{\Σ}}}\underset{m=0}{\overset{M}{\mathrm{\Σ}}}{c}_{i,\mathrm{km}}{F}_L^k{G}_L^m,$ i=1,2

其中参数i表示能谱索引，q_i,L是沿着给定X射线路径L所探测到的多色投影，表示为基函数

的一个线性组合，函数C_i的多项式系数c_i,km由前述的标定方法确定。

进一步地，步骤S33中由双基材料密度图像估计值的误差求基材料密度图像修正值的公式为：

$\left\{\begin{array}{c}{e}_{F_L}^{\left(n\right)}=\frac{\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂F_L}}{{\left(\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂F_L}\right)}^2+{\left(\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂G_L}\right)}^2}(q_{i,L}-q_{i,L}^{\left(n\right)}),\\ e_{G_L}^{\left(n\right)}=\frac{\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂G_L}}{{\left(\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂F_L}\right)}^2+{\left(\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂G_L}\right)}^2}(q_{i,L}-q_{i,L}^{\left(n\right)}),\end{array}\right.i=\mathrm{1,2}.$

其中为第n次迭代中多色投影的估计值，与测量得到的多色投影真实值q_i,L比较得到多色投影误差，沿着当前射线路径L直接分配给两种基材料，分别得到对两种基材料图像的修正值

和

进一步地，步骤S34中利用第n次迭代计算得到的双基材料密度图像估计值的误差，更新当前重建的双基材料密度图像的公式为：

其中参数N是所重建的基材料密度图像的像素总数，

和

分别是两幅基材料密度图像的第j个像素值，

和

是更新后的像素值，r_L,j表示基材料密度图像的第j个像素对沿着射线L的投影的贡献权值，参数λ是松弛因子，一旦一个投影Q_φ中所有射线对应的图像修正值都得到后，同时修正双基材料密度图像，获得更新后的像素值。

进一步地，被测物体的双基材料密度图像重建过程采用并行方法计算，并基于硬件并行计算平台加速实现。

与现有技术相比，本发明方法具有以下有益效果：

(1)无需预先测定双能CT系统的高、低能谱以及基材料在不同X射线能量下的衰减系数；

(2)适用性强，同时适用于高、低能谱几何参数一致和几何参数不一致扫描的情况，以及适用于常用的各种扫描模式，还可以推广应用于多能谱CT和慢电压切换CT系统；

(3)算法并行度高，适用于图形处理器(GPU)等并行硬件加速实现；

(4)图像重建精度高，仅通过少量迭代，即可重建出高质量的双基材料密度图像，并能有效去除所合成的衰减系数图像中的硬化伪影。

附图说明

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

图1为本发明一个实施例的双能谱CT迭代成像方法流程图；

图2为用作测试模体的标准的胸腔模型图像；

图3为用于标定双能谱CT多色正投影方程的实验用标定模体图像；

图4为实验用X射线高、低能谱示意图；

图5a为本发明一个实施例的在140kV管电压下采集到的测试模体的多色投影数据图像；

图5b为本发明一个实施例的在80kV管电压下采集到的测试模体的多色投影数据图像；

图6a为本发明一个实施例的水基材料密度图像重建的结果图像；

图6b为本发明一个实施例的骨基材料密度图像重建的结果图像；

图6c为本发明一个实施例的能量为70keV的光子的线性衰减系数图像；

图7a为采用传统图像域分解法重建的骨基材料密度图像的结果图像；

图7b为采用传统图像域分解法重建的水基材料密度图像的结果图像；

图7c为用图7a的骨基图像和图7b的水基图像合成的能量为70keV的光子的线性衰减系数图像；

图8为本发明一个实施例的沿着测试模体、本发明方法重建图像、传统图像域分解法重建图像的中心列线上的像素值之间的比较示意图；

图9a为本发明一个实施例的在所测多色投影数据中有噪声的情况下水基材料密度图像重建的结果图像；

图9b为本发明一个实施例的在所测多色投影数据中有噪声的情况下骨基材料密度图像重建的结果图像；

图9c为本发明一个实施例的在所测多色投影数据中有噪声的情况下关于能量为70keV的光子的线性衰减系数图像。

具体实施方式

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

在忽略X射线散射效应的情况下，双能谱CT图像重建问题可表示为如下的数学模型：

$q_{i,L}=-\ln \∫S_i\left(E\right)e^{-{\∫}_L\mathrm{\μ}(x.E)\mathrm{dl}}\mathrm{dE},L\∈{\mathrm{\Ω}}_i,i=\mathrm{1,2}.---\left(1\right)$

其中，参数i表示能谱索引，q_i,L表示沿着给定X射线路径L所探测到的多色投影值，S_i(E)表示第i个归一化的X射线能谱，μ(x,E)为点x处关于能量E的线性衰减系数，Ω_i是对应于第i个能谱的X射线扫描路径集合。在理想的情况下，Ω₁=Ω₂，意味着对应不同X射线能谱的多色投影在几何上相匹配。在一般的情况下，Ω₁≠Ω₂，意味着不是所有的射线都能被两个不同的X射线能谱测量到。例如在高、低电压快速切换的单源单探测器的双能谱CT系统中，高、低能谱在不同角度交替测量，测得的多色投影在几何上完全不匹配，即

。双能谱CT图像重建问题可归结为利用测得的多色投影q_i,L的集合求解μ(x,E)。

为了求解方程(1)，将μ(x,E)分解为与物质特性相关量的线性组合：

μ(x,E)=φ(E)f(x)+θ(E)g(x)            (2)

这里f(x)和g(x)分别表示被测物体的两种基材料的密度空间分布，即基材料密度图像。φ(E)和θ(E)分别表示这两种基材料的质量衰减系数。将式(2)代入式(1)得到：

$q_{i,L}=-\ln \∫S_i\left(E\right)e^{-\mathrm{\φ}\left(E\right)F_L-\mathrm{\θ}\left(E\right)G_L}\mathrm{dE},L\∈{\mathrm{\Ω}}_i,i=\mathrm{1,2}---\left(3\right)$

并且：

$\left\{\begin{array}{c}{F}_L=R\left(f\left(x\right)\right)={\∫}_{L}f\left(x\right)\mathrm{dl},\\ G_L=R\left(g\left(x\right)\right)={\∫}_{L}g\left(x\right)\mathrm{dl},\end{array}\right.---\left(4\right)$

这里F_L和G_L分别表示沿着射线L穿过f(x)和g(x)的线积分；R表示Radon变换（拉东变换）。

在数值实现时，被重建的基材料密度图像被离散化为像素。设f=(f₀,f₁,…,f_J-1)^τ和g=(g₀,g₁,…,g_J-1)^τ分别表示f(x)和g(x)的离散形式，这里f_j和g_j分别是在f(x)和g(x)的第j个像素上的采样值，参数J是每幅图像的像素总数，τ表示向量的转置。设R_L=(r_L,0,r_L,1,…,r_L,J-1)是一个投影向量，这里r_L,j表示f或g的第j个像素对沿着射线L的投影的贡献权值。公式(4)能重写为如下向量形式：

$\left\{\begin{array}{c}{F}_L=R_{L}f=\underset{j=0}{\overset{J-1}{\mathrm{\Σ}}}{r}_{L,j}{f}_j\\ G_L=R_{L}g=\underset{j=0}{\overset{J-1}{\mathrm{\Σ}}}{r}_{L,j}{g}_j\end{array}\right.---\left(5\right)$

本发明方法是双能谱CT图像的一种迭代重建方法，因此需要模拟双能谱CT的多色正投影过程。通常用多项式方程来近似公式(3)所示的多色正投影方程。

$q_{i,L}=C_i(F_L,G_L)=\underset{k=0}{\overset{K}{\mathrm{\Σ}}}\underset{m=0}{\overset{M}{\mathrm{\Σ}}}{c}_{i,\mathrm{km}}{F}_L^k{G}_L^m,i=\mathrm{1,2},---\left(6\right)$

公式中将多色投影表示为基函数的一个线性组合。函数C_i的多项式系数c_i,km可用标定方法确定。

设f⁽ⁿ⁾和g⁽ⁿ⁾分别是真实双基材料图像f和g在第n次迭代时的估计图像，

和

分别是估计图像和真实图像之间的误差，则:

$\left\{\begin{array}{c}f=f^{\left(n\right)}+e_f^{\left(n\right)},\\ g=g^{\left(n\right)}+e_g^{\left(n\right)}.\end{array}\right.---\left(7\right)$

利用Radon变换R_L的线性性质，在第n次迭代时双基材料图像的真实线积分F_L和G_L和估计线积分

和

之间的关系是：

$\left\{\begin{array}{c}{F}_L=R_{L}f=R_L{f}^{\left(n\right)}+R_L{e}_f^{\left(n\right)}=F_L^{\left(n\right)}+e_{F_L}^{\left(v\right)},\\ G_L=R_{L}g=R_L{g}^{\left(n\right)}+R_L{e}_g^{\left(n\right)}=G_L^{\left(n\right)}+e_{G_L}^{\left(n\right)}.\end{array}\right.---\left(8\right)$

将公式(8)代入公式(6)中，则得到：

$q_{i,L}=C_i(F_L^{\left(n\right)}+e_{F_L}^{\left(n\right)},G_L^{\left(n\right)}+e_{G_L}^{\left(n\right)}),i=\mathrm{1,2}.---\left(9\right)$

公式(9)用一阶泰勒展开近似：

$q_{i,L}\≈C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})+\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂F_i}{e}_{F_L}^{\left(n\right)}+\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂G_L}{e}_{G_L}^{\left(n\right)},i=\mathrm{1,2}.---\left(10\right)$

即：

$q_{i,L}-q_{i,L}^{\left(n\right)}\≈\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂F_L}{e}_{F_L}^{\left(n\right)}+\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂G_L}{e}_{G_L}^{\left(n\right)},i=\mathrm{1,2}.---\left(11\right)$

这里i=1,2,是在第n次迭代时多色投影的估计值，所以公式(11)的左式是测量得到的真实多色投影和估计的多色投影之间的误差。

本发明的方法不仅适用于高、低能谱几何参数一致扫描情况，而且适用于几何参数不一致扫描的情况。所谓高、低能谱几何参数一致扫描，是指能够获得每条X射线路径高、低能谱的投影数据；所谓几何参数不一致扫描，即采集高、低能谱投影的X射线路径部分或全部不一致。

在高、低能谱几何参数一致扫描情况下，公式(11)可用经典的牛顿-拉普森方法求解。在高、低能谱几何参数不一致扫描情况下，由于沿着每一条X射线，不一定均可以采集到高、低能两种多色投影，对于CT重建问题，公式(11)涉及到沿着所有射线同时求解。这种情况下公式(11)对应的雅可比矩阵太大，难以求逆，因此不能用经典的牛顿-拉普森方法求解。本发明使用逐条射线求解的方法，避免了对雅可比矩阵求逆。如公式(12)所示，多色投影残差

沿着当前射线直接分配给两种基材料，分别得到对两种基材料图像的修正值

和

：

$\left\{\begin{array}{c}{e}_{F_L}^{\left(n\right)}=\frac{\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂F_L}}{{\left(\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂F_L}\right)}^2+{\left(\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂G_L}\right)}^2}(q_{i,L}-q_{i,L}^{\left(n\right)}),\\ e_{G_L}^{\left(n\right)}=\frac{\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂G_L}}{{\left(\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂F_L}\right)}^2+{\left(\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂G_L}\right)}^2}(q_{i,L}-q_{i,L}^{\left(n\right)}),\end{array}\right.i=\mathrm{1,2}.---\left(12\right)$

根据公式(6)，可得到：

$\left\{\begin{array}{c}\frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂F_L}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{m=0}{\overset{M}{\mathrm{\Σ}}}{\mathrm{kc}}_{i,\mathrm{km}}{\left(F_L^{\left(\text{n}\right)}\right)}^{k-1}{\left(G_L^{\left(n\right)}\right)}^m\\ \frac{\∂C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\∂G_L}=\underset{k=0}{\overset{K}{\mathrm{\Σ}}}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{mc}}_{i,\mathrm{km}}{\left(F_L^{\left(n\right)}\right)}^k{\left(G_L^{\left(n\right)}\right)}^{m-1}\end{array}\right.,i=\mathrm{1,2}.---\left(13\right)$

设

和

分别是f⁽ⁿ⁾和g⁽ⁿ⁾第j个像素值，

和

是更新后的像素值。一旦一个投影Q_φ中所有射线对应的图像修正值都得到后，用下式同时修正双基材料密度图像，获得更新后的像素值。

这里参数λ是一个松弛因子，取值范围为0<λ<1，具体值由经验确定，通常取值为0.9，且随着迭代次数增加而减小。整个图像更新过程逐角度重复。当前的投影角度的选择很重要，应该按照连续投影之间尽量不相关的原则选择用于修正图像的投影。本发明的方法避免了求雅可比矩阵的逆，适用于求解大规模的稀疏的非线性方程组。以上所有公式计算在数值实现时，均可并行计算，利用并行硬件如图形处理器(GPU)等加速实现。

本发明方法中，关于双基材料密度图像迭代重建的具体实施步骤如下：

1）为待重建的双基材料密度图像赋初值f⁽⁰⁾和g⁽⁰⁾，设迭代次数n=0；

2）沿着给定的X射线路径L利用公式(5)分别计算当前双基材料密度图像f⁽ⁿ⁾和g⁽ⁿ⁾的线积分

和，然后用公式(6)计算高能和低能多色投影

；

3）由公式(12)求解出图像修正值

和；

4）一旦得到对应当前投影中所有射线的图像修正值，根据公式(14)更新当前双基材料密度图像f⁽ⁿ⁾和g⁽ⁿ⁾；

5）如果收敛标准满足，则停止迭代；否则设n=n+1，重复步骤2）-4）在下一个投影角度继续计算直到双基材料密度图像收敛。

实验表明，本发明方法仅通过少量迭代，就可能够重建出高质量的基材料密度图像。

图1为本发明一个实施例的双能谱CT图像迭代重建方法流程图；如图所示，其中左侧的流程图为标定双能谱CT多色正投影的步骤，右侧的流程图为利用标定好多色正投影方程，构造迭代格式，进行基材料密度图像迭代重建的步骤。图1实施例具体包括以下步骤：

S1：根据多色投影值与基材料厚度组合之间的曲线关系获取双能谱多色正投影方程；

在具体实现时，步骤S1可以包括以下步骤：

S11：根据被测物体的尺寸和材料，选定一个包含两种不同的基材料的标定模体；

其中步骤S11中标定模体由两种基材料组合而成，这两种基材料在不同CT扫描方向的厚度组合能够涵盖被测物体中两种基材料所有可能的厚度组合。

S12：利用双能谱CT系统扫描标定模体，采集标定模体在各个角度下的双能谱多色投影数据；

其中，双能谱CT系统包括双源双探测器的双能CT、快速能量切换的双能CT、重复两次扫描的双能CT等。双能谱CT系统所能适用的扫描模式包括圆周扫描、螺旋扫描、有限角度扫描、稀疏角度扫描及其他一般扫描模式。

步骤S12中的双能谱多色正投影方程可以用多项式方程、分段多项式方程或样条曲线方程表示；在用多项式方程表示多色正投影方程的情况下：

$q_{i,L}=C_i(F_L,G_L)=\underset{k=0}{\overset{K}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{M}{\mathrm{\&Sigma;}}}{c}_{i,\mathrm{km}}{F}_L^k{G}_L^m,$ i=1,2

其中参数i表示能谱索引，q_i,L是沿着给定X射线路径L所探测到的多色投影，表示为基函数

的一个线性组合，函数C_i的多项式系数c_i,km由权利要求2中的标定方法确定。

S13：通过最小化标定模体中不同双基材料厚度组合与其双能谱多色投影值之间的距离，标定双能谱CT系统的双能谱多色正投影方程。

S2：利用双能谱CT系统扫描被测物体，采集被测物体的高、低能谱多色投影数据；

S3：根据所述双能谱多色正投影方程，构造迭代格式，利用所测得的高、低能谱多色投影数据，迭代重建被测物体的双基材料密度图像。

其中，步骤S3中的迭代格式可以是每个角度修正一次，也可以是用有序子集方法进行修正。

在具体实现时，步骤S3可以包括以下步骤：

S31：根据所述双能谱多色正投影方程，构造求双基材料密度图像的迭代格式；

S32：为待重建的被测物体的双基材料密度图像赋初始值；

S33：利用所述双能谱多色正投影方程，获取当前双基材料密度图像的多色投影的估计值，并与被测物体实测的多色投影的真实值相比较，得到所重建的双基材料密度图像估计值的误差；

其中，步骤S33中由双基材料密度图像估计值的误差求基材料密度图像修正值的公式为：

$\left\{\begin{array}{c}{e}_{F_L}^{\left(n\right)}=\frac{\frac{\&PartialD;C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\&PartialD;F_L}}{{\left(\frac{\&PartialD;C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\&PartialD;F_L}\right)}^2+{\left(\frac{\&PartialD;C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\&PartialD;G_L}\right)}^2}(q_{i,L}-q_{i,L}^{\left(n\right)}),\\ e_{G_L}^{\left(n\right)}=\frac{\frac{\&PartialD;C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\&PartialD;G_L}}{{\left(\frac{\&PartialD;C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\&PartialD;F_L}\right)}^2+{\left(\frac{\&PartialD;C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\&PartialD;G_L}\right)}^2}(q_{i,L}-q_{i,L}^{\left(n\right)}),\end{array}\right.i=\mathrm{1,2}.$

其中

为第n次迭代中多色投影的估计值，与测量得到的多色投影真实值q_i,L比较得到多色投影误差，沿着当前射线路径L直接分配给两种基材料，分别得到对两种基材料图像的修正值

和

。

S34：利用计算得到的双基材料密度图像估计值的误差，更新当前重建的双基材料密度图像；

其中，步骤S34中利用第n次迭代计算得到的双基材料密度图像估计值的误差，更新当前重建的双基材料密度图像的公式为：

其中参数N是所重建的基材料密度图像的像素总数，

和

分别是两幅基材料密度图像的第j个像素值，

和

是更新后的像素值，r_L,j表示基材料密度图像的第j个像素对沿着射线L的投影的贡献权值，参数λ是松弛因子，一旦一个投影Q_φ中所有射线对应的图像修正值都得到后，同时修正双基材料密度图像，获得更新后的像素值。

S35：加入物理约束条件，使当前重建的双基材料密度图像满足所述物理约束条件；

其中，这里的物理约束条件包括像素值非负约束和/或分段光滑约束。

S36：重复步骤S33-S35，直到满足迭代结束条件，重建出双基材料密度图像。

此外，为提高重建的速度，图1实施例的被测物体的双基材料密度图像重建过程还可以采用并行方法计算，并基于硬件并行计算平台加速实现。

以下通过一个具体实施例来说明本发明的双能谱CT图像迭代重建方法的具体实现过程。首先生成双能谱CT的多色投影数据。图2为本实施例中所用的测试模体，该测试模体是一个标准的胸腔模型，尺寸为400mm×200mm。选定水和骨组织为双基材料，水的密度为1.0g/cm³，骨组织的密度为1.92g/cm³。图3为本实施例中所用的标定模体。该模体包含一个较大的圆盘和一个较小的圆盘，其中较大圆盘的直径为240mm，较小圆盘的直径为60mm。大圆盘中填充了水，而小圆盘为骨组织。双能谱CT系统发出高、低X射线能谱如图4所示。这两个能谱是X光管的管电压分别为80kV和140kV（1mm铜滤波片）下发出的。在数值实现时，能谱的采样间隔设为1keV。在双能谱CT扫描中，采用扇束扫描。系统的扫描参数为：射线源到探测器的距离为1300mm，射线源到转台中心的距离为1000mm，探测器由512个1.345mm的探测单元组成。射线源旋转一周，交错均匀采集720个角度的高能和低能投影数据各一组，即相邻的高能投影和低能投影有0.25度的夹角，这样得到了几何上不匹配的高、低能多色投影数据。重建图像大小为512×512。在80kV和140kV管电压下分别采集到的测试模体的多色投影数据如图5所示，其中图5a为140kV管电压下采集的投影数据；图5b为80kV管电压下采集的投影数据。

图6为采用本发明方法实施例图像初值置0迭代12次后重建的结果图像,其中图6a是骨基材料密度图像，图6b是水基材料密度图像，图6c是关于能量为70keV的光子的线性衰减系数图像。这里线性衰减系数图像是通过将重建的水基图像与水对指定能量的X射线光子质量衰减系数的乘积以及重建的骨基图像与骨对指定能量的X射线光子质量衰减系数的乘积求和得到的。为了加速迭代收敛，可以以牛顿-拉普森方法的重建图像作为初始值，则3次迭代即可重建出和图6相同效果的图像。图7为采用传统图像域分解法重建的结果图像，用于和本发明方法重建的结果图像作比较。其中7a是用该方法重建的骨基材料密度图像，图7b是水基材料密度图像，图7c是用骨基图像水基图像合成的能量为70keV的光子的线性衰减系数图像。为了凸显所重建图像的特征，图6和图7中用了窄显示灰度窗，其中图6a和图7a的显示灰度窗为[0.61.1]，图6b和图7b的显示灰度窗为[0.91.1]，图6(c)和图7c的显示灰度窗为[-100HU100HU]。

比较图6和图7，可以看出本发明实施例的方法可以很好地分解出双基材料密度图像，基材料密度图像中仅有一些轻微的采样伪影，在线性衰减系数图像中基本消除了射线硬化伪影，图像质量明显优于传统图像域分解法所重建的图像。

图8表示沿着图6和图7所示的单能量光子线性衰减系数图像的中心列线所取得的像素值的曲线，与从测试模体图像中心列线所取得的像素真实值之间的比较。从图8中可以看到本方法的重建结果非常接近于真实值，精度远高于传统图像域分解法重建的结果图像。

进而用归一化均方距离（NMSB）和归一化平均绝对距离（NMAD）作量化比较。传统图像域分解法重建的单能量光子线性衰减系数图像的NMSB和NMAD值分别为0.205690和0.096764，而本发明方法重建图像的相应NMSB和NMAD值分别为0.030976和0.014484，由此可知本发明方法的重建的图像误差远小传统图像域分解法所重建的图像误差，进一步验证了本发明方法的优点。

图9显示了在所测多色投影数据中有噪声的情况下，本发明方法的重建效果。在多色投影数据中所添加的噪声为对应每条射线1,000,000个光子的泊松噪声。同样图9a是用本方法重建的骨基材料密度图像，图9b是水基材料密度图像，图9c是线性衰减系数图像。从图9中可以看到本发明方法即使在所测投影数据被噪声污染情况下，仍然能够重建出高质量的基材料图像和线性衰减系数图像。

本实施例说明本发明方法适用于投影数据几何完全不一致的情况，不需要预先标定能谱信息，且收敛速度快，抗噪性能强，通过少量迭代，便可以重建出高质量的基材料图像和线性衰减系数图像。

综上所述，与现有技术相比，本发明方法具有以下有益效果：

(1)无需预先测定双能CT系统的高、低能谱以及基材料在不同X射线能量下的衰减系数；

(2)适用性强，同时适用于高、低能谱几何参数一致和几何参数不一致扫描的情况，以及适用于常用的各种扫描模式，还可以推广应用于多能谱CT和慢电压切换CT系统；

(3)算法并行度高，适用于图形处理器(GPU)等并行硬件加速实现；

(4)图像重建精度高，仅通过少量迭代，即可重建出高质量的双基材料密度图像，并能有效去除所合成的衰减系数图像中的硬化伪影。

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

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

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

Claims (10)

1.一种双能谱CT图像迭代重建方法，其特征在于，包括以下步骤：

S1：根据多色投影值与基材料厚度组合之间的曲线关系获取双能谱多色正投影方程；

S2：利用双能谱CT系统扫描被测物体，采集被测物体的高、低能谱多色投影数据；

S3：根据所述双能谱多色正投影方程，构造迭代格式，利用所测得的高、低能谱多色投影数据，迭代重建被测物体的双基材料密度图像。

2.根据权利要求1所述的双能谱CT图像迭代重建方法，其特征在于，其中所述步骤S1包括：

S11：根据被测物体的尺寸和材料，选定一个包含两种不同的基材料的标定模体；

S12：利用双能谱CT系统扫描标定模体，采集标定模体在各个角度下的双能谱多色投影数据；

S13：通过最小化标定模体中不同双基材料厚度组合与其双能谱多色投影值之间的距离，标定双能谱CT系统的双能谱多色正投影方程。

3.根据权利要求1所述的双能谱CT图像迭代重建方法，其特征在于，其中所述步骤S3包括：

S31：根据所述双能谱多色正投影方程，构造求双基材料密度图像的迭代格式；

S32：为待重建的被测物体的双基材料密度图像赋初始值；

S33：利用所述双能谱多色正投影方程，获取当前双基材料密度图像的多色投影的估计值，并与被测物体实测的多色投影的真实值相比较，得到所重建的双基材料密度图像估计值的误差；

S34：利用计算得到的双基材料密度图像估计值的误差，更新当前重建的双基材料密度图像；

S35：加入物理约束条件，使当前重建的双基材料密度图像满足所述物理约束条件；

S36：重复步骤S33-S35，直到满足迭代结束条件，重建出双基材料密度图像。

4.根据权利要求1所述的双能谱CT图像迭代重建方法，其特征在于，所述双能谱CT系统的扫描模式包括圆周扫描、螺旋扫描、有限角度扫描、稀疏角度扫描及一般扫描模式中的至少一种。

5.根据权利要求1所述的双能谱CT图像迭代重建方法，其特征在于，步骤S3中所述迭代格式是每个角度修正一次或是用有序子集方法修正。

6.根据权利要求3所述的双能谱CT图像迭代重建方法，其特征在于，所述物理约束条件包括像素值非负约束和/或分段光滑约束。

7.根据权利要求2所述的双能谱CT图像迭代重建方法，其特征在于，步骤S12中所述双能谱多色正投影方程用多项式方程、分段多项式方程或样条曲线方程表示；在用多项式方程表示多色正投影方程的情况下：

$q_{i,L}=C_i(F_L,G_L)=\underset{k=0}{\overset{K}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{M}{\mathrm{\&Sigma;}}}{c}_{i,\mathrm{km}}{F}_L^k{G}_L^m,$ i=1,2

其中参数i表示能谱索引；q_i,L是沿着给定X射线路径L所探测到的多色投影，表示为基函数

的一个线性组合；K和M分别是多项式的阶数，即多项式共包含(K+1)×(M+1)个基函数

函数C_i的多项式系数c_i,km由权利要求2中的标定方法确定。

8.根据权利要求3所述的双能谱CT图像迭代重建方法，其特征在于，步骤S33中由双基材料密度图像估计值的误差求基材料密度图像修正值的公式为：

$\left\{\begin{array}{c}{e}_{F_L}^{\left(n\right)}=\frac{\frac{\&PartialD;C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\&PartialD;F_L}}{{\left(\frac{\&PartialD;C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\&PartialD;F_L}\right)}^2+{\left(\frac{\&PartialD;C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\&PartialD;G_L}\right)}^2}(q_{i,L}-q_{i,L}^{\left(n\right)}),\\ e_{G_L}^{\left(n\right)}=\frac{\frac{\&PartialD;C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\&PartialD;G_L}}{{\left(\frac{\&PartialD;C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\&PartialD;F_L}\right)}^2+{\left(\frac{\&PartialD;C_i(F_L^{\left(n\right)},G_L^{\left(n\right)})}{\&PartialD;G_L}\right)}^2}(q_{i,L}-q_{i,L}^{\left(n\right)}),\end{array}\right.i=\mathrm{1,2}.$

其中

为第n次迭代中多色投影的估计值，与测量得到的多色投影真实值q_i,L比较得到多色投影误差，沿着当前射线路径L直接分配给两种基材料，分别得到对两种基材料图像的修正值

和

9.根据权利要求3所述的双能谱CT图像迭代重建方法，其特征在于，步骤S34中利用第n次迭代计算得到的双基材料密度图像估计值的误差，更新当前重建的双基材料密度图像的公式为：

其中参数N是所重建的基材料密度图像的像素总数，

和

分别是两幅基材料密度图像的第j个像素值，

和

是更新后的像素值，r_L,j表示基材料密度图像的第j个像素对沿着射线L的投影的贡献权值，参数λ是松弛因子，一旦一个投影Q_φ中所有射线对应的图像修正值都得到后，同时修正双基材料密度图像，获得更新后的像素值。

10.根据权利要求1所述的双能谱CT图像迭代重建方法，其特征在于，被测物体的双基材料密度图像重建过程采用并行方法计算，并基于硬件并行计算平台加速实现。

Images