同步辐射CT投影数据中的带状伪影及重建误差分析*

孔慧华, 杨玉双

(1. 中北大学 理学院， 山西 太原 030051；

2. CSIRO澳大利亚联邦科学与工业研究组织， 澳大利亚 Victoria 3169；

3. 山西大学 理论物理研究所， 山西 太原 030006)



同步辐射CT投影数据中的带状伪影及重建误差分析*

孔慧华1, 杨玉双2,3

(1. 中北大学 理学院， 山西 太原 030051；

2. CSIRO澳大利亚联邦科学与工业研究组织， 澳大利亚 Victoria 3169；

3. 山西大学 理论物理研究所， 山西 太原 030006)

摘要:同步辐射计算机断层成像技术集合了同步辐射光源X射线和无损检测的优越特性. 针对同步辐射光源在高光子能量下光源强度具有时间波动性和空间非均匀性的问题， 分析了光源不稳定性对投影造成的伪影， 进而研究了伪影导致的重建误差. 研究结果表明： 光源的时间波动性会造成投影数据中出现横向带状伪影， 降低重建图像的分辨率； 光源的空间非均匀性会造成投影数据中出现纵向带状伪影， 不仅会降低图像的分辨率， 还会使得重建图像中出现环状伪影. 数值模拟和实际实验结果验证了分析的正确性.

关键词:同步辐射CT； 光源强度的不稳定性； 带状伪影； 重建误差

0引言

同步辐射计算机断层成像技术(Synchrotron Radiation Computed Tomography， SR-CT)将同步辐射光源结合到CT技术中， 是一种用于检测物体内部结构的方法， 具有高分辨率、 无损性等独特的优越特性， 该技术已经广泛应用于生物医学、 材料科学等许多领域[1-5].

SR-CT与传统微焦点CT在原理上是完全一致的， 基本思想都是根据投影重建图像. 它们在组成装置上的主要差别在于光源的不同， 同步辐射光源利用电子在圆形轨道中的运动， 发出单色光， 具有高强度、 高分辨、 高准直、 强穿透等特性， 其光源可近似为平行光束， 而传统CT装置中的X射线是具有一定发散角的扇形光束. 另外同步辐射光源的频谱范围也很宽泛， 包含红外、 可见光直到X光等各种波长的光， 实验可根据需要选择合适的波长， 而传统CT光源为多色光， 难以实现波长的选择[6].

同步辐射光源潜在的试验问题是强度的稳定性不好[7-8]， 这与同步辐射光源的短暂性能有关， 如储存环中电子流的变化和轨道漂移明显影响入射线的强度. 由于同步辐射光源的不稳定性使得测量数据存在误差， 导致投影出现伪影， 进而影响重建图像的分辨率， 使得研究人员无法对研究对象的微细结构做出正确的判断. 通常来说光子能量越大， 光源的不稳定性越大， 投影中的伪影越明显. 本文首先分析高光子能量下同步辐射光源的不稳定性导致投影数据中出现的伪影问题， 进而分析此问题引起的重建误差.

1计算机断层成像原理

设待重建物体的衰减系数为f(x,y)， 强度为I0的X射线在某一方向l， 行进L距离后衰减至I， 则由比尔(Beer)定理[9]可知总衰减

(1)

式中：I0为入射光强度；I为透射光强度；p为射线投影. 计算机断层成像技术就是从一系列的投影p重建物体的衰减系数f(x,y).

图 1 为计算机断层成像原理示意图[10],Oxy代表坐标系，Oxryr为原坐标系顺时针旋转θ角度后的坐标系. 设待重建图像f(x,y)在视角θ下与yr轴距离xr处的投影为p(θ,xr)， 则根据滤波反投影的重建公式[11]有

(2)

图 1　反投影过程Fig.1　Back projection process

2数据的获取及伪影分析

2.1数据的获取

由前面的讨论可知， 测得入射光强度I0与透射光强度I， 即可得到投影p. 下面介绍同步辐射过程中I0及I的获取过程.

图 2 给出了同步辐射X射线显微CT示意图[12].

图 2　同步辐射X射线显微CT示意图Fig.2　Schematic diagram for synchrotron radiation X-ray micro CT

如图 2 所示， X射线保持不动， 样品在旋转台上做旋转运动， CCD探测器在样品后面连续采集数据. 投影成像的光源即为同步辐射光源. 在SR-CT技术中， CCD探测器在样品后面连续采集的数据， 即透射光强度I， 通常称为透射图像(tomography image). 实际中很难直接获得同步辐射光源入射物体前的强度I0， 一般通过背景法来获取， 即在未放置物体前存储一幅同步辐射光源的背景图像， 称为亮场图像(flat-field image), 以此获取入射光强度. 为了减小误差， 通常是间隔一定的角度就拍摄一次亮场图像. 另外为避免系统噪声对重建图像的影响， 还会在数据采集结束前， 拍摄不开光源时的暗场图像(dark-field image). 有了这三种数据就可以通过计算得到各个层面上的投影(sinogram image).

