锥束工业CT射束硬化校正方法

张益海，张 催，潘小东，商宏杰，骆岩红，马思汉，李公平

(兰州大学 核科学与技术学院，兰州 730000)

锥束工业CT射束硬化校正方法

张益海，张 催，潘小东，商宏杰，骆岩红，马思汉，李公平

(兰州大学 核科学与技术学院，兰州 730000)

锥束X射线CT系统中，X射线束具有宽能谱的特性，会导致X射线硬化效应的产生，从而使得重建图像中出现杯状和条状伪影等，故有必要研究硬化校正方法以减少硬化效应带来的影响。提出了基于试验能谱和蒙特卡罗模拟的多项式拟合硬化校正方法，基于试验测量X射线初始能谱，利用蒙特卡罗模拟得到透射厚度与多能投影的关系，再利用多项式拟合得到多能投影与等效单能投影的转换关系，从而实现硬化校正。分别采用该方法和传统的滤波片法对一铝制标准件进行射束硬化校正试验。结果表明，所提出的方法比滤波片法能更好地消除锥束CT的硬化伪影。

锥束CT；硬化校正；多项式拟合法；滤波片法

工业CT(Industrial Computed Tomography，ICT) 能够获取被检测物体内部断层的图像信息，具有其他无损检测方法无法替代的特点，在工业无损检测领域具有广阔的应用前景[1-6]。工业CT通常使用X射线管作为X射线源，其出射的X射线具有宽能谱的特性。由于低能光子比高能光子在透射物体过程中的衰减更多，透射物体后的X射线能谱发生变化，X射线谱峰向高能方向偏移，高能光子所占比例增大，射线变得更难衰减，探测器所探测到的高能光子所占比例增大，从而导致投影值偏小。在CT图像重建常用的算法中，如滤波反投影算法(FBP)、代数迭代算法(ART)等，都是基于X射线能量单一的假设[7]。而实际上锥束工业CT的X射线是多色的，因此只能得到多能投影Pm(Multi-Energy Projection)数据，若直接由Pm数据进行重建图像，射束硬化效应会导致本来密度均匀的被测物体在CT图像上表现为亮度不一，灰度值不相等，通常直观显示为杯状或条状伪影等[8]。硬化伪影会降低CT图像的分辨率，影响检测的精度，因此有必要对其进行校正。多项式拟合硬化校正方法是CT系统中常用的校正方法，该方法是利用多项式拟合得到Pm和透射厚度L的关系曲线，再取曲线原点处的切线作为等效单能投影Ps(Single-Energy Projection)与L的关系曲线，最后通过多项式拟合得到Pm与Ps的转换关系式。该方法的主要缺点是需要根据要求专门制备特定的吸收片，但吸收片的厚度很难做到足够小，这就使得在原点处的切线变得不够精确，并且对于由多种材料混合组成的样品，很难找到与其材料相同的吸收片。因此，笔者采用蒙特卡罗方法模拟射线穿过吸收片的试验，得到相应结果，相较于传统的试验测量，其优点主要在于可以方便地定义吸收片的厚度和材料，并且减少了大量的测量时间和花费。但由于在利用蒙特卡罗方法模拟前没有得到X射线管的初始能谱而使得结果误差较大，因此，笔者基于前期对X射线管初始能谱测量的工作，提出了基于试验能谱和蒙特卡罗模拟的多项式拟合硬化校正方法；并且将该方法应用于一铝制标准件开展试验，并和传统的滤波片法[9]做对比，试验结果表明，相较于滤波片方法，提出的方法能给出更好的校正效果，且对射线利用率更高。

1 射线硬化校正原理分析

对于单能X射线，并且样品材料均匀的情况，X射线的衰减量可由Lambert-Beer[10]定律来描述,见式(1)。

式中:I0为X射线的初始强度;I为X射线穿过样品后的透射强度;L为射线透射过样品的厚度，简称透射厚度;μ为样品的线性吸收系数。

由式(1)变换得到ln(I0/I)=μL，通常令Ps=ln(I0/I)。则Beer定律可写成：

