雷阳，连茂森，罗守华

东南大学 生物科学与医学工程学院，江苏 南京 210096

引言

在X线计算机断层扫描（Computed Tomography，CT）中，不同元素组成的物质可能在CT图像中表现为相同的像素值，区分不同的元素和组织存在困难。双能谱CT通过获取同一被测物体在两个不同能谱下的X线衰减数据，从而获得被测物体更多的信息[1]。物体能量相关的线性衰减系数可以近似表示为两种基物质衰减系数的线性组合[2]，从而实现不同组成物质成像。双能CT基图像重建算法主要分为图像域近似重建算法[3-5]和迭代重建算法[6-8]两大类。图像域近似重建算法是一种经验性的方法，先通过校准过程获取基物质衰减系数，再对重建后的图像在图像域进行基物质分解，相对容易实现且算法速度较快，但校准过程中的错误将导致硬化伪影和基物质分解错误。迭代重建算法遵照正投影过程中多色X光比尔-朗伯吸收定律的非线性积分模型来匹配投影数据，通过迭代方法重建出基物质图像，该类方法能够有效去除重建图像中的硬化伪影，且基物质分解准确性更高。Chen等[9]提出的基于优化的能谱CT重建算法（Adaptive Steepest Descent-Non Convex-Projection Onto Convex Sets，ASD-NC-POCS）将物质衰减系数分解成能量相关和能量无关两部分，通过正投影和反投影迭代过程匹配不同能谱的投影数据，该算法具有迭代重建算法的优点，并遵照多色X光比尔-朗伯吸收定律的非线性积分模型，具有理论上的优势。但是，由于算法迭代中多色X光前投影和反投影过程计算量大，应用于显微CT双能成像实际系统中存在困难。本文基于Chen等[9]方法，考虑到两个能谱下数据之间的相关性，在迭代过程中对其进行改进，并采用有序子集联立代数重建方法[10]替换Chen等[9]方法中代数重建方法，从而加快了算法收敛速度，使其能更好的应用在工程实践中。

双能CT在医学成像中的应用包括虚拟单能成像[11]，骨去除成像[12]、动脉粥样硬化斑块成像[13]和造影剂成像[14]等。近年来，特发性肺纤维化是临床上难以治愈的健康杀手[15]，科学研究发现基于干细胞的免疫疗法可帮助肺的再生[16]，对干细胞的显影和示踪具有重大价值。CT作为肺部检测的重要医学成像手段，具有较高的空间分辨率，在肺部成像方面有着显著优势[17]，可结合碘造影剂及双能成像算法的应用，以达到通过识别造影剂来间接示踪干细胞的目的。碘造影剂是CT成像和临床应用中最常用的造影材料。本文选择模型小鼠肺部碘造影双能成像，来验证算法在显微CT上的可行性。

相比于临床CT，显微CT系统采用微焦点X射线源和高分辨X射线探测器，其系统空间分辨率可达10 μm至100 μm[18-19]。由于射线源焦点和探测器像元尺寸较小，导致成像过程X射线总光通量较少[20]，为提高图像质量常选择较长的探测器曝光时间，并提高扫描角度数。现有的双能CT数据采集方法主要有两次连续扫描[21]、X射线球管电压快速转换[22]、多层探测器采集[23]、双源CT[14,24]和光子计数探测器采集等[25]。考虑显微CT采样特点，有效减少显微CT中两次扫描间隔时间以实现同源、同时、同向的图像采集条件实现困难，因此选择两次连续扫描方式来获得两个不同能谱的投影数据，以减少双能显微CT成像对硬件的需求。

综上所述，本文对Chen等[9]提出的ASD-NC-POCS方法进行改进，基于小鼠肺部碘造影剂识别，对显微CT双能成像展开研究。对于算法所需的能谱信息，为降低能谱重建问题的复杂性，本文利用单纯形模拟退火算法[26]实现基于参数模型的X射线源能谱重建[27]。

1 方法

1.1 ASD-NC-POCS算法

Chen等[9]提出了基于优化的双能CT重建算法。物质的线性衰减系数可以近似表示为两种其他物质衰减系数的线性组合：

