双能DR物质识别算法在CT成像系统中的应用

作者：汤昕烨， 张丽， 陈志强， 高河伟， 张国伟

【摘要】  本文介绍一种适用于计算机断层成像（CT）的双能透视（DR）物质识别算法，此算法基于单能CT重建图像，通过对CT图像的分割及对断层几何信息的提取，分块重建材料有效原子序数及电子密度的分布。结合各种扫描轨迹的CT成像系统，可以实现有效的物质材料识别。同时，对比于传统的双能CT方法，本方法结合单能CT重建图像，改善了双能DR识别效果，能够实现较为精确的物质识别，在安全检查等应用领域有着现实意义。

【关键词】  双能DR； CT成像； 物质识别

    Abstract: A dualenergy digital radiography (DR) method assisted by segmentation of single energy reconstruction image is proposed for material recognition in Xray CT inspection systems. The effective atomic number and equivalent electron density distribution of the scanned objects can be reconstructed with this method. Compared to the conventional dual energy computed tomography (CT) technique, this method markedly reduces the amount of dualenergy detectors in the system of CT scanning, which shows its significance in application in the field of security inspection.

    Key words: dualenergy DR; CT image system; material recognition

    引言

    自“九一一”事件以来，在全世界范围内公众对恐怖事件的关注日益增加。出于公共安全因素的考虑，对乘客行李中炸药及其他危险品或违禁品的检查正越发得到重视。目前，机场和火车站入口使用的安全检查系统主要通过双能X射线方法来对乘客行李内的物质进行区分及检查［1］。

    Alvarez和Macovski ［2］于1976年最早提出了双能X射线断层扫描的基本重建方法，通过解决非线性方程式重建物质原子序数及密度的分布。近年来提出的一种双能曲线方法［3,4］，则通过对探测器得到的高、低能透明度进行曲线拟合，实现对未知均匀材料有效原子序数和质量厚度的精确重建。

    计算机断层成像（CT）技术解决了一般透视成像中的物体重叠问题，如果投影数据完备，就可以得到精确的断层图像。至今，CT成像系统已经发展到第五代。英国工程师Hounsfield于1967年设计了第一代CT机，所采用的是平行束平移—旋转方式。第二代CT机采用扇束平移—旋转方式。第三代CT机采用的是目前常用的扇形束探测器与源旋转扫描方式。第四代CT机的扫描方式和第三代的相似，但是扫描过程中只有X射线源旋转运动。第五代CT机是电子束扫描机（Electronbeam Scanner），它通过电子束的旋转来实现X光源的旋转运动，从而实现毫秒级的快速扫描［5］。

    结合上述技术的特点，开发基于双能技术的CT成像系统用于机场、海关、车站等行李物品的检查，不但能得到断层图像，辨认叠放物品，还能够获得被检物品的材料信息，实现较为准确的物质识别。

    本文提出了一种基于CT重建图像的双能DR重建方法，能够较为准确地重建材料的有效原子序数和电子密度分布。与常规双能DR算法相比，此方法可以解决材料重叠的问题。与基于求解非线性方程组的常规双能CT算法相比，该方法通过图像分割技术，采用分块重建和求解线性方程组的方法。与CT成像系统相结合，能够通过使用少量双能探测器实现较为准确的快速物质识别，适用于安全检查系统对高通关率的要求。

    1  双能断层成像重建

    传统的双能CT算法利用双能量投影数据信息重建待测物体的物质特性。通过对线性吸收系数进行分解，可将其表示为与物质特性相关量的线性组合，基本分解函数可分为不同的两类: 由康普顿散射和光电吸收所引起衰减的线性组合：μ(E)=c1μp(E)+c2μC(E),(1)其中μp(E)是光电效应引起的线性衰减，μC(E)是由于康普顿散射引起的线性衰减。常数c1和c2是所测材料的未知比例系数。

    两种所选基础材料的衰减系数的线性组合：μ(E)=a1μA(E)+a2μB(E),(2)其中μA(E)和μB(E)是两种所选基础材料的线性衰减系数。

    根据Beer法则，对一个实际X射线成像系统来说，当X射线通过均匀物质发生衰减时，探测器的高、低能透明度输出可分别表示为: TH=∫EHN(E)exp -∫μ(E)dlEPd(E)dE，(3a)

    TL=∫ELN(E)exp -∫μ(E)dlEPd(E)dE，(3b)其中N(E)是X射线能量谱，Pd(E)是探测器的能量响应函数，μ是线性衰减系数。

    因此，材料的物质特性可以通过已知的双能量投影数据以及对其中线性衰减系数的基本分解进行重建。分别基于重建图像和投影数据正弦图的postreconstruction方法及prereconstruction方法是两种经典的双能量CT重建算法。

    2  双能DR物质识别算法

    在传统的双能CT重建算法中，物质的有效原子序数和密度需要逐个像素进行重建。为了得到密度ρ和原子序数Z的分布信息，我们必须分别得到高能和低能的CT投影数据，而这需要在扫描过程中使用大量的双能探测器。但是在实际应用中，我们常常并不需要逐个像素的重建ρ和Z值。例如对于断层几何结构简单，物理特性均匀的物体，由于各均匀部分的ρ和Z值近似，因此只需要逐块重建即可达到物质识别的要求。因此我们提出由单能CT重建图像辅助的双能DR方法，通过分割单能CT重建图像提取物体断层的几何结构并逐块重建ρ和Z值。在直线轨迹扫描中，通过在某一方向使用双能量探测器，可同时完成单能CT及双能DR的扫描过程。

    对于固定的双能量探测器而言，待测物体与系统的相对直线运动可等效为对物体的DR扫描，如图1所示。假设对单能CT重建图像的分割将断层分割为相对均匀的N块进行重建，而由双能探测器获得了M组DR投射数据。则我们用TH(i)和 TL(i)表示第i组高、低能投影数据，lj(i)表示第i组投影数据所对应的射束经过第j块的长度。

    通过对线性衰减系数的分解：μ(E)=ρNAA(aσel+aσcohSC+aσincohSC)

    =ρNAA(aσel+aσSC)，（4）其中光电效应反应截面aσel=42Znα4m0c2E7/20，散射反应截面

    aσSC≈ZKN(E)，

    式中，α=1/137为精细结构常数，m0为电子质量，0=83π(e2/mc2)2，KN(E)表示康普顿散射的KleinNishina截面。

    由此可以定义物质的特征系数a1、a2：a1=ρZnA,(5a)

    a2=ρZA,(5b)其中指数n由光电效应与物质原子序数的关系取值在4～5之间，本文中取n=4.5。则衰减系数可以简化表示为μ(E)=ρZA(Zn-1F(E,Z)+G(E,Z))，其中F(E,Z)，G(E,Z)均与物质的原子序数Z弱相关。

    由此我们可以将高、低能透明度表示为：TH=∫EHN(E)exp-∫μ(E)dlEPd(E)dE

    =∫EHN(E)exp -A11E3-A2fKN(E)

    ·EPd(E)dE,(6a)

    TL=∫ELN(E)exp -∫μ(E)dlEPd(E)dE

    =∫ELN(E)exp -A11E3-A2fKN(E)

    ·EPd(E)dE，(6b)A1和A2是系数a1和a2沿射线方向的线积分，即

    A1=∫a1dl，  A2=∫a2dl.

    对于每一组高低能透明度(TH,TL)，通过双能曲线拟合方法，可以唯一确定等效均匀物质的原子序数Z和质量厚度tm。因此与之相应的一组系数(A1,A2)可以通过计算下式得到，A1=tmZn/A,(7a)

    A2=tmZ/A.(7b)由于各分块中系数a1和a2相同，其线积分可表示为各分块系数的加权和，A1(i)=∫a1dl=nj=1a1(i,j)=ka1klk(i),(8a)

    A2(i)=∫a2dl=nj=1a2(i,j)=ka2klk(i).(8b)由此我们可以得到系数a1={a1,1,a1,2,…,a1,N}及a2={a2,1,a2,2,…,a2,N}的线性方程组：l1(1)l2(1)…lN(1)

    l1(2)………

    

    l1(M)……lN(M)a1,1

    a1,2

    

    a1,N

    =A1(1)

    A1(2)

    

    A1(M),          (9a)

    l1(1)l2(1)…lN(1)

    l1(2)………

    

    l1(M)……lN(M)a2,1

    a2,2

    

    a2,N

    =A2(1)

    A2(2)

    

    A2(M).          (9b)根据各块系数(a1,k,a2,k)与相应物质的原子序数及密度的关系，a1,k=ρkZnkAk,

    a2,k=ρkZkAk,(10)可得，a1,ka2,k=Zn-1k，  k∈1,2,…,N.(11)因此通过求解方程组(9a)和(9b)，我们可以对物质的有效原子序数及电子密度分布进行重建。

