基于MCNP提取X射線透射圖像雙能值的數(shù)值方法

項 安 陳瑞燾

（同濟大學(xué) 電氣工程系 上海 201804）

基于MCNP提取X射線透射圖像雙能值的數(shù)值方法

項 安 陳瑞燾

（同濟大學(xué) 電氣工程系 上海 201804）

傳統(tǒng)的雙能X射線檢測技術(shù)使用與物質(zhì)原子序數(shù)相關(guān)的特征值進行物質(zhì)分類，但由于該特征值與被檢物體的厚度有關(guān)（厚度效應(yīng)），影響了雙能檢測方法的可靠性。為此提出了一種利用辛普森公式以及序列二次規(guī)劃法的數(shù)值算法，以消除雙能值的厚度效應(yīng)。用基于MCNP的蒙特卡羅方法模擬雙能X射線檢測過程，結(jié)合概率密度分類法證明該算法能顯著降低物質(zhì)分類的錯誤概率，且得到了不同物質(zhì)的分類邊界。用X射線安檢機進行測試，結(jié)果表明該算法顯著消除了物質(zhì)的厚度效應(yīng)，且分類邊界有效。說明本文算法的有效性與實用性。

X射線，特征值，厚度效應(yīng)，數(shù)值算法，MCNP

雙能X射線透射檢測技術(shù)廣泛應(yīng)用于行李安檢領(lǐng)域[1]。該技術(shù)通過高低能X射線透射信號TH和TL計算與物質(zhì)原子序數(shù)相關(guān)的雙能值，以此來區(qū)分物質(zhì)的類別，如式(1)所示：

式中，T0H和T0L分別是高低能X射線通過無障礙物的自由空間時的透射信號；R是與原子序數(shù)相關(guān)的特征值，稱為雙能值；μ(Z,E)為能量為E的X射線在原子序數(shù)為Z的物質(zhì)中透射的衰減系數(shù)；σ(Z,E1)與σ(Z,E2)分別是有效原子序數(shù)為Z的物質(zhì)與能量為E1和E2的反應(yīng)總截面，與物質(zhì)的原子序數(shù)及光子能量E有確定的關(guān)系[2-3]。

式(1)是根據(jù)單色譜X射線的透射規(guī)律得到[2]，但實際使用的X光源是通過軔致輻射產(chǎn)生的，特點是具有連續(xù)的能譜分布[4]。實際X射線穿過物體的透射信號如式(2)所示：

式中，N(E)為能量為E的光子數(shù)量；Pd(E)為X光探測器對能量為E的光子的檢測效率。由于總的透射信號是不同能量透射信號的積分，將式(2)代入式(1)后，得到的R值不僅與物質(zhì)的原子序數(shù)相關(guān)，還與X射線透射厚度t有關(guān)。如式(3)所示：

式中，t為X射線經(jīng)過物體的厚度。

X射線安檢設(shè)備實用過程中，行李中的物體擺放位置隨意，X射線穿透厚度t也不同。圖1為一種行李中物品可能的擺放方式。當(dāng)t值不同時，相同材料可能測得不同的R值，而兩種不同的材料i、j當(dāng)滿足μiti=μjtj時可能會測得近似的R值，導(dǎo)致雙能值的厚度效應(yīng)，使得根據(jù)R值進行分類的方法可靠性降低。為了解決該問題，提出了一種消除雙能檢測技術(shù)厚度效應(yīng)的數(shù)值方法，建立了雙能X射線檢測的蒙特卡羅模型并進行仿真，仿真結(jié)果的處理中使用概率分類方法，得到了三種實驗樣品的分類邊界。最后通過實際X射線異物檢測機驗證了其有效性。

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

1 解決厚度效應(yīng)的數(shù)值算法

1.1 利用辛普森公式的數(shù)值積分

如果能根據(jù)式(2)得到μ(Z,E)，則可直接利用式(1)求得準(zhǔn)確的雙能值R，并且與t無關(guān)。下面提出一種用于求解μ(Z,E)的數(shù)值算法。

根據(jù)式(2)，多色X射線的透射信號是不同能量透射信號的積分。首先利用數(shù)值積分的辛普森公式將式(2)進行轉(zhuǎn)換，積分區(qū)間0-Ein可以分為M個等間距的區(qū)間，M+1個分割點定義為：

令式(2)中

且令

則式(2)可寫為：

對式(6)中的定積分使用復(fù)合辛普森公式[5]得：M為偶數(shù)時

式中，h=Ein/M，為等分的能量區(qū)間長度。根據(jù)辛普森公式的誤差估計式，式(7)的誤差為：

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

其中：

式中，ψm是能量為Em的X射線在透射厚度為t時的衰減值。

1.2 利用序列二次規(guī)劃方法求最優(yōu)解

