能谱CT图像重建算法研究

赵 星， 张 朋

(首都师范大学 数学科学学院， 北京 100048)

0 引言

计算机断层成像(Computed Tomography，CT)具有扫描速度快、图像清晰等特点，已成为当今医学诊断和工业检测常用的影像技术。医学或工业CT中常用的X射线源通过加速电子撞击钨靶(或钼靶等)产生不同能量的X光子，这些不同能量X光子的数目服从一定概率分布曲线，称为X光子能谱，其中包含了靶材料的特征光子和轫致辐射光子。X光管在100 kV电压下轰击钨靶产生的X射线能谱如图1所示。

图1 X光管在100 kV电压下轰击钨靶产生的X射线能谱，特征辐射光子表现为X射线能谱曲线上的阶跃式尖峰[1]

常规CT图像重建算法没有考虑X射线能谱，而是基于单能X射线与物质的作用模型进行计算，其重建图像的CT值是各个像素处对不同能量X光子的衰减系数的某种非线性平均。这种平均可能会导致“低密度的高原子序数物质”与“高密度的低原子序数物质”的CT值相近甚至相同，使得通过CT值难以区分某些组织和造影剂。这种平均并非全局一致的，因而会导致CT图像带有某些固有的图像畸变，如硬化伪影、放射伪影等，影响CT图像临床诊断的准确性。

为了减少图像伪影，提高CT对物体成分识别和定量化分析的能力，先后有多种双能谱CT(Dual-spectral CT)和多能谱CT(Multi-spectral CT)成像模式被提出，本文统称为能谱CT。能谱CT以某种方法采集穿过被测物体的多个能谱(或能段)的X射线测量数据，并将组成材料关于X射线能量的衰减系数近似为一组基函数的线性组合，通过计算实现被测物体基材料成分的分解和图像重建。根据基材料图像，或者根据由此计算出的被测物体化学成分的等效电子密度和等效原子序数，可实现对物体成分的识别和定量化分析。同时由基材料图像合成的虚拟单能图像可以有效减少图像伪影。利用射线源分别在高电压为110 kV和低电压为60 kV时产生的服从一定能谱分布的多能X射线，扫描同一个上肢模体[图2(a)]，由高低能谱相应的投影数据分别用FBP算法直接重建的其中一个断层的高能图像[图2(b)]、低能图像[图2(c)]，使用基材料分解算法重建的骨基材料图像[图2(d)]和软组织基材料图像[图2(e)]，以及用两个基材料图像合成的虚拟单能图像[图2(f)]，可看出虚拟单能图像中伪影显著减少。

图2 针对上肢模体一个断层的双能谱CT成像(a)上枝模体；(b)高能谱重建CT图像；(c)低能谱重建CT图像；(d)骨组织基材料图像；(e)软组织基材料图像；(f)110 keV虚拟单能图像

双能谱CT经过近20年的发展，已在医学诊疗中大量应用，有效改善了对一些病灶的辨识度，提升了对这些病症的诊断能力[2-4]。例如用于确定肾结石尿结石成分[5-6]、肝脏铁浓度、骨髓成分、斑块特征[7-8]等；用于钙含量定量化，以评估骨骼和牙齿健康；评估乳腺癌风险[9-10]等。另外，通过基材料分解后合成的虚拟单能图像，能谱CT还可以显著减轻医源性高吸收植入物(如金属固定物、起搏器、栓塞弹簧、支架、种植牙等)在常规CT图像中导致的严重金属伪影等，提高骨折金属固定术后结构信息的观察水平[11]；能谱CT还能比常规CT血管造影法(CTA)更好地区分血管与血管壁钙化、血管与支架，提高血管堵塞鉴别诊断能力，提高对血管尤其是静脉栓塞的成分定性诊断准确性；此外，还可以诊断尿酸盐结晶沉积，有助于痛风鉴别诊断[12-13]。

在下面的章节，首先介绍能谱CT的数据采集模式；然后聚焦于能谱CT图像重建算法的综述，包括基材料分解数学模型，基材料分解算法分类以及各类算法的优缺点、适用性分析，基于深度学习的基材料分解算法等，并进一步探讨相关的研究问题。

1 能谱CT的数据采集模式

