赵金龙， 刘 祎， 桂志国， 杨一鸣

(1. 中北大学 山西省生物医学成像与影像大数据重点实验室， 山西 太原 030051；2. 中北大学 信息与通信工程学院， 山西 太原 030051)

0 引 言

X射线计算机断层扫描(Computed Tomography， CT)成像是目前较为先进的无损检测技术， 其通过在不同角度下采集被测物体的投影信息并以一定重建算法最终重建物体内部结构， 现已被广泛应用于工业、 医疗诊断等领域[1]. 在生产生活中， 很多工件具有结构复杂、 材质密度差异大的特点， 对这类待测物体进行CT成像时， 常常由于探测器感应范围有限， 单一能量无法采集全部的投影信息而影响重建结果. 目前常用的改进方法主要是提高成像器件的动态范围， 但该方法对探测器系统、 采集传输系统有较高要求， 研制难度大， 成本高， 难以普遍推广[2]. 为此提出了变能量成像方法， 通过调整电压采集不同能量下的投影， 最终得到完整工件的投影信息以达到重建目的. 变能量成像方法最早由韩焱等在2008年提出， 其根据被测物体厚度调整X射线能量， 从而完成对被测物体的扫描[3]. 2010年， 段彦杰[4]验证了在利用变能量对工件扫描时， 会存在最佳灰度带. 陈平等[5]根据最佳灰度带进一步研究了变能量 X射线多谱成像技术， 成功地对复杂结构进行高动态成像， 证明了变电压CT成像的可行性. 对于递变能量 CT图像重建， 赵晋利[6]将递变电压下的投影融合后重建， 陈平、 李权等[7-8]以低能加权投影重建图像为初值对高能投影进行重建. 张雪英[9]对基于先验结构信息的变电压全变分-代数重建法(TV-ART)叠加迭代重建算法进行了研究， 该方法利用代数重建法(ART)依次对各能量下的投影进行重建直至最高能量， 随后以全变分(TV)算法优化图像. 为了弥补高能重建图像边缘结构的缺失， 将低能量下重建图像的边缘部分作为先验结构信息叠加到高能量 CT重建结果上， 得到物体的完整结构.

本文研究了基于先验图像压缩感知多能重建方法， 首先从低到高依次采集多个能量下的投影数据， 对最低能量下的投影采用凸集投影-全变分最小化(POCS-TVM)算法进行重建， 重建后的图像能较好地呈现厚度较薄区域的结构特征； 然后， 将重建图像作为先验信息， 利用先验图像约束压缩感知算法(PICCS)对下一个递加能量下的投影进行重建， 重建后的图像再次作为先验信息， 依次类推直到最高能量， 完成工件图像重建.

1 变电压投影序列

为了模拟密度变化较大的工件， 采用银、 铁、 铝3种材质， 建立如图 1 所示的仿真模型.

仿真模型为三层圆环结构， 外层材质为铝， 椭圆形， 长轴140 mm， 短轴88 mm； 中间层材质为铁， 长轴64 mm， 短轴48 mm； 内层材质为银， 圆形， 半径为8 mm. 分别采用30, 40, 60, 80, 100, 130, 150, 170, 200, 230, 250, 270, 300, 350, 400, 450, 500, 600 kV电压照射工件. 为了简化仿真实验， 采用平行束扫描方式并做如下假设： 源电压为80 kV时， 背景处的投影值达到探测器的上限， 设上限值为4 095， 即探测器感应范围为[0,4 095]， 其余投影等比折算. 设定探测器个数为200， 投影角度为0°～180°之间均匀分布的180个角度. 按上述条件得到部分电压投影图如图 2 所示.

图 1 仿真模型Fig.1 Simulation model

图 2 部分电压投影序列Fig.2 Projection sequence of partial voltage

由图 2 可知， 某一电压下的投影会出现欠曝光或过曝光或欠曝光与过曝光共存的现象， 导致单一电压照射下无法获得工件的完整投影， 影响图像重建的准确度. 电压较低时， 射线可以穿透工件较薄区域， 并根据厚度不同呈现不同的灰度等级， 从而采集到这一部分的投影信息. 而对于较厚或密度较大区域则由于X射线能量不足而无法穿透， 探测器无法接受到该部分的射线， 投影图像的灰度值为0， 会在投影图上呈现黑色区域； 随着电压逐渐增大， 较厚或密度较大的区域逐渐被穿透， 投影图中的黑色区域逐渐消失， 较厚或密度较大区域的投影信息逐渐被采集到， 而较薄区域的投影则由于过曝光而失去价值.

