基于投影匹配的X射线双能计算机层析成像投影分解算法_李保磊|层析成像

第31卷　第3期2011年3月

光　学　学　报V ol . 31, N o . 3March , 2011

基于投影匹配的X 射线双能计算机层析成像

投影分解算法

李保磊

1, 2

　张耀军

1

1

(公安部第一研究所, 北京100048; 2北京大学数学科学学院, 北京100871)

摘要　X 射线双能计算机层析成像(CT ) 技术是安全检查领域一种重要的材料探测与识别手段。双能C T 投影分解是双能CT 预处理重建算法的核心内容和关键步骤。针对现有投影分解算法的不足, 提出了一种基于投影匹配的双能CT 投影分解算法。依据系统能谱和基材料线性衰减系数曲线, 通过求解投影积分方程组建立高低能投影查找表。对于给定的高低能投影, 在查找表中寻找最佳匹配点, 进而获取基材料分解投影。该算法避免了现有算法复杂的迭代优化过程, 简化了系统的标定过程, 分解精度取决于查找表的设定步长。相对现有算法该算法有实现过程简单、易于并行计算的优点。仿真实验结果验证了算法的有效性。关键词　X 射线光学; 双能计算机层析成像; 投影匹配; 投影分解; 预处理重建中图分类号　T P391　　　文献标识码　A 　　　do i :10. 3788/AOS 201131. 0311002

Projection De composition Algorithm of X -Ray Dual -Ene rgy Compute d Tomography Base d on Projection Matching

Li Baolei 1, 2　Zhang Yaojun 1

1

The First R esearch Instit ut e of Minist ry of Public Security , Beijing 100048, China

2

School of Mathematical Sciences , Peking University , Beijing 100871, China

Abstract 　X -ray dual -energy computed tomography imaging technique is an important material detec tion and recognition method in the field of security inspection . The projection decomposition is the nuclear content and key technique in the pre -rec onstruction algorithm of dual -energy c om puted tom ography . Ac cording to the disadvantages of the c urrent algorithms , a projection decomposition algorithm based on projection matching is proposed . Firstly , based on energy spectrum of the system and the linear attenuation coefficient curve of the ba sic materials , the high -and low -energy projection lookup table can be got by solving the projection integral equations set . For a given dual -energy projection , find the best match point in the lookup table and then obtain the decomposition projection of basic m ateria ls . The proposed algorithm avoids the process of complex iteration and optimization and simplifies the proc ess of system calipation . The dec omposition ac curacy depends on the step setting of lookup table . Compared to the c urrent algorithms , the proposed algorithm ′s realization is more simple and easy for parallel computation . The feasibility of the a lgorithm is va lidated by the results of sim ulation experiment . Key wo rds 　X -ray optic s ; dual -energy computed tomography ; projection matching ; projection decomposition ; pre -reconstruction

OCIS co des 　110. 0110; 340. 0340; 020. 0020; 250. 0250

1　引　　言

X 射线双能计算机断层成像(C T ) 作为一种重要的成像技术, 在安全检查领域起着越来越重要的作用

[1～3]

包括单能X 射线透视、双能X 射线透视、多视角、背散射和X 射线CT 技术等) , 双能C T 技术具有较高

的探测能力, 它能够较准确的获取材料的电子密度、有效原子序数二维信息, 而利用此二维信息能够实

