CN110533734A - 基于传统单能ct的多能谱分段稀疏扫描迭代重建方法 - Google Patents

Abstract

基于传统单能CT的多能谱分段稀疏扫描迭代重建方法，包括的步骤有：步骤一、在传统单能CT上采用分段稀疏的扫描模式，应用多个能量的X射线分别在相应的多个有限角度范围内扫描被测物体，获得被测物体在多个能量下的多个分段稀疏投影数据。步骤二、将步骤一获得的多个分段稀疏投影数据，应用联合全变分正则化与鲁棒性主成分分析约束的图像迭代重建方法，分别重建被测物体在多个能量下的CT图像。本发明利用不同能量下被测物体结构的一致性，设计分段稀疏扫描模式，大大降低对球管电压切换频率的技术要求。同时本发明为解决投影数据的不完备性，提出联合全变分正则化与鲁棒性主成分分析约束的图像迭代重建方法，有效地去除噪声与伪影，得到高质量的重建结果。

Description

基于传统单能CT的多能谱分段稀疏扫描迭代重建方法

技术领域

本发明涉及X射线CT成像领域，特别涉及基于传统单能CT的多能谱分段稀疏扫描迭代重建方法。

背景技术

计算机断层成像技术(computed tomography，简称CT)可以重建被扫描物体的断层结构信息，已被广泛应用于医学诊断、放射治疗、工业检测、航空航天和考古等领域。然而目前常用的X射线球管所发出的X射线并不是理想的单色X线束，而是服从一定谱分布的多色X线束，因此会导致射束硬化与低组织对比度等问题。此外，不同成份构成的不同组织可由同一CT值表示，导致了不同组织类型的分类变得困难。为解决该问题，Alvarez和Macovski提出了双能谱CT，利用两个能谱衰减系数的线性组合可以很好地对组织进行分类。

近十年，双能谱CT已被证实是一种可以实现物质分解、化学组成分析、伪单色谱成像等等功能的医学成像方法。同理，多能谱CT在双能谱CT的基础上增加了更多的能谱信息，能对物质进行更深入的分析，因此具有广泛的应用前景。多能谱CT由于采用多个能谱信息，想要获取不同的能谱信息就需要采集多个不同能量下的投影数据。然而若直接使用现有四种主流的能谱CT扫描模式(单源双次扫描，双源扫描，双层探测器扫描，快速管电压切换扫描)实现多能谱投影数据采集，将导致复杂的技术要求与高昂的系统花费，如传统CT系统难以满足快速管电压切换扫描对球管电压切换频率的高要求。因此降低多能谱CT的技术要求与系统成本，同时精确重建多能图像具有重要的实用价值。

目前有两类实现多能谱CT投影数据采集与图像重建的方式：第一类基于光子计数探测器识别多能信号，以采集不同能量下的投影数据；第二类基于传统单能CT系统获取多个不同能量下的投影数据，同时配合相应重建算法重建多能图像。光子计数探测器基于脉冲高度分析，可辨别不同能量。但由于脉冲响应时间、电荷共享效应及硬件等方面的约束，光子计数探测器型的能谱CT仍处于研究阶段。因此，许多基于传统单能CT系统实现多能谱CT的方法也被提出。Kim等人于2015年提出了基于管电压切换的多能谱CT系统(K.Kim,etal.,Sparse-view spectral CT reconstruction using spectral patch-based low-rank penalty,2015,IEEE Trans.Med.Imaging)。该系统使用传统能量积分型探测器采集投影数据，扫描过程中管电压在三个能量下迅速切换；同时提出极大似然估计约束的方法，利用基于能谱块的低秩约束来重建图像。然而该系统要求球管的电压切换频率超过100Hz,较难实现于传统的CT系统。

为了降低该基于管电压切换的多能谱CT系统的技术要求，Shen等人提出了分段扫描的多能谱CT系统(L.Shen and Y.Xing,Multienergy CT acquisition andreconstruction with a stepped tube potential scan,2015,Medical Physics)。该系统在投影数据采集时将扫描范围分为几个分段，球管在各分段均采用不同的能量进行扫描，以此将球管电压切换频率降低到Hz级别；同时该系统采用基于熵全变差的迭代重建方法进行图像重建。然而该方法重建结果中仍残留有限角度伪影，在某些情况下重建质量仅与传统方法相近。