其中，μk(E)是物质k能量相关的分解系数，即其衰减系数；bk是其空间位置相关的基图像。根据Lambert-Beer衰减定律，对于能谱为s的射线j，穿过物体的衰减系数积分(b)可表示为：

其中，表示归一化能谱，满足，m为能谱离散化后区间个数编号，i为投影物体三维空间体素号，为射线j与体素i的相交长度。

双能CT需要获取两个能谱投影数据，根据已知的物质分解系数μkm，反解式(2)，求解组成物质k的基图像bk。求解过程即最小化正投影模型计算数据(b)和实际采集数据之差的二范数D(g(b),gM)，并使得基图像的全变分（TV）最小来限制噪声，表示为：

将物质k的衰减系数μkm分解成能量相关和能量无关的两部分：

由X射线能谱与衰减系数μkm加权平均求得，和表示为：

从而式(2)中(b)可表示为线性(b)和非线性两部分：

获得之后，通过式(8)求得。重复上述迭代过程重建出基图像。

1.2 算法改进

Chen等[9]通过式(9)求得基图像，但根据实验结果，收敛较慢。为了减少迭代次数，提高收敛速度，对式(9)的迭代过程进行如下改进。双能CT图像域近似重建算法先用两个能谱下的投影数据直接重建得到重建图像，通过校准过程获取基物质衰减系数，再根据式(10)对重建后图像进行基物质分解求得基图像bk。

其中，μs1和μs2分别是两个能谱下重建图像像素CT值，μ1,s1、μ2,s1、μ1,s2和μ2,s2是校准得到的两种基物质衰减系数。受双能CT图像域重建算法启发，我们在凸集投影过程中考虑高低两个能谱下数据之间的相关性，联立两个基图像更新方程来更新基图像。首先按式(11)分别得到两个能谱下基图像更新残差值

和即式(10)中μs1和μs2，式(10)中两种基物质衰减系数矩阵用ASD-NC-POCS方法中基物质能量无关的衰减系数矩阵替代，通过反解式(10)更新基图像，即：

基图像的更新公式可表示为(13)所示：

为了进一步加快算法速度，本文用有序子集联立代数重建方法替换Chen等人方法中代数重建方法部分，将投影数据按投影角度划分成Nos个子集，考虑每个子集数据后对基图像更新一次，st表示投影数据的第t个子集。基图像的更新残差如式(14)所示。

伪代码如下，其中方框内为算法改进部分。

重复迭代

POCS 更新

TV 下降更新

for t = 1 to NTVdo

for k = 1 to K do

end for

end for

非线性项更新

for s = 1 to S do

for j = 0 to J[s]-1 do

鉴于党内歧见纷呈，蒋介石试图在思想上统合各派。1931年1月10日，蒋介石在日记中提到，“近日甚思研究哲学，觅一系统而研究之，以为收其放心之本也”⑪。1月15日，蒋又提到，“对军校训话，思有系统之讲演，如何能振起已死之人心”⑫。蒋介石认识到思想统一的重要性，认为，“现在一般青年的沉闷，社会的纷扰，统统是由于中国人思想不确定、不统一，由此可知思想之统一，比什么事情都要紧!”蒋进而提出以三民主义作为统一国民思想的唯一基础，“我们中国要在二十世纪的世界谋生存，没有第二个适合的主义，只有依照总理的遗教，拿三民主义来做中心思想，才能统一中国”⑬。

end for

end for

直到收敛条件满足

1.3 X射线能谱测量方法

目前，工程应用中主要通过选择合适的数学算法，根据衰减数据与入射能谱的对应关系实现能谱重建[28]。基于参数的钨靶能谱模型主要由阳极靶面产生的韧致辐射、球管对射线束的衰减及钨靶的特征辐射三部分构成[26]。基于参数的钨靶X射线能谱模型可表示为式(15)。

其中，a0～a3、h及c1、c2为 7个待求解的未知参数。在某个特定电压电流下，X射线源能谱重建转化为未知参数的求解。多色X射线穿过均匀物质的衰减值P如式(16)所示。