根據(jù)式(11)，只要分別求得高能、低能光源能量為Ej和Ei的X射線透射衰減值ψh(Z,Ej)和ψl(Z,Ei)，即可根據(jù)式(12)求得精確的R值，且消去了厚度t的影響。

顯然無法利用式(9)直接求出ψm(m=0,1,2,…,M)。但通過序列二次規(guī)劃(Sequential Quadratic Programming, SQP)方法可以得到最優(yōu)解，SQP是一種求解非線性約束優(yōu)化問題的有效方法。主要求解過程分為三步：(1) 拉格朗日函數(shù)Hessian矩陣的更新；(2) 二次規(guī)劃問題求解；(3) 一維搜索和目標(biāo)函數(shù)的計算。使用SQP方法首先要定義優(yōu)化函數(shù)與約束條件[6]，對于式(9)求解問題可以定義式(13)所示的優(yōu)化函數(shù)。

又根據(jù)式(1)和式(11)，結(jié)合σ(Z,E)與原子序數(shù)Z和光子能量E的關(guān)系[7-8]可知：

式中，m,n≥0。E取值不同時，X射線光子與物質(zhì)原子的主要作用機制有所不同，m與n的取值相應(yīng)有所不同[7-8]。

當(dāng)E＜EM1時，主要為光電效應(yīng)，這時n=3；當(dāng)EM1＜E＜EM2時，主要作用機制為相干散射，這時n=2；當(dāng)EM2＜E＜EM時，主要作用機制為非相干散射（康普頓散射），這時n=0[7-8]。由此可以定義式(15)所示的約束條件：

結(jié)合式(12)和式(14)，即可利用SQP方法求解能夠使式(13)中優(yōu)化函數(shù)取得最小值的ψm(m=0,1,2,…,M)，并且求得的ψm滿足式(15)中的約束條件。Matlab中優(yōu)化工具箱提供的fmincon函數(shù)可以用于解決該問題。

綜上所述，該數(shù)值方法的求解思路如下：

(a) 計算式(10)所示的am及相關(guān)參數(shù)。其中(m=0,1,2,…,M)，構(gòu)建式(13)所示的優(yōu)化函數(shù)和式(15)所示的約束條件；

(b) 測量被檢物體的高低能X射線透射信號TH和TL；

(c) 令T=TL，結(jié)合式(13)用SQP方法計算ψl向量，其中l(wèi)=0,1,2,…,M；

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

(e) 取ψl中的ψi與ψh中的ψj計算雙能值R=ψi/ψj，ψi和ψj是向量ψl和ψh中的第i個和第j個元素，i和j取值根據(jù)實驗確定。選取原則是使得不同厚度的同種材料的R值盡量接近，根據(jù)求得的R值即可對物體進行分類。事先要測定不同樣品的R值，確定物質(zhì)的分類邊界。本文基于MCNP構(gòu)建了雙能X射線透射檢測模型，根據(jù)該模型測得不同材料樣品的高低能透射數(shù)據(jù)，其優(yōu)勢是避免了用實際X光機測量時樣本種類不夠的問題。

2 基于MCNP構(gòu)建蒙特卡羅仿真模型

蒙特卡羅仿真方法，也稱隨機模擬方法，其基本思想是：為了求解數(shù)學(xué)、物理、工程技術(shù)以及生產(chǎn)管理等方面的問題，首先建立一個概率模型或隨機過程，使他的參數(shù)等于問題的解；然后通過對模型或過程的觀察或抽樣實驗來計算所求參數(shù)的統(tǒng)計特征，最后給出所求解的近似值，而解的精確度也可以用估計值得標(biāo)準(zhǔn)誤差來表示。

MCNP是美國阿拉姆斯(Los Alamos)國家實驗室研制開發(fā)的一個大型多功能的蒙特卡羅程序。能解決中子、光子、電子或者偶合中子、光子、電子的輸運（適用的能量范圍為：中子10-11-20 MeV，光子1 keV-100 GeV，電子1 keV-1 GeV），以及計算核臨界系統(tǒng)（包括次臨界和超臨界系統(tǒng)）的特征值。本文使用的MCNP程序版本為MCNP5，用以模擬雙能透射中X射線光子的運輸過程。

