项 安 陈瑞焘

（同济大学 电气工程系 上海 201804）

基于MCNP提取X射线透射图像双能值的数值方法

项 安 陈瑞焘

（同济大学 电气工程系 上海 201804）

传统的双能X射线检测技术使用与物质原子序数相关的特征值进行物质分类，但由于该特征值与被检物体的厚度有关（厚度效应），影响了双能检测方法的可靠性。为此提出了一种利用辛普森公式以及序列二次规划法的数值算法，以消除双能值的厚度效应。用基于MCNP的蒙特卡罗方法模拟双能X射线检测过程，结合概率密度分类法证明该算法能显著降低物质分类的错误概率，且得到了不同物质的分类边界。用X射线安检机进行测试，结果表明该算法显著消除了物质的厚度效应，且分类边界有效。说明本文算法的有效性与实用性。

X射线，特征值，厚度效应，数值算法，MCNP

双能X射线透射检测技术广泛应用于行李安检领域[1]。该技术通过高低能X射线透射信号TH和TL计算与物质原子序数相关的双能值，以此来区分物质的类别，如式(1)所示：

式中，T0H和T0L分别是高低能X射线通过无障碍物的自由空间时的透射信号；R是与原子序数相关的特征值，称为双能值；μ(Z,E)为能量为E的X射线在原子序数为Z的物质中透射的衰减系数；σ(Z,E1)与σ(Z,E2)分别是有效原子序数为Z的物质与能量为E1和E2的反应总截面，与物质的原子序数及光子能量E有确定的关系[2-3]。

式(1)是根据单色谱X射线的透射规律得到[2]，但实际使用的X光源是通过轫致辐射产生的，特点是具有连续的能谱分布[4]。实际X射线穿过物体的透射信号如式(2)所示：

式中，N(E)为能量为E的光子数量；Pd(E)为X光探测器对能量为E的光子的检测效率。由于总的透射信号是不同能量透射信号的积分，将式(2)代入式(1)后，得到的R值不仅与物质的原子序数相关，还与X射线透射厚度t有关。如式(3)所示：

式中，t为X射线经过物体的厚度。

X射线安检设备实用过程中，行李中的物体摆放位置随意，X射线穿透厚度t也不同。图1为一种行李中物品可能的摆放方式。当t值不同时，相同材料可能测得不同的R值，而两种不同的材料i、j当满足μiti=μjtj时可能会测得近似的R值，导致双能值的厚度效应，使得根据R值进行分类的方法可靠性降低。为了解决该问题，提出了一种消除双能检测技术厚度效应的数值方法，建立了双能X射线检测的蒙特卡罗模型并进行仿真，仿真结果的处理中使用概率分类方法，得到了三种实验样品的分类边界。最后通过实际X射线异物检测机验证了其有效性。

图1 行李中物品可能摆放方式Fig.1 Possible position of objects in the luggage.

1 解决厚度效应的数值算法

1.1 利用辛普森公式的数值积分

如果能根据式(2)得到μ(Z,E)，则可直接利用式(1)求得准确的双能值R，并且与t无关。下面提出一种用于求解μ(Z,E)的数值算法。

根据式(2)，多色X射线的透射信号是不同能量透射信号的积分。首先利用数值积分的辛普森公式将式(2)进行转换，积分区间0-Ein可以分为M个等间距的区间，M+1个分割点定义为：

令式(2)中

且令

则式(2)可写为：

对式(6)中的定积分使用复合辛普森公式[5]得：M为偶数时

式中，h=Ein/M，为等分的能量区间长度。根据辛普森公式的误差估计式，式(7)的误差为：

将式(5)代入式(7)，可得：

其中：

式中，ψm是能量为Em的X射线在透射厚度为t时的衰减值。

1.2 利用序列二次规划方法求最优解

根据式(11)，只要分别求得高能、低能光源能量为Ej和Ei的X射线透射衰减值ψh(Z,Ej)和ψl(Z,Ei)，即可根据式(12)求得精确的R值，且消去了厚度t的影响。