图1给出了单能投影Ps和多能投影Pm的投影值与透射厚度L的关系示意。由图1可看出：对于单能X射线，投影值与L呈线性关系；对于多能X射线，投影值与L呈非线性关系，投影值随着L的增加，增大的趋势逐渐趋于平缓。这是由于出射X射线光强不但与L有关，还与射线的能量相关，如式(3)所示。

式中:E为射线的能量；S(E)为射线的能谱分布；μ(E)为样品对应能谱的线性衰减系数。

则Beer定律可写成多能投影Beer定律形式，如式(4)所示。

式中：Pm为过原点的曲线(见图1)。

由于X射线管发出的X射线具有连续能谱分布，当X射线透射样品时，射线束中低能光子比高能光子吸收得更快而导致射线硬化效应，使得X射线吸收不均匀，产生高频信号；如将此线衰减系数μ(E)[11]代入滤波反投影重建算法经过计算后，会产生条状或杯状伪影，这就降低了CT图像的质量，影响检测的精度。因此，必须进行相应的校正：一方面可以通过找到Pm与Ps的转换关系以实现硬化校正；另一方面可以降低X射线的多色性，使其近似为“单色光”以减小硬化效应。

图1 多能投影、单能投影的投影值与透射厚度的关系

2 试验过程

2.1 基于试验能谱和蒙特卡罗模拟的多项式拟合硬化校正方法

2.1.1 X射线管初始能谱测量

在用蒙特卡罗方法进行模拟试验前，需要得到X射线管出射X射线能谱，而这可通过试验测量或者计算机模拟得到。为了得到更准确的结果，笔者采用试验测量的方法，采用瑞士COMET公司生产的MXR225/22 X射线管进行。主要试验参数为：管电压100 kV，管电流0.2 mA，小焦点d为1.0 mm，靶材为钨，靶角20°，辐射范围为40°，固有滤波为0.8 mm铍窗，出射口有2 mm铝滤波片。利用美国ORTEC公司GLP-10180/07P4型高纯锗探测器测量能谱。试验结果如图2所示，X射线由连续谱和特征谱构成。特征峰主要包括钨和银元素产生的特征峰，这是由于靶材为钨，而银具有极佳的导热性，可作为钨靶与铜衬底的钎焊材料[12]。

图2 试验测量X射线管小焦点初始能谱(100 keV)

2.1.2 蒙特卡罗方法的软件模拟计算

蒙特卡罗方法是一种随机过程方法，其并不严格求解一些物理方程，而是通过模拟单粒子在介质中的随机运动过程，然后将大量粒子的平均结果作为系统结果的近似，给出粒子系统在介质中运输的最后结果。笔者采用MCNP-4C软件工具包建立X射线与不同厚度吸收片作用的模拟试验模型(见图3)。X射线源为点源，所出射光子能量服从如图2所示的能谱分布。吸收片材料为Al，其厚度从0.1 mm增加到20 mm，每次增加的厚度为0.1 mm，共200组模拟试验。探测器置于样品后方，三者的几何位置根据试验测量所得，其中源-片距离为51 mm，源-探距离为815 mm。由于空气对射线的吸收可忽略不计，为简单起见，模拟环境设置为真空。

图3 模拟试验模型

所编写的MCNP软件输入文件代码，包括曲面块设计，栅元块设计和数据块设计，模拟了射线穿过吸收片后射线强度发生变化的过程。在给定射线初始强度I0的情况下，通过探测器探测射线穿过不同吸收片后的透射强度I，得到200组试验数据，利用投影值计算公式P= ln(I0/I)计算得到Pm，以L为横坐标，Pm为纵坐标作图，即可得到两者的关系曲线，如图4所示。

图4 多能投影值与透射厚度的关系

2.1.3 校正公式计算

吸收片厚度很薄时，硬化效应很弱，投影值与L近似成线性关系，因此在图4曲线原点处作其切线，该切线即代表Ps与L的关系，结果如图5所示。由图5可看到：Ps与Pm在L较小时，相差很小，但随着L的增加，两者的差距逐渐增大。这表明随着样品厚度的增大，硬化效应带来的影响也将逐渐增大。

图5 多能投影值、单能投影值与透射厚度的关系

在相同透射厚度下，对应得到200组(Ps,Pm)值，并以多能投影值与单能投影值的关系Pm为横坐标，Ps为纵坐标，得到两者的关系曲线，如图6所示。从图6可看出：Ps随Pm的增大而增大，且两者逐渐趋于线性关系。