因此针对现有技术不足，提供基于传统单能CT的多能谱分段稀疏扫描迭代重建方法以解决现有技术不足甚为必要。

发明内容

本发明的目的在于避免现有技术的不足之处而提供基于传统单能CT的多能谱分段稀疏扫描迭代重建方法。该基于传统单能CT的多能谱分段稀疏扫描迭代重建方法可在传统单能CT系统上实现，且能有效地去除噪声与伪影，得到高质量的多能图像。

本发明的上述目的通过以下技术措施实现：

提供一种基于传统单能CT的多能谱分段稀疏扫描迭代重建方法,包括的步骤有：

步骤一、在传统单能CT上采用分段稀疏的扫描模式，应用多个能量的X射线分别在相应的多个有限角度范围内扫描被测物体，获得被测物体在多个能量下的多个分段稀疏投影数据；

步骤二、将步骤一获得的多个分段稀疏投影数据，应用联合全变分正则化与鲁棒性主成分分析约束的图像迭代重建方法，分别重建被测物体在多个能量下的CT图像。

优选的，上述图像迭代重建方法的模型如式(Ⅰ)：

其中为保真项，A为系统矩阵，X为目标图像且X＝{X₁,...,X_k，k≤E_N},X_k为E_N个能量中第k个能量下的目标图像，k为正整数且k≤E_N,为测得的投影矩阵，为E_N个能量中第k个能量下测得的投影矩阵，R(X,X_P)为联合全变分正则化与鲁棒性主成分分析约束，其中先验图像为 为E_N个能量中第k个能量的先验图像。

优选的，上述联合全变分正则化与鲁棒性主成分分析约束R(X,X_P)通过式(Ⅱ)表达：

其中X_L为低秩矩阵，X_S为稀疏矩阵，X_M为目标图像X与先验图像X_P的拼接矩阵，T表示对矩阵的稀疏转换可以为全变差或者紧框架，||·||₁和||·||_*分别表示矩阵的L₁范数和核范数，λ_F、λ_L和λ_S分别表示相应矩阵范数的权重因子。

优选的，上述低秩矩阵X_L和稀疏矩阵X_S由X_M分解得到。

优选的，上述步骤二具体包括有，

步骤2.1、将步骤一得到的多个分段稀疏投影数据全部合并并进行正则化，得到估算出的完整投影数据，进入步骤2.2；

步骤2.2、将步骤2.1得到的完整投影数据进行滤波反投影，得到先验图像，进入步骤2.3；

步骤2.3、对待重建的多个能量下的CT图像赋予初始值，初始值为与待重建图像矩阵大小相同的全零矩阵，进入步骤2.4；

步骤2.4、q为迭代次数，令q＝1进入步骤2.6；

步骤2.5、令q＝q+1进入步骤2.6；

步骤2.6、将初始值和步骤一得到的多个分段稀疏投影数据带入共轭梯度最小二乘算法，得到初始图像；

步骤2.7、将步骤2.2得到的先验图像和步骤2.6得到的初始图像带入联合全变分正则化与鲁棒性主成分分析约束重建得到重建结果，然后对重建结果进行迭代更新，得到迭代终值；

步骤2.8、判断q的值，当q＜L时则令步骤2.7得到的迭代终值为初始值，进入步骤2.5；当q＝L时，则进入步骤2.9，L≥2且为正整数；

步骤2.9、将迭代终值定义为被测物体在多个能量下的CT图像。

优选的，上述L的值范围为50～200。

优选的，上述步骤一中多个能量的X射线通过传统单能CT切换不同的X射线源球管电压值获得。

优选的，上述步骤一中应用不同能量的X射线分别在相应的多个有限角度范围内扫描被测物体即以分段稀疏的扫描模式进行多能谱投影数据的采集。

优选的，上述扫描模式的扫描范围以一定的分段角度被平均分为多个分段，球管电压在扫描每个分段时保持不变，分段结束后切换至另一电压，多个电压交替切换，扫描一周即完成多能谱投影数据的采集。

优选的，上述分段角度通过式(Ⅲ)定义如下：

θ＝2π/(Ω*E_N) 式(Ⅲ)，

其中，θ为分段角度，Ω为正整数，E_N为能谱数量。

