滕月阳, 郑孙易, 卢子鹏, 康雁
东北大学 中荷生物医学与信息工程学院, 辽宁 沈阳 110169
收稿日期：2016-06-30
基金项目：国家自然科学基金资助项目(61372014, 61302013)；辽宁省科技厅工业攻关及成果产业化项目(2014305001)。
作者简介：滕月阳(1979-), 男, 辽宁沈阳人, 东北大学讲师, 博士;
康雁(1964-), 男, 辽宁沈阳人, 东北大学教授, 博士生导师。

摘要：基物质分解是双能CT重建的重要步骤, 其中双物质分解是常用的分解模型之一, 该模型的核心关键是计算分解系数投影.为了更快计算它, 提出了基于误差反馈梯度下降的双能CT双物质分解算法和基于Armijo-Goldstein梯度下降的双能CT双物质分解算法.由于计算了梯度下降步长, 这两种方法能快速迭代求解基物质分解系数投影.同时他们有效地解决了双能CT重建的非线性问题.仿真实验结果显示, 与传统查表匹配法相比, 这两种算法稳定收敛, 计算速度快, 重建精度高, 对临床应用有重要的意义.在重建结果精度近似的情况下, 基于Armijo-Goldstein梯度下降的算法采用不精确线性搜索步长, 因此它的运行速度更快.
关键词：双能CT双物质分解梯度下降图像重建计算步长
Two-Material Decomposition Algorithm of Dual-Energy CT Based on Gradient Descent Method
TENG Yue-yang, ZHENG Sun-yi, LU Zi-peng, KANG Yan
School of Sino-Dutch Biomedical and Information Engineering, Northeastern University, Shenyang 110169, China
Corresponding author: KANG Yan, E-mail: kangyan@bmie.neu.edu.cn
Abstract: asis material decomposition is a very essential step in dual-energy CT (DECT) reconstruction and two-material decomposition is one of the most common model whose key point is to obtain the projections of decomposition coefficient. To improve the speed of it, two-material decomposition algorithms were proposed, which are the dual-energy CT based on the error feedback gradient descent method and the Armijo-Goldstein rule gradient descent method, respectively. These two methods were able to get the projections of decomposition coefficient quickly because of the computed step size in gradient descent. Moreover, the nonlinear problem in dual-energy CT reconstruction was also effectively and efficiently solved by using the proposed methods. Simulation results indicated that compared to the projection matching method, the two proposed methods can get stable convergence and high reconstruction precision with a short span of time, which has an important significance to the clinical application. With the same reconstruction precision, the algorithm based on the Armijo-Goldstein rule gradient descent is faster, using the inexact linear search step size.
Key Words: dual-energy CTtwo-material decompositiongradient descentimage reconstructioncomputed step size
计算机断层成像技术(computed tomography, CT)在医学诊断中起着重要的作用, 双能CT是新一代的CT系统.传统CT只利用X射线的衰减信息而忽略了能谱信息, 只能得到扫描组织的形态学结构图像^[1-2].而双能CT充分考虑了X射线的衰减信息和能谱信息, 能同时得到扫描组织的形态学结构图像和化学组织成分图像, 从而进行化学组织成分分析.
双能CT重建算法是双能成像技术的关键.现有的双能CT重建算法为了解决重建过程中的非线性问题, 需要复杂的优化处理过程, 计算时间长且精度低, 因此临床应用性并不强.本文提出的基于梯度下降法的双能CT双物质分解算法可以根据高低能投影数据循环迭代进而计算分解系数投影.该算法沿着一条X射线经过的路径测得高低能投影数据, 并使用梯度误差反馈和Armijo-Goldstein法则计算梯度下降步长，从而进行迭代求解对应的分解系数投影.Armijo-Goldstein法则属于不精确线性搜索方法, 因为每一次下降步长不要求达到精确最小, 所以迭代次数增加, 但是整体收敛的速度依然很快.在得到分解系数投影以后, 对应的分解系数可以通过传统的CT重建算法求出.最后实验重建各个能量下的单能量图像, 拟合得到局部区域物质的线性衰减曲线.该算法只需简单的梯度计算和步长计算即可加快迭代收敛, 提高算法的运行速度.并且该算法在保证运行速度的同时, 随着迭代次数的增加, 计算的精度也被不断地提高.
1 方法1.1 双能CT重建算法的原理双能CT重建模型主要为基物质分解模型^[3-5]:

(1)

式中，u(E)代表物质在能量E下的线性衰减系数.式(1) 假设任何一种物质的衰减系数都可由两种基物质的衰减系数线性表示.其中, u₁(E), u₂(E)是两种基物质在能量E下对应的衰减值, 该数据可以在NIST库^[6]中查询得到, b₁, b₂是两种基物质对应的分解系数.根据基物质分解模型, 沿着某一条X射线路径下有:

(2)

(3)

其中，B₁, B₂分别表示两种基物质的分解系数投影.根据混合能谱X射线条件下的朗伯比尔定律可得

(4)

(5)

其中：P_L, P_H代表沿着某条X射线路径下低能和高能投影值；S_L, S_H为X射线的低能和高能能谱数据.由式(4) 和式(5) 可知, 只要求出基物质分解系数投影B₁, B₂, 再通过传统CT重建算法即可得到基物质分解系数图像.最后根据式(1) 就可以重建出各个能量下对应的单能量图像.
1.2 基于梯度下降法^[7]的双能CT双物质分解算法该算法目标是求解基物质分解系数投影.其思想是寻找一个凸函数作为目标函数, 目标函数包括两个变量B₁, B₂, 然后计算出关于点(B₁, B₂)的梯度值, 使用梯度下降法来迭代更新B₁, B₂的值, 最后收敛得到最终解.假设X射线有效能量范围是40~140 keV, 把式(2) 和式(3) 离散化以后可以得到^[8]

(6)

(7)

然后构造一个目标函数, 其中P_L, P_H表示式(6) 和式(7) 的低高能投影, P_L, P_H代表实测低高能投影数据, 则得到以下目标函数:

(8)

该目标函数是一个凸函数, 如果(B₁, B₂)是最终解, 则式(9) 成立:

(9)

凸函数的局部最优解就是全局最优解, 而局部最优解可以通过梯度下降法求得.首先对目标函数f(B₁, B₂)关于B₁, B₂两个变量分别求偏导:

(10)

(11)

再对P_L, P_H关于变量B₁, B₂求偏导:

(12)

(13)

(14)

(15)

至此, 可以得到迭代计算过程中所需的梯度向量.
1.3 基于误差反馈计算梯度下降步长要最小化目标函数, 需要在搜索最优解的时候, 沿着梯度的负方向下降.这时, 收敛下降速度才是最快的.则

(16)

(17)

(18)

(19)

(20)

(21)

其中:k代表迭代次数; ?f(B^(k))为第k次迭代的负方向梯度向量; ?^(k)为第k次负方向梯度单位向量; err(B^(k))为第k次迭代误差; C为某个常数, 该值影响收敛的速度, 因此通常根据经验设定其值; ρ^(k)为第k次迭代的梯度下降步长.对于每一次迭代, 都要计算一次负方向的梯度向量, 进而求出负方向的梯度单位向量, 用ρ^(k)和?^(k)的乘积作为误差反馈, 来更新分解系数投影变量B₁, B₂.ε表示误差阈值, 当ε＜err(B^(k))时, 迭代终止.因此误差阈值ε的不同将会导致结果精度的差异.经过多次迭代更新, 分解系数投影变量B₁, B₂会收敛直到最终解.
1.4 Armijo-Goldstein法则^[9]计算梯度下降步长误差反馈计算得到的梯度下降步长可以使基于梯度下降的双能重建算法收敛到精确解, 但其下降步长并不是最适合的, 该方法存在收敛速度过慢的问题.因此, 本文提出了基于Armijo-Goldstein法则计算梯度下降步长的方法.
该方法需要设置一些变量, 设α表示所求的梯度下降步长, α_k表示第k次迭代值, α_max为所求下降步长的最大值; 给定常数; 根据式(10) 和式(11) 可以求出第k次迭代的梯度向量, 设d_k=?^(k)为单位梯度负方向法向量; 设s_k=α_kd_k, x_k=(B₁^(k), B₂^(k)); 设f(x_k)表示式(8).根据Armijo-Goldstein法则得出以下两个式子:

(22)

(23)