利用不同能量的X射线扫描被测物体进行CT成像的概念早在20世纪70年代就提出来了[14]，至今发展出如图3所示的大致6种能谱CT数据采集模式[15]，包括两次扫描模式[图3(a)]、双源双探模式[图3(b)]、电压快速切换模式[图3(c)]、双层探测器模式[图3(d)]、滤波片预调制模式[图3(e)]、光子计数探测器(photon-counting detector, PCD)模式[图3(f)]等。其中两次扫描模式就是射线源在较高电压和较低电压下对同一个被测物体分别扫描两次，分别获得高低能谱的投影数据。两次扫描模式的优点是能谱CT系统硬件成本低，缺点是延长了扫描时间，增大了辐射剂量，且存在配准问题。典型的代表如Canon Aquilion Precision 能谱CT[16]。

图3 能谱CT的数据采集模式[15](a) 两次扫描模式;(b) 双源双探模式;(c) 电压快速切换模式;(d)双层探测器模式;(e) 滤波片预调制模式;(f) 光子计数探测器模式

2005年，Siemens公司推出了国际上首款实用化医用双能谱CT系统Somatom Force。其中采用了垂直放置的两套“射线源和探测器”系统。工作时，两套系统分别以不同电压扫描患者，采集相应高低能谱的两套数据。其优点是扫描时间短，不足之处是需要两套射线源探测器系统，硬件成本高，且双源采集技术会受到交叉散射的影响。该能谱CT已获得了一系列成功的临床应用[17-19]。

2009年，GE推出了一种红宝石能谱CT系统，其中通过射线源电压瞬变技术，来实现采集相应高低能谱的两套数据，因此,该CT也被称为电压快速切换双能CT[20]。电压快速切换模式具有结构紧凑、简单，成像速度快，视野大等优点，但其对电源技术要求较高，且电压瞬变的延迟效应会导致能谱估计误差，进而导致图像CT值畸变。此外，难以通过滤波片优化高低电压相应的射线能谱，且难以进行变电流优化扫描。还有一种射线源电压慢变的能谱CT成像技术[21]，扫描时相对较慢地变化射线源的电压，降低了对电源技术的要求，且可离散为多个电压变换段的扫描，但对于每个电压变换段是个稀疏角成像问题。

2013年，PHILIPS推出了IQon能谱CT，其中采用了单射线源、双层探测器(也称“三明治”探测器)模式[22]。双层探测器一次扫描即可获得高低能谱两套数据。单射线源、双层探测器模式的主要优点是硬件成本相对较低，扫描剂量也较低。“三明治”探测器采用了上下两层闪烁晶体的设计，两类闪烁体对高低能光子都会吸收。因此，两套数据相应的高低能等效能谱具有较大的重叠区域，高能谱和低能谱之间的相关性较高，会影响图像对双基材料的区分能力。

滤波片调制模式利用在射线源端加装特别设计的滤波片对X射线能谱进行分区或分时进行调制，使射线源在电压不变情况下产生不同能谱，以获取多能谱CT扫描数据[23-24]。这种方法硬件实现简单，只需要在射线源端加装一个滤波器，设备成本增加不明显。通过精心设计滤波片和标定几何参数，会改善CT图像质量。缺点是滤波片本身结构会引入一些伪影，同时射线源的发出的光子经过滤波片被衰减掉一些，不能充分利用射线源发出的光子，增加了统计噪声。另外，通过滤波片调制出的子能谱之间的区分度较低，增加了能谱CT基材料图像重建的难度。已有多种滤波片设计方案被提出，如基于一维栅条状滤波片调制X射线能谱数据获取方法[25]；如在不同旋转角度下对应不同的滤波片厚度的动态滤波片方法[1]；如Siemens的添加了分裂滤波器(Split filter)的Somatom Edge Plus能谱CT，分别用铝和金的组合层以及铝和锡的组合层来对同一个射线源发出的能谱调制出两个子能谱，有效提高了对头颈部扫描的成像质量[24]。