优选的，上述步骤2.1中，通过式(Ⅳ)将步骤一得到的多个分段稀疏投影数据全部合并并进行正则化，得到被测物体在2π范围内估算的完整投影数据；

其中w_k为正则化算子且w_k表示第k个能量下投影的L₁范数的倒数，为E_N个能量中第k个能量下的投影，Y_P包含E_N个不同能量下估算的完整投影，为E_N个能量中第k个能量下估算的完整投影，′为转置。

优选的，上述步骤2.2具体为，将步骤2.1得到的完整投影数据，通过式(Ⅴ)滤波反投影算法得到多个能量下的先验图像；

其中X_P为先验图像，为第k个能量下的先验图像。

优选的，上述步骤2.6具体为，将初始值和步骤一得到的多个分段稀疏投影数据带入式(Ⅵ)共轭梯度最小二乘算法，得到初始图像：

其中X_C为初始图像,为第k个能量下的初始图像。

优选的，上述步骤2.7具体为，将步骤2.2得到的先验图像和步骤2.6得到的初始图像通过联合全变分正则化与鲁棒性主成分分析约束重建得到重建结果，其中采用奇异值阈值法以最小化低秩矩阵X_L的核范数，采用TV转换以稀疏目标图像X与稀疏矩阵X_S得到最优解，然后对重建结果进行迭代更新，得到迭代终值。

在求解稀疏矩阵X_S与目标图像X时，||TX||₁通过式(Ⅶ)计算得到，||TX_S||₁式(Ⅷ)计算得到，

其中r和c为像素点所在的行和列；

求解低秩矩阵X_L时，由奇异值分解XL＝U∑V^*，X_L通过奇异值阈值法式(Ⅸ)计算得到，

SVT(X_L，b)＝UΛ(∑)V^* 式(Ⅸ)，

其中b为阈值,∑＝diag{a_i}为对角矩阵，Λ(∑)＝diag{max(a_i-b，0)}为奇异值收缩算子。

本发明的基于传统单能CT的多能谱分段稀疏扫描迭代重建方法，包括的步骤有：步骤一、在传统单能CT上采用分段稀疏的扫描模式，应用多个能量的X射线分别在相应的多个有限角度范围内扫描被测物体，获得被测物体在多个能量下的多个分段稀疏投影数据。步骤二、将步骤一获得的多个分段稀疏投影数据，应用联合全变分正则化与鲁棒性主成分分析约束的图像迭代重建方法，分别重建被测物体在多个能量下的CT图像。本发明利用不同能量下被测物体结构的一致性，设计分段稀疏扫描模式，大大降低基于快速管电压切换的多能谱CT对球管电压切换频率的技术要求。同时本发明为解决投影数据的不完备性，提出联合全变分正则化与鲁棒性主成分分析约束的图像迭代重建方法，有效地去除噪声与伪影，得到高质量的重建结果。

附图说明

利用附图对本发明作进一步的说明，但附图中的内容不构成对本发明的任何限制。

图1为本发明基于传统单能CT的多能谱分段稀疏扫描迭代重建方法的流程图。

图2为实施例2所用的测试体模。

图3为传统单能CT系统在球管电压分别为80kVp、100kVp和120kVp时发出的能谱曲线图。

图4为分段稀疏扫描示意图,其中E₁,E₂,E₃表示球管电压值分别为80kVp、100kVp和120kVp。

图5为X射线球管电压随扫描角度变化的曲线图。

图6a为球管电压80kVp时采集到的投影数据。

图6b为球管电压100kVp时采集到的投影数据。

图6c为球管电压120kVp时采集到的投影数据。

图7a为基于图6的投影数据采用本发明方法合并正则化处理后，在球管电压为80kVp时估算出的2π范围内的完整投影。

图7b为基于图6的投影数据采用本发明方法合并正则化处理后，在球管电压为100kVp时估算出的2π范围内的完整投影。

图7c为基于图6的投影数据采用本发明方法合并正则化处理后，在球管电压为120kVp时估算出的2π范围内的完整投影。

图8a基于图7a的投影数据，采用FBP算法得到在球管电压为80kVp时的先验图像。

图8b基于图7b的投影数据，采用FBP算法得到在球管电压为100kVp时的先验图像。