文献[4]证明了在不同电压下的投影存在一最佳灰度带， 该灰度带上的成像效果最佳. 文献[6]表明最佳灰度带不是唯一的， 在决定最佳灰度带时只要确保相邻投影有一定的重复率， 即可达到完整重建的目的. 因此需要确定一最佳灰度带并对各电压下的投影进行截取， 截取后的投影应既不包含欠曝光也不包含过曝光， 同时投影应涵盖整个工件的结构信息. 本次仿真实验设定最佳灰度带为[200，3 500]. 为了更好地观察各电压下的有效投影， 当投影的灰度小于200或大于3 500时， 设其值为探测器感应上限值， 即4 095[10]， 此时无效投影区域为白色背景， 部分电压投影图如图 3 所示.

图 3 部分电压有效投影序列Fig.3 Effective projection sequence of partial voltage

随着电压的逐步增强， 模型的有效投影信息(灰黑色部分)逐渐完备， 最终包含整个工件的结构信息.

2 基于先验图像压缩感知多能CT重建

2.1 基于先验图像的压缩感知算法

2004年， D．Donoho等提出了革命性的理论 “压缩感知”(CS)[11]. 该理论指出， 当图像是稀疏的， 即使投影是高度欠采样， 仍可以较好地重建图像. 很多图像本身并不具备稀疏特性， 因此在图像重建前对 CT图像进行稀疏变换， 如有限差分变换， 将有限差分的 L1范数即全变分(TV)作为最优化问题的目标函数.

将图像fi,j排列为一维向量F， 并以非负条件对图像进行约束， 使像素值都为正值， 得到

min‖F‖TV,s.t.P=AF,F≥0,

(1)

式中：P为投影值；A为投影系数矩阵.

2008年， Chen等[12-13]在CS的基础上提出了先验图像压缩感知算法(PICCS)， 该算法引入先验图像F对重建图像进行约束. 先验图像与重建图像在结构上具有较高的相似性， 当重建图像与先验图像相减， 显然得到的图像(F-Fp)更加稀疏[14-15]， 于是在式(1)的基础上得到新的目标函数

min‖F-FP‖TV,s.t.P=AF,F≥0.

(2)

由于先验图像可能存在条纹伪影， 如果不进行仔细的处理， 这些伪影将出现在最终的图像中. 因此， 在式(2)的目标函数基础上增加一个附加项min‖F‖TV， 该附加项可以有效消除伪影， 并使用权重系数α平衡两者的关系， 一般取值0.5[16]. 最终目标函数为

min[α‖F-FP‖TV+(1-α)‖F‖TV]

s.t.P=AF,F≥0.

(3)

对于上式这类最优化问题求解采用梯度下降法[17]， 其中约束条件用凸集投影(POCS)算法实现， 将凸集投影和全变分最小化相结合的方法称为 POCS-TVM算法. POCS-TVM算法主要分为非负约束条件下的 ART算法(POCS)以及梯度下降法求解全变分最小化过程(TVM)， 全变分梯度公式为

fi+1,j-fi,j(fi+1,j-fi,j)2+(fi+1,j-fi+1,j-1)2.

一般地， POCS-TVM算法分为2个步骤：

步骤 1： POCS过程

2) 用ART算法完成一次运算

3) 非负约束

步骤 2： TVM 过程

2) 计算梯度下降搜索增量因子

3) 计算全变分梯度

Gn(k)=

及梯度方向

式中：Fp为先验图像;n为TVM迭代次数，n=1,2,…,N(N为TVM最大迭代次数).

4) 以梯度反方向修正图像

n=n+1,

式中：φ为调节因子. 反复更新上式并修改图像， 直到n=N为止.

5) 判定是否满足迭代终止条件

式中：T为一极小的正实数. 若上式不成立， 则将结果作为POCS初值进行下一次算法循环

k=k+1.

2.2 多能CT重建

任意一单电压下的有效投影都只包含工件的一部分信息， 利用该电压的有效投影单独重建时也只能重建该工件的一部分. 为了能够完整地重建整个工件， 采用多个能量投影进行重建. 图 4 为200 kV和230 kV两个相邻能量的重建图像.

图 4 相邻电压重建图像Fig.4 Reconstruction image of adjacent voltage

由图 4 可以看出， 两个重建图像在相对较薄的部分具有较高的相似性. 因此本文对传统的ART-TV变电压重建进行改进. 根据PICCS理论， 若将图4(a)作为先验图像Fp， 并以式(1)作为目标函数对图4(b)进行重建， 将有效改善重建图像质量. 同时， 先验图像也包含了高能重建缺失的部分低密度边缘信息， 将先验图像纳入重建过程将对边缘起到一定保护作用. 为此， 将该方法扩展到多个能量， 即本文的基于先验图像压缩感知多能重建方法， 具体过程如图 5 所示.

图 5 重建示意图Fig.5 Reconstruction diagram

利用POCS-TVM对最低电压下的有效投影进行重建， 得到工件的轮廓图像， 将该图像作为先验信息， 利用PICCS对相邻高电压有效投影进行重建， 以融合两个能量下的投影信息. 将重建结果再次作为先验图像， 重复上述步骤直至最高电压. 重建完毕后， 整个工件的信息将融合在同一CT图像中. 重建算法具体步骤如下：