在确保断层分块数量不超过所获取的双能量数据组数，即M>N的情况下，我们采用最优化方法减小误差的影响，令Si=Nj=1lj(i),(12)并选择(N-1)个不相关方程使得Si尽可能大。如果存在两个方程(编号p和q，p,q∈1,2,…,M)使得Sp=Sq，且A1(p) < p>

    a*1,2=f2(a1,1),

    

    a*1,N=fN(a1,1).(14)余下(M-N+1)个方程：

    l1(N)l2(N)…lN(N)

    l1(N+1)………

    

    l1(M)……lN(M)a1,1

    a1,2

    

    a1,N

    =A1(N)

    A1(N+1)

    

    A1(M),(15)

    也可表示为，

    l1(N)a1,1+l2(N)a1,2+…

    +lN(N)a1,N-A1(N)=0,

    l1(N+1)a1,1+l2(N+1)a1,2+…

    +lN(N+1)a1,N-A1(N+1)=0,

    

    l1(M)a1,1+l2(M)a1,2+…

    +lN(M)a1,N-A1(M)=0.(16)

    由于误差的存在，我们定义函数fe来刻画由a*1,1,a*1,2,…,a*1,N引起的差别。因此，

    fe(N)=l1(N)a*1,1+l2(N)a*1,2+…

    +lN(N)a*1,N-A1(N),

    fe(N+1)=l1(N+1)a*1,1+l2(N+1)a*1,2+…

    +lN(N+1)a*1,N-A1(N+1),        (17)

    

    fe(M)=l1(M)a*1,1+l2(M)a*1,2+…

    +lN(M)a*1,N-A1(M).

    结合方程(15)有，

    fe(N)=l1(N)a1,1+l2(N)f2(a1,1)+…

    +lN(N)fN(a1,1)-A1(N),

    fe(N+1)=l1(N+1)a1,1+l2(N+1)f2(a1,1)

    +…+lN(N+1)fN(a1,1)-A1(N+1),

    

    fe(M)=l1(M)a1,1+l2(M)f2(a1,1)+…

    +lN(M)fN(a1,1)-A1(M),       (18)

    即

    fe(N)=gN(a1,1),

    fe(N+1)=gN+1(a1,1),

    

    fe(M)=gM(a1,1),(19)

    我们得到fe(a1,1)=fe(N)+fe(N+1)       

    +…+fe(M)

    =gN(a1,1)+gN+1(a1,1)

    +…+gM(a1,1).(20)目标是得到使得fe(a1,1)最小的a1,1，即a1,1=arg min (fe(a1,1)).(21)结合方程组(15)我们得到a*1,1,a*1,2,…,a*1,N以及方程组(9)的最优解a*1={a*1,1,a*1,2,…,a*1,N}。

    综上所述，双能DR重建物质识别方法可以总结如下：

     对单能CT重建图像按吸收系数进行分割，提取几何量以构建线性方程组(9a)和(9b)的系数矩阵；

     通过双能曲线方法或者多项式拟合的方法由一组双能DR投影数据(TH,TL)得到相应的等效均匀物质的原子序数和质量厚度(Z,tm)；

     通过公式(7a)和(7b)计算相应的系数A1和A2，建立线性方程组(9a)和(9b)；

     用最优化方法求解线性方程组(9a)和(9b)，计算各分块系数a1，a2。通过(10)计算有效原子序数Zeff和电子密度ρZ/A。

    3  数值模拟计算验证

    模拟计算采用140 keV及160 keV的X射线能谱，对5种均匀材料组成的水浴模型进行前向投影，得到双能DR投影数据。表1为各种材料的有效原子序数及密度等参数。图2  水浴模型（a）及所使用的X射线能谱（b）示意图表1  模型材料及参数模拟计算中，单能CT重建图像由对带有高斯噪音(均值0，方差0.01)的模型进行平行束FBP重建图像模拟。重建图像大小为256×256，双能投影数据组数为256。重建结果如表2所示。表2  模拟计算结果及误差

    模拟结果显示，各材料的有效原子序数相对误差均小于5%，电子密度的相对误差在1%左右。从重建结果可以看出，双能DR物质识别算法对于结构简单的均匀材料的识别具有较高的准确性。然而，由于重建过程依赖于对单能CT重建图像的分割，因此单能CT图像的重建质量以及模型断层的几何复杂性是影响该算法准确性的重要因素。

    4  结语

    本文介绍了一种采用双能X射线DR的物质识别方法，通过单能CT图像的辅助，重建被扫描物体的有效原子序数及电子密度的分布。结合具有拟平行束性质的直线轨迹CT成像系统，一方面能够得到断层图像，解决物体重叠问题；同时能够实现较为精确的物质识别。因此能够满足安检系统所需的快速通关、能够剥离重叠物体并识别材料等要求。由于单能CT重建图像对重建物质ρ和Z值的影响较大，今后研究的方向将着重于对直线轨迹CT重建图像分割效果的改善。

 