图6 多能投影值与单能投影值的关系

再对该试验数据点进行多项式拟合。由于4阶多项式拟合就能得到很好的拟合效果，由此得到Ps与Pm的转换关系，如式(5)所示。

成像过程中，只需将试验采集得到的Pm数据先代入式(5)转换为Ps数据，再进行图像重建，即可实现射束硬化校正的目的。

2.2 滤波片硬化校正方法

滤波片硬化校正方法，就是在CT系统的X射线源出束口处加装滤波片，以对X射线进行预硬化，滤除低能X射线，使尽量多的高能射线通过。图7为初始透射能谱与滤波后能谱的对比，其中两个谱的总计数均为400 000±2 000。笔者选择0.46 mm铜片作为滤波片，可看出射线的低能部分被很好地过滤掉了，X射线的平均能量得到了增大，从而降低了硬化效应。该方法最大的不足在于，其硬化校正质量与滤波片厚度成正相关，但随着滤波片的增厚，射线强度损失严重，因此必须增加辐照时间或者增大X射线管的管电流以保证探测器拥有足够的曝光。

图7 初始能谱与过滤能谱的对比图

3 校正结果及对比

3.1 多项式拟合硬化校正结果 试验使用的是COMET公司的X射线管，型号为MXR225/22。试验主要参数为：管电压100 kV，管电流0.2 mA，采用小焦点(d=1.0 mm)。使用美国VARIAN公司的非晶硅平板探测器，型号为PaxScan1313DX，该探测器的探测面积为13 cm×13 cm，像元矩阵总数为1 024×1 024，像元尺寸为127 μm，极限分辨率为3.94 lp·mm-1，动态范围为16位。对一铝制标准件进行锥束X射线CT试验，得到360个角度下的多能投影数据，选择第512层数据进行断层图像重建，重建图像像素为400像素×400像素，采用FBP重建算法。

图8 多项式拟合硬化校正前后重建切片效果对比

得到基于试验能谱和蒙特卡罗模拟的多项式拟合硬化校正前后切片重建效果的对比，如图8所示。从图8(a)、(b)可看到经过校正以后，投影正弦图变得更加清晰。从图8(c)、(d)可看到：校正前，样品周围比较模糊，同时样品中间相对较暗；校正后，样品更为清晰，也更为均匀；从图8(e)、(f)中能更加直观地看出，未校正的重建图像存在明显的杯状伪影且噪声起伏较大，而经过校正后的重建图像中的杯状伪影得到了有效消除，切片的灰度值几乎均匀，同时噪声起伏减小，重建图像质量得到明显提高。

图9 滤波片硬化校正前后重建切片效果对比

3.2 滤波片法硬化校正结果

采用滤波片硬化校正方法的结果如图9所示。试验条件同3.1节，此外加装0.46 mm铜片，管电流0.5 mA。通过校正前后结果的对比可看出，该校正方法能在一定程度上减少硬化伪影，使图像质量得到改善，但图像仍存在硬化伪影。

图10 三维CT图像硬化校正前后对比

3.3 两种校正方法对比

由3.1节和3.2节的结果可看到，传统的滤波片硬化校正方法只能在一定程度上抑制硬化伪影现象，不能完全解决射束硬化带来的问题，并且由于滤波片对X射线的吸收强烈，探测器曝光减少，因此不得不增大管电流以保证足够的曝光。笔者提出的基于试验能谱和蒙特卡罗模拟的多项式拟合硬化校正方法，能在不降低射线利用率的情况下，得到更好的硬化校正结果。

此外，对于由较高原子序数组成的样品，例如铁和铜样品，滤波片硬化校正方法开始不太适用，而基于试验能谱和蒙特卡罗模拟的多项式拟合硬化校正方法仍有较好的效果，图10给出了将该方法应用于锥束工业CT系统三维CT图像的结果，样品为钢卷尺。由图10可看到，经过校正后的三维CT图像在细节处更为清晰，伪影更少。试验条件同3.1节，对于铁样品的硬化校正，Ps与Pm的转换关系如式(6)所示。

4 结论