sinogramimage=

(3)

2.2伪影分析

在给出数据获取过程之后， 下面分析同步辐射光源的不稳定性对投影数据的影响.

图 3　同步辐射X射线显微CT坐标系统示意图Fig.3　Schematic diagram for synchrotron radiation X-ray micro CT coordinate system

实验时， 由于同步辐射光源在时间上的波动性及在空间上的非均匀性， 会使得用于计算的入射光强度不等于实际的入射光强度， 从而使得投影数据存在伪影， 继而影响重建结果. 同步辐射光源的时间波动性是指同一区域上的光强在不同时刻是不一样的， 其空间非均匀是指同步辐射光源在同一时刻， 不同区域上的光强是不同的， 体现为CCD探测器在不同位置(a,b)接收的光源强度不同.

将入射光强度和透射光强度分别记作I0(t,a,b),I(t,a,b)； 变量t表示不同时刻； 变量a表示探测器平面的横向位置； 变量b表示探测器平面的纵向位置， 这里代表被检测物体的不同层面. 由于光源随位置b而引起的变化主要体现在不同层接收到的光源强度不同， 从而导致物体某一层面的投影图像整体偏亮或者偏暗. 本文主要讨论时间t及位置a的变化引起的投影伪影问题， 即被检测物体某一层面的投影中的伪影， 设b=b0.

2.2.1时间不稳定性引起的伪影

首先忽略光源强度的空间非均匀性， 讨论光源的时间波动性， 即假设光源强度在同一时间不随a的变化而变化. 设未放置物体时的光强为I0(t0,a,b0)， 对于透射光强度， 由于不同角度下的测量数据需要按顺序先后依次测量， 设测量θ1,θ2,…,θn角度下透射光强度的时间分别为t1,t2,…,tn， 将其分别记为I(t1,a,b0),I(t2,a,b0),…,I(tn,a,b0)， 于是相应的不同角度下的投影可以表示为

而实际的投影应该为

当I0(t0,a,b0)≠I0(ti,a,b0)时， 势必会给投影带来误差

(4)

2.2.2空间不稳定性引起的伪影

下面忽略光源强度的时间波动性， 考虑空间非均匀性， 讨论由位置a的变化引起的伪影. 设ti时刻相邻位置a1和a2的光源强度分别为I0(ti,a1,b0)和I0(ti,a2,b0). 为了便于分析， 假设物体的衰减系数为常数μ(即f(x,y)≡μ)， 则在a1和a2位置的投影分别为

(5)

(6)

从而a1,a2位置处的投影差异为

(7)

在f(x,y)≡μ的假设下， Δp1(θi,a1,b0)是由射线穿过物体的长度所造成的， 是正常差异. 在理想状态下， 同一时刻的光源是均匀稳定的， 即I0(ti,a1,b0)=I0(ti,a2,b0)， 从而Δp2(θi,a1,b0)=0. 当辐射光源不满足空间均匀性， 即I0(ti,a1,b0)≠I0(ti,a2,b0)时， Δp2(θi,a1,b0)≠0， 该项误差是由光源强度的空间不稳定性造成的， 当θ取遍所有角度θ1,θ2,…,θn时， 在投影的纵向a1位置处就形成了带状伪影， 带宽与投影的空间变化程度有关. 当然随着同步辐射技术的提高， 光源的空间不稳定性越来越小， 但在实验中CCD探测器某个部位的不灵敏性或空气中灰尘依附在探测器上造成的影响等同于光源的空间不稳定性.

从以上分析可以看出， 同步辐射光源的时间波动性和空间非均匀性会分别造成投影数据中的横向带状伪影和纵向带状伪影. 下面分析这两种伪影导致的重建误差.

3重建误差分析

这里仍然只考虑某一层面投影的变化， 即b=b0. 在第2节的假设下， 滤波反投影的重建公式(2)还可以表示为

(8)

其离散化形式为

(9)

(10)

于是有

(11)

由滤波反投影重建公式可知， 对每一个θi， 把[Δp(θi,a,b0)*h(a)] 反投影于满足a=xcosθi+ysinθi的直线上的各点， 如图 1 中虚线所示， 滤波投影的值沿着整个直线路径被“涂抹”. 当在θi角度出现误差时， 误差被沿着a=xcosθi+ysinθi的直线“反抹”， 从而使经过这些直线的点的值出现增加或减小， 即图像在某处的值增加或减小， 从而降低图像的分辨率.