上述几种能谱CT数据获取方法所使用的均是能量积分型探测器(Energy-Integrating Detector，EID)。EID先将射线通过闪烁体转换成可见光，再将可见光通过光敏器件转换成电信号，电信号积分后由A/D转换成数字信号，因此对光子能量没有区分能力[图4(a)]。近年来，光子计数探测器(photon-counting detector, PCD)发展迅速。PCD将半导体材料联结在专用集成电路上，单个光子撞击在半导体敏感器上被转换为电荷信号，电荷信号经过ASIC放大并与给定的阈值比较，当脉冲大于阈值时计数为1，否则计数为0，从而达到对光子的能量区分和计数[图4(b)]。PCD与EID相比，具有以下优势[26]：① PCD具有光子能量区分能力，可按像素设定光子能量阈值，一次扫描即可获取多套能谱段的投影数据，不但可以有效降低辐射剂量[27]，而且比现有EID型双能谱CT获取的高低能谱数据有更好的能量区分性，提高了材料定量分解的精度[28-29]；② PCD可以选择能量阈值，因此，可利用高Z元素对比剂(如钆、金或铋)线性衰减系数在诊断能量处的不连续性，可以方便地针对造影剂的K-edge进行能谱CT成像；③ PCD还可利用能量阈值可调节的特性，设定合适的阈值消除探测器内部元件的暗电流噪声造成的假计数[30]。实验证明，光子计数探测器进行多能谱成像可以有效地提高图像分辨能力，增强不同组织之间的对比，并且有效降低辐射剂量[31-32]。围绕PCD-CT的研究与开发是近年来国际CT的重点发展方向。光子计数探测器自2001年以来得到了迅速的发展，目前已有不少公司有相对成熟且已面市的PCD产品[33-34]。在探测器芯片方面，最著名的是CERN(欧洲核子研究中心)率先研发的Medipix3产品[26]、荷兰ASI公司(与CERN合作)的Timepix系列产品[35]。在探测器材料方面，典型的光子敏感半导体材料有硅(Si)、砷化镓(GaAs)、碲化镉(CdTe)、碲锌镉(CZT)等。目前基于碲化镉的PCD发展最快。在适用于小型系统的PCD产品方面，著名的产品包括瑞典XCounter公司推出的XC-FLITE CdTe-CMOS系列产品，瑞士Dectris公司推出的Pilatus3 X CdTe系列产品。近期多家国际厂商也公布了其全身光子计数CT新产品，其中德国西门子公司的采用CdTe材料，美国GE公司采用深硅 (Deep Silicon)材料，荷兰飞利浦公司和日本佳能公司采用CZT材料。其中西门子医疗的NAEOTOM Alpha光子计数CT于2021年9月30日率先获得美国食品药品监督管理局(FDA)许可。

图4 能量积分型(EID)探测器和光子计数型(PCD)探测器的原理及对比(a)能量积分型(EID)探测器原理；(b)光子计数型(PCD)探测器原理

2 基于基材料分解的能谱CT数学模型

在考虑X光子能谱的情况下，第k个探测器沿着第L条射线路径扫描物体的测量数据相应的多色投影pk(L)，可用如下数学模型表示：

L∈Πk,k=1,2,…,K

(1)

式中，Sk(E,L)为第k个归一化等效能谱；μ(x,E)为在x处的物体对能量E的光子的线衰减系数；Πk为第k个能谱的扫描射线路径集合。Sk(E,L)与射线源、滤波片、探测器晶体、探测器电路等有关，可进一步写为:

Sk(E,L)=Sk(E)·Q(E,L)

(2)

这里Q(E,L)是能量依赖的探测器响应函数。在能谱CT图像重建时，待求的是未知函数μ(x,E)，由所采集的投影数据来直接求解未知函数是欠定的。为此将物质的线衰减系数近似表示成M个基函数的线性组合：

(3)

式中，fm(E)为预先标定的基函数；am(x)为基函数的组合系数。这样未知函数由求解μ(x,E)变为求解am(x)。

目前近似线衰减系数的模型主要有两种：一种称为双效应分解模型，由Alvarez等于1976年提出[13]。该方法的原理是根据低能段X光子衰减的主要原因是光电效应和康普顿散射，将光电效应截面和康普顿散射截面随X光子能量变化的函数作为两个基函数，即f1(E)=E-3,f2(E)=fKN(E)，这里fKN(E)是Klein-Nishina函数。分别求解相应的光电效应系数图像a1(x)和康普顿散射系数图像a2(x)，进而计算每个像素处的等效原子序数和平均电子密度。在低能段虽然光电效应和康普顿散射是光子衰减的主因，但还存在瑞利散射的贡献。该模型忽略瑞利散射，所以为近似模型。

另一种称为基材料分解模型，由Kalender等于1986年提出[36]。在低能段，如医学CT常用的 20～140 keV能段，对人体中常见的元素(Z≤20)随X光子能量的线衰减系数曲线做奇异值SVD分解，会发现前两个特征值远远大于后面的特征值，表明前20 号元素的衰减系数所组成的空间，使用两个线性无关的向量基本就可以近似表达。为此，可以将基函数设为第m种材料的质量衰减系数或者标准密度下的线衰减系数，即fm(E)=μm(E)，求解相应基材料的等效密度或者体积分数分布图像am(x)。相比于双效应分解模型，基材料分解模型更加灵活，可方便地转化为虚拟单能量图像序列以及算出双效应图像。因此，该模型的使用更为广泛。将(3)式代入(1)得：

