CN111968060A - 一种基于倾斜投影修正技术的多能谱ct快速迭代重建方法 - Google Patents

Abstract

本发明公开了一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法，该方法包括：步骤1，获取被测物体在多个能谱下的多色投影数据；步骤2，获取各种基材料密度图像估计值；步骤3，对基材料密度图像估计值进行正投影，得到各种基材料的投影估计值及被测物体的多色投影估计值；步骤4，计算多色投影估计值与多色投影数据之间的误差，利用倾斜投影技术对基材料投影进行修正，得到各种基材料的投影残差；步骤5，将投影残差进行反投影，得到各种基材料的残差图像，更新被测物体的基材料密度图像估计值；步骤6，判断终止条件是否满足，若满足则终止迭代，否则转向步骤3。本发明能够由采集到的多个能谱的多色投影数据快速重建出被测物体的多种基材料密度图像。

Description

一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法

技术领域

本发明涉及X射线CT成像技术领域，特别是关于一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法。

背景技术

X射线计算机断层成像技术(X ray Computed Tomography,简称X射线CT)可以在不破坏或损伤物体的情况下呈现物体的内部细节，已广泛应用于医学、生物、工业、材料、古化石和航天等众多领域。传统CT成像理论假设X射线由单一能量的光子组成，忽略了X射线的多色性，因此利用传统单能量CT重建算法重建实采数据时会产生射束硬化伪影，如杯状伪影和条状伪影，严重影响成像质量。

多能谱CT成像系统利用多个不同能谱下的X射线扫描被测物体，获得被测物体在多个不同能谱下的投影数据。利用这些投影数据可以重建被测物体的基材料密度图像，或者是等效原子序数和电子密度图像。与单能谱CT相比，多能谱CT获取了更多的被测物体的信息，具有更好的物质区分能力，在硬化伪影去除、骨密度测量、PET衰减校正和伪单能图像计算等方面有广泛的应用前景。

现有的多能谱CT多色投影数据获取方式大致有两类，第一类方式是利用X射线源获取两个或多个X射线能谱，分别对被测物体进行扫描，获取两组或多组多色投影数据。代表性的技术有“全扫描”模式、双源双探扫描模式和快速电压切换扫描模式等。其中，“全扫描”模式是利用传统CT设备，在不同的管电压和管电流下对被测物体分别进行多次扫描，这种方法可在传统CT设备上完成，不需要额外添加硬件设备。第二类方式是只使用一个X射线能谱，使用三明治式探测器或者光子计数探测器获取多组多色投影数据。通常情况下，第一类方式获取的多组多色投影数据是射线路径几何不一致的，第二类方式获取的多组多色投影数据是射线路径几何一致的。

重建算法大致可以分为四类：图像域重建方法、投影域重建方法、基于深度学习的重建方法和迭代重建方法。图像域重建方法是对采集到的多组多色投影数据分别用传统CT重建算法进行重建，然后将重建图像进行线性组合得到基材料的密度图像。该方法重建的图像被认为是对真实图像的一阶近似，无法准确描述基材料密度图像与多色投影之间的非线性关系，通常在重建图像中还会存在硬化伪影，且受噪声的影响较大。为了提高图像域重建算法的图像质量，研究者提出了一些基于传统低通滤波器或者基于统计先验滤波器的改进算法，这些算法可以在一定程度上抑制噪声对重建结果的影响，然而对于图像分解精度的改进有限。投影域的重建方法首先对多组多色投影进行组合，得到多种基材料密度图像的投影数据，然后对多种基材料密度图像的投影数据分别利用传统CT重建方法进行重建，得到基材料密度图像。通常，利用投影域重建方法重建的图像要优于利用图像域重建方法重建的图像。但是该方法要求不同能谱下采集的投影数据是几何一致的，即要求每条X射线路径下的所有能谱的投影数据都要采集到。然而实际的多能谱CT系统采集的数据，并不能保证这一要求，比如双源双探的能谱CT扫描模式。近年来，利用深度学习技术，研究者在有完备训练集的情况下，在多能谱CT重建领域得到了高质量的图像重建结果。然而很多情况下，本实施例无法获得充足的训练样本。考虑到多能谱CT问题的非线性性，理论上迭代重建算法更适合这类问题求解，其利用数值方法或者优化方法构造迭代结构，通过对图像重建结果逐步修正，可以得到高精度的被测物体的各基材料密度图像信息。由于能谱CT重建问题的非线性性以及病态性，使得现有的多能谱CT重建算法无法快速重建出高质量的基材料密度图像。