当投影中出现纵向带状伪影Δp2(θi,a1,b0),θi∈{θ1,θ2,…,θn}时，

(12)

由滤波反投影公式可知， 对每一个θi， 把Δp2(θi,a1,b0)*h(a1) 反投影于满足a1=xcosθi+ysinθi的直线上的各点. 由此可知， 0°～180°内， 误差被“涂抹”于所有与原点距离为a1的直线上的点， 从而会形成一个以原点为圆心，a1为半径的半圆环， 圆环的宽度反映了纵向伪影的宽度.

4实验

4.1模拟实验

基于同步辐射光源的不稳定性引起的重建误差的理论分析， 对上述两种误差形式分别进行数值模拟运算. 测试模型为Shepp-Logan模型， 大小为256×256像素， 投影采样间隔为1°， 采样区域为0°～180°， 对投影数据分别加入横向带状伪影和纵向带状伪影， 横向带宽设为20个像素， 误差分别为-5, -3, -1, 1, 3, 5, 7, 9, 11, 13, 15, 13, 9, 7, 5, 3, 1, -1； 纵向伪影位置为91～93像素， 带宽为3像素， 误差为10. 图 4 给出了加入伪影的投影数据. 图 5 是滤波反投影算法的重建结果图. 结合重建误差分析可知， 纵向带状伪影会将误差反抹于与原点距离为35～37像素的直线上， 误差较大的值出现在以原点为圆心， 半径为35～37像素长的半圆环上， 圆环的宽度为3像素， 反映了纵向伪影的宽度， 如图5(c)所示. 相比于纵向带状伪影， 横向带状伪影的误差主要是降低了图像的分辨率， 如图5(b)所示.

图 4　投影数据Fig.4　Projection data

图 5　重建结果图Fig.5　Reconstruction result picture

4.2实际实验

为进一步验证分析的正确性， 下面将通过实际采集的实验数据说明上述两种情况.

实验在上海同步辐射光源的X射线成像及生物医学光束线(BL13W)开展. 该光束线采用扭摆器(Wiggler)光源， 通过双晶单色器将白光转化为单色光， 光子能量可调范围为8～72 keV. 实验所用的岩石样品高20 mm， 直径3 mm， 取自中石油鄂尔多斯姚店油田的一块高200 mm， 直径100 mm 的致密砂岩样本. 实验中采用高分辨率X射线CCD， 像素尺寸为3.7 μm， 转台精度为0.000 1°.