L∈Πk,k=1,2,…,K

(4)

定义Am(L)是对am(x)的线积分值，即：

(5)

代入式(4)则为：

L∈Πk,k=1,2,…,K

(6)

在数学上，基于基材料分解的能谱CT问题就是由各个能谱的多色投影数据pk(L),k=1,2,…,K，重建各个基材料的CT图像am(x),m=1,2,…,M，以及由式(3)得到μ(x,E)。

最常用的是以双能谱X射线扫描实现双基材料的CT图像重建，这时M=K，是个正定方程组求解。能谱CT还常常需要对多种基材料成像，一方面用于细分被测物体的成分，如在乳腺能谱CT成像中对脂肪、纤维腺、钙化点等三种基材料成像；另一方面对于含有I、Gd、Au 等元素的造影剂，其线衰减系数曲线的K-edge能量明显高于1-20号元素的K-edge，使得造影剂的衰减系数曲线与原基函数有很大不同，需要增加针对这些造影剂的基函数，提高对组织成分和造影剂的区分能力。当基材料数目大于能谱数时，即M>K，这时为欠定方程组求解。虽然可以通过体积分数之和为1的约束，实现逐像素三基材料图像分解，但体积分数之和为1的约束常常不满足。对于如基于PCD探测器的能谱CT、慢电压切换等扫描模式，可以实现多个能段的扫描，例如4个能段数据的采集。当基材料数目小于能谱数时，即M

3 能谱CT基材料图像重建算法的类型及优缺点

基于基材料分解的能谱CT图像重建算法大致可分为四类，分别是图像域直接分解算法、投影域直接分解算法、迭代分解算法和深度学习分解算法，其中图像域直接映射分解算法和投影域直接映射分解算法是迭代分解算法和深度学习分解算法的基础。下面分别论述各类算法的实现及特点。

3.1 图像域直接分解算法

在假设K个单能射线扫描的情况下，式(6)可整理为：

k=1,2,…,K

(7)

式中,fm,k为第m种基材料在第k个单能时的衰减系数。

在K个能谱射线扫描的情况下，忽略探测器单元的探测效率的非一致性，沿着射线L的多色投影pk(L)也可近似为：

(8)

L∈Πk,k=1,2,…,K

(9)

式中,R-1为Radon逆变换，可由传统重建算法实现。设：

bk=R-1(pk(L)),L∈Πk,k=1,2,…，K

(10)

bk为对投影数据pk直接进行图像重建产生的图像；设：

am=R-1(Am(L)),L∈Πk，m=1,2,…,M

(11)

am为待求得第m种基材料图像，则式(8) 写成矩阵形式为：

(12)

如图5(a)所示，图像域直接分解算法首先对投影数据pk进行图像重建产生图像bk，然后通过线性组合图像bk得到基材料图像am，可表示为式(13)。只要将式(13)中的分解矩阵的逆矩阵标定出，就可由图像bk直接分解为基材料图像aM，该算法也被称为Brooks算法[37]。

图5 两种基材料分解算法的主要步骤对比(a) 图像域直接分解算法;(b) 投影域直接分解算法

(13)

图像域直接分解算法的优点是仅对重建后的CT图像处理，可以直接使用商业CT设备产生的重建图像，易于临床实践。该算法的缺点是基材料分解通过像素之间的直接映射实现，这种映射难以消除由于投影域数据缺陷导致的重建图像中的结构性伪影[38-39]。为此图像域分解算法一般需要先对初始多色投影重建图像进行如水预校正、骨预校正等预处理以减轻硬化伪影，或者通过多色投影高次方重建图像的多项式近似组合减轻硬化伪影[40-41]，然后进行分解。

图像域基材料分解也可以根据式(12)作为线性方程组求解进行迭代计算，在给定初始基材料图像的基础上，迭代更新基材料图像直至收敛。为了更好地抑制噪声对深入分析噪声传递特性，可以将噪声关系作为权重因子加入到分解数据模型中去加权迭代分解，并针对图像域分解仅靠像素映射难以对图像的边缘准确保持，通过结合边缘加权正则化项，较好地抑制了基材料图像噪声[3, 38]。图像域分解这种逐像素的分解方法，通过结合同一像素内3种基材料体积分数之和为1的约束，双能谱CT还可以实现多基材料分解[42]。