发明内容

本发明的目的在于提供一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法，其能够由采集到的多个能谱的多色投影数据快速重建出被测物体的多种基材料密度图像。

为实现上述目的，本发明提供一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法，该方法包括：步骤1，利用多能谱CT系统扫描被测物体，获得被测物体在多个能谱下的真实多色投影数据；步骤2，为被测物体的各种基材料密度图像赋初值，作为各种基材料密度图像估计值；步骤3，对各种基材料密度图像估计值进行正投影，得到各种基材料的投影估计值，进而利用X射线能谱信息和基材料的质量衰减系数信息得到被测物体在各个能谱下的多色投影估计值；步骤4，计算多色投影估计值与真实多色投影数据之间的误差，并利用倾斜投影技术对基材料投影进行修正，得到各种基材料的投影残差A；步骤5，将各种基材料的投影残差进行反投影操作，得到各种基材料的残差图像，进而更新被测物体的各种基材料密度图像估计值；步骤6，判断终止条件是否满足，若满足则终止迭代，否则转向步骤3。

本发明由于采取以上技术方案，其具有以下优点：。

本发明方法用于多能谱CT快速迭代重建，适用于多种常用的多能谱CT扫描模式。与现有方法相比，本发明方法能够在保证多能谱CT图像重建质量的同时，显著加快重建图像的收敛速度。

附图说明

图1为本发明提出的一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法的流程图。

图2a为本发明一个实施例用作测试模体的FORBILD胸腔模型图像。

图2b为本发明一个实施例用作测试模体的水基材料密度图像。

图2c为本发明一个实施例用作测试模体的骨基材料密度图像。

图3为本发明一个实施例所用的X射线源分别在140kV和80kV管电压下发出的高能谱和低能谱示意图。

图4a为本发明一个实施例的140kV能谱下采集到的测试模体添加初始光子数为10⁶的泊松噪声的多色投影数据图像。

图4b为本发明一个实施例的80kV能谱下采集到的测试模体添加初始光子数为10⁶的泊松噪声的多色投影数据图像。

图5a为一个实施例采用本发明方法的骨基材料密度图像重建的结果图像。

图5b为一个实施例采用本发明方法的水基材料密度图像重建的结果图像。

图5c为本发明一个实施例所重建的关于能量为70keV的光子的线性衰减系数图像。

图6为一个实施例，用以说明本发明方法与已公开的相关专利文献方法的不同。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

参阅图1，为本发明实施例提供的基于倾斜投影修正技术的多能谱CT快速迭代重建方法包括以下步骤：

步骤1，利用多能谱CT系统扫描被测物体，获得被测物体在N个能谱下的N个真实多色投影数据，其中，一个能谱下相应获得一个多色投影数据。

步骤2，为被测物体的各种基材料密度图像赋初值，作为各种基材料密度图像估计值。

步骤3，对各种基材料密度图像估计值进行正投影，得到各种基材料的投影估计值，进而利用X射线能谱信息和基材料的质量衰减系数信息得到被测物体在各个能谱下的多色投影估计值。

步骤4，计算多色投影估计值与真实多色投影数据之间的误差，并利用倾斜投影技术对基材料投影进行修正，得到各种基材料的投影残差。

步骤5，将各种基材料的投影残差进行反投影操作，得到各种基材料的残差图像，进而更新被测物体的各种基材料密度图像估计值。

步骤6，重复步骤3至步骤5，直到满足终止条件。

在一个实施例中，步骤1中的“多能谱CT系统”大致有两类。

第一类方式是利用X射线源获取两个或多个X射线能谱，分别对被测物体进行扫描，获取两组或多组多色投影数据。代表性的技术有“全扫描”模式、双源双探扫描模式，快速电压切换扫描模式等。其中，“全扫描”模式是利用传统CT设备，在不同的管电压和管电流下对被测物体分别进行多次扫描，这种方法可在传统CT设备上完成，不需要额外添加硬件设备。