【参考文献】
  ［1］王琪，陈志强，邬小平，等.X 射线安全检查技术综述［J］. CT 理论与应用研究，2004，13(1)：32-37

　　［2］Robert E. Alvarez and Albert Macovski, Energyselective Reconstructions in Xray Computerized Tomography ［J］. Phys. Med. Biol., 1976, 21(5): 733-744.

　　［3］Guowei Zhang, Li Zhang, Zhiqiang Chen, An HL Curve Method for Material Discrimination of Dual Energy Xray Inspection Systems ［A］. Conf. Record of 2005 IEEE Nuclear Science Symposium ［C］. 2005， 326-328.

　　［4］Guowei Zhang, Zhiqiang Chen, Li Zhang, and Jianping Cheng, Exact Reconstruction for Dual Energy Computed Tomography Using an HL Curve Method ［A］.Conf. Record of 2006 IEEE Medical Imaging Conference［C］.2006， M14462.

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

　　［6］Sidky E Y, Zou Y, and Pan X. Volume image reconstruction from a straightline source trajectory, IEEE Symposium on NSSMIC ［C］, 2005.

　　［7］Deborah J. Walter, Eric J. Tkaczyk, Xiaoye Wu. Accuracy and precision of dual energy CT imaging for the quantification of tissue fat content ［J］. Proceedings of SPIE ［C］. 2006, 6142.

　　［8］Ohnoa Y, Torikoshia M, Tsunooa T, et al. Dualenergy Xray CT with CdTe array and its extension ［J］. Nuclear Instruments and Methods in Physics Research, 2005, A (548): 72-77.