图8c基于图7c的投影数据，采用FBP算法得到在球管电压为120kVp时的先验图像。

图9a基于图8a的投影数据，采用CGLS算法迭代50次得到在球管电压为80kVp时的初始图像。

图9b基于图8b的投影数据，采用CGLS算法迭代50次得到在球管电压为100kVp时的初始图像。

图9c基于图8c的投影数据，采用CGLS算法迭代50次得到在球管电压为120kVp时的初始图像。

图10a基于图6a的投影数据，采用CGLS算法迭代200次在球管电压为80kVp时的重建结果。

图10b基于图6b的投影数据，采用CGLS算法迭代200次在球管电压为100kVp时的重建结果。

图10c基于图6c的投影数据，采用CGLS算法迭代200次在球管电压为120kVp时的重建结果。

图11a基于图8a的投影数据，并结合图9a与图10a所示的先验图像与初始图像，采用本发明方法在球管电压为80kVp时的重建结果。

图11b基于图8b的投影数据，并结合图9b与图10b所示的先验图像与初始图像，采用本发明方法在球管电压为100kVp时的重建结果；

图11c基于图8c的投影数据，并结合图9c与图10c所示的先验图像与初始图像，采用本发明方法在球管电压为120kVp时的重建结果。

具体实施方式

结合以下实施例对本发明的技术方案作进一步说明。

实施例1。

基于传统单能CT的多能谱分段稀疏扫描迭代重建方法，如图1所示，包括的步骤有：

步骤一、在传统单能CT上采用分段稀疏的扫描模式，应用多个能量的X射线分别在相应的多个有限角度范围内扫描被测物体，获得被测物体在多个能量下的多个分段稀疏投影数据；

步骤二、将步骤一获得的多个分段稀疏投影数据，应用联合全变分正则化与鲁棒性主成分分析约束的图像迭代重建方法，分别重建被测物体在多个能量下的CT图像。

本发明的多个能量的X射线是指具有不同能谱曲线的多个X射线。

图像迭代重建方法的模型如式(Ⅰ)：

其中为保真项，A为系统矩阵，X为目标图像且X＝{X₁,...,X_k，k≤E_N},X_k为E_N个能量中第k个能量下的目标图像，k为正整数且k≤E_N,为测得的投影矩阵，为E_N个能量中第k个能量下测得的投影矩阵，R(X,X_P)为联合全变分正则化与鲁棒性主成分分析约束，其中先验图像为 为E_N个能量中第k个能量的先验图像。

联合全变分正则化与鲁棒性主成分分析约束R(X,X_P)通过式(Ⅱ)表达：

其中X_L为低秩矩阵，X_S为稀疏矩阵，X_M为目标图像X与先验图像X_P的拼接矩阵，T表示对矩阵的稀疏转换可以为全变差或者紧框架，||·||₁和||·||_*分别表示矩阵的L₁范数和核范数，λ_F、λ_L和λ_S分别表示相应矩阵范数的权重因子。

低秩矩阵X_L和稀疏矩阵X_S由X_M分解得到。

本发明的方法不仅充分利用了不同能谱间的结构相似信息，同时充分利用了梯度信息和主成分分析以约束由分段稀疏扫描导致投影数据不完备引起的噪声和伪影并保持结构细节。其中联合全变分正则化与鲁棒性主成分分析约束不仅通过传统的鲁棒性主成分分析对图像的低秩部分及稀疏部分进行约束以保留结构细节信息，同时还加入了先验信息与梯度信息，在引入不同能谱间的结构相似信息的同时去除噪声提高平滑度。

步骤二具体包括有，

步骤2.1、将步骤一得到的多个分段稀疏投影数据全部合并并进行正则化，得到估算出的完整投影数据，进入步骤2.2；

步骤2.2、将步骤2.1得到的完整投影数据进行滤波反投影，得到先验图像，进入步骤2.3；

步骤2.3、对待重建的多个能量下的CT图像赋予初始值，初始值为与待重建图像矩阵大小相同的全零矩阵，进入步骤2.4；

步骤2.4、q为迭代次数，令q＝1进入步骤2.6；

步骤2.5、令q＝q+1进入步骤2.6；

步骤2.6、将初始值和步骤一得到的多个分段稀疏投影数据带入共轭梯度最小二乘算法，得到初始图像；