第二类方式是只使用一个X射线能谱，使用三明治式探测器或者光子计数探测器获取多组多色投影数据。通常情况下，第一类方式获取的多组多色投影数据是射线路径几何不一致的，第二类方式获取的多组多色投影数据是射线路径几何一致的。

在一个实施例中，步骤1中，根据多能谱CT系统中所用探测器不同，得到的X射线多色投影数据维度也不同。

当探测器为包含W个探测器单元的线阵性探测器时，在圆轨迹扇束扫描模式下，在V个扫描角度共采集到W×V个X射线投影数据。其中，每个角度的投影包括W个X射线投影数据，将投影数据按照角度顺序进行排序，组合得到一幅二维X射线多色投影图。

当探测器为包含W×H个探测器单元的面阵性探测器时，在圆轨迹扇束扫描模式下，在V个扫描角度共采集到W×H×V个X射线投影数据。其中，每个角度的投影包括W×H个X射线投影数据，将投影数据按照角度顺序进行排序，组合得到一幅三维X射线多色投影图。

在一个实施例中，步骤2中的“为被测物体的各种基材料密度图像赋初值”是以0作为初始图像，即图像中的每个像素均为0。

在一个实施例中，以

表示第n次迭代得到的第m种基材料的投影估计值、

表示第n次迭代得到的第m种基材料密度图像估计值、M表示待分解基材料的个数(被测物体包含M种基材料)，则步骤3所述的对各基材料密度图像估计值进行正投影，得到各种基材料的投影估计值的公式为：

式中，l表示一条X射线扫描路径；R_l表示沿射线路径l的正投影算子。

在一个实施例中，以p_k，l，k＝1，2，…，N表示扫描值得到的第k个能谱对应的多色投影估计值，N表示能谱的个数，假定基材料数目M＝N，则步骤3所述的由各种基材料的投影估计值得到各个能谱下的多色投影估计值p_k，l的公式为：

式中，S_k，m为第k个归一化能量谱的离散形式；θ_i，m为第i种基材料的质量衰减系数的离散形式，m＝1，2，…，M_k；ξ_k表示对应于第k个等效能谱的X射线扫描路径的集合；δ表示X射线能谱和物质的质量衰减系数的采样间距，通常为1kev；f_i为第i种基材料的密度图像的估计值。

通常，若多能谱CT系统获取的多色投影数据是几何一致的，那么在ξ₁∩ξ₂∩...∩ξ_N射线路径下的多色投影数据用于多能谱CT重建是完备的。若多能谱CT系统获取的多色投影数据是几何不一致的，那么在ξ₁∩ξ₂∩...∩ξ_N射线路径下的多色投影数据用于多能谱CT重建是不完备的。l表示ξ_k中的一条X射线扫描路径。

在一个实施例中，以多色投影射线路径几何不一致为例，说明步骤4中的“利用倾斜投影技术对步骤3中的基材料的投影估计值

进行修正，得到各种基材料的投影残差”的方法：

步骤4.1，将式(2)表示的各个能谱下的多色投影估计值p_k，l对应的非线性投影方程线性化：

假设经过n次迭代后，已有基材料密度图像估计值(f₁ ⁽ⁿ⁾，f₂ ⁽ⁿ⁾，…，f_N ⁽ⁿ⁾)，N≥M，将式(2)在(f₁ ⁽ⁿ⁾，f₂ ⁽ⁿ⁾，…，f_N ⁽ⁿ⁾)处做一阶泰勒展开，得到线性化的投影方程(3)：

其中，

式中，

和

是用于简化方程的中间参数，无实质性含义；f_i ⁽ⁿ⁾表示第i种基材料密度图像第n次迭代后的结果，i表示基材料的序号；

表示第n次迭代得到的第k个能谱对应的多色投影估计值；m表示将第k个等效能谱的有效能量范围等分为M_k个小能量区间，m＝1，2，…，M_k，M_k的数值X射线能谱的最大能量有关，每个小能量区间的长度为δ，δ表示X射线能谱和物质的质量衰减系数的采样间距，通常为1kev，在每个小能量区间内对S_k(E)，k＝1，2，…，N和θ_i(E)，i＝1，2，…，N进行采样，于是S_k(E)和θ_i(E)被离散为如下向量形式：

