双能CT系统中物质密度的识别方法*

王 辉，康南生，杨 波

(1.中国科学院合肥智能机械研究所，安徽合肥230031;2.中国科学技术大学 自动化系，安徽 合肥230027)

0 引言

计算机断层扫描成像技术(computerized tomography，CT)是计算机技术与放射学相结合产生的一门新的成像技术，它是在无损状态下获得被检测断面的灰度信息，以图像的形式清晰准确直观地展现被检测物体内部结构特征情况、材料的密度、有无缺陷等，广泛地应用于医学、工业等领域。

毒品、爆炸物对个人和社会都有极为严重的危害，如何在火车站，海港和飞机场等重要通道处快速准确地识别出毒品、爆炸物非常重要。利用双能量算法进行毒品、爆炸物检测获得广泛的应用，根据不同材料，不同密度的物质对不同能级的X射线衰减程度不同，对双能量系统获取的2种图像数据进行综合考虑计算，可以得到与被识别物体密度相关的值，可以提高系统的物质分辨能力，在毒品、爆炸物检测领域中起到重要作用。

国内外学者对双能量重建算法做了相关的研究，但基本都是建立在单色谱双能量的基础上的，所得到的透射图像与X射线穿过物质的厚度无关，但是在实际应用中，所用的系统基本都是连续谱X射线，这样就必须考虑待测样品的厚度变化对探测结果的影响。因此，常用的单色谱双能量重建算法并不能在连续谱双能量系统中得到应用［1］。

Alvarez RE和Macovsky A于1976年首次提出经典的双能量X射线断层扫描基本重建算法，它采用牛顿—拉夫逊迭代法去迭代2个非线性双能量方程组来重建原子序数和物质密度的分布，算法简单，但该方法耗时严重，而且，数字迭代过程中可能是发散的，因此，本文提出了一种新颖的基于双能量曲线的重建算法，通过对探测器得到的高、低能图像灰度值进行曲线拟合，实现对物质电子密度的仿真重建。实验结果表明:该算法效果很好，相对误差较小［2］。

1 投影重建与FBP算法

投影重建的理论基础是基于Lambert-Beers定律

式中 I0为入射X射线强度，I为出射X射线强度，Δx为厚度，μ为材料的线性衰减系数，通常μ随着X射线的能量改变，并因所选材料的不同而不同。

滤波反投影(FBP)算法的基础是傅立叶切片理论，物体在某个角度上得到的平行投影的傅立叶变换，等于该物体在同一角度下进行的二维傅里叶变换的一条直线，如图1所示。

图1 傅立叶切片理论示意图Fig 1 Diagram of Fourier slice theory

基于傅立叶切片理论［3］，只需要获取各个方面的投影数据的一维傅立叶变换，并将它们拼在一起得到整个图像的二维傅立叶变换，然后，做二维傅立叶反变换就可以得到重建图像。但由于傅立叶空间中产生的采样模式不是笛卡儿坐标的，所以，投影采样是在极坐标栅格上，为了执行二维傅立叶反变换，这些采样不得不被栅格化到一个笛卡儿坐标系中，然而，在频率域空间的差值会给重建图像造成漂移、伪影等干扰，所以，直接的差值傅立叶变换是不可取的。

物体的二维傅立叶变换是由许多一维傅立叶变换拼接起来的。由于通过投影采样得到的一维傅立叶变换是长条状，所以，拼接的过程必须加权，这个加权的过程就是加权滤波，而整个重建过程，一般称为滤波反投影，如图2所示［4］。

图2 滤波反投影示意图Fig 2 Diagram of filtered back-projection

2 双能量重建算法

2.1 双能量重建算法的基本原理

Alvarez R E和Macovsky A［5］在1976年发表的文章中指出给定物质对光子能的线性衰减系数μ在30～100kV范围内，可以表达为康普顿散射和光电效应衰减的线性组合，光电效应衰减成分很大程度地取决于被探测物质的原子序数，而康普顿散射给出的是物质的电子密度分布，它正比于物质的密度

式中 ρ为物质密度，g/cm3;A为原子质量，n=4.5。