算法的步骤:
步骤1?在初始搜索区间[0, α_max]中确定一个初始点α₀, 令a₀=0, b₀=α_max, k=0;
步骤2?计算判断式(22) 是否成立, 成立则转步骤3;否则令a_k+1=a_k, b_k+1=α_k, 再转步骤4;
步骤3?计算判断式(23) 是否成立, 成立则输出α_k; 否则令a_k+1=α_k, b_k+1=b_k, 再转步骤2;
步骤4?令α_k+1=(a_k+1+b_k+1)/2, 并且令k=k+1, 再转步骤2.
按照以上步骤, 在梯度下降法中每迭代一次, 计算一次最合适步长α, 如此循环迭代计算, 最后收敛到精确值.
2 实验结果2.1 实验数据的仿真本实验的双能重建数据是通过GATE仿真软件模拟获得.GATE仿真软件是一个基于Monte Carlo仿真方法构建的核医学仿真平台, 能够真实地模拟双能CT系统的扫描过程.
实验使用的体模是放置在同一平面上的6个小球, 每个小球填充不同的单质或者化合物，如H₂O, Ca, Al, S等物质.然后实验通过GATE软件设置相关参数模拟双能CT扫描过程, 从而获取高低能投影数据.实验采用的高低管电压分别是140和80 kV, 扫描系统采用等距扇束扫描方式, 放射源到旋转中心距离为550 mm, 探测器中心到旋转中心距离为90 mm, 探测器通道数为128, 一共360个旋转角度.图 1为通过GATE仿真的双能投影数据.
图 1(Fig. 1)

图 1 投影数据Fig.1 Projection data (a)—高能；(b)—低能.

2.2 实验结果与比较本实验采用Matlab语言编程, 并在联想笔记本Y560(CPU:i3主频:2.1 GHz)上运行.一共使用了三种方法获取基物质分解系数投影, 分别为查表匹配算法、误差反馈梯度下降算法和Armijo-Goldstein梯度下降算法, 查表匹配算法采用的查找表大小为8 001×8 001, 设定的两种基物质分别为C和Al.
表 1比较了三种算法计算基物质分解系数投影的时间.
表 1(Table 1)

表 1 算法运行时间Table 1 Running time of algorithms s 算法 查表匹配 误差反馈 Armijo-Goldstein 时间 236 480.32 2 329.53 433.37

表 1 算法运行时间 Table 1 Running time of algorithms

由表 1可得, 相对于查表匹配算法, 误差反馈梯度下降算法使得运行时间减少了约100倍, Armijo-Goldstein梯度下降算法使得运行时间减少了约540倍, 算法运行效率大大提高.因为Armijo-Goldstein梯度下降算法采用不精确线性搜索步长, 所以相对于误差反馈梯度下降算法, Armijo-Goldstein梯度下降算法所花的时间更少.
经过以上三种算法计算得到的基物质分解系数投影, 实验使用FBP滤波反投影算法或者ART迭代算法^[10]即可得到基物质分解系数图像, 再根据式(1) 计算可得各个能量下的单能量图像.本实验设定的能量范围是40~140 keV.由于扫描体模内填充的物质及其衰减系数是已知的, 因此可以生成扫描体模在各个能量下理论的真实单能量图像, 然后将其和重建单能量图像进行比较从而验证算法的优良性.图 2以70 keV单能量为例显示了重建结果.
图 2(Fig. 2)

图 2 70 keV下的重建结果对比Fig.2 Comparison of reconstructed images with 70 keV (a)—真实单能量图像; (b)—查表匹配算法的单能量图像; (c)—误差反馈梯度下降算法的单能量图像；(d)—Armijo-Goldstein梯度下降算法重建的单能量图像.

从图中可以看出, 三种算法重建出的单能量图像和真实单能量图像基本相同.
在得到各个能量下的单能量图像以后, 可以在某个单能量图像上的每个小圆内随机选择一个像素点, 求出该像素点8邻域内每个像素点值的和, 再求平均值, 从而计算出每个小圆内填充的物质在该能量下对应的线性衰减系数.每个能量下都按此方法进行计算和存储, 即可描绘出该物质对应的线性衰减曲线.从NIST库中可以得到每种物质对应理论的真实线性衰减曲线.
本文以小圆填充物质之一的Ca和S为例进行比较.
根据图 3和图 4, 无论是Ca还是S, 误差反馈梯度下降算法的重建曲线和Armijo-Goldstein梯度下降算法的重建曲线都比查表匹配算法的重建曲线更加接近于真实的线性衰减曲线.这说明误差反馈梯度下降算法和Armijo-Goldstein梯度下降算法的重建结果精度比查表匹配算法的重建结果精度更高.
图 3(Fig. 3)