1) 逐步提高X射线源电压并依次对工件进行照射， 采集到N组投影数据In,n=1,2,…,N；

2) 观察步骤1)采集到的投影数据灰度直方图， 确定最佳灰度范围并对投影数据截取， 将截取后的投影做如下对数变换[18]

式中：I0为相应的入射能量.

3) 选取最低电压下的投影， 以POCS-TVM算法对其进行重建得到重建图像Fp；

4) 将Fp作为先验图像， 其对应的投影为第n组， 利用PICCS对n+1组投影数据进行重建得到重建图像F；

5) 判断是否达到最高能量， 若达到最高能量输出重建图像F， 否则令Fp=F， 重复步骤4).

3 实验仿真及分析

3.1 实验结果及分析

根据上文所述方法进行实验， 得到部分结果如图 6 所示.

图 6 部分电压重建结果Fig.6 Reconstruction results of partial voltage

起始电压为30 kV， 上限电压为60 kV时， 如图 6(a) 所示， 上限电压较低， 密度较大和较厚的区域X射线未能穿透， 无法采集投影信息， 因而不能对该部分进行重建. 因此， 此处像素的灰度值仍为初始图像的默认值0， 在图中呈现黑色. 而工件较薄的部分可以被X射线穿透并采集到有效投影， 利用这些投影数据可以有效地重建工件薄边缘. 随着能量的增大， X射线的穿透能力增强， 较厚和密度较大区域逐渐被穿透， 相应的这一部分的投影信息也逐渐被采集完全， 高密度和较厚的区域得以被重建， 如图6(b)～(f)所示. 当上限电压达到最大值600 kV时， 如图6(f)所示， 工件的所有结构信息都已被采集到， 此时的重建图像可以较好地展现工件内部结构. 重建过程将相邻低电压作为先验图像， 有效地改善了重建图像质量并保护了低密度边缘投影信息.

3.2 实验对比分析

在相同的实验条件下， 利用文献[7]和[9]所述的重建方法进行重建， 并与本文所研究的算法进行比较， 结果如图 7 所示.

图 7 不同重建算法CT图像对比Fig.7 Comparisons of CT images with different reconstruction algorithms

图7(a)为文献[7]提出的灰度加权算法重建图像， 随着电压的逐步增大， 灰度加权算法虽然完整地展现了工件的内部结构， 但是由于高能量投影缺少低电压边缘信息， 随着迭代算法的反复进行， 出现边缘与背景融合的现象， 难以观察到被测物体的边缘结构. 图7(c)算法中由于将相邻低电压重建图像作为先验图像， 为高能重建提供了部分丢失的低密度边缘信息， 使边缘得到了有效的保护， 最后的重建结果中低密度边缘仍较为清晰. 图7(b)为文献[9]提出的变电压迭代算法， 由于每个电压下的投影都是不完备的， 射线的缺失导致重建后的图像出现较多的伪影. 虽然该算法在迭代至最高电压时采用TV最小化来抑制伪影， 但由于这些伪影是多个电压重建累积下来的， 很难有效地消除， 另外先验结构的截取也是影响本方案的一个重要因素. 基于先验图像压缩感知多能重建算法将先验图像与重建图像相减， 增加了图像的稀疏性， 即使采集到的投影不完备依然可以较好地重建图像， 有效地避免了伪影的产生和累积， 从图7(c)的重建图结果可以看出， 因射线缺失造成的伪影得到了有效的抑制.

为了更加客观地对比3种算法， 从归一化均方距离、 最坏情况距离[19]、 结构相似度三方面进行比较， 结果如表 1 所示. 归一化均方距离、 最坏情况距离值越小， 结构相似度值越大， 说明误差越小， 重建效果越理想.

表 1 评价参数对比

本文算法相对于灰度加权算法在归一化均方距离和最坏情况距离方面分别减少了50.97%和20.94%， 相对于变电压迭代算法分别减少2.7%和37%， 结构相似度相对于灰度加权算法提升了17%， 比变电压迭代算法略低. 但综合考虑3个评价参数以及图像视觉效果， 本文算法仍使图像质量得到了较大的改善. 图 8 为原图以及各算法重建图像第70行像素灰度变化情况.

图 8 灰度曲线对比图Fig.8 Gray curve comparison diagram

由图 8 可见， 本文研究的重建算法图像灰度曲线更加平滑平稳， 且不同物质灰度对比较明显， 特别是工件区域(箭头所指)更接近原图.

4 结 论

针对结构复杂、 内部材质密度相差较大的工件， 单一能量无法准确重建的问题， 研究了基于先验图像压缩感知多能重建算法， 将低能量的重建图像作为先验图像运用到高能量CT重建中. 通过仿真实验表明， 本文算法可以有效地保护低密度边缘， 并减少因密度差异大、 单一能量投影信息不完全造成的伪影， 重建图像质量得到了改善.