步骤2.7、将步骤2.2得到的先验图像和步骤2.6得到的初始图像带入联合全变分正则化与鲁棒性主成分分析约束得到重建结果，然后对重建结果进行迭代更新，得到迭代终值；

步骤2.8、判断q的值，当q＜L时则令步骤2.7得到的迭代终值为初始值，进入步骤2.5；当q＝L时，则进入步骤2.9，L≥2且为正整数；

步骤2.9、将迭代终值定义为被测物体在多个能量下的CT图像。

需说明的是，本发明的步骤2.7得到的迭代终值后再重新进入步骤2.5至步骤2.7时，这时迭代终值可以理解为约束更新重建图像。

本发明L的值建议范围为50～200，也可以为其他的任意的范围，具体的L的值根据实际情况而定。本实施例具体为150。

步骤一中多个能量的X射线通过传统单能CT切换不同的X射线源球管电压值获得。

步骤一中应用不同能量的X射线分别在相应的多个有限角度范围内扫描被测物体即以分段稀疏的扫描模式进行多能谱投影数据的采集。

扫描模式的扫描范围以一定的分段角度被平均分为多个分段，球管电压在扫描每个分段时保持不变，分段结束后切换至另一电压，多个电压交替切换，扫描一周即完成多能谱投影数据的采集。

分段角度通过式(Ⅲ)定义如下：

θ＝2π/(Ω*E_N) 式(Ⅲ)，

其中，θ为分段角度，Ω为正整数，E_N为能谱数量。

需说明的是，本发明在每次完整的分段稀疏扫描中，共有E_N个多能子集。每个子集有360/E_N个投影，且每个子集由Ω个角度为θ的有限角度投影组成。该扫描模式只需少次切换球管电压，由此可将球管电压切换频率降低至Hz级别。球管相对旋转中心扫描一周即完成投影数据的采集，因此可在传统单能CT系统上实现且不引入额外辐射剂量。

本发明不同能谱的X射线扫描范围彼此不同，直接将所有能谱的投影数据组合可得到被测物体在2π范围内的完整投影。但是由于各个能量的X射线投影数据幅度不同，直接组合获得的2π范围内的完整投影不能真实反映各能量的幅度信息。为获得各个能量在2π范围内的完整投影，需将所有投影数据进行归一化，然后对归一化后的投影数据整体赋予相应能量下的幅度信息。该操作可得到不同能量在2π范围内估算的完整投影Y_P。

步骤2.1中，通过式(Ⅳ)将步骤一得到的多个分段稀疏投影数据全部合并并进行正则化，得到被测物体在2π范围内估算的完整投影数据。

其中w_k为正则化算子且w_k表示第k个能量下投影的L₁范数的倒数，为E_N个能量中第k个能量下的投影，Y_P包含E_N个不同能量下估算的完整投影，为E_N个能量中第k个能量下估算的完整投影，′为转置。

步骤2.2具体为，将步骤2.1得到的完整投影数据，通过式(Ⅴ)滤波反投影算法得到多个能量下的先验图像。

其中X_P为先验图像，为第k个能量下的先验图像。

步骤2.6具体为，将初始值和步骤一得到的多个分段稀疏投影数据带入式(Ⅵ)共轭梯度最小二乘算法，得到初始图像：

其中X_C为初始图像,为第k个能量下的初始图像。

步骤2.7具体为，将步骤2.2得到的先验图像和步骤2.6得到的初始图像通过联合全变分正则化与鲁棒性主成分分析约束重建得到重建结果，其中采用奇异值阈值法以最小化低秩矩阵X_L的核范数，采用TV转换以稀疏目标图像X与稀疏矩阵X_S得到最优解，然后对重建结果进行迭代更新，得到迭代终值。

在求解稀疏矩阵X_S与目标图像X时，||TX||₁通过式(Ⅶ)计算得到，||TX_S||₁式(Ⅷ)计算得到，

其中r和c为像素点所在的行和列。

求解低秩矩阵X_L时，由奇异值分解X_L＝U∑V^*，X_L通过奇异值阈值法式(Ⅸ)计算得到，

SVT(X_L，b)＝UΛ(∑)V^* 式(Ⅸ)。