显然无法利用式(9)直接求出ψm(m=0,1,2,…,M)。但通过序列二次规划(Sequential Quadratic Programming, SQP)方法可以得到最优解，SQP是一种求解非线性约束优化问题的有效方法。主要求解过程分为三步：(1) 拉格朗日函数Hessian矩阵的更新；(2) 二次规划问题求解；(3) 一维搜索和目标函数的计算。使用SQP方法首先要定义优化函数与约束条件[6]，对于式(9)求解问题可以定义式(13)所示的优化函数。

又根据式(1)和式(11)，结合σ(Z,E)与原子序数Z和光子能量E的关系[7-8]可知：

式中，m,n≥0。E取值不同时，X射线光子与物质原子的主要作用机制有所不同，m与n的取值相应有所不同[7-8]。

当E＜EM1时，主要为光电效应，这时n=3；当EM1＜E＜EM2时，主要作用机制为相干散射，这时n=2；当EM2＜E＜EM时，主要作用机制为非相干散射（康普顿散射），这时n=0[7-8]。由此可以定义式(15)所示的约束条件：

结合式(12)和式(14)，即可利用SQP方法求解能够使式(13)中优化函数取得最小值的ψm(m=0,1,2,…,M)，并且求得的ψm满足式(15)中的约束条件。Matlab中优化工具箱提供的fmincon函数可以用于解决该问题。

综上所述，该数值方法的求解思路如下：

(a) 计算式(10)所示的am及相关参数。其中(m=0,1,2,…,M)，构建式(13)所示的优化函数和式(15)所示的约束条件；

(b) 测量被检物体的高低能X射线透射信号TH和TL；

(c) 令T=TL，结合式(13)用SQP方法计算ψl向量，其中l=0,1,2,…,M；

(d) 令T=TH，结合式(13)用SQP方法计算ψh向量，其中h=0,1,2,…,M；

(e) 取ψl中的ψi与ψh中的ψj计算双能值R=ψi/ψj，ψi和ψj是向量ψl和ψh中的第i个和第j个元素，i和j取值根据实验确定。选取原则是使得不同厚度的同种材料的R值尽量接近，根据求得的R值即可对物体进行分类。事先要测定不同样品的R值，确定物质的分类边界。本文基于MCNP构建了双能X射线透射检测模型，根据该模型测得不同材料样品的高低能透射数据，其优势是避免了用实际X光机测量时样本种类不够的问题。

2 基于MCNP构建蒙特卡罗仿真模型

蒙特卡罗仿真方法，也称随机模拟方法，其基本思想是：为了求解数学、物理、工程技术以及生产管理等方面的问题，首先建立一个概率模型或随机过程，使他的参数等于问题的解；然后通过对模型或过程的观察或抽样实验来计算所求参数的统计特征，最后给出所求解的近似值，而解的精确度也可以用估计值得标准误差来表示。

MCNP是美国阿拉姆斯(Los Alamos)国家实验室研制开发的一个大型多功能的蒙特卡罗程序。能解决中子、光子、电子或者偶合中子、光子、电子的输运（适用的能量范围为：中子10-11-20 MeV，光子1 keV-100 GeV，电子1 keV-1 GeV），以及计算核临界系统（包括次临界和超临界系统）的特征值。本文使用的MCNP程序版本为MCNP5，用以模拟双能透射中X射线光子的运输过程。

在MCNP中建立的双能X射线检测仿真模型如图2所示。X射线源使用厚度为1 cm扁平面源，光子能量分布根据实际光源能谱曲线设定[8-9]，距离铜片15 cm，金属薄片距离待测样品15 cm，待测样品距离闪烁晶体15 cm。铜片厚度为1 mm，且只用于高能X射线的检测，低能X射线检测时不使用，其目的是使得减少高低能谱的重叠。闪烁晶体为CsI（碘化铯），密度为4.51 g·cm-3。闪烁晶体后端通过导光材料接到光电倍增管，闪烁体其他表面与包装材料直接接触，要求包装材料对可见光有较高的反射率。入射到晶体的X射线光子与晶体原子相互作用，产生次级带电粒子，这些带电粒子又与晶体原子相互作用，发射光子。带电粒子沉积在闪烁晶体中的能量与发射的荧光强度成一定的比例，故在仿真模型中，只要用*F6命令记录闪烁晶体栅元内的电子沉积能量。再根据式(16)即可算得检测得到的光强[10]。