其中，S_k，m为S_k(E)在第m个子能量区间内的采样值；θ_i，m为θ_i(E)θ_i，m在第m个子能量区间内的采样值。

在几何不一致的情况下，在某些X射线路径下只已知一个多色投影数据，假定当前已知的第1个能谱对应的多色投影估计值p_1，l是在能谱S₁(E)，X射线路径属于ξ₁的情形下得到的，那么此时在能谱为S_k(E)，k＝2，3，…，N，X射线路径属于ξ₁下的p_k，l，k＝2，3，…，N是未知的。在当前给定射线路径下，投影方程(3)转化为线性多色投影方程组(6.1)至(6.N)：

其中：

i表示第i个x射线能谱对应的数据，i＝1…N；j表示第j种基材料密度图像对应的数据，j＝1…N；f_n表示待求的第n种基材料密度图像；x_n表示待求的第n种基材料的投影估计值，即步骤3中基材料到投影估计值的下一次迭代点；p_1，l表示第1个能谱对应的多色投影估计值；

表示第n次迭代得到的第1个能谱对应的多色投影估计值；f_k ⁽ⁿ⁾与步骤2中得到的

物理意义相同；

n＝2，3，…，N表示在路径l下，第n种基材料密度图像在第n次迭代的结果的多色投影值；

都为用于简化方程的中间参数，表达形式由式(4)中的

和

类推得到。

步骤4.2，获取当前迭代点向超平面做正交投影的投影方向：利用线性多色投影方程(6.1)，计算当前迭代点(R_lf₁ ⁽ⁿ⁾，R_lf₂ ⁽ⁿ⁾，…，R_lf_N ⁽ⁿ⁾)向线性多色投影方程(6.1)对应的超平面做正交投影的投影方向dir2，其表示为式(8)：

式中，dir21表示当前迭代点(R_lf₁ ⁽ⁿ⁾，R_lf₂ ⁽ⁿ⁾，…，R_lf_N ⁽ⁿ⁾)向线性多色投影方程(6.1)对应的超平面做正交投影的投影方向；a_1，k表示线性多色投影方程(6.1)中的参数，x_k表示当前计算的第k种基材料的投影数据。

取式(8)的单位向量作为当前射线路径下已知能谱下的线性多色投影方程(6.1)对应的超平面做正交投影的投影方向dir2即如式(9)：

其中，||·||为向量取模运算符；dir2表示当前射线路径下下已知能谱对应的投影方程(6.1)所表示的超平面的法方向；dir21表示式(8)的结果。

步骤4.3，计算由线性多色投影方程组(6.2)～(6.N)的法向量在N维空间中张成的超平面的法方向：

由线性多色投影方程(6.2)～(6.N)可知，向量(a_i1，a_i2，…，a_iN)，i＝2，3，…，N分别为方程(6.2)～(6.N)所表示的超平面的法方向。由(a_i1，a_i2，…，a_iN)，i＝2，3，…，n…，N这N-1个向量在N维空间中所张成的超平面的法方向表示为式(10)：

其中：|·|为行列式；e₁，e₂，…，e_N为N维空间下的一组标准正交基。那么，e_n的系数D_n表示为：

由线性多色投影方程(6.2)～(6.N)的法向量在N维空间中张成的超平面的法方向可以表示为式(11)：

dir11＝(D₁，D₂，…，D_n) or dir12＝-dir11 (11)

dir11表和dir12均表示由线性多色投影方程(6.2)～(6.N)的法向量在N维空间中张成的超平面的法方向，dir11和dir12为求解过程中的中间变量，最终使用的是dir1。

取式(11)所表示的两个方向中与dir2成锐角的方向，作为由线性多色投影方程(6.2)～(6.N)的法向量在N维空间中张成的超平面的法方向dir1，即式(12)：

其中，<，>表示两个向量的内积。