其中b为阈值,∑＝diag{a_i}为对角矩阵，Λ(∑)＝diag{max(a_i-b，0)}为奇异值收缩算子。

在奇异值阈值法中，奇异值逐步变小或为零，最终低秩矩阵得到最优解。

该基于传统单能CT的多能谱分段稀疏扫描迭代重建方法，包括的步骤有：步骤一、在传统单能CT上采用分段稀疏的扫描模式，应用多个能量的X射线分别在相应的多个有限角度范围内扫描被测物体，获得被测物体在多个能量下的多个分段稀疏投影数据。步骤二、将步骤一获得的多个分段稀疏投影数据，应用联合全变分正则化与鲁棒性主成分分析约束的图像迭代重建方法，分别重建被测物体在多个能量下的CT图像。本发明利用不同能量下被测物体结构的一致性，设计分段稀疏扫描模式，大大降低对球管电压切换频率的技术要求。同时本发明为解决投影数据的不完备性，提出联合全变分正则化与鲁棒性主成分分析约束的图像迭代重建方法，有效地去除噪声与伪影，得到高质量的重建结果。

实施例2。

基于传统单能CT的多能谱分段稀疏扫描迭代重建方法，在传统单能CT上采用分段稀疏扫描模式，应用不同能量的X射线分别在相应的多个有限角度范围内扫描被测物体，获得被测物体在各个能量下的投影数据，然后利用所测得的多个能量下的投影数据，使用联合全变分正则化与鲁棒性主成分分析约束的图像迭代重建方法分别重建各个能量下被测物体的CT图像。

本实施例中所用的测试体模如图2所示，为NCAT体模(一个具有丰富信息的胸腹部体模)对本方法进行仿真分析，尺寸为512mm×512mm。

系统的扫描参数为：射线源至平板探测器距离为800mm，旋转中心至平板探测器距离为400mm，探测器由512个1mm的探测单元组成。

本实施例选用三个能量进行分段稀疏扫描，图3所示为本实施例中传统单能CT系统在球管电压分别为80kVp、100kVp和120kVp时发出的能谱曲线图，其中能谱的采样间隔为1keV。

本实施例的分段稀疏扫描示意图，如图4所示，其中能谱数量E_N为3，三个能量E₁,E₂,E₃分别为80kVp、100kVp和120kVp，分段角度θ为2π/15，此时正整数Ω为5。在以分段稀疏扫描模式采集投影数据时，2π范围内每隔1度采集一个投影，扫描一周即完成所有多能投影数据的采集。依次切换三个能量80kVp、100kVp、120kVp，以分段角度为2π/15进行投影数据采集。

图5为本实施例分段稀疏扫描模式下，球管电压随扫描角度变化的曲线图。

当扫描角度为时，球管电压为80kVp；

当扫描角度为时，球管电压为100kVp；

当扫描角度为时，球管电压为120kVp；其中i＝1，2，3，4。

图6为本实施例分段稀疏扫描模式下采集到的不同能量下分段稀疏的多能投影数据Y_E。其中图6a为球管电压为80kVp时的分段稀疏的投影数据图6b为球管电压为100kVp时的分段稀疏的投影数据图6c为球管电压为120kVp时的分段稀疏的投影数据

采用式(Ⅳ)对图6所示的分段稀疏投影的数据进行合并正则化，估算出各个能量下2π范围内的完整投影数据Y_P分别如图7所示，其中图7a为球管电压为80kVp时估算出的2π范围内的完整投影数据图7b为球管电压为100kVp时估算出的2π范围内的完整投影数据图7c为球管电压为120kVp时估算出的2π范围内的完整投影数据

将估算出的完整投影数据Y_P带入式(Ⅴ)，通过滤波反投影算法，获得各个能量下的先验图像X_P分别如图8所示，其中图8a为球管电压为80kVp时的先验图像图8b为球管电压为100kVp时的先验图像图8c为球管电压为120kVp时的先验图像

将各能量下分段稀疏的投影数据Y_E带入式(Ⅵ)，通过共轭梯度最小二乘算法迭代50次，得到各个能量下的初始图像X_C分别如图9所示，其中图9a为球管电压为80kVp时的初始图像图9b为球管电压为100kVp时的初始图像图9c为球管电压为120kVp时的初始图像