3.2 投影域直接分解算法

式(8)写成矩阵形式为：

(14)

或者为：

(15)

如图5(b)所示，在将式(15)中的分解矩阵的逆矩阵标定后，投影域直接分解算法可由采集的多能量投影pk(L)组合计算出各基材料的单色投影Am(L)，即基材料图像的线积分，然后由Am(L)通过FBP算法等重建出基材料图像。需要注意，投影域直接分解算法要求不同能量下的扫描几何参数必须一致，这样扫描射线路径才能重合，即Π1=Π2=…=ΠK=Π，才能实现式(15)所示的在投影域的分解计算。

如果考虑pk(L)与Am(L)之间的非线性关系，可以表示为下式：

Am(L)=Dm(p1(L),p2(L),…pK(L)),

L∈Π,m=1,2,…,M

(16)

式中,Dm为将L射线上测得的投影pk(L)映射为第m个基材料单色投影的非线性函数。对这个非线性函数的准确建模涉及能谱、探测器响应函数、噪声模型、物质在不同能量射线下的衰减系数等参数。在有噪声和信号串扰的情况下，Dm不再对应一个经过原点的单调上升光滑曲面。为了简化这个映射标定过程，常用简单函数近似实现Dm，如采用查找表[43-44]、多项式拟合函数[45-46]、分片多项式[47]、分片查找表等。这些函数可以通过阶梯状模体、多个圆柱模体、阴阳模体等已知简单模体标定[46]。这个方法的优点是简单易用，考虑了能谱的多色性，且具有较强的抗噪能力。缺点是这几种近似函数的表达能力有限，且通常标定模体中包含的信息难以完全涵盖实际被测物体的信息，限制了投影域分解的精度。

3.3 迭代分解算法

迭代分解算法是能谱CT基材料图像分解算法的研究重点，已有多种迭代分解算法被提出。这些迭代分解算法又可进而分为两步迭代法和一步迭代法[48-51]。

两步迭代法的具体步骤可以归纳为：① 由基材料当前图像求本次迭代的多色投影残差，然后由多色投影残差求各个基材料单色投影残差；② 用单色投影残差计算基材料修正图像，并用之更新基材料当前图像；③ 判断收敛条件，满足则结束，不满足则返回第(1)步继续迭代。两步法中间有求单色投影残差的步骤，在给定式(6)后，可以采用Newton-Raphson迭代法及Newton-Gauss迭代法等求解[13, 52, 53]。

一步迭代法的具体步骤可以归纳为：① 由基材料当前图像求本次迭代的多色投影残差，采用共轭梯度法、Chambolle-Pock算法、可分离二次替代函数(separable quadratic surrogate，SQS)等方法直接求解基材料修正图像，并用之更新基材料当前图像[54-55]；② 判断收敛条件，满足则结束，不满足则返回第①步继续迭代。

在考虑射线能谱的情况下，由多色投影到基材料图像的映射关系是非线性的，适合用迭代算法求解。同时迭代算法容易结合成像的物理模型，在过程中还可以结合正则化函数约束求解，通常能获得更高质量的重建结果。迭代算法的缺点是计算量较大，特别是一步迭代法计算量更高，但具有较好的抗噪性能。两步迭代法增加了计算基材料图像单色投影的步骤，将含有能谱的非线性方程部分和CT图像重建的线性部分分离，可以更加方便进行并行计算，且增加了一个变换域，可以在多个域中对图像进行约束和处理。