。在常见的X 射线安全检查技术中(主要

　　收稿日期:2010-06-10; 收到修改稿日期:2010-11-30

基金项目:国家科技支撑计划(2009BA K 64B02) 资助课题。

作者简介:李保磊(1980—), 男, 博士, 主要从事X 射线计算机断层成像技术及其应用等方面的研究。E -mail :bao lei li @yahoo . com . cn

光　　　学　　　学　　　报

现对绝大部分材料的探测识别[4]。同时CT 型安检设备也是美国FAA 唯一认定的EDS 型安检设备。可见双能C T 技术在整个X 射线安检技术中的地位。

双能CT 的概念在20世纪70年代由Robert E . Alvarez 等首次提出并进行了研究。最初该项技术主要用于医学射线图像的骨骼与软组织的分离和病灶诊断, 之后其应用逐渐扩展到工业检测和安全检查领域, 如利用双能技术对CT 图像进行射束硬化校正等。在安全检查领域, 国外Reveal 和Analogic 公司均已推出自己的双能CT 行包检查系统[6, 7]。

双能CT 安检系统的核心技术是双能CT 重建算法。双能CT 重建算法大致可以分为三类。1) 迭代重建算法。这类方法通过迭代的过程, 得到与能谱无关或者近似无关的CT 图像, 进而抑制由射束硬化效应带来的伪影, 改善图像质量; 或者通过双能迭代, 得到材料分解图像。不过此类方法计算速度慢, 不易实际应用。2) 后处理重建算法[10～12]。该方法首先利用高低能投影重建出高低能CT 图像, 然后再对高低能CT 图像进行处理, 获取物体断层的物理参数分布图像。由于其自身固有的特点, 此类方法不能从根本上消除射束硬化效应, 因此精度难以保证。3) 基于投影分解的双能CT 重建算法[13～16]。对于这类方法, 基于两种常见的物质线性衰减系数分解模型———双效应分解模型和基材料分解模型[13, 17], 对应着两种投影分解过程。这类方法通过对投影数据的预处理分解, 利用分解的参数完成双能C T 重建进而对材料进行识别探测, 因此又称为预处理重建算法。这种算法理论上可以得到不受能谱影响的物质有效原子序数和电子密度图像, 是目前双能CT 技术领域的主流重建算法。

在预处理重建算法中, 投影分解是关键技术和核心内容, 现有文献中已有多种投影分解方法, 如多项式间接拟合法、多项式直接拟合法、分段多项式拟合法、等值线法和高次曲面函数拟合法等。这些方法在一定程度上促进了双能CT 技术的发展,

[16, 18]

[8, 9][5]

但是它们求解过程复杂, 需要复杂的标定过程, 计算速度慢, 对噪声十分敏感, 分解精度还有待提高。2006年, Zhengro ng Ying 等提出了双能CT 的优化投影分解方法, 提高了分解算法的稳定度和收敛性能, 但是仍然存在计算速度与初值的选取问题。总之, 现有投影分解方法还不能满足安检实际工作需求, 需要研究简便实用的投影分解算法。

[19, 20]

2　双能CT 重建原理

物质的线性衰减系数有两种常见的分解模型:

μ(E )=a c f K N (E ) +a p f p (E ) , μ(E )=b 1μ1(E ) +b 2μ2(E ) ,

(1) (2)

式中f p (E ) , f K N (E ) 为只与能量E 有关而与材质无关的分解系数。且有

