计算机断层成像(computed tomography, CT)技术应用十分广泛^[1].但尚存在一些不足, 比如金属伪影^[2].双能CT(dual-energy CT, DECT)能够显著降低金属伪影, 并且能够分辨被检测目标的组成成分, 大大改善CT的检测质量^[3-6].

DECT概念是1976年由Alvarez和Macovski提出的.物质对X射线的吸收系数是随X射线的能量变化的, 对于大多数物质而言, 这种变化的衰减系数可以由光电子效应和康普顿散射的线性组合精确近似.DECT通过两束能谱不同的X射线扫描目标物体, 得到两套投影数据.通过相应的重建算法, 得到光电子效应和康普顿散射的系数图.自DECT的概念提出后, 国内外研究人员对其重建算法进行了大量研究, 随能级变化吸收系数的线性组合也由原来的光电效应和康普顿散射转换为基物质分解^[7].目前国内的双能重建算法主要是查表法(look up table, LUT)^[3].查表法原理简单, 易于实现.但是, 为建立查找表, 需要事先知道目标物体的组成成分.此外, 对于不同的检测对象, 需要建立不同的查找表.因而, 实用性不太强.为此, 本文提出一种基于最大似然-可分离抛物面替代函数(maximum likelihood-separable paraboloidal surrogate, ML-SPS)的迭代重建算法.

本文首先介绍ML-SPS算法, 然后利用该算法对GATE得到的模拟数据投影进行重建实验, 并将结果和查表法的重建结果对比, 最后根据实验结果得出结论.

1 模型和算法 1.1 双能CT的物理模型和统计模型

被扫描物体对X射线的衰减系数随X射线能级大小的变化而改变.据此, Alvarez等提出了精确X射线衰减物理模型的数学表达式^[2]:

(1)

其中:I₀表示X射线的入射强度; I表示经过目标物体衰减后X射线的出射强度; μ_i表示射线穿过的目标物体的第i个体素的衰减系数; x_i表示射线穿过第i个体素的长度; s(E)表示入射的X射线能谱.

DECT重建中, 需要两个具有不同能谱的X射线扫描目标物体, 得到两套投影数据.一束X射线的能谱所包含的能级数远大于2.为了能够用两套投影数据重建DECT, 需要对衰减系数进行分解.常用的分解方法是双效应分解和双物质分解.本文采用双物质分解, 并取空气和碘作为基物质.如此, 任意一种物质的衰减系数都可以表示为空气和碘的衰减系数的线性组合, 即

(2)

其中：μ_I(E)和μ_A(E)分别表示碘和空气在能级E的衰减系数；c和d表示其线性组合的系数.将式(2)代入式(1)可以得到DECT重建的物理模型:

(3)

当探测器为光子计数型探测器时, CT的测量结果通常满足泊松分布^[5]:

(4)

此泊松分布的均值为

(5)

1.2 ML目标函数

假设DECT的扫描过程中, 目标物体包含有n个体素, 发射了m条X射线.记扫描系统的投影矩阵为

(6)

其中a_ij表示第i条射线穿过第j个体素的长度.

第i条射线穿过目标物体后测量得到投影值y_i的概率为

(7)

因此, 最大似然函数ML为

(8)

对式(8)取负对数并略去常数项:

(9)

综合式(5)、式(6)、式(9)可以得到最大似然函数:

(10)

其中:

(11)

(12)

C，D表示碘和水的线性组合系数，其余各变量与前文提到的相同.式(10)即为目标函数.

1.3 凸优化方法简化目标函数

式(10)所示的目标函数难以直接求解.为简化目标函数的优化过程, 需要对目标函数进行优化转换.目标函数中存在负指数和负对数的形式, 因此考虑采用凸优化的方法对目标函数进行简化^[8].

首先需要构建凸组合系数.凸组合系数需要满足两个条件^[8]: ①所有系数之和为1；②每个系数都是非负的.

目标函数(10)中既含有指数项又含有对数项, 因此, 对两个部分分别构建两个不同的凸组合系数.对于指数项部分, 构建凸组合系数:

(13)

其中, 上标k表示第k次迭代.对于对数项部分, 构建凸组合系数:

(14)

由式(13)和式(14)能够得出: 和，因此, 所构建的凸组合系数满足上述第一个条件.能谱s(E)、投影矩阵元素a_ij和衰减系数c_j^kμ_I(E)+d_j^kμ_A(E)都是非负的, 因此，α_ijE^k≥0, ∀j=1, 2, …, n; β_iE^k≥0, ∀E∈[E_min, E_max]，所构建的凸组合系数满足第二个条件.

根据凸函数的性质可得

(15)