能谱CT数据采集时，常常有不同能谱的射线扫描路径几何不一致的情况[56]，即高能谱射线路径集合Π高能谱与低能谱射线路径集合Π低能谱不相等，如SIEMENS的双源双探数据采集模式、GE的快速电压切换模式，以及电压缓变下的慢电压切换扫描模式等。这种情况下由多色投影计算基材料单色投影时，不适合用传统Newton-Raphson等方法求解。为此作者团队提出了一种求解几何不一致情况下迭代重建算法即E-ART算法[57]。该算法利用ART算法实现了多色投影到基材料单色投影修正值的计算，然后利用基材料图像域的一致性，更新基材料图像。该算法的优点是重建图像质量高、形式灵活，适于各种常用的CT扫描模式，也可由电压缓变的能谱CT测量数据重建相应基图像和虚拟单能图像[58]。该算法的缺点是收敛速度慢、计算量大。为此，作者团队又提出了多种算法以改进该算法的收敛速度[59-61]。这些算法的核心思想都是估算出一致投影，同时通过参数控制，减少估算误差和噪声对图像质量重建造成的影响。为了每次迭代都同时利用高低能谱两组数据，从而获得更为准确的残差投影，提高算法的收敛速度，作者团队又提出一种适用于几何不一致情况下的IFBP算法[62]。该算法首先在投影域求得所有角度的逐射线加权的多色投影残差，然后对这些加权多色投影残差用FBP算法重建出加权残差图像，最后在图像域利用加权残差图像组合求得基材料修正图像。该算法巧妙地利用了在投影域虽然几何不一致，但在图像域具有像素一一对应的性质。在扫描角度完全情况下，具有抗噪性强、收敛速度快的特点，一般只要2～3次迭代IFBP算法即可得到较好的图像重建结果。在扫描角度稀疏的情况下，重建图像中会有较多稀疏角伪影，且计算基材料修正图像过程属于病态方程组求解，会放大噪声，进而影响基材料分解图像的质量。

在单能CT图像的迭代重建中，已有多种图像先验模型被用于去求CT图像的稀疏解。这些模型如全变分(Total Variation ,TV)、向量化TV[63]、非局部TV[64]、相对TV[65]、方向TV[66]，非局部均值滤波, 紧框架(tight frame)[67]、字典学习模型[68]、基于先验图像模型等，在不完备数据下的CT图像重建中获得了良好的效果。能谱CT成像与单能CT成像一样，也有在测量数据不完备时成像的情况[69]。与单能CT不同的是能谱CT增加了能谱维度，这样既具有了更多种的测量数据不完备的情况，又具有了更大的图像稀疏化空间。为此，传统先验模型被研究人员推广以充分利用能谱CT图像在能谱维度上的结构相似性或低秩性[70]，如总变分核范数(Total Nuclear Variation,TNV)最小化方法[71]、张量字典学习(Tensor Dictionary Learning,TDL)[72]、PRISM(Prior Rank, Intensity and Sparsity Model)[73]、tensor PRISM模型[32]、集合差分相似先验(Correlative multi-channel prior)[74]，以及将全能谱图像作为先验图像(spectral PICCS)[75]等。另外避免单种稀疏模型的局限性，一些正则化模型的相组合方法也被提出，如L0范数与TDL方法相结合模型[76]、非局部思想与低秩约束相结合[77]、分块与低秩及TV相结合[78-79]等。不同正则化项的作用顺序对图像噪声的去除和结构的恢复也很重要，为此，还提出了有限角度下的基于序列正则化的基材料分解方法[80]。迭代分解算法涉及多个变换域，可以增加在不同变换域中对图像的正则化约束，提高基材料图像重建质量[81-82]。以上图像先验模型与优化方法的迭代重建框架一起，可以衍生出很多材料分解算法，相对于非正则化方法获得了显著的效果[81，83-84]。但受限于对CT数据采集过程和图像数据的更深入理解和认知，以及成像条件的多变，面对实际应用，基于以上事先定义好的先验模型的方法仍然有很强的局限性，缺乏随实际数据变化而调整的自适应性。

3.4 深度学习分解算法

近年来深度学习发展迅猛，在CT图像降噪、结构化伪影抑制等方面已经展现了巨大的潜力[85-87]。基于深度学习的CT图像重建方法基本可以分为两类，第一类将深度学习模块用作图像后处理模块，与传统的解析重建算法、迭代重建算法、优化框架等结合，作为正则化模块通过提取图像特征实现去除伪影和噪声。不同于压缩感知需要人为定义图像的稀疏变换，深度学习可以通过数据驱动得到更具适应性的图像先验。因此，可获得比压缩感知方法更显著的效果。第二类将结合物理模型，用深度学习模块不仅实现图像后处理，而且实现重建部分。其主要挑战是CT系统矩阵巨大，虽然有通过神经网络直接映射的方法[88-89]，但网络结构庞大，难以处理大尺寸CT图像重建。因此，目前的方法主要集中于第一类。在各种CNN配置中，U-Net在图像特征提取方面表现出色[90]，是用于CT图像特征提取[91]和伪影消除[92]的最主要主干网络构型。