f p (E 3,

E

(3)

-ln (1+2α+f KN (α2

1+2αα

ln (1+2α) -2, 2α(1+2α)

(4)

式中α=E /510. 975keV , a p 表示光电效应系数, a c

为康普顿散射效应系数, a p , a c 是独立于能量只与材质有关的物理量, 且有

, 　a c =l 2, n ≈4～5, a p =l 1A A

n

(5)

式中l 1, l 2为两常数, ρ为物质密度, Z 为原子序数, A 为原子量; μ1(E ) , μ2(E ) 分别为两种基材料的线性衰减系数。b 1, b 2为对应两种基材料的分解系数, 对

于某一固定的物质, b 1, b 2是两个常数。(1) 式表示在一定的射线能量范围内, 物质的衰减可由光电效应和康普顿散射两种作用共同组成。(2) 式表示任何一种物质的线性衰减系数都可由两种基材料的线性衰减系数线性叠加而成[13, 17]。

依据这两种衰减系数分解模型, 记

A c =a c d l , A p =a p d l ; B 1=b 1d l , B 2=b 2d l , 根据宽能谱射线条件下的Beer 定律可得

∫∫∫∫

P L =-ln {S L (E ) ex p [-A c f K N (E ) -A p f p (E ) ]d E }+ln S L (E ) d E P H P L P H

H

c

KN

p

p

H

∫∫,

=-ln {S (exp [-A f (E ) -A f (E ) ]d E }+ln S (d E

∫E ) ∫E ) =-ln {S (exp [-B μ(E ) -B μ(E ) ]d E }+ln S (d E

∫E ) ∫E )

,

=-ln {S (ex p [-B μ(E ) -B μ(E ) ]d E }+ln S (d E

∫E ) ∫E )

L

1

1

2

2

L

H

1

1

2

2

H

(6)

(7)

李保磊等:　基于投影匹配的X 射线双能计算机层析成像投影分解算法

式中S L (E ) , S H (E ) , P L , P H 分别为高低能谱和高低能投影。基于投影分解的双能C T 预处理重建算法的核心即为(6) 式或者(7) 式的求解, 即根据(6) ,

(7) 式求解A c , A p ; B 1, B 2, 我们称这个过程为投影分解过程。由于A c , A p , B 1, B 2为a c , a p , b 1, b 2的线积分投影值, 求解出A c =a c d l , A p =a p d l ; B 1=b d l , B =b d l 后, 根据CT 原理, 利用滤波反投∫∫

1

2

2

2) 投影匹配

给定实际得到的高低能投影P L , P H , 寻找i =i 0, j =j 0, 使得[P L (i , j ) -P L ]2+[P H (i , j ) -P H ]2的值最小, 即

～～

[P L (i 0, j 0)-P L ]2+[P H (i 0, j 0)-P H ]2=[P L (i , j ) -P L ]

～

～

～

2

～

～

～

～

∫∫

+[P H (i , j )-P H ]

～

2

, (13)

也就是说投影P L , P H 在查找表中找到了最佳匹配点。

3) 计算P L , P H 对应的分解投影B 1(i ) , B 2(j )

B 1(i )=B 1min +(i 0-1) ·d ; B 2(j )=B 2min +(j 0-1) ·d .

～～

～

～

～

～

影重建算法, 便可计算出a c , a p , b 1, b 2, 由此可以计算材质的等效原子序数和电子密度信息, 以完成材料的探测识别。计算公式为

p

Z eff =K 1

a c

[5, 21]

1/n

(14)

, (8) (9)

1/n

e =ρK 2a c ,

式中K 1, K 2为两常数, n ≈3～4, 且

　　可以看出, 算法没有根据P L , P H 通过复杂的过程逆向求解, 而采用了正向投影匹配的方法实现了B 1, B 2的求解。

使得投影分解过程简单且易于实现。值得指出, B 1, B 2的步长越小, 算法的精度就越高。而且, 可以将利用文中方法得到的分解结果作为非线性最小二乘方法的初值, 由于初值已逼近理论值, 这样就可以得到稳定收敛的解。从这个角度分析, 本文方法也有其理论与现实意义。算法的流程图如图1所示。

Z eff

n n

1e112e22b 1ρe1+b 2ρe2

, (10)

ρe =b 1ρe1+b 2ρe2, (11)

式中Z 1, Z 2分别为两种基材料的原子序数; ρe1, ρe2分别为两种基材料的电子密度[21]。

3　基于投影匹配的双能CT 投影分解

算法

依据(6) , (7) 式的相似性, 不失一般性, 投影分解算法以基材料模型为例, 即对(7) 式进行求解。现有投影分解算法多对(7) 式做正向优化求解, 过程复杂。这里提出的分解方法和现有方法的思路相反。算法的思想是:首先依据C T 系统能谱信息生成合理物理意义范围内的各种B 1, B 2数值对应的P L , P H , 建立数据查找表。而后再根据实际得到的P L , P H 去匹配查找表中与之吻合的高低能投影, 同时获取对应的B 1, B 2值。该投影分解算法的步骤如下:

1) 查找表生成

设B 1, B 2范围为B 1∈[B 1min 　B 1max ], B 2∈[B 2min 　B 2max ]。设定一定的步长d , 生成B 1, B 2数值序列B 1(i ) , B 2(j ) , i =1, 2, …M , j =1, 2, …N , 依据系统能谱和(7) 式生成高低能投影序列:

P L (i , j ) i =1, 2, …, M ; j =1, 2, …, N , P H (i , j ) i =1, 2, …, M ; j =1, 2, …, N , (12) 建立高低能投影P L (i , j ) , P H (i , j ) 和B 1(i ) , B 2(j ) 数据查找表, 一对P L (i , j ) , P H (i , j ) 对应着一对B 1(i ) , B 2(j ) 。

图1双能计算机层析成像投影分解算法流程图Fig . 1Flow cha rt of pr ojectio n decomposition alg orithm

of X -ray dual -energ y computed to mog raphy

4　仿真实验结果

4. 1　系统能谱及查找表建立