基于图6所示的本实施例分段稀疏扫描模式下采集到投影数据Y_E，直接采用传统共轭梯度最小二乘算法迭代200次后的重建结果如图10所示。其中图10a为球管电压为80kVp时传统迭代重建方法的重建结果；图10b为球管电压为100kVp时传统迭代重建方法的重建结果；图10c为球管电压为120kVp时传统迭代重建方法的重建结果。由于存在有限角度和稀疏角度的投影数据且数据量不足，导致传统迭代重建方法得到的重建结果受噪声及伪影干扰，严重影响了图像诊断信息的提取。

基于图6所示的本实施例分段稀疏扫描模式下采集到投影数据Y_E，并结合图8及图9所示的先验图像和初始图像，带入式(Ⅰ)，即采用本发明图像扫描模式与联合全变分正则化与鲁棒性主成分分析约束的图像迭代方法在迭代200次后的重建结果如图11，其中图11a为球管电压为80kVp时本发明方法的重建结果X₁；图11b为球管电压为100kVp时本发明方法的重建结果X₂；图11c为球管电压为120kVp时本发明方法的重建结果X₃。可以看出，相较于传统重建方法重建结果图10，本发明实施例的方法可以更有效地抑制噪声及伪影，边缘等结构信息保存更完整，可重建更高质量的能谱图像。

通过采用结构相似度(SSIM)、均方误差(MSE)和对比度噪声比(CNR)对重建图像进行定量分析,可以发现：本发明的方法重建结果与参考图像SSIM平均达0.98，MSE平均值为0.0487，CNR平均值为14.48。由此可知，本发明方法重建结果误差很小，进一步说明了本发明方法可以准确重建优质的多能图像。

本实施例说明本方法适用于传统单能CT系统下的多能投影数据采集与多能图像精确重建。该方法不仅易于实现，且可以重建优质图像，能够大大降低多能谱CT系统的技术要求与系统花费。

最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。

Claims (10)

1.基于传统单能CT的多能谱分段稀疏扫描迭代重建方法,其特征在于，具有如下步骤：

步骤一、在传统单能CT上采用分段稀疏的扫描模式，应用多个能量的X射线分别在相应的有限角度范围内扫描被测物体，获得被测物体在多个能量下的多个分段稀疏投影数据；

步骤二、将步骤一获得的多个分段稀疏投影数据，应用联合全变分正则化与鲁棒性主成分分析约束的图像迭代重建方法，分别重建被测物体在多个能量下的CT图像。

2.根据权利要求1所述的基于传统单能CT的多能谱分段稀疏扫描迭代重建方法,其特征在于:所述图像迭代重建方法的模型如式(Ⅰ)：

其中为保真项，A为系统矩阵，X为目标图像且X＝{X₁,...,X_k，k≤E_N},X_k为E_N个能量中第k个能量下的目标图像，k为正整数且k≤E_N,为测得的投影矩阵，为E_N个能量中第k个能量下测得的投影矩阵，R(X,X_P)为联合全变分正则化与鲁棒性主成分分析约束，其中先验图像为 为E_N个能量中第k个能量的先验图像。

3.根据权利要求2所述的基于传统单能CT的多能谱分段稀疏扫描迭代重建方法,其特征在于：所述联合全变分正则化与鲁棒性主成分分析约束R(X,X_P)通过式(Ⅱ)表达：

其中X_L为低秩矩阵，X_S为稀疏矩阵，X_M为目标图像X与先验图像X_P的拼接矩阵，T表示对矩阵的稀疏转换可以为全变差或者紧框架，||·||₁和||·||_*分别表示矩阵的L₁范数和核范数，λ_F、λ_L和λ_S分别表示相应矩阵范数的权重因子；

所述低秩矩阵X_L和稀疏矩阵X_S由X_M分解得到。

4.根据权利要求3所述的基于传统单能CT的多能谱分段稀疏扫描迭代重建方法,其特征在于：所述步骤二具体包括有，

步骤2.1、将步骤一得到的多个分段稀疏投影数据全部合并并进行正则化，得到估算出的完整投影数据，进入步骤2.2；

步骤2.2、将步骤2.1得到的完整投影数据进行滤波反投影，得到先验图像，进入步骤2.3；