在MCNP中建立的雙能X射線檢測仿真模型如圖2所示。X射線源使用厚度為1 cm扁平面源，光子能量分布根據(jù)實際光源能譜曲線設(shè)定[8-9]，距離銅片15 cm，金屬薄片距離待測樣品15 cm，待測樣品距離閃爍晶體15 cm。銅片厚度為1 mm，且只用于高能X射線的檢測，低能X射線檢測時不使用，其目的是使得減少高低能譜的重疊。閃爍晶體為CsI（碘化銫），密度為4.51 g·cm-3。閃爍晶體后端通過導(dǎo)光材料接到光電倍增管，閃爍體其他表面與包裝材料直接接觸，要求包裝材料對可見光有較高的反射率。入射到晶體的X射線光子與晶體原子相互作用，產(chǎn)生次級帶電粒子，這些帶電粒子又與晶體原子相互作用，發(fā)射光子。帶電粒子沉積在閃爍晶體中的能量與發(fā)射的熒光強度成一定的比例，故在仿真模型中，只要用*F6命令記錄閃爍晶體柵元內(nèi)的電子沉積能量。再根據(jù)式(16)即可算得檢測得到的光強[10]。

式中，Qpe是光電倍增管陽極輸出的電子數(shù)；Np是晶體沉積能量后產(chǎn)生的閃爍光的光子數(shù)；η是閃爍光子在晶體中傳輸直到被光電倍增管陰極收集的傳輸效率；QE是光電倍增管陰極將光子轉(zhuǎn)化為電子的量子效率；M為光電倍增管放大倍數(shù)，對R7224倍增管，M=3.3×106。根據(jù)實驗測定，閃爍晶體中可見光約為540 nm，對應(yīng)QE=8.7。對于本文使用的150 keV及以下的X射線，Np≈10%。η值參照文獻[11]測定和計算，可以得到不同沉積能量時的不同η取值。

閃爍晶體除光電倍增管讀出端以外，均使用白色Teflon材料（FR104-1，上海三愛富新材料股份有限公司生產(chǎn)）包裝。晶體與光電倍增管之間有1mm的硅脂作為導(dǎo)光材料，晶體為10 mm× 10mm×50 mm的長方體。

圖2 雙能X射線檢測的MCNP仿真模型Fig.2 MCNP simulation model of dual-energy X-ray detection.

仿真中待測樣品選取了三種典型的材料：鋁、聚乙烯、有機玻璃。分別測得三種材料不同厚度下的高低能透射信號，代入本文提出的數(shù)值算法后，即可得到消除了厚度效應(yīng)的雙能值。算法中計算am需要N(Em)與Pd(Em)，N(Em)是光源發(fā)出的不同能量X射線光子數(shù)，根據(jù)實際X光源特性設(shè)定。Pd(Em)是探測器對能量為Em的X射線的檢測效率，可以通過實驗測定[11-12]。根據(jù)經(jīng)驗，式(15)中的M值取31，M1取5，M2取15。算法中的i取18，j取24，此時式(12)中的Ei=E18=42 keV，Ej=E24=107 keV，雙能值為：

3 仿真結(jié)果分析

3.1 雙能值的比較

通過仿真獲得鋁、聚乙烯、有機玻璃不同厚度的高能透射信號THt-Al、THt-PE、THt-PMMA和低能透射信號TLt-Al、TLt-PE、TLt-PMMA，并且測得沒有待測樣品的自由空間透射信號T0H和T0L，即可根據(jù)式(1)用傳統(tǒng)方法計算雙能值RTt-Al、RTt-PE、RTt-PMMA，該方法沒有消除厚度效應(yīng)對雙能值的影響。另一方面，將THt-Al、THt-PE、THt-PMMA和TLt-Al、TLt-PE、TLt-PMMA分別代入本文提出的數(shù)值算法中，結(jié)合式(12)，可以算得三種材料不同厚度下的雙能值RSt-Al、RSt-PE、RSt-PMMA，該方法消除了厚度效應(yīng)的影響，理論上雙能值只與材料原子序數(shù)有關(guān)，而與材料厚度無關(guān)。

計算得到的結(jié)果如表1所示，且算得了每種材料雙能值的平均值與標(biāo)準(zhǔn)差。

表1 鋁、聚乙烯和有機玻璃高低能透射信號與雙能值Table 1 High and low energy transmission signal and dual energy value of aluminum, polyethylene and polymethyl methacrylate.

（續(xù)表1）

3.2 概率密度分類方法

得到物體不同厚度的雙能值后，希望據(jù)此找到不同物質(zhì)之間的分類邊界。因為雙能值與材料厚度有關(guān)，可以定義雙能值隨厚度變化的概率密度函數(shù)。又因?qū)嶋H安檢過程中，被檢測物體的材料、厚度等參數(shù)是隨機變化的，所以該概率密度函數(shù)可以定義為正態(tài)分布函數(shù)，如式(18)所示。

式中，Ci表示第i類物質(zhì)材料（i為Al、PE、PMMA）；μi和σi是Ci材料雙能值的平均值和標(biāo)準(zhǔn)差；P(r/Ci)為Ci物質(zhì)雙能值的概率密度分布。據(jù)此定義式(19)所示的決策函數(shù)gi(r)為：

式中，P(Ci)是Ci材料的先驗概率。決策規(guī)則如滿足：