仿真实验采用的高低能峰值管电压分别为80kV 和140kV , 并假设探测器响应理想线性, 系统能谱由开源能谱生成程序SpectrumGUI 生成如图2所示。

设定基材料为碳和铝, 且设定

B 1min =0, B 1max =10. 0cm

B 2min =0, B 2max =10. 0cm , d =0. 001cm (15) 依据(7) 式生成数据查找表。如图3所示。

[22]

,

光　　　学　　　学　　　

报

图2仿真系统能谱

Fig . 2Simulated ener gy spectrum of sy

stem

图3数据查找表F ig . 3Da ta loo kup

table

图4实验1的实验结果F ig . 4Re sults o f experiment 1

李保磊等:　基于投影匹配的X 射线双能计算机层析成像投影分解算法

4. 2　仿真实验1

依据上述双能实验条件, 对外层为1. 28m 厚度Fe 、内层为直径76. 8m m 水的圆柱体作双能扫

描。其中探测器256pixel , 像素间隔1m m ; 焦距590mm 。根据(7) 式生成高低能投影, 其中水和Fe 的线性衰减系数可由NIS T 数据库得到。选择碳C 和铁Fe 作为基材料, 设d =0. 01cm , 利用本文提出的算法, 得到如图4所示实验结果。实验过程中, 首先进行投影分解得到分解C 投影和Fe 投影, 然后再根据(7) 式计算高低能投影值, 将计算的高低能投影值与原始高低能投影值作对比, 以检验分解

[23]

算法的有效性。4. 3　仿真实验2

为了考察算法对材料探测的有效性, 依据仿真实验1实验条件对外层为0. 64cm 厚度Al 、内层为直径5. 12cm M g 柱的圆柱体作双能C T 断层扫描, 以C 和Al 作为基材料, 设d =0. 001cm 。利用本文提出的算法, 得到如图5所示实验结果。计算过程中(10) 式中n 取3. 5, 其中图5(i ) , 图5(j ) 中的曲线为图5(g ) , 图5(h ) 中所示的横线处有效原子序数和电子密度与真实值的对比

。

图5实验2的实验结果F ig . 5Re sults o f experiment 2

　　从上面的仿真实验可以看出, 利用所提出的投影分解算法进行双能投影分解, 通过分解投影计算的高低能投影与原高低能投影高度吻合, 这也证明了算法的有效性。同时, 利用提出的投影分解算法进行双能C T 重建, 获取的材料电子密度和有效原子序数和理想真实值十分吻合, 经计算可知误差在0. 5%以内, 满足材料探测需求。另外, 在寻找最佳匹配点时, 由于算法自身的特点和结构, 十分适合并行计算, 这也为算法速度的提高奠定了基础。相反迭代算法

[20]

况, 始终可以得到稳定的解。

5　结　　论

提出了一种基于投影匹配的双能CT 投影分解算法, 该算法相对现有算法, 实现过程简单, 易于并行运算, 实际仿真实验结果验证了算法的有效性。另外, 在实际双能C T 系统中, 无须每次扫描均用本文中的算法进行投影分解, 而可以事先用所建立高低能投影到分解投影的数据查找表, 此为系统标定过程。显然, 相比现有的算法, 该算法简化了系统的标定过程。而对于实际获取的高低能投影, 基于查找表利用线性插值的方法即可得到对应的分解投影, 这样可以保证双能CT 系统在实际进行物体材

则难以实现并行计算。同时在迭代算

法的计算过程中, 收敛结果对迭代初值的依赖性很大, 如果初值设置不当, 会导致算法不收敛或者出现错误的结果。而本文算法则不可能出现这样的情

光　　　学　　　学　　　报

料探测时的高效性。

参

考

文

献

1Li Baolei , Yang M in , Fu Jian et al . . T wo ring artifact cor rection m ethods for compu ted tom og raphy [J ]. Acta Op tica S inica , 2009, 29(7) :1849～1853

　李保磊, 杨　民, 傅　健等. 两种CT 成像环状伪影抑制方法[J ]. 光学学报, 2009, 29(7) :1849～1853

2Fang Zheng , S un Xiaomin , Luo Qingming . New sim ple m ethod of X -ray beam h ardening co rrection for CT recon struction [J ]. Acta Optica S in ica , 2007, 27(2) :302～306　方　正, 孙小敏, 骆清铭. 一种简便的计算机层折系统X 射线硬化校证方法[J ]. 光学学报, 2007, 27(2) :302～306