设S(E)为入射能量谱，则连续谱X射线CT投影公式［6］

设SH(E)，SL(E)为入射的高、低能量谱，则

双能量算法实质就是通过探测到的双能量灰度值TH，TL值［7］，由以上两式得到A1，A2的值，然后，利用 FBP重建算法重建得到 a1，a2，最后，由 a1=ρZn/A，a2=ρZ/A，可以得到被探测物体的原子序数和电子密度ρZ/A，从而可以得到物体的相对密度。

2.2 基于双能量曲线的重建算法

在双能量CT系统中，当探测有效原子序数为Z的某样品时，当样品的厚度tm变化时，可得到一系列点集(TH，TL)，每一组(Z，tm)都对应一组灰度值(TH，TL)，形成一条TH-TL曲线，不同的Z值对应不同的曲线。因此，双能量曲线重建算法要首先标定出各样品的双能量曲线图TH-TL。

算法的具体实现步骤:

1)曲线标定:标定出若干违禁物品的双能量曲线图，在给定的高，低能量值下，对5种厚度的待探测样品进行分别透射，由得到的5种厚度下的高，低双能量值，绘出待探测样品的双能量曲线图TH-TL。

2)反向查找:当对某样品进行探测时，由探测到的双能量值(TH，TL)，根据已标定出的双能量曲线图TH-TL，反向去查找得到该样品的有效原子序数和厚度信息，即(Z，tm)。

3)计算投影数据集:由得到的Z，tm值计算投影数据集

4)FBP重建:由得到的投影数据集A1，A2，利用FBP重建算法重建出 a1，a2:a1=ρZn/A，a2=ρZ/A。

5)计算待测样品密度:a1，a2是与待测样品有效原子序数和电子密度相关的量，由此可以得到待测样品的密度。

3 实验与结果

取样品氯化钠(Zeff=14.88)，糖(Zeff=8.876)，TNT(Zeff=7.104)并分别装在自制的直径为6，8，10，12，14 mm的圆柱形泡沫器皿中，分别做双能量透射，实验设备为skyscan1076 CT机。

本文中双能量系统选取40，100 keV，10 W，扫描角度范围为0°～360°，旋转角度为1°，空间分辨率为18μm。由得到的投影数据，绘制4种物质的TH-TL曲线图，如图3所示。

图3 双能量曲线图Fig 3 Dual-energy curves

标定出物质的TH-TL曲线图后，就可以对该样品的任意厚度进行透射，通过曲线图可以反向查找得到样品的任意探测点的厚度，进而可以计算出样品的密度，得到的数据如表1。

由表1结果可知，本文提出的算法可以很好地得到被探测物质的密度，精度较高，相对误差小于3%。

4 结论

本文针对连续谱双能量X射线CT系统，提出了一种新颖的基于双能量曲线的重建算法，以得到被探测物质密度，并结合实验证明了该算法具有较高的精确度，相对误差小于3%，但现实系统中存在噪声，因此，降低这些噪声，可减小相对误差。

表1 实验数据Tab 1 Experimental data

［1］孙丽娜.X射线安检设备中物质分类识别方法的研究［D］.沈阳:东北大学，2003:2.

［2］Chuang K S，Huang H K.A fast dual-energy computational method using isotransmission lines and table look up［J］.Med Phys，1987，14(2):186－192.

［3］Jiang Hsieh.计算机断层成像技术原理、设计、伪像和进展［M］.张朝宗，译.北京:科学出版社，2006.

［4］庄天戈.CT原理与算法［M］上海:上海交通大学出版社，1992:7.

［5］Alvarez E，Macovski A.Energy-selective reconstructionsin X-ray computerized tomogrphy［J］.Phys Med Biol，1976，21(5):733－744.

［6］Tsunoo T，Torikoshi M，ohno Y.Measurement of electron density and effective atomic number using dual-energy X-ray CT［J］.Medical Physics，2008，35(1):4924－4932.

［7］Zhang Guowei，Chen Zhiqiang，Zhang Li.Exact reconstruction for dual energy computed tomography using an H—L curve method［J］.IEEE Nuclear Science Symposium Conference Record，2006，6:3485－3488.