式中，Qpe是光电倍增管阳极输出的电子数；Np是晶体沉积能量后产生的闪烁光的光子数；η是闪烁光子在晶体中传输直到被光电倍增管阴极收集的传输效率；QE是光电倍增管阴极将光子转化为电子的量子效率；M为光电倍增管放大倍数，对R7224倍增管，M=3.3×106。根据实验测定，闪烁晶体中可见光约为540 nm，对应QE=8.7。对于本文使用的150 keV及以下的X射线，Np≈10%。η值参照文献[11]测定和计算，可以得到不同沉积能量时的不同η取值。

闪烁晶体除光电倍增管读出端以外，均使用白色Teflon材料（FR104-1，上海三爱富新材料股份有限公司生产）包装。晶体与光电倍增管之间有1mm的硅脂作为导光材料，晶体为10 mm× 10mm×50 mm的长方体。

图2 双能X射线检测的MCNP仿真模型Fig.2 MCNP simulation model of dual-energy X-ray detection.

仿真中待测样品选取了三种典型的材料：铝、聚乙烯、有机玻璃。分别测得三种材料不同厚度下的高低能透射信号，代入本文提出的数值算法后，即可得到消除了厚度效应的双能值。算法中计算am需要N(Em)与Pd(Em)，N(Em)是光源发出的不同能量X射线光子数，根据实际X光源特性设定。Pd(Em)是探测器对能量为Em的X射线的检测效率，可以通过实验测定[11-12]。根据经验，式(15)中的M值取31，M1取5，M2取15。算法中的i取18，j取24，此时式(12)中的Ei=E18=42 keV，Ej=E24=107 keV，双能值为：

3 仿真结果分析

3.1 双能值的比较

通过仿真获得铝、聚乙烯、有机玻璃不同厚度的高能透射信号THt-Al、THt-PE、THt-PMMA和低能透射信号TLt-Al、TLt-PE、TLt-PMMA，并且测得没有待测样品的自由空间透射信号T0H和T0L，即可根据式(1)用传统方法计算双能值RTt-Al、RTt-PE、RTt-PMMA，该方法没有消除厚度效应对双能值的影响。另一方面，将THt-Al、THt-PE、THt-PMMA和TLt-Al、TLt-PE、TLt-PMMA分别代入本文提出的数值算法中，结合式(12)，可以算得三种材料不同厚度下的双能值RSt-Al、RSt-PE、RSt-PMMA，该方法消除了厚度效应的影响，理论上双能值只与材料原子序数有关，而与材料厚度无关。

计算得到的结果如表1所示，且算得了每种材料双能值的平均值与标准差。

表1 铝、聚乙烯和有机玻璃高低能透射信号与双能值Table 1 High and low energy transmission signal and dual energy value of aluminum, polyethylene and polymethyl methacrylate.

（续表1）

3.2 概率密度分类方法

得到物体不同厚度的双能值后，希望据此找到不同物质之间的分类边界。因为双能值与材料厚度有关，可以定义双能值随厚度变化的概率密度函数。又因实际安检过程中，被检测物体的材料、厚度等参数是随机变化的，所以该概率密度函数可以定义为正态分布函数，如式(18)所示。

式中，Ci表示第i类物质材料（i为Al、PE、PMMA）；μi和σi是Ci材料双能值的平均值和标准差；P(r/Ci)为Ci物质双能值的概率密度分布。据此定义式(19)所示的决策函数gi(r)为：

式中，P(Ci)是Ci材料的先验概率。决策规则如满足：

则决策为Ci类材料。

如果各种材料的先验概率P(r/Ci)是相等的，则决策函数变为：