由式(15)可见, 虽然通过凸组合系数的优化转换, 使得原来目标函数(10)中的对数项转换为一次线性函数, 但其指数项部分仍然存在指数项, 且该指数项中, 两个变量c_i和d_i仍然未得到分离, 依然难以求出其最小值点.为此, 再次利用凸组合系数进行变量分离.

构建凸组合系数:

(16)

(17)

显然γ₁^k≥0和γ₂^k≥0, 并且γ₁^k+γ₂^k=1.

利用凸函数的性质可以得到指数项的替代函数:

(18)

其中,

(19)

由式(18)可以看出, 所有的自变量都被分离开了，但目标函数中仍然存在指数项.想通过对目标函数求偏导数, 并令偏导数等于零的方法求解目标函数最小值仍然非常困难.本文提出了采用泰勒级数展开的方法构建替代函数.

1.4 泰勒级数展开进一步简化目标函数

替代函数相对于原来的函数需要满足三条性质, 才能保证利用替代函数求得的极值点能够保证原函数收敛.对于函数h(x), 假设在x^k处构建其替代函数h′(x, x^k), 这三条性质可表示为

(20)

(21)

(22)

记, 根据泰勒展开式得到的替代函数为

(23)

为了式(23)满足三条性质, 构建函数:g(x)=Q₃(x, x^k)-e^-x.容易求得:g(x^k)=0, 且.由于, 因此, ∀x∈X, .根据 , x^k为函数g(x)的最小值点, 且最小值为0, 即g(x)=Q₃(x, x^k)-e^-x≥0. Q₃(x, x^k)满足替代函数的所有条件.

把替代函数Q₂和Q₃代入到Q₁可以得到可分离的抛物面型替代函数(SPS).对替代函数各个自变量求偏导数, 并令偏导数等于零得到迭代公式:

(24)

(25)

其中，

(26)

1.5 ML-SPS算法伪代码

1) C:= 1; D:=1;

2) repeat

For j=1, n do:

式(24);

c′_j=t*c_j;

式(25)；

d′_j=t*d_j;

End for

3) C:=C’; D:= D’;

4) until{stopping criteria}

2 模拟实验结果与分析

为了验证ML-SPS算法的重建效果, 利用GATE模拟光子计数型探测器CT扫描的过程, 得到两组投影数据, 并利用这两组数据进行了重建.

2.1 模拟实验参数

GATE模拟中, 两束入射X射线的电压分别为80 kVp和140 kVp, 其能谱范围分别为40~80 keV和40~140 keV.被扫描物体是一个有6个等大的圆形区域的方形模体, 如图 1所示.

模体的6个圆形区域自左至右, 从上到下, 依次填充的是水、碳、钙、铝、硫和氯化钠.背景区域填充空气.模体边长为128 mm.

GATE模拟CT扫描过程的几何参数:射线源到模体中心的距离为550 mm, 探测器为直线型光子计数探测器, 通道数为128, 探测器到模体中心距离为90 mm, 共扫描360个角度, 角度间隔为1°.

2.2 重建结果分析

针对GATE模拟得到的投影数据分别进行了ML-SPS算法重建和LUT算法重建.其中, LUT算法建立的查找表的精度为10^-5, ML-SPS算法中投影矩阵的计算采用Siddon的算法^[9], 迭代次数设置为200次.图 2展示了两种方法重建得到的10 keV能级下的图像.

对比图 1和图 2可以看出，ML-SPS和LUT算法得到的重建结果与原始图像差别不大, 且两种算法的重建结果相似.

为了更为客观地对比分析ML-SPS算法的重建结果, 分别计算了ML-SPS和LUT算法相对于原始图像的相关系数(correlation coefficient, CC)和信噪比(signal noise ratio, SNR)^[10].由于DECT重建结果中包含多个能级的图像(比如文中能谱包含40~140 keV共101个能级的图像), 因此, 需要对每个能级的图像分别求取CC和SNR.为了便于分析所有能级的CC和SNR值, 以能级为横坐标, CC值和SNR值为纵坐标, 画出了CC值和SNR值曲线，如图 3和图 4所示.

由图 3和图 4可见, 对于所有能级, ML-SPS算法重建图像的信噪比和相关系数均略大于LUT算法.ML-SPS算法重建质量略优于LUT算法.

3 结论

1) 根据DECT物理模型和统计模型的特征, 得出ML目标函数, 并利用其凸函数性质, 通过凸优化和拟泰勒级数展开的方法构建了SPS, 将原目标函数的优化问题转换为线性问题, 并证明了其收敛性.

2) 所提出的ML-SPS算法模拟实验结果略优于传统的LUT算法.

3) 所提出的ML-SPS算法弥补了LUT算法需要事先根据被扫描物体的成分建立查找表的不足, 更具有实用性.