步骤4.4，利用式(13)，将dir1和dir2两个方向加权求和，得到从当前迭代点(R_lf₁ ⁽ⁿ⁾，R_lf₂ ⁽ⁿ⁾，…，R_lf_N ⁽ⁿ⁾)向当前射线路径下已知的多色投影方程做倾斜投影的方向dir，令

dir＝λ₁dir1+λ₂dir2 (13)

为从当前迭代点(R_lf₁ ⁽ⁿ⁾，R_lf₂ ⁽ⁿ⁾，…，R_lf_N ⁽ⁿ⁾)向当前射线路径下已知的多色投影方程做倾斜投影的方向，其中λ₁和λ₂为控制两个方向权重的系数。适当增大λ₁/λ₂的值，可以加快本发明实施例方法的收敛速度。由线性多色投影方程(6.1)～(6.N)组成的线性方程组的条件数越大，通常多能谱CT求解问题病态性越高，此时应减少λ₁/λ₂的值，使求解过程不会过度放大噪声。通常情况下，取λ₁/λ₂∈[0.5，1.5]，可在保证基材料密度图像重建质量的同时，显著加快基材料密度图像收敛速度。

沿该方向dir在当前射线路径下已知的多色投影方程上的投影点(R_lf₁ ⁽ⁿ⁺¹⁾...R_lf_n ⁽ⁿ⁺¹⁾...R_lf_N ⁽ⁿ⁺¹⁾)^T表示为方程(14)：

(R_lf₁ ⁽ⁿ⁺¹⁾ ... R_lf_n ⁽ⁿ⁺¹⁾ ... R_lf_N ⁽ⁿ⁺¹⁾)^T＝(R_lf₁ ⁽ⁿ⁾ ... R_lf_n ⁽ⁿ⁾ ... R_lf_N ⁽ⁿ⁾)^T+A(14)

式中，A为各种基材料的投影残差，即当前迭代点(R_lf₁ ⁽ⁿ⁾，R_lf₂ ⁽ⁿ⁾，…，R_lf_N ⁽ⁿ⁾)与该点沿着倾斜投影方向dir在已知的线性多色投影方程(6.1)上的投影点之间的差。

在一个实施例中，步骤5中，根据式(16)获得被测物体的各基材料密度图像的估计值：

(f₁ ⁽ⁿ⁺¹⁾ ... f_n ⁽ⁿ⁺¹⁾ ... f_N ⁽ⁿ⁺¹⁾)^T＝(f₁ ⁽ⁿ⁾ ... f_n ⁽ⁿ⁾ ... f_N ⁽ⁿ⁾)^T+R_l ^-1A (16)

式中，T表示转置。

得到本发明实施例更新后的基材料密度图像估计值f₁ ⁽ⁿ⁺¹⁾，…，f_n ⁽ⁿ⁺¹⁾，…，f_N ⁽ⁿ⁺¹⁾，完成一次迭代。重复步骤4.1-步骤4.4以及步骤5，直到遍历完所有能量的多色投影数据。

需要说明的是，在多能谱CT重建问题中，通常多色投影射线路径几何不一致时有效的方法，可以直接推广到几何一致时使用，所以，几何一致下的步骤和几何不一致下的步骤是一样的。

在一个实施例中，步骤6所述的终止条件，包括但不限于最大迭代次数，得到的被测物体的各种基材料密度图像估计值是否符合主观判断或所需要求，两次迭代得到的各种基材料密度图像估计值的差小于设定的阈值。

在一个实施例中，无论是步骤3所述将各种基材料密度图像估计值进行正投影得到各种基材料的投影估计值，还是步骤4所述由多色投影估计值与真实多色投影数据之间的误差利用倾斜投影修正技术得到各种基材料的投影残差，或是步骤5所述将各种基材料的投影残差进行反投影操作得到各种基材料的图像残差，均可采用并行方法计算，并基于硬件并行计算平台加速实现。

以下通过一个具体实施例来说明本发明提出的一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法的具体实现过程。