深度学习本质上只是提供了一种功能更加强大的数学建模和特征提取的方法，基于传统数学模型拓展网络结构才能更加有效率的解决问题。能谱CT基材料图像重建比传统单能CT图像重建涉及更多类型的变换域和图像，因此，能发展出比传统图像重建更多的网络构型和重建算法。例如，文献[93-94]分别提出了使用深度学习在投影域拟合从高低能投影值到材料分解系数积分之间的非线性关系。正如对单能CT图像的重建主要集中于将深度学习模块用作图像后处理模块一样，提出的更多方法用于生成和处理各个变换域的图像。例如，利用多域变换过程结构的一个蝴蝶型网络被用来学习从重建图像到材料分解图像之间的映射过程[95]，或者使用U-net、GAN网络来从多能重建图像生成材料分解图像和虚拟单能图像[96-98]，或者使用ResNet来学习从单能谱下的衰减图像生成两个能量下的虚拟单能图像，进一步由两张虚拟单能图像计算出材料分解图像[99]等。

4 进一步的研究问题

能谱CT的基材料分解问题在数学上是需要一个大规模非线性积分方程组的求解问题。同时，该方程组有高度的病态性，这是因为基材料的质量衰减系数随着射线能量(keV)的变化曲线非常接近[图6(a)]，且扫描物体的不同能谱之间较多重叠[图6(b)]，使得材料分解矩阵的条件数变大。即使如PCD探测器，虽然能对不同能量段的X光子分别计数，理论上能谱互不重叠，但结合探测器的响应函数，等效能段窗不是矩形的，不同能段的等效能谱同样有较多重叠。此外，相对于高能段光子，低能段的光子的测量数据统计误差较大，特别当射线经过高吸收物质时，后者测量数据的统计误差更大，甚至于数据会完全失效。方程组的病态性使得基材料分解过程是一个噪声放大过程，容易受噪声和各类伪影的影响。

图6 骨和软组织的质量衰减系数随着射线能量的变化曲线及有较多重叠的高低能谱(a) 质量衰减系数随着射线能量的变化曲线;(b) 有较多重叠的高低能谱

医学CT对剂量的限制要求严格，主动降低剂量使得穿过病人的光子数较少，或者为了提高时间分辨率，降低采样角度数，被动降低剂量。工业CT受到射线源功率，被测物体材质、尺寸和形状等客观条件限制，很多情况下沿着某些射线路径穿过被测物体的光子数很少，特别是对于像叶片等异型结构。在穿过物体光子数少的X射线路径上，所对应探测器单元的测量数据中统计噪声明显，甚至会完全淹没有掉用的信号。这使得基材料分解图像的噪声水平和伪影程度较能谱投影直接重建图像的更高。怎样通过数据采集方法的改善和算法的改进，降低方程组的病态性，提高能谱基材料分解算法的抗噪性；在较低的辐射剂量下，获得临床可接受的图像，是能谱CT成像领域的一个重要问题。对于PCD探测器，虽然允许应用低能段阈值去除本底电子噪声，但由于对各个能段的光子分别探测和计数，在能段划分过多的情况下，每个能段的光子数显著减少，统计噪声显著增大[55]。因此，不仅要对能段划分数目和每个能段的统计噪声之间寻求平衡，而且同样需要提高算法的抗噪能力。另外和普通CT一样，能谱CT也可以通过稀疏角采集和降低每个角度的射线源能量值降低扫描剂量。设置采集角度数和每个角度的剂量之间需要平衡，这个关系是与被测物体相关的，如果被测物体结构复杂，则需要提高扫描角度数。

投影域迭代重建法依托式(4)定义的较为精确的成像物理模型，该模型依赖于X射线能谱、材料衰减系数以及探测器响应函数等的精确标定，需要标定的参数众多且在标定这些参数的过程中本身又会引入误差。等效能谱的标定虽然已经有很多研究成果，但精度和稳定性仍然有较大的提高空间。探测器的响应不一致，会造成重建图像中的高频环状伪影和低频条带状伪影。特别是对于PCD探测器，半导体将X光子转换为电子的过程中不同探测器单元之间有一定差别，且逐个计数，探测器响应不一致现象比EID探测器更加明显，严重影响基材料图像的分解精度。能谱CT的环状伪影的去除仍然是一个重要的研究问题。