由式(15)(16)可计算出基于模型的多色X射线穿过均匀长度铝棒的理论衰减值。

实验选择图1所示均匀圆柱体铝棒作为能谱测量模体，其中红线部分表示X射线穿过模体长度。在某个特定电压电流下，对模体进行CT扫描并三维重建获得断层图像，阈值分割去掉图像噪声，使图像中铝棒所在区域像素值为1，其余区域为0。对阈值分割后的铝棒三维体数据正投影，从而获得若干组X射线穿过铝棒长度I和探测器测量衰减值的对应数据。

图1 用于能谱参数计算的测量模型示意图

为了拟合理论计算衰减值和实际测量衰减值，我们利用最小二乘原则构造式(17)所示的优化函数，选择单纯形模拟退火算法求解能谱模型中7个未知参数来最小化目标函数值，实现能谱重建。

2 仿真

2.1 能谱和基图像

本文选取水基图像和碘基图像，碘和水能量相关的物质衰减系数数据来自于NIST网站X射线衰减数据库。能谱数据选用Hiscan VM型显微CT系统40 kVp和80 kVp电压下测得的能谱。以1 keV为步长对碘和水的物质衰减系数采样插值，得到两种基物质与能谱对应的1～80 kVp能量下离散衰减值。

2.2 仿真模型

我们构造一个包含9个圆柱体的三维仿真模型，其材料密度如表1所示。图2a为仿真模型中平面断层，图像大小为700×700像素，图2b～c分别为水、碘基图像，准确地反映了模体中水和不同密度碘的分布。

表1 仿真模体材料密度

图2 仿真模体及重建图

材料 水 碘 碘 碘 碘 碘 碘 碘 碘密度 (g/mL) 1 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01

2.3 双能投影数据获取

假定某个特定电压电流下射线源发出的每条X射线能谱分布一致，在能谱范围内每个能量（keV）下，对仿真模体进行正投影，再将所有能量下正投影数据求和，得到一个角度下射线源40 kVp和80 kVp电压下正投影数据Ps，如式(16)所示。我们选取360°下的720个角度进行正投影，角度采样间隔0.5°。

2.4 仿真验证结果

通过FDK重建算法分别获得40 kVp和80 kVp电压下的重建图像，如图2d～e所示，红色箭头所指处存在由于CT成像过程中多色X光与FDK重建算法假设不一致导致的硬化伪影。采用Chen等[9]方法迭代8000次和本文改进方法迭代20次，来比较收敛速度及双能重建基物质成像效果，其中，γ取1，NTV取10，αk取0.2，Nos取90。基物质重建图像如图3a～b所示，为了便于显示，对基图像值乘以比例系数1000，并绘制了像素值沿图像中水平虚线的变化曲线，横坐标为像素号，纵坐标为像素值。从图像结果可知，Chen等[9]方法和本文改进方法都能有效去除硬化伪影。改进方法迭代20次后，水基图像中圆柱体1的像素值近似为1000，碘基图像中圆柱体4和8的像素值分别近似为60和20，与仿真模体材料密度一致。相比于Chen等[9]方法迭代8000次的结果，本文改进方法明显提高了收敛速度，显著减少了算法的迭代次数。进一步，选取收敛性指标(b(n))来衡量算法收敛性，选取基图像相似性指标(b(n))来衡量迭代过程中求得的基图像与真实基图像之间的误差。(b(n))和(b(n))分别如式(18)(19)所示，ε取10-8。图4为Chen等[9]方法和改进方法(b(n))和(b(n))关于迭代次数n的曲线，其结果与图3结果一致，进一步表明改进方法能有效提高收敛速度。本文算法在Visual Studio 2010环境实现，为了提高运算速度，采用CUDA 6.0编程实现算法的GPU并行加速。

图3 Chen等方法及改进方法结果图

图4 Chen等方法迭代8000次和改进方法迭代20次收敛性指标(b(n))和基图像相似性指标(b(n))曲线

3 实验验证

为了进一步验证该方法在实际显微CT系统中的可行性，基于Hiscan VM型显微CT，测量了该系统X射线球管能谱，并进行了小鼠肺部碘造影双能成像实验。