图2a为本具体实施例中用作测试模体的FORBILD胸腔模型图像。图2b为本具体实施例中用作测试模体的水基材料密度图像。图2c为本具体实施例中用作测试模体的骨基材料密度图像。在本实施例中使用的高低能X射线能谱用开源软件使用开源软件SpectrumGUI，仿真了GE Maxiray 125球管在80kVp和140kVp下的能谱，并对能谱做了归一化处理。实验管电压分别为80kV和140kV。在模拟实验中，模体假定为不同密度的水和骨的组合。水和骨也是重建图像的基材料。材料的质量衰减系数从美国国家标准技术研究院(NIST)网站获得，并根据能谱的取值对其做了相应的插值。在模拟多色投影时，X射线能谱和材料的质量衰减系数的采样间隔均为1kV。

实验扫描参数设置如下：射线源到转台中心的距离为490mm，射线源到探测器的距离为880mm。线探测器由512个探测器单元组成，每个探测器单元的尺寸为0.2mm。

根据方程(2)模拟多色投影，某一射线路径下X射线与基材料密度图像的交线长由光线投射算法计算。无噪声模拟实验和含噪声模拟实验均以快速电压切换的形式采集得到，即偶数多色投影是在某一种电压下模拟扫描得到，奇数多色投影是在另一种电压下模拟扫描得到，因此得到的投影为几何不一致的投影。两个电压下的多色投影，均在CT系统旋转一周的情况下采集了720个角度的投影数据。对投影数据添加初始光子数为10⁶的泊松噪声，采集到的测试模体添加泊松噪声的高、低能多色投影数据如图4a和图4b所示。

具体实施步骤如下：

步骤1：为基材料密度图像f₁ ⁽⁰⁾，…，f_n ⁽⁰⁾，…，f_N ⁽⁰⁾赋初值；

步骤2：假定经过n次迭代后，有基材料密度图像估计值f₁ ⁽ⁿ⁾，…，f_n ⁽ⁿ⁾，…，f_N ⁽ⁿ⁾。对给定的一条X射线路径，设其对应的多色投影数据为p_1，l，根据公式(4)和公式(7)计算a_ij，(i＝1…N，j＝1…N)和

步骤3：首先根据公式(8)和(9)计算出dir2；然后根据公式(10)、(11)和(12)计算出dir1；最后根据公式(13)计算出本文算法的倾斜投影方向dir；

步骤4：根据公式(14)计算当前路径下更新后的基材料密度图像投影R_lf₁ ⁽ⁿ⁺¹⁾，…，R_lf_n ⁽ⁿ⁺¹⁾，…，R_lf_N ⁽ⁿ⁺¹⁾，进而利用公式(15)得到更新的基材料密度图像f₁ ⁽ⁿ⁺¹⁾，…，f_n ⁽ⁿ⁺¹⁾，…，f_N ⁽ⁿ⁺¹⁾；

步骤5：判断终止条件是否满足，若满足则终止迭代，否则转向步骤2。

选取参数λ₁＝1，λ₂＝1，进行10次迭代得到的基材料密度图像如图5所示，其中图5a是骨基材料密度图像(灰度窗为[0，1.92])，图5b是水基材料密度图像(灰度窗为[0，1.18])，图5c是关于能量为70keV时的线性衰减系数图像(灰度窗为[0，0.03])。这里线性衰减系数图像是通过将重建的水基图像与水对指定能量的X射线光子质量衰减系数的乘积以及重建的骨基图像与骨对指定能量的X射线光子质量衰减系数的乘积求和得到的。由图5结果可看出，本发明方法可以由采集到的多个能量谱的多色投影数据重建出被测物体的多种基材料密度图像。与现有方法相比，本发明方法具有重建图像质量高、收敛速度快的优点。

以下通过一个具体实施例来说明本发明方法与已公开的相关专利文献的不同。

使用开源软件Spectrum GUI生成，仿真了GE Maxiray 125球管140kVp下的能谱，并对能谱做了归一化处理。将140kV下的能谱分为1-70kV和71-140kV两段，能谱的采样间隔δ为1kV，并分别归一化，分别作为本实施例中使用的X射线能谱

和

选用的基材料分别为水和骨，线性衰减系数从美国国家标准技术研究院(NIST)网站获得，并根据能谱的取值对其做了相应的插值。

根据(2)式，本实施例构建在一条特定射线路径下的二元非线性方程组