光子散射问题一直是CT图像重建的一个重要研究问题。能谱CT在数据采集过程中，对于双源、双探等扫描模式，还有高低能谱数据采集交叉散射问题。对于PCD能谱CT，由于散射光子的能量比初始光子的能量较低，被PCD探测器低能段错误计数的概率更大。低能光子不仅穿透性差且散射影响更严重，因此PCD低能段数据中不仅统计噪声高，而且散射光子多，需要特别处理。此外PCD探测器的电荷共享、逃逸峰等也都会严重影响它的响应函数。环境对X光子的散射和CT图像也有一定影响。这些信号串扰因素难以准确建模和标定，对基材料准确分解造成了困难。目前一个方法是先对能谱数据采集正过程用多项式投影近似公式建模，然后再用两步法进行迭代求解[13, 100]。该方法简化了正过程的建模和标定，通过迭代中增加正则化约束，减少图像重建误差。另外一个方法是在基材料分解与能谱估计修正交替进行[101]，以及在基材料分解的同时去除散射伪影[102]。这些过程互相促进，可有效提高基材料分解精度。怎样提高抗信号串扰模型的精度，设计更合理的优化条件，仍然是一个需要努力攻克难题。另外，如同传统CT研究中需要面对的稀疏角、有限角等情况下的图像重建，在探测器尺寸受限情况下的扫描视野扩大等问题，在能谱CT图像重建中仍然是需要进一步研究。

深度神经网络在CT图像重建中的成功依赖于大规模带注释的数据集，而标签采集工作量大且在很多情况下难以获得。为此，基于无监督学习方法的CT图像处理方法近几年发展起来[103-104]，典型的如Noise2Noise[105]、Noise2Self[106]、Noise2Void[107]、Noisier2noise[108]、Noise2Inverse[109]等。这些方法的基本思路是对于同一场景的不同图像，或者同一图像的不同图像块之间，如果这些图像或图像块中所含有噪声相互独立且均值为零，则可以彼此作为标签进行网络训练，实现降噪。对于去除无偏噪声这类任务，无监督学习就可以取得很好的效果。对于带有自身结构伪影的抑制，当前的无监督学习有较大局限性，还需要发展新的深度学习方法。

迁移学习(Transfer Learning，TL)是又一种弥补大规模标签数据集不足的方法。TL通过学习源任务和目标任务之间的共同特征或因素，将在先前领域/任务中学到的知识和技能推广到新任务或新领域[110-111]；最常见的是在大型注释数据集上预先训练模型，然后在更具体、更小的数据集上微调模型[112-113]。TL方法对于CT图像标签难以获得的这种情况非常适用。对于扫描对象基本固定的医学CT，在设备调试时本来就会用到一些与人体尺寸和组织类型相近的物理模体，用于标定CT设备，可以充分利用这些模体的实采数据去训练网络[114]。另外，可以用计算机仿真数据来进行迁移学习，通过将三维CT图像分割为特定材料(软组织和骨骼)、投影材料体积和模拟光子计数数据，并考虑扫描仪的能量响应，计算数值数据，然后用数值数据作为标签进行网络训练[89]。迁移学习中，数值数据必须接近于真实数据且具有代表性，否则效果不理想，还需要进一步研究数据来增强策略。

5 总结

提高图像空间分辨率、密度分辨率、时间分辨率、物质分辨率，降低X射线剂量，抑制噪声与伪影是CT研究者追求的目标与方向。相比于传统单能CT，能谱CT通过分解与重建被测物体的基材料图像，可提高对物体成分的识别和定量化分析能力，同时合成的虚拟单能图像可以有效减少图像伪影。能谱CT基材料分解算法已有数十年的发展，特别是随着近十余年来各种实用型能谱CT的出现，研究成果大量涌现，资料浩如星辰，难以一一列举。本文在介绍能谱CT的数据采集模式、基材料分解数学模型的基础上，分类介绍了能谱CT基材料分解算法及其优缺点、适用性等，并探讨了进一步的研究问题。

随着目前PCD能谱CT硬件系统的加速研制，以及深度神经网络方法的革命性应用，能谱CT仍可能是未来较长时间内的研究热点，高性能且对训练数据依赖少的算法仍然值得深入研究，希望本文内容对促进能谱CT的交流与研究有所裨益。

后记

在邱佩璋先生逝世一周年之际，我们仅以此文表达对先生的深切缅怀。邱佩璋先生是我国CT理论与技术研究的元老、开拓者，是《CT理论与应用研究》期刊的主要创始人，为我国CT理论与技术发展和专业人才的培养做出了巨大的贡献。首都师范大学CT团队创建人张朋曾师从于邱先生学习和研究，深受先生的教诲和恩泽。该团队的建设发展也得到了先生的长期关怀和指导，团队成员深受先生为人风范和治学精神的潜移默化的影响。