3.1 能谱测量结果

图5a为显微CT系统Hiscan VM型。该系统X射线球管为滨松微焦点X射线源L9421-02，管电压范围20～90 kVp，管电流范围0～200 μA，最大输出功率8 W，X射线微焦点尺寸为5 μm。探测器为Dexela的CMOS平板X射线探测器，阵列为1944×1536，动态范围0～16383，像元大小75 μm。在80 kVp电压下，我们选取纯度99%，直径为2 cm的铝棒模体，由于40 kVp电压下射线强度较弱，选取1 cm的铝棒模体。图5b为球管和80 kVp能谱测量曲线，横坐标为能量，纵坐标为归一化X射线光通量。图5c为40 kVp和80 kVp下不同长度铝棒实际测量X射线衰减值和基于能谱模型计算衰减值的拟合曲线，横坐标为铝棒长度，纵坐标为探测器采样值。模体实际测量X射线衰减值和基于能谱模型计算衰减值曲线基本重合，证明基于单纯形的模拟退火算法可获得目标函数一个较好的解。

3.2 小鼠肺部碘造影双能成像实验

（1） 实验动物及分组。8周龄雄性c57正常小鼠（22～25 g）12只，购自江苏省南京市青龙山动物繁殖场，采用随机数法随机抽取其中6只作为肺部碘造影成像模型小鼠。

（2）实验试剂及仪器。实验的主要试剂包括：1%戊巴比妥胺（麻醉剂）；70%无水乙醇；碘海醇溶液200 mg/mL。主要仪器包括：Hiscan VM型显微CT；眼科镊子-直/0.8 mm宽/10 cm；眼科镊-弯/0.5 mm/10 cm；Fine精细剪-直/尖头/9.5 cm。

（3）成像实验。本文采用气管滴注[29]的方式将造影剂注射到小鼠肺部。用1%戊巴比妥胺腹腔内注射麻醉，纵向剪开小鼠颈部，并撕开肌肉组织，从白色气管注射200 mg/mL的碘海醇0.15 mL；由于两次连续扫描间隔时间长，为了避免小老鼠的呼吸运动对结果的影响，立即用颈椎脱臼法处死，然后分别进行40 kVp和80 kVp电压下的两次连续扫描。

图6 小鼠肺部碘造影成像实验结果图

（4）小鼠肺部碘造影成像实验结果。图6为小鼠肺部碘造影成像实验结果。

小鼠肺部碘造影双能成像实验结果表明，改进的ASDNC-POCS方法能够应用在显微CT双能成像上，有效地从小动物体内分离出碘造影剂。而且，相比于Chen等[9]显著减少了算法迭代次数，提高了算法收敛速度。

图5 显微CT系统及能谱曲线图

4 讨论与总结

本文基于Chen等[9]提出的ASD-NC-POCS算法，在凸集投影过程中考虑高低两个能谱下数据之间的相关性，对其方法进行改进。为了验证改进方法在显微CT双能成像中的有效性，对其进行了仿真，验证了改进方法能够从含有水和不同密度碘材料的混合模体中分离出碘和水两种基物质，并有效去除了硬化伪影。同时，相比于Chen等[9]的方法有效减少了迭代次数，明显提高了收敛速度，从而能更好地应用于工程实践中。基于实际显微CT系统Hiscan VM，采用两次连续扫描的数据获取方法进行了小鼠肺部碘造影双能成像实验。结果表明，改进方法成功从小鼠肺部分离出碘造影剂，并有效地提高算法收敛速度，能够有效应用于显微CT双能成像上。

由于两次扫描间隔时间相对较长且小鼠肺部形变较快，同时CT系统在连续两次扫描中位置的不一致，会造成高低能投影数据不完全匹配，从而导致在基图像中物体结构的边缘处分离效果不太理想。另一方面，得到的碘和水的基图像的噪声水平相对较大。我们考虑在后续的工作中提高目前CT机械系统的角度精度，保证两次连续扫描获得图像的角度位置一致性，使两个能谱下同一角度图像相匹配；同时提高图像质量和降低图像噪声水平。