将表1中得到的三类材料的平均值μAl、μPE、μPMMA与方差σAl、σPE、σPMMA代入式(21)，可得到三种材料双能值的概率密度曲线。图3(a)所示为用基于式(1)的传统方法得到的RTt-Al、RTt-PE、RTt-PMMA的概率密度曲线。图3(b)为本文提出的优化数值算法得到的RSt-Al、RSt-PE、RSt-PMMA的概率密度曲线。其中C1为聚乙烯，C2为有机玻璃，C3为铝。

图3 用传统方法(a)和优化算法(b)得到的三种材料双能值概率密度曲线Fig.3 Probability density curve of three materials’ dual energy value obtained by traditional methods (a) and optimization algorithms (b).

图3 同时展示了两种方法分别得到的三种材料之间的分类边界，如图3(a)中C1与C2的分类边界r=1.4121，C2与C3的分类边界r=1.6140。所以若算得的双能值r＜1.4121，则决策为C1；若r＞1.4121并且r＜1.6140，则决策为C2；若r＞1.6140，则决策为C3。结合图3，使用传统方法计算双能值进行物质分类的错误概率为式(22)所示：

使用本文提出的数值算法计算双能值进行物质分类的错误概率为：

由式(22)、(23)可见，使用本文提出的数值算法后，双能X射线检测技术的物质分类错误概率明显降低，说明了该方法的有效性。

4 实验测试

通过以上的仿真得到了三类材料的分类边界，为了验证该边界有效性，下面用上海太易公司的双能X射线异物检测机进行测试。检测对象为t0=2 cm厚的铝板与有机玻璃板，图4为待测物体与X射线源以及X射线检测器的相对位置示意图。X射线发射扇形射线束，经过厚度为t0的待测物体时穿透厚度不同，如图4中t0、t1、t2和t3所示，导致双能透射信号TH和TL随厚度变化。图5为实验得到的两种材料的双能透射图像，图5中行号与图4中行号对应。由图5，不同行号对应不同X射线穿透厚度，透射灰度值随行号逐渐变化，且透射厚度越短时，灰度值越高，图像越亮。

图4 待测物体与检测机的相对位置Fig.4 Relative position between the object to be detected and inspection machine.

图5 铝(a)和有机玻璃(b)双能透射图像Fig.5 Dual energy transmission images of aluminum (a) and polymethyl methacrylate (b).

根据图5得到的透射信号是图像灰度值，由于本文提出的优化算法中需要提供的透射信号为归一化的X射线强度，所以对于算法中的式(12)修正如下：

式中，A为常数；I为图像灰度值。

根据图5得到的灰度值，用基于式(1)的传统方法计算铝和有机玻璃的双能值R，并绘出R值随行号（即穿透厚度）的变化，如图6(a)所示。由图6(a)，两种材料的双能值随行号的增加有明显变化，且均有超过分类边界r=1.6140的部分，所以利用传统方法计算双能值进行物质分类就可能出现错误。用本文提出的优化数值算法计算两种材料双能值，绘出R值随行号的变化，如图6(b)所示。由图6(b)，优化的数值算法使得双能值基本不随X射线穿透厚度而变化，且两种材料可以用分类边界r=1.1021有效分类。说明通过MCNP仿真得到的分类边界有效。

图6 用传统方法(a)和优化算法(b)计算的双能值随行号的变化Fig.6 Variation of dual-energy values according to the line number calculated by traditional methods (a) and optimization algorithms according (b).

5 结语

仿真以及实际测试的结果都说明本文提出的优化数值算法能够有效消除双能值随X射线透射厚度不同而产生的变化，即消除了厚度效应，显著降低了物质分类的错误概率。利用MCNP构建仿真模型，模拟双能X射线检测，可以获得各种材料分类边界的精确值。通过实际的双能X射线异物检测机进行测试，验证了使用本文提出的优化数值算法得到的双能值和分类边界能够更有效区分物质，再次证明了该算法的实用性与有效性。

1 王琪, 陈志强, 邬小平, 等. X射线安全检查技术综述[J]. CT理论与应用研究, 2004, 13(1): 32-37

WANG Qi, CHEN Zhiqiang, WU Xiaoping, et al. Review of X-ray security inspect ion technology[J]. CT Theroy and Applications, 2004, 13(1): 32-37