提出了基于试验能谱和蒙特卡罗模拟的多项式拟合硬化校正方法，并将其应用于锥束CT系统的硬化校正工作中。结果表明，与传统的滤波片硬化校正方法相比，所提出的方法不但有更好的硬化校正效果，还能提高射线的利用率，可成为锥束工业CT射束硬化校正的一种有效方法。由于研究有限，此方法暂时仅应用于单材料样品的硬化校正，还未应用于多材料混合样品的硬化校正工作，后续将做进一步的工作。

[1] KASTNER J, HARRER B, REQUENA G, et al. A comparative study of high resolution cone beam X-ray tomography and synchrotron tomography applied to Fe- and Al-alloys[J]. Ndt& E International, 2010, 43(7):599-605.

[2] AGRAWAL A K, SARKAR P S, KASHYAP Y S, et al. Application of X-Ray CT for non-destructive characterization of graphite fuel-tube[J]. Journal of Nondestructive Evaluation, 2016, 35(3):1-8.

[3] LIU T. Cone-beam CT reconstruction for planarobject[J]. Ndt& E International, 2012, 45(1):9-15.

[4] IZUMI S,KAMATA S, SATOH K, et al. High energy X-ray computed tomography for industrial applications[J]. IEEE Transactions on Nuclear Science, 1993, 40(2):158-161.

[5] 张丰收,张定华,邓效忠,等. 工业CT射束硬化研究[J]. 机械强度,2007,29(6):1022-1025.

[6] 吴志宏, 丛鹏, 刘锡明. 基于重投影的CT图像硬化伪影校正[J]. 原子能科学技术, 2015, 49(5):935-938.

[7] JIANG H. Computed tomography principles, design, artifacts, and recent advances[M]. Bellingham:SPIE Press, 2009.

[8] 曾钢, 刘力, 阙介民,等. 工业CT图像均匀性校正[J]. 高能物理与核物理, 2004, 28(10):1131-1134.

[9] 吴世法. 近代成像技术与图像处理[M].北京:国防工业出版社,1997:143-144.

[10] HANKE R, BOBEL F. Determination of material flaw size by intensity evaluation of polychromatic X-ray transmission[J]. NDT & E International,1992,25(2):87-93.

[11] KROHN B R, SILVER M D. Optimizing contrast resolution in industrial computed tomographic scanners[J]. Materials Evaluation,1990,49(10):1296-1331.

[12] 张催,李公平,潘小东,等. 工业CT用X射线管射线强度分布研究[J]. 核技术, 2015, 38(9): 8-12.

The Methods for Beam Hardening Correction of Cone-beam Industrial Computed Tomography

ZHANG Yihai, ZHANG Cui, PAN Xiaodong, SHANG Hongjie, LUO Yanhong, MA Sihan, LI Gongping

(School of Nuclear Science and Technology, Lanzhou University, Lanzhou 730000，China)

Since the X-ray beam has the characteristics of wide spectrum in the cone beam X-ray computed tomography imaging system, this leads to the X-ray hardening effect, so as to make the reconstruction image in the goblet artifact and strip artifact.It is necessary to study methods of the hardening correction in order to reduce the effects of hardening effect.A method of beam hardening correction is put forward, which is based on the experimental spectrum and Monte Carlo simulation of polynomial fitting.Based on the initial X-ray energy spectrum by experimental measurement, the relationship between transmission thickness and pluripotent projection can be obtained by exploiting Monte Carlo simulation. And then using polynomial fitting can obtain relationship between pluripotent projection and equivalent single projection, so as to realize hardening correction. This method and the traditional filter method for X-ray hardening correction on the aluminum workpiece were respectively used. The experimental results show that the proposed method is better than the filter method which can decrease cone-beam CT beam hardening artifacts effectively.It provides an effective correction method for cone beam industrial computed tomography beam hardening correction.

cone-beam CT; beam hardening correction; polynomial fitting method; filtering method

2016-10-12

兰州大学中央高校基本科研业务费专项资金资助项目(lzujbky-2016-32, lzujbky-2016-208)

张益海(1988-)，男，硕士，主要从事X射线成像研究

李公平， ligp@lzu.edu.cn

10.11973/wsjc201706002

TL99;TG115.28

A

1000-6656(2017)06-0008-05