3Chen Xin , S un Yi , Zhu Peiping . Recon struction algorithm of cone -beam p hase X -ray com puter -tom ography based on g rating im aging [J ]. Acta Optica S inica , 2008, 28(6) :1079～1084　陈　欣, 孙　怡, 朱佩平. 光栅成像的锥束相位CT 重建算法[J ]. 光学学报, 2008, 28(6) :1079～10844R . F . Eilbert , K . D . Krug . Aspects of image recognition in vivid technology ′s dual -energy X -ray sys tem for explos ive detection [C ]. S P IE , 1992, 1824:127～1435R . Alvarez , A . M acovski . Energy -selective recon structions in X -ray com puterized tomography [J ]. Phys . Med . Bio l . , 1976, 21(5) :733～7446h ttp ://w w w . revealimaging . com /, 2009. 10. 87h ttp ://w w w . analogic . com /, 2009. 10. 8

8J . A . Fessler , I . Elb ak ri , P . S ukovic et al . . M aximum -likel ihood dual -energy tom og raphic image recons tru ction [C ]. S P IE , 2002, 4684:38～499P . S ukovic , N . H . C linthorne . Penalized w eigh ted least -s qu ares im age recons tru ction for dual energy X -ray transmission tomography [J ]. IE EE Tr ans . M ed . Imag . , 2000, 19(11) :1075～108110G . Ch rist . Exact treatment of the dual -energy method in C T u sing polyen ergetic X -ray spectra [J ]. Phy s . Med . Biol . , 1984, 29(12) :1501～151011J . C . M . Steenbeek , C . van Kuijk , J . L . Gras huis . S election of

fat -equivalent materials in pos tp roces s ing dual -energy quantitative

CT [J ]. Med . Ph ys . , 1992, 19(4) :1051～105612M . M . Goodsitt . Beam h ardening errors in pos t -processing dual en ergy quantitative computed tomography [J ]. Med . Ph ys . , 1995, 22(7) :1039～104713J . Peter S tones trom , R . Alvarez , A . M acovski . A framew ork for spectral artifact corrections in X -ray C T [J ]. IEE E Transactions on B iomed ical E ng ineer in g , 1981, 28(2) :128～14114K . S . Chuang , H . K . Hu ang . A fas t dual -energy compu tational method using iso -tran smission lines and table lookup [J ]. Med . Phys . , 1987, 14(2) :186～19215R . E . Alvarez , J . A . Seibert , S . K . Th om pson . C om parison of dual energy detector sys tem perfo rm ance [J ]. Med . Phy s . , 2004, 31(3) :556～565

16H . N . Cardinal , A . Fenster . An accu rate meth od for direct dual -energy calib ration and decomposition [J ]. Med . Ph ys . , 1990, 17(3) :327～341

17W . A . Kalen der , W . H . Perman , J . R . Vetter et a l . . Evaluation of a prototype dual -energy com puted tomographic apparatus . I . Phantom studies [J ]. 3) :Med . Phys . , 1986, 13(334～33918K . S . Chuang , H . K . Huang . Com paris on of fou r dual energy image decom position m ethods [J ]. Ph ys . Med . Biol . , 1988, 33(4) :455～466

19Z . R . Ying , R . Naidu , K . Guilbert et a l . . Dual energy volumetric X -ray tom og raphic s ens or for luggage screening [C ]. Proceedings of IEEE Sensors Applications Sym posium , 2007, (6-8) :1～620Z . R . Ying , R . Naidu , C . R . C raw fo rd . Dual energy computed tomograp hy for exp losive detection [J ]. J . X -r ay S cience an d

Tech no logy , 2006, 14(4) :235～256

21Gu ow ei Zhang , Jianping Ch eng , Li Zh ang et a l . . A practical reconstruction method for dual energy compu ted tom og raphy [J ]. J . X -Ray Science an d Technolog y , 2008, 16(2) :67～8822http ://sourcefo rge . net /p rojects /spectru mgui /, 2009. 12. 1823http ://physics . nis t . gov /Phys RefData /XrayM assCoef /Cover . html , 2009. 12. 18