2 Xie W. Simulation of X-ray imaging systems for luggage inspection[D]. Blacksburg Virginia, USA: Virginia Polytechnic Institute and State University, 1995

3 Abbott L A, Conners, Wei X, et al. Simulation of X-ray imaging for luggage inspection[J]. Proceedings of the Second Explosives Detection Symposium & Aviation Security Conference, 1996, (11): 248-253

4 Fewell T R, Shuping R E. Photon energy distribution of some typical diagnostic X-ray beams[J]. Medical Physics, 1977, 4(3): 187-197

5 肖筱南. 现代数值计算方法[M]. 北京: 北京大学出版社, 2003: 203-206

XIAO Xiaonan. The modern numerical calculation method[M]. Beijing: Peking University Press, 2003: 203-206

6 张德丰. Matlab数值分析与应用[M]. 北京: 国防工业出版社, 2009: 282-287

ZHANG Defeng. Matlab numerical analysis and application[M]. Beijing: National Defense Industry Press, 2009: 282-287

7 Lu Q. The utility of X-ray dual-energy transmission and scatter technologies for illicit material detection[D]. Virginia: Virginia Polytechnic Institute and State University, 1999

8 阮书州. 高能X射线物质识别的蒙特卡罗模拟[D]. 吉林: 吉林大学, 2011

RUAN Shuzhou. Monte Carlo simulations of high energy X-ray in material identification[D]. Jilin: Jilin University, 2011

9 Ay M R, Shahriari M, Sarkar S, et al. Monte Carlo simulation of X-ray spectra in diagnostic radiology and mammography using MCNP4C[J]. Physics in Medicine and Biology, 2004, (49): 4897-4917

10 岳珂, 徐瑚珊, 梁晋洁, 等. 基于Geant4的CsI(Tl)闪烁体探测器的Monte Carlo模拟[J]. 原子核物理评论, 2010, 27(4): 85-89

YUE Ke, XU Hushan, LIANG Jinjie, et al. Monte Carlo simulation of CsI scintillator detector based on Geant4[J]. Nuclear Physics Review, 2010, 27(4): 85-89

11 Giulia Hull, Choong Woon-Seng, Moses William W, et al. Measurements of NaI(TI) electron response: comparison of different samples[J]. IEEE Transactions on Nuclear Science, 2009, 56(1): 331-336

12 Aitken D W, Beron B L, Yenicay G, et al. The fluorescent response of NaI(Tl), CsI(Tl), CsI(Na) and CaF2(Eu) to X-rays and low energy gamma rays[J]. IEEE Transactions on Nuclear Science, 1967, 14(1): 468-477

CLC TL99

A numerical method for extracting dual-energy value of X-ray transmission image based on MCNP

XIANG An CHEN Ruitao
(Department of Electrical Engineering, Tongji University, Shanghai 201804, China)

Background: Traditional dual-energy X-ray inspection technology serves to classify different materials based on the proper value related with effective atomic number. Purpose: This paper aims to solve the problem of the relevance between the proper value and material’s thickness (thickness effect), which reduces the reliability of dual-energy inspection technology. Methods: A numerical method which combines Simpson formula and sequential quadratic programming (SQP) method is proposed to eliminate the thickness effect of dual-energy values. The process of X-ray inspection is simulated by Monte Carlo method based on MCNP. Probability-classification method is used to evaluate the effectiveness of diminishing error probability of material-classification, and classification boundaries of different materials. A practical X-ray inspection machine is used to verify this method. Results: The experimental results indicate that not only the thickness effect is diminished obviously by this method, but also the classification boundaries are effectively achieved. Conclusion: The dual-energy value of X-ray transmission image can be extracted by our proposed method with effectiveness and practical applicability.

X-ray, Proper value, Thickness effect, Numerical method, MCNP

TL99

10.11889/j.0253-3219.2014.hjs.37.070203

项安，男，1965年出生，2001年于南昌大学获博士学位，研究方向为X射线安检相关技术

陈瑞焘，E-mail: 082663@tongji.edu.cn

2014-03-27，

2014-04-16