图 3 钙的线性衰减曲线Fig.3 Linear attenuation curves of Ca

图 4(Fig. 4)

图 4 硫的线性衰减曲线Fig.4 Linear attenuation curves of S

为了更加形象地反映重建精度的情况, 实验将每个能量下的重建值减去真实值, 之后取绝对值, 再除以真实值得到相对误差值, 进而得到物质线性衰减系数的相对误差曲线.
从图 5和图 6可以看出, 误差反馈梯度下降算法和Armijo-Goldstein梯度下降算法的重建结果相对误差和重建精度比较相近, 并且在大部分能量下, 他们的相对误差都比查表匹配算法的重建结果相对误差更小, 精度更高.
图 5(Fig. 5)

图 5 钙的线性衰减系数相对误差曲线Fig.5 Relative error of linear attenuation of Ca

图 6(Fig. 6)

图 6 硫的线性衰减系数相对误差曲线Fig.6 Relative error of linear attenuation of S

3 结论本文提出的基于梯度下降算法的双能CT双物质分解算法, 使用了误差反馈和Armijo-Goldstein法则来计算下降步长, 然后迭代求解基物质分解系数投影, 接着通过传统CT重建算法得到基物质的分解系数图像, 最后准确地重建出各个能量下的单能量图像.误差反馈梯度下降算法和Armijo-Goldstein梯度下降算法比查表匹配算法运行时间更快, 重建精度更高.Armijo-Goldstein梯度下降算法和误差反馈梯度下降算法相比, 它们的重建精度都很高且相差不大, 但是Armijo-Goldstein梯度下降算法采用不精确线性搜索步长的方法, 因此它的运行速度更快.
参考文献

[1]	Xing Y X, Zhang L, Duan X H, et al. A reconstruction method for dual high-energy CT with MeV X-rays[J].IEEE Transactions on Nuclear Science, 2011, 58(2): 537–546.DOI:10.1109/TNS.2011.2112779
[2]	An D C, Casselman J, Hoof T V, et al. Analysis of metal artifact reduction tools for dental hardware in CT scans of the oral cavity:peak kilovoltage, iterative reconstruction, dual-energy CT, metal artifact reduction software:does it make a difference[J].Neuroradiology, 2015, 57(8): 841–849.DOI:10.1007/s00234-015-1537-1
[3]	Petrongolo M, Dong X, Zhu L. A general framework of noise suppression in material decomposition for dual-energy CT[J].Medical Physics, 2015, 42(8): 4848–4862.DOI:10.1118/1.4926780
[4]	Zhao Y, Zhao X, Zhang P. An extended algebraic reconstruction technique (E-ART) for dual spectral CT[J].IEEE Transactions on Medical Imaging, 2014, 34(3): 761–768.
[5]	Jiang Y, Wang J, Wang J, et al. Clinical application of dual energy spectral CT method in patients with hepatic alveolar echinococcosis[J].Radiology of Infectious Diseases, 2016, 3(3): 120–127.DOI:10.1016/j.jrid.2016.07.003
[6]	Hubbell J H, Seltzer S M.Note on NIST X-ray attenuation databases[DB/OL].[2016-02-08].http://www.nist.gov/pml/data/xraycoef/.
[7]	Tripathi A, Mcnulty I, Shpyrko O G. Ptychographic overlap constraint errors and the limits of their numerical recovery using conjugate gradient descent methods[J].Optics Express, 2014, 22(2): 1452–1466.DOI:10.1364/OE.22.001452
[8]	赵云松, 张慧滔, 赵星, 等. 双能谱CT的迭代重建模型及重建方法[J].电子学报, 2014, 42(2): 1–2. ( Zhao Yun-song, Zhang Hui-tao, Zhao Xing, et al. Iterative reconstruction model and reconstruction method for dual energy computed tomography[J].Acta Electronica Sinica, 2014, 42(2): 1–2.)
[9]	Livieris I E, Pintelas P. A modified Perry conjugate gradient method and its global convergence[J].Optimization Letters, 2015, 9(5): 999–1015.DOI:10.1007/s11590-014-0820-0
[10]	Liu J, Ding H, Molloi S, et al. TICMR:total image constrained material reconstruction via nonlocal total variation regularization for spectral CT[J].IEEE Transactions on Medical Imaging, 2016, 35(12): 2578–2586.DOI:10.1109/TMI.2016.2587661