为保证X射线对样品具有较好的穿透性， 根据实验经验， 选取光子能量分别为25 keV, 35 keV, 45 keV， 实验中选取样品到探测器的距离为4 cm， 0～180°内的角度采样数为720. 不同能量下的曝光时间分别为1 s， 1.5 s和8 s. 采集图像(亮背景、 暗背景和透射图像的尺寸为1 377×721 像素.

所有采集到的数据通过X-TRACT 软件[12]前处理(暗背景和亮背景的校正、 图像归一化、 相位恢复和伪影校正)后， 得到每层图像的投影数据， 图 6 给出了第103层(第一行)和第360层(第二行)的投影数据.

由图6可以清楚地看到， 同步辐射光源的时间不一致性所导致的投影数据中的横向带状伪影， 电压越大， 伪影越细密明显. 同时还可以清楚地看到， 在第103层的投影数据中存在纵向带状伪影. 另外， 第103层切片的投影的横向带状伪影明显比第360层的伪影严重， 这说明同步辐射光源强度随空间位置的变化而变化， 且中心位置的变化较边缘位置的变化小.

图 6　样品的投影Fig.6　Projection of sample

利用X-TRACT 软件对投影数据进行重建， 图 7 给出了与图6所对应的投影的重建结果图， 第一行是不同电压下第103层的重建结果， 第2行是第360层的重建结果， 由此可以看出， 25 keV时， 投影中的纵向带状误差造成了明显的环状伪迹， 而在35 keV和45 keV时， 由于横纵向带状伪迹同时存在， 降低了图像的分辨率， 使得环状伪迹有些模糊.

图 7　样品的重建结果Fig.7　Reconstruction result picture of sample

5结论

本文主要分析了同步辐射光源的时间和空间不稳定性引起的投影数据和重建结果的误差， 并通过数值模拟给出了这两种情形引起的误差， 最后通过实际的实验数据验证了分析的正确性. 目前作者正在研究如何降低重建结果误差.

参考文献:

[1]Lareida A, Beckmann F, Schrott-Fischer A, et al. High-resolution X-ray tomography of the human inner ear: synchrotron radiation-based study of nerve fibre bundles, membranes and ganglion cells[J]. Journal of Microscopy, 2009, 234 (1): 95-102.

[2]Reiche I, Müller K, Staude A, et al. Synchrotron radiation and laboratory micro X-ray computed tomography-useful tools for the material identification of prehistoric objects made of ivory, bone or antler[J]. Journal of Analytical Atomic Spectrometry, 2011, 26: 1802-1812.

[3]Makoto A, Kentaro U, Masato H, et al. In vitro assessments of white-spot lesions treated with NaF plus tricalcium phosphate(TCP) toothpastes using synchrotron radiation micro computed tomography (SR micro-CT)[J]. Journal of Dentistry and Oral Hygiene, 2014, 6 (1): 10-21.

[4]Cao Yong, Hu Jianzhong, Wu Tianding, et al. 3D angioarchitecture of spinal cord in a rat model detected by synchrotron radiation micro-computed tomography[J]. Journal of Orthopaedic Translation, 2014, 2(4): 241-242.

[5]Borstnar G, Mavrogordato M N, Helfen L, et al. Interlaminar fracture micro-mechanisms in toughened carbon fibre reinforced plastics investigated via synchrotron radiation computed tomography and laminography[J]. Composites: Part A. 2017, 71:176-183.

[6]汪敏. 同步辐射CT技术研究及应用[D]. 合肥： 中国科学技术大学, 2006.

[7]汪敏, 胡小方, 伍小平. 同步辐射计算机断层技术衬度误差机理分析[J]. 物理学报, 2006, 55 (8): 4065-4069.

Wang Min, Hu Xiaofang, Wu Xiaoping. Analysis of contrast error mechanism for symchrotron radiation computed tomography technique[J]. Acta Physica Sinica, 2006, 55 (8): 4065-4069. (in Chinese)

[8]汪敏, 岑豫皖, 胡小方, 等. 同步辐射计算机断层技术光源误差机理分析[J]. 物理学报, 2008, 57 (10): 6202-6206.

Wang Min, Cen Yuwan, Hu Xiaofang, et al. Error mechanism of light source for synchrotron radiation computed tomography technique[J]. Acta Physica Sinica, 2008, 57 (10): 6202-6206. (in Chinese)

[9]Yu Hengyong, Xu Qiong, He Peng, et al. 基于Medipix的谱微-CT[J]. CT理论与应用研究, 2012, 21(4): 583-596.

Yu Hengyong, Xu Qiong, He Peng, et al. Medipix-based spectral micro-CT[J]. Computed Tomography Theory and Application, 2012, 21(4): 583-596. (in Chinese)

[10]谢强. 计算机断层成像技术: 原理、 设计、 伪像和进展[M]. 北京： 科学出版社， 2006.

[11]Kong Huihua, Yu Hengyong. Analytic reconstruction approach for parallel translational computed tomography[J]. Journal of X-ray Science and Technology, 2015, 23: 213-228.

[12]王玉丹, 彭冠云, 佟亚军, 等. 影响同步辐射X射线螺旋显微CT的若干因素研究[J]. 物理学报, 2012, 61(5): 054205-1-9.

Wang Yudan, Peng Guanyun, Tong Yajun, et al. Effects of some factors on X-ray spiral micro-computed tomography at synchrotron radiation[J]. Acta Physica Sinica, 2012, 61(5): 054205-1-9. (in Chinese)

Analysis on Reconstructed Errors and Band-Artifacts in Projections for Synchrotron Radiation Computed Tomography

KONG Hui-hua1, YANG Yu-shuang2,3

(1. School of Science, North University of China, Taiyuan 030051, China;2. CSIRO Materials Science and Engineering, Victoria 3169, Australia;3. Theoretical Physics Institute, Shanxi University, Taiyuan 030006, China)

Abstract:The synchrotron radiation computed tomography has the advantages of synchrotron radiation X ray and nondestructive examination. In view of the problem of time fluctuation and spatial non-uniformity of synchrotron light source under the high photon energy, the artifacts caused by the instability of light source in the projection was analyzed, and the reconstruction error caused by the artifacts was studied. The result shows that temporal fluctuation of the beam causes the horizontal band artifacts in the projection data and reduces the resolution of the reconstruction image. The spatial non-uniformity of the light source causes the longitudinal banding artifacts in the projection data which not only reduces the resolution of the image but also causes the ring artifacts of reconstruction image. Numerical simulation and experimental results verify the correctness of the analysis.

Key words:synchrotron radiation CT; instability of light intensity; banding artifacts; reconstructed errors

中图分类号:TP391

文献标识码:A

doi:10.3969/j.issn.1673-3193.2016.01.013

作者简介：孔慧华(1977-)， 女， 副教授， 博士， 主要从事图像重建领域的研究.

基金项目：山西省自然科学基金资助项目(2012021011-2)； 山西省高等学校优秀创新团队支持计划和山西省百人计划资助

*收稿日期:2015-05-12

文章编号:1673-3193(2016)01-0061-06