則決策為Ci類材料。

如果各種材料的先驗概率P(r/Ci)是相等的，則決策函數(shù)變?yōu)椋?/p>

將表1中得到的三類材料的平均值μAl、μPE、μPMMA與方差σAl、σPE、σPMMA代入式(21)，可得到三種材料雙能值的概率密度曲線。圖3(a)所示為用基于式(1)的傳統(tǒng)方法得到的RTt-Al、RTt-PE、RTt-PMMA的概率密度曲線。圖3(b)為本文提出的優(yōu)化數(shù)值算法得到的RSt-Al、RSt-PE、RSt-PMMA的概率密度曲線。其中C1為聚乙烯，C2為有機玻璃，C3為鋁。

圖3 用傳統(tǒng)方法(a)和優(yōu)化算法(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。結(jié)合圖3，使用傳統(tǒng)方法計算雙能值進行物質(zhì)分類的錯誤概率為式(22)所示：

使用本文提出的數(shù)值算法計算雙能值進行物質(zhì)分類的錯誤概率為：

由式(22)、(23)可見，使用本文提出的數(shù)值算法后，雙能X射線檢測技術(shù)的物質(zhì)分類錯誤概率明顯降低，說明了該方法的有效性。

4 實驗測試

通過以上的仿真得到了三類材料的分類邊界，為了驗證該邊界有效性，下面用上海太易公司的雙能X射線異物檢測機進行測試。檢測對象為t0=2 cm厚的鋁板與有機玻璃板，圖4為待測物體與X射線源以及X射線檢測器的相對位置示意圖。X射線發(fā)射扇形射線束，經(jīng)過厚度為t0的待測物體時穿透厚度不同，如圖4中t0、t1、t2和t3所示，導(dǎo)致雙能透射信號TH和TL隨厚度變化。圖5為實驗得到的兩種材料的雙能透射圖像，圖5中行號與圖4中行號對應(yīng)。由圖5，不同行號對應(yīng)不同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).

根據(jù)圖5得到的透射信號是圖像灰度值，由于本文提出的優(yōu)化算法中需要提供的透射信號為歸一化的X射線強度，所以對于算法中的式(12)修正如下：

式中，A為常數(shù)；I為圖像灰度值。

根據(jù)圖5得到的灰度值，用基于式(1)的傳統(tǒng)方法計算鋁和有機玻璃的雙能值R，并繪出R值隨行號（即穿透厚度）的變化，如圖6(a)所示。由圖6(a)，兩種材料的雙能值隨行號的增加有明顯變化，且均有超過分類邊界r=1.6140的部分，所以利用傳統(tǒng)方法計算雙能值進行物質(zhì)分類就可能出現(xiàn)錯誤。用本文提出的優(yōu)化數(shù)值算法計算兩種材料雙能值，繪出R值隨行號的變化，如圖6(b)所示。由圖6(b)，優(yōu)化的數(shù)值算法使得雙能值基本不隨X射線穿透厚度而變化，且兩種材料可以用分類邊界r=1.1021有效分類。說明通過MCNP仿真得到的分類邊界有效。

圖6 用傳統(tǒng)方法(a)和優(yōu)化算法(b)計算的雙能值隨行號的變化Fig.6 Variation of dual-energy values according to the line number calculated by traditional methods (a) and optimization algorithms according (b).

5 結(jié)語

仿真以及實際測試的結(jié)果都說明本文提出的優(yōu)化數(shù)值算法能夠有效消除雙能值隨X射線透射厚度不同而產(chǎn)生的變化，即消除了厚度效應(yīng)，顯著降低了物質(zhì)分類的錯誤概率。利用MCNP構(gòu)建仿真模型，模擬雙能X射線檢測，可以獲得各種材料分類邊界的精確值。通過實際的雙能X射線異物檢測機進行測試，驗證了使用本文提出的優(yōu)化數(shù)值算法得到的雙能值和分類邊界能夠更有效區(qū)分物質(zhì)，再次證明了該算法的實用性與有效性。

1 王琪, 陳志強, 鄔小平, 等. X射線安全檢查技術(shù)綜述[J]. CT理論與應(yīng)用研究, 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 肖筱南. 現(xiàn)代數(shù)值計算方法[M]. 北京: 北京大學(xué)出版社, 2003: 203-206

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

6 張德豐. Matlab數(shù)值分析與應(yīng)用[M]. 北京: 國防工業(yè)出版社, 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射線物質(zhì)識別的蒙特卡羅模擬[D]. 吉林: 吉林大學(xué), 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年于南昌大學(xué)獲博士學(xué)位，研究方向為X射線安檢相關(guān)技術(shù)

陳瑞燾，E-mail: 082663@#edu.cn

2014-03-27，

2014-04-16