步骤2.3、对待重建的多个能量下的CT图像赋予初始值，初始值为与待重建图像矩阵大小相同的全零矩阵，进入步骤2.4；

步骤2.4、q为迭代次数，令q＝1进入步骤2.6；

步骤2.5、令q＝q+1进入步骤2.6；

步骤2.6、将初始值和步骤一得到的多个分段稀疏投影数据带入共轭梯度最小二乘算法，得到初始图像；

步骤2.7、将步骤2.2得到的先验图像和步骤2.6得到的初始图像带入联合全变分正则化与鲁棒性主成分分析约束重建得到重建结果，然后对重建结果进行迭代更新，得到迭代终值；

步骤2.8、判断q的值，当q＜L时则令步骤2.7得到的迭代终值为初始值，返回步骤2.5；当q＝L时，则进入步骤2.9，L≥2且为正整数；

步骤2.9、将迭代终值定义为被测物体在多个能量下的CT图像。

5.根据权利要求4所述的基于传统单能CT的多能谱分段稀疏扫描迭代重建方法,其特征在于：所述L的值范围为50～200；

所述步骤一中多个能量的X射线通过传统单能CT切换不同的X射线源球管电压值获得。

6.根据权利要求5所述的基于传统单能CT的多能谱分段稀疏扫描迭代重建方法,其特征在于：所述步骤一中应用不同能量的X射线分别在相应的多个有限角度范围内扫描被测物体即以分段稀疏的扫描模式进行多能谱投影数据的采集；

所述扫描模式的扫描范围以一定的分段角度被平均分为多个分段，球管电压在扫描每个分段时保持不变，分段结束后切换至另一电压，多个电压交替切换，扫描一周即完成多能谱投影数据的采集；

分段角度通过式(Ⅲ)定义如下：

θ＝2π/(Ω*E_N) 式(Ⅲ)，

其中，θ为分段角度，Ω为正整数，E_N为能谱数量。

7.根据权利要求6所述的基于传统单能CT的多能谱分段稀疏扫描迭代重建方法,其特征在于：所述步骤2.1中，通过式(Ⅳ)将步骤一得到的多个分段稀疏投影数据全部合并并进行正则化，得到被测物体在2π范围内估算的完整投影数据；

其中w_k为正则化算子且w_k表示第k个能量下投影的L₁范数的倒数，为E_N个能量中第k个能量下的投影，Y_P包含E_N个不同能量下估算的完整投影，为E_N个能量中第k个能量下估算的完整投影，为转置。

8.根据权利要求7所述的基于传统单能CT的多能谱分段稀疏扫描迭代重建方法,其特征在于：所述步骤2.2具体为，将步骤2.1得到的完整投影数据，通过式(Ⅴ)滤波反投影算法得到多个能量下的先验图像；

其中X_P为先验图像，为第k个能量下的先验图像。

9.根据权利要求8所述的基于传统单能CT的多能谱分段稀疏扫描迭代重建方法,其特征在于：所述步骤2.6具体为，将初始值和步骤一得到的多个分段稀疏投影数据带入式(Ⅵ)共轭梯度最小二乘算法，得到初始图像：

其中X_C为初始图像,为第k个能量下的初始图像。

10.根据权利要求9所述的基于传统单能CT的多能谱分段稀疏扫描迭代重建方法,其特征在于：所述步骤2.7具体为，将步骤2.2得到的先验图像和步骤2.6得到的初始图像通过联合全变分正则化与鲁棒性主成分分析约束重建得到重建结果，其中采用奇异值阈值法以最小化低秩矩阵X_L的核范数，采用TV转换以稀疏目标图像X与稀疏矩阵X_S得到最优解，然后对重建结果进行迭代更新，得到迭代终值；

在求解稀疏矩阵X_S与目标图像X时，||TX||₁通过式(Ⅶ)计算得到，||TX_S||₁式(Ⅷ)计算得到，

其中r和c为像素点所在的行和列；

求解低秩矩阵X_L时，由奇异值分解X_L＝U∑V^*，X_L通过奇异值阈值法式(Ⅸ)计算得到，

SVT(X_L，b)＝UΛ(∑)V^* 式(Ⅸ)，

其中b为阈值,∑＝diag{a_i}为对角矩阵，Λ(∑)＝diag{max(a_i-b，0)}为奇异值收缩算子。