其中，

和

分别为水和骨的质量衰减系数；p_1，l和p_2，l分别是在S₁(E)和S₂(E)能谱下获得的多色投影值，其对应的水基材料投影R_lf₁＝4，骨基材料投影R_lf₂＝1。本实施例假定S₁(E)、S₂(E)、θ₁(E)和θ₂(E)已知，R_lf₁和R_lf₂为待求值，那么式上述方程组构成了关于未知数(R_lf₁，R_lf₂)的二元非线性方程组。实验中本实施例E-ART方法法和本专利方法求解该方程组，每次迭代只使用其中一个方程，假定另一个方程的多色投影值未知。初值设为R₁f₁ ⁽⁰⁾＝0，R_lf₂ ⁽⁰⁾＝0，所有方法均迭代50次。

图6为用E-ART方法和本发明方法求解该方程组时的迭代轨迹图，黑色实线和黑色虚线分别为S₁(E)和S₂(E)能谱下的多色投影曲线。A、B和C构成的折线和D和E折线分别为E-ART方法和本发明方法求解上述问题的迭代路径。为了在图中清楚地显示本发明方法迭代路径，此实施例中本发明方法参数选取为λ₁＝1，λ₂＝1。

可以看出，本发明方法与上述已公开的相关专利文献方法的不同之处在于E-ART方法以正交投影的方式更新迭代点，迭代点收敛到真解的速度很慢，迭代50次时仍不能收敛至真解附近。而本发明方法采用倾斜投影技术，大大加快了迭代点到真解的收敛速度，20次左右即收敛到真解。

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

本领域普通技术人员可以理解：本实施例采用快速电压切换模式获得所需数据，但是本发明保护范围并不限于此扫描模式。

本领域普通技术人员可以理解：本实施例获得的两组数据是几何不一致的，但是本发明可以对几何一致的数据进行，获得相应的基材料密度图像。

本领域普通技术人员可以理解：本发明可以对三个或更多个能谱下的扫描数据进行，获得三个或更多种基材料的密度图像。

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

Claims (7)

1.一种基于倾斜投影修正技术的多能谱CT快速迭代重建方法，其特征在于，包括：

步骤1，利用多能谱CT系统扫描被测物体，获得被测物体在多个能谱下的真实多色投影数据；

步骤2，为被测物体的各种基材料密度图像赋初值，作为各种基材料密度图像估计值；

步骤3，对各种基材料密度图像估计值进行正投影，得到各种基材料的投影估计值，进而利用X射线能谱信息和基材料的质量衰减系数信息得到被测物体在各个能谱下的多色投影估计值；

步骤4，计算多色投影估计值与真实多色投影数据之间的误差，并利用倾斜投影技术对基材料投影进行修正，得到各种基材料的投影残差A；

步骤5，将各种基材料的投影残差进行反投影操作，得到各种基材料的残差图像，进而更新被测物体的各种基材料密度图像估计值；

步骤6，判断终止条件是否满足，若满足则终止迭代，否则转向步骤3。

2.如权利要求1所述的基于倾斜投影修正技术的多能谱CT快速迭代重建方法，其特征在于，步骤4中的投影残差A表示为式(15)：

式中：

N表示能谱的个数；

M表示基材料的数目，M＝N；

p_1,l表示第1个能谱对应的多色投影估计值；

表示第n次迭代得到的第1个能谱对应的多色投影估计值；

分别表示第n次迭代得到的第k、n个能谱对应的多色投影估计值；

dir表示从当前迭代点(R_lf₁ ⁽ⁿ⁾,R_lf₂ ⁽ⁿ⁾,…,R_lf_N ⁽ⁿ⁾)向当前射线路径下已知的多色投影方程做倾斜投影的方向；

dir2为当前迭代点(R_lf₁ ⁽ⁿ⁾,R_lf₂ ⁽ⁿ⁾,…,R_lf_N ⁽ⁿ⁾)向线性多色投影方程(6.1)对应的超平面做正交投影的投影方向；

dir1、dir11表和dir12表示线性多色投影方程组(6.2)～(6.N)的法向量在N维空间中张成的超平面的法方向，dir11和dir12为求解过程中的中间量；

dir21表示当前迭代点(R_lf₁ ⁽ⁿ⁾,R_lf₂ ⁽ⁿ⁾,…,R_lf_N ⁽ⁿ⁾)向线性多色投影方程(6.1)对应的超平面做正交投影的投影方向；

||·||为向量取模运算符；

λ₁表示控制dir1方向权重的系数；

λ₂表示控制dir2方向权重的系数；

<,>表示两个向量的内积；

x_n表示第n种基材料的投影数据；

表示用于简化线性多色投影方程(6.1)的中间参数；

a_ij、

表示用于简化线性多色投影方程组(6.1)～(6.N)的中间参数，i表示第i个X射线能谱对应的数据，i＝1…N；j表示第j种基材料密度图像对应的数据，j＝1…N；

ξ_k、ξ_n分别表示对应于第k、n个等效能谱的X射线扫描路径的集合；

S_k,m为S_k(E)在第m个子能量区间内的采样值，第k个归一化能量谱的离散形式；S_k(E)为表示第k个归一化的等效能谱，∫S_k(E)dE＝1；

δ表示X射线能谱和物质的质量衰减系数的采样间距；

θ_i,m为θ_i(E)在第m个子能量区间内的采样值，第i种基材料的质量衰减系数的离散形式；θ_i(E)为第i种基材料的质量衰减系数；

R_l表示沿射线路径l的正投影算子；

f_i ⁽ⁿ⁾、f_k ⁽ⁿ⁾分别表示第n次迭代得到的第i、k种基材料密度图像估计值；

m表示将第k个等效能谱的有效能量范围等分为M_k个小能量区间，m＝1,2,…,M_k；

f_n表示待求的第n种基材料密度图像；

表示在路径l下，第n种基材料密度图像在第n次迭代的结果的多色投影估计值；

都为用于简化方程的中间参数，表达形式由式(4)中的

和

类推得到。

3.如权利要求2所述的基于倾斜投影修正技术的多能谱CT快速迭代重建方法，其特征在于，步骤3中为：利用式(1)，对各种基材料密度图像估计值进行正投影，得到各种基材料的投影估计值；

式中：

表示第n次迭代得到的第m种基材料的投影估计值；

表示第n次迭代得到的第m种基材料密度图像估计值。

4.如权利要求3所述的基于倾斜投影修正技术的多能谱CT快速迭代重建方法，其特征在于，步骤3中为：利用式(2)，得到被测物体在各个能谱下的多色投影估计值；

式中：

p_k,l表示第k个能谱对应的多色投影估计值；

f_i表示待求的第i种基材料密度图像。

5.如权利要求4所述的基于倾斜投影修正技术的多能谱CT快速迭代重建方法，其特征在于，步骤4中的线性多色投影方程组(6.2)～(6.N)由式(2)在(f₁ ⁽ⁿ⁾,f₂ ⁽ⁿ⁾,…,f_N ⁽ⁿ⁾)处做一阶泰勒展开得到线性化的投影方程(3)，再由投影方程(3)转化得到：

6.如权利要求1至5中任一项所述的基于倾斜投影修正技术的多能谱CT快速迭代重建方法，其特征在于，步骤2中的“为被测物体的各种基材料密度图像赋初值”是以0作为初始图像，即图像中的每个像素均为0。

7.如权利要求6所述的基于倾斜投影修正技术的多能谱CT快速迭代重建方法，其特征在于，步骤1中，根据多能谱CT系统中所用探测器不同，得到的X射线多色投影数据维度也不同：

当探测器为包含W个探测器单元的线阵性探测器时，在圆轨迹扇束扫描模式下，在V个扫描角度共采集到W×V个X射线投影数据。其中，每个角度的投影包括W个X射线投影数据，将投影数据按照角度顺序进行排序，组合得到一幅二维X射线多色投影图。

当探测器为包含W×H个探测器单元的面阵性探测器时，在圆轨迹扇束扫描模式下，在V个扫描角度共采集到W×H×V个X射线投影数据。其中，每个角度的投影包括W×H个X射线投影数据，将投影数据按照角度顺序进行排序，组合得到一幅三维X射线多色投影图。

Images