基于蒙特卡羅平臺的多能CT閃爍探測器研究

孫永剛，孔慧華，張海嬌

(中北大學 信息探測與處理山西省重點實驗室， 太原　030051)

隨著X射線CT成像技術的發展日益成熟， CT成像技術在國民生活各個領域得到廣泛的應用。目前CT成像的研究主要針對傳統單能CT和多能譜CT。傳統單能CT容易出現射束硬化現象，且圖像質量不高，因此多能CT應用而生。多能CT成像的仿真模擬成本高昂，利用計算機搭建平臺模擬CT成像是一種更好的方法。陳平等[1]使用計算機平臺對能譜進行濾波實現了CT成像的仿真。楊曉飛等[2]在計算機平臺上構建雙層探測器實現了對雙能CT的模擬。但這些模擬只是針對雙能傳統CT探測器，無法同時獲得多能譜下的投影。多層閃爍探測器可以探測到X射線中每一個光子的能量并判斷其所屬能量區間，進而對在不同能譜下的光子進行統計處理，并生成投影圖像。該探測器與傳統CT探測器的區別在于其能同時判斷出所屬能量區間的光子數，繼而得到每個能量段下的投影圖像，比傳統探測器效率更高且投影效果更好。本文所搭建的蒙特卡羅模擬仿真平臺利用Geant4及GATE數據包來模擬CT成像過程。

1　Geant4蒙特卡羅模擬平臺的搭建

CT成像系統主要包括X 射線源、掃描模體、X射線探測器、計算機操作系統等，其主要部分X射線探測器有線陣探測器和面陣探測器。平臺搭建時需分別對射線源、模體及探測器進行相應的設計，以最大程度地完成對系統的模擬。根據上述原理，構建出的錐束CT系統模型如圖1所示。

1.1　Geant4模擬軟件的安裝

蒙特卡羅模擬平臺所使用的Geant4是模擬粒子穿過物體的物理過程的一個數據包，它包含了X射線穿過物體時所有可能發生的物理現象。因此，相對其他模擬方法，Geant4數據包能更真實有效地進行CT模擬仿真。

本文主要在Windows系統上使用Cygwin虛擬環境軟件安裝Geant4模擬仿真數據包。Cygwin是Red Hat公司開發的Linux虛擬平臺軟件，可在Windows平臺上運行Linux系統和搭建好的蒙特卡羅模擬平臺。在安裝時需要Clhep以及VC++2010對Geant4數據包進行環境變量的設置，并在安裝完成后運行數據包自帶的實例來檢測數據包是否安裝完成[3]。

1.2　GATE及數據處理軟件ROOT的安裝

圖2　GATE結構

在Geant4數據包安裝完成后，需安裝CT掃描模擬軟件GATE以及數據處理軟件ROOT,對Geant4模擬得出的能譜數據進行數據分析并生成投影文件。同時，在GATE軟件中定義入射光子數(即光子強度I)模擬在一定能譜下的CT投影過程。

GATE是模擬CT成像的軟件，其結構包含用戶層、應用層、核心層及Geant4數據包。在用戶層可以將Gean4模擬得出的能譜編譯成mac數據文件，通過設置射線源、模體及探測器模擬投影過程，GATE結構如圖2所示。

ROOT是對數據進行統計分析的軟件包，利用ROOT可以對Geant4生成的能譜數據進行處理，得到可輸入到GATE的能譜文件。在Cygwin中安裝GATE軟件包，在設置好環境變量后，再安裝ROOT數據處理軟件。

1.3　CT成像模擬在平臺上的實現

在利用Geant4模擬X射線粒子穿過物體時，首先需要在Geant4程序中設置相關參數：在Primary Generator Action文件中定義入射粒子的種類、入射能量以及入射方向；然后在Detector Messenger文件中根據所需要模擬的物體材料及形狀進行定義[4]；最后在Stepping Action文件中編寫輸出文件，從而得到能譜數據，輸入到basic文件中，至此完成粒子在Geant4平臺上的參數設置。在此過程中，可以在Event Action中加入所需的物理反應過程，比如康普頓散射、瑞利散射及X射線的衰減，因而相比其他模擬仿真手段，Geant4模擬仿真能更加客觀、真實地模擬對X射線的探測過程。

利用Geant4得到basic文件后，需要用ROOT數據分析軟件對數據進行分析統計，得到在一定的入射管電壓下的能譜數據Source.mac文件，在該文件中設置點源的位置及對應入射電壓下的能譜數據。在GATE中，可設置γ光子、正負電子等相關粒子,包括一些放射性核素的衰變時間、放射粒子的放射過程。電磁過程用來模擬粒子間的電磁相互作用。λ稱為粒子作用的路徑長度，由式(1)表示。

(1)

其中σ(Zi,E)表示射線所穿過含有原子i材料的橫截面。

Geant4支持的粒子分布主要有線性分布、均勻分布、高斯分布以及用戶自定義的柱狀分布。在GATE中通過設置軸向及徑向夾角來定義二維和三維的發射粒子源。對于X放射源集合形狀的定義，GATE提供的射線源分布類型有“Point”“Plant”“Volume”“Surface”等，每一種分布類型又包含不同的放射源幾何體類型。這里為了便于數據的統計，一般采用的粒子分布為高斯分布，高斯模型為

(2)

Phantom.mac可以對掃描物體進行設置，設置的參數包含物體的材料、大小、形狀以及放置的坐標。對于物體材料的設置，可以在GATE提供的材料定義模塊Material.xml中定義，由物質的密度、原子數來定義需要的材料。在GATE平臺可以設置長方體(box)、球體(sphere)、圓柱體(cylinder)等幾種不同的模體形狀。如果想在一個模體中設置其他不同的模體，定義復雜的幾何結構，可在GATE中通過設置母體與子體的繼承關系來實現。這里設置一個母體為整個平臺(world)的圓柱體，并設置填充材料、幾何體大小及位置。類似地，可以繼續以圓柱體為母體設置其子體，方法與之前相同。這樣就完成了對CT系統中模體的設置。通常在GATE平臺上，模體在模擬的過程會自動圍繞中心軸旋轉，也可設置模體靜止不動，讓點源圍繞模體自動旋轉，從而得到投影數據。通過在Acquisition.mac文件中設置每次旋轉的時間間隔以及總時間決定旋轉次數，即投影數據文件的個數。最后，在Output.mac文件中設定輸出投影文件的數據格式，再利用Matlab程序對得到的投影數據進行轉換，從而得到投影圖像。

在Geant4模擬X射線穿過物體的過程中，會發生康普頓散射、瑞利散射及X射線衰減等物理現象。由光電效應可知，電子在原子中的束縛能只相當于紫外光子的能量，比X光子的能量小得多。因此，康普頓效應可看作X光子與自由電子的散射，電子在散射前靜止，其散射公式如下：

Δλ=λ-λ0=λe(1-cosθ)

(3)

(4)

式(3)稱為康普頓方程，Δλ稱為康普頓位移。式(4)稱為電子的康普頓波長[5]。

1.4　探測器對光子的識別及投影的生成

對于閃爍探測器的材料，一般采用碘化鈉(NaI)、硅酸釓(GSO)和碘化銫(CSI)來區分光子所屬的能量區間。針對探測器厚度的選擇，可根據楊曉飛等的結論：在設定每個晶體單元的厚度為1.5 mm時，探測器的探測效率最高。除了探測器厚度及材料外，還需要確定探測器的晶體設置對圖像的影響。首先設定所需的晶體探測器；然后利用Geant4模擬計數器對探測器接收到的光子數進行統計；最后按照GATE所支持的3種格式進行輸出，分別為ROOT、ASCII及RAW投影數據。在設定好所截取的閾值T1、T2后，Geant4會記錄在截取的能量范圍內探測器所接收到的光子數，一個光子數對應一個事件，即光子從射線源穿過物體到達探測器的過程，最終輸出模擬數據并記錄在Singles.data的二進制文件里。所有在模擬過程中產生的粒子記錄在Hits.data的文件里。因此，可根據所截取的能量段計算光子的強度及投影值：I(T1,T2)和P(T1,T2)。對每個事件都要調用輸出管理器，即使衰變的粒子沒有到達探測器。這里還需要注意的是，由于粒子的衰變受事件的隨機特性支配，所以處理過的衰變粒子的數量可能與預期的數字有輕微的不同，表示為：

N=A×Δt

(5)

其中：A表示記錄到的粒子數；Δt表示衰變時間[6]。

對于多層閃爍探測器來說，當探測器吸收光子時，將產生一個光譜，它的電壓與能量沉積成正比。所得到的光譜包含多個能量通道，探測器所能截取的能譜范圍為截取閾值。當接受的光譜能量高于閾值時，將記錄下一個光子數。這樣可以得到每個能量通道上所對應的光子數，設定多個能量閾值即可得到多個能譜下的光子數。對于給定一個能量閾值T1，探測器所接收到的光子強度為

(6)

其中：E∈[0,∞)表示能量；I0(E)表示從X射線源發射的光子強度；D(E)表示探測器的探測效率；I(T1)表示探測器接收到的光子強度；l表示射線路徑；μ(E,l)表示在能量E下的衰減系數。

當設定2個能量閾值0

(7)

與傳統CT相同，利用對數運算可以近似表示能量在T1

(8)

2　多層閃爍探測器的蒙特卡羅模擬實驗

2.1　多層探測器對多能CT的模擬

利用單個板狀探測器一次模擬得到的是一段能譜下的投影結果，而多層閃爍探測器區別于傳統探測器，可以同時得到多個能譜段下的投影數據，因此探測效率大幅提高。多層閃爍探測器需要多個不同材料的板狀晶體探測器來同時截取多個能譜，以此實現對多能的模擬。由于模擬所需的時間及其他條件較多，所以這里構建3層閃爍探測器為模擬對象，每層探測器分別由像素相同的3種晶體材料組成，利用晶體在不同條件下對能譜的識別來同時獲得多能譜下的CT投影。

2.2　模擬實驗的射線源及模體的相關參數

圖3　模擬實驗參數

點源是一個可發射光子數為15×106bq的X射線源。所掃描的模體是一個充滿水的圓柱體，其內部放置了兩個材料為硅(Si)和硫(S)且半徑均為1 mm的小球。模型的物距和像距都是150 mm,可繞中軸線逆時針旋轉，旋轉次數由每次旋轉的時間和所選取的最大旋轉角度決定，一般是1 s旋轉1°。模擬仿真模型及模體橫截面如圖3所示。模體中的具體參數如表1所示。

2.3　探測器對能譜的截取范圍

由于構成閃爍探測器的晶體材料、種類及厚度對探測器截取能譜會有影響，因此需要對多層閃爍探測器驗證所能截取的能譜范圍。首先確定所選單個晶體單元的形狀及晶體探測面積大小，這里的晶體探測面積可理解為晶體單元接收光子的面積。根據劉延等的研究理論，設定探測器是一個包含256×256個單元的板狀探測器，晶體探測面積的大小以及形狀都會影響探測器對光子的探測效率，所以每個探測單元的表面積不宜過大，以獲得更高的探測效率。在蒙特卡羅平臺設置由厚度為1.5 mm且探測面積均相同的正方體、圓柱和三棱柱探測晶體所組成的單層板狀探測器，分別模擬3種形狀的晶體對光子的探測效率，以此確定晶體的形狀[7]。探測效率如表2所示。

表1　模體中設置的參數

表2　不同形狀下的晶體探測效率

對比蒙特卡羅模擬結果可以看出:正方形晶體的探測效率比其他形狀的晶體探測效率要高，由此確定探測器的晶體單元為1.5 mm×1.5 mm×1.5 mm的正方形晶體。

之后即可驗證由NaI、GSO及CSI晶體構成的探測器對能譜的識別范圍。利用蒙特卡羅模擬平臺分別模擬由3種大小為1.5 mm×1.5 mm×1.5 mm的正方體晶體材料所構成的單層閃爍探測器，其像素為256×256。分別做出3個單層探測器的能譜圖像，并對結果進行高斯擬合，結果如圖4所示。

在對3個能譜圖像進行擬合后可得出探測器的探測效率會在其中一段能譜范圍內達到最高，從而確定探測器可識別的能譜范圍[8]。以NaI晶體組成的探測器所能識別的能譜范圍大致為20～50 keV，GSO晶體探測器的能譜識別范圍為60～90 keV，CSI晶體探測器的能譜識別范圍為90～120 keV。

圖4　單層探測器能譜圖像

2.4　模擬實驗

由于探測器所需的參數已知，因此多層閃爍探測器設計時在原NaI晶體探測器的基礎上再增加2個不同晶體材料的探測器，分別為GSO晶體和CSI晶體。3個探測器像素保持一致，都是由256×256個厚度為1.5 mm的探測單元組成，選取的管電壓在125 kV以下、入射光子數為15×106bq[9-10]。為了便于分析統計，物體旋轉360°度獲取360個投影數據，選取投影角度為10°時的投影數據。根據之前驗證的探測器能譜識別的范圍，設定截取能譜分別是20～50 keV、60～90 keV及90～120 keV，像素為256×256。得到的投影結果如圖5所示。

圖5　截取能譜20～50 keV、60～90 keV及90～120 keV時的投影

通過圖5可以看出：利用多層閃爍探測器同時得到了多個能譜下的投影圖，截取能量過高或者過低都會影響投影。能譜范圍在60～90 keV時探測器對光子吸收效率最高，所以對物質的識別效果最好。入射光子在不同材料中的衰減系數不同，因而在投影不同區域的對比度差別較大[11-12]。

2.5　與傳統單能CT的對比試驗

為了驗證平臺對多能譜CT成像的模擬效果，本文利用傳統單能CT在相同點源和入射能量下，對相同模體做出單能譜下的投影，并對模擬結果進行對比，具體結果如圖6所示。

圖6　單能譜下的投影模擬結果

可以看出，在使用3層閃爍探測器時截取的多能譜投影的效果遠優于傳統單能閃爍探測器得到的投影效果。傳統單能譜的投影比多能譜噪聲大，物質間的成像差異較小，不利于對物質的識別。因此，在蒙特卡羅模擬平臺上的模擬結果與實際情況吻合，結果可靠，可以對多能CT進行準確的模擬實驗。

2.6　實驗結果驗證

這里截取投影圖第69行的切片圖像，并利用FBP算法重建圖像做出相應灰度圖像。可以看出：重建的結果符合之前的模擬，驗證了實驗的合理性。實驗結果如圖7所示。

圖7　實驗結果

在利用蒙特卡羅平臺模擬CT成像的過程中，構建多層閃爍探測器同時截取到多個能譜所獲得的投影圖相對單個能譜下的圖像噪聲要小，而且對于物質的識別效果也更加顯著[13-15]。從投影圖可以看出：不同物質對光子的吸收是不同的，對光子吸收較多的物質在投影圖上亮度較高，反之則較暗。

3　結束語

在Linux系統下搭建蒙特卡羅模擬仿真平臺。在確定好模體的形狀、材料后，通過對多層閃爍探測器的晶體探測面積、晶體形狀以及能譜截取參數的確定，完成對多能譜CT模擬的準備工作。并在平臺上編寫相關的模擬程序，模擬出相關圖像，驗證了在蒙特卡羅模擬平臺的可用性。

通過在GATE中設置多層閃爍探測器，同時模擬出在多個能譜下的投影圖像，并利用FBP算法重建切片圖像和灰度圖像進行對比。結果表明：多能譜下的投影對物體的識別效果更好，實現了在蒙特卡羅仿真平臺上對多能譜CT成像過程的模擬。

參考文獻：

[1]郭蓉，潘晉孝，陳平，等.X射線譜的Geant4模擬及其濾波效應研究[J].核電子學與探測技術，2015，35(8):767-769.

[2]楊曉飛.基于GATE平臺的雙能CT成像系統模擬設計[D].沈陽：東北大學，2014.

[3]廖玉婷，王永波，許可欣，等.基于光子計數探測器的X線CT蒙特卡羅仿真[J].中國醫學物理學雜志,2016,33(2):122-127.

[4]Geant4 Collaboration Geant4 User’s Guide for Application Developers[Z].CERN,2011:159-163.

[5]VERHAEGHE J，READER A J.Lower variance FBP image reconstruction via new filter families[C]//Nuclear Science Symposium Conference Record.USA:IEEE，2010：2277-2281.

[6]郝佳，張麗，陳志強,等.多能譜X射線成像技術及其在CT中的應用[J].CT理論與應用研究，201l，20(1)：141-150.

[7]何小安，李朝光，易榮清，等.基于單光子計數的NaI(Tl)閃爍探測器靈敏度標定方法[J].強激光與粒子束,2015,27(3):176-180.

[8]GIMENEZEN，BALLABRIGA R，CAMPBELL M，et al.Characterization of Medipix3 with synchrotron radiati[J].IEEE Transactions on Nuclear Science，2011，58(1)：323-332.

[9]徐佳，宋玉收，孫光智，等.X射線能譜過濾的MC模擬[J].核電子學與探測技術，2014，48(7)：1281-1285.

[10] 陳平，郭蓉，潘晉孝,等.基于Geant4 的能譜濾波分離虛擬平臺研究[J].中國測試，2016，42(6)：65-69.

[11] XU Q,YU H,BENNETT J,et al.Image reconstruction for hybrid true-color micro-CT[J].IEEE Transactions on Biomedical Engineering,2012,59：1711-1719.

[12] HE P,YU H,THAYER P,et al.Preliminary Experimental Results from a MARS Micro-CT System[J].Journal of X-ray Science and Technology,2012,20：199-211.

[13] SIDKY E Y,KAO C M,PAN X.Accurate image reconstruction from few views and limited angle data in divergent-beam CT[J].Journalof X ray Scienceand Technology,2006,14:119-139.

[14] XU Q,MOU X,WANG G,et al.Statistical Interior Tomography[J].IEEE Transactions on Medical Imaging,2011,30(5):1116-1128.

[15] 王璽，肖沙里，張流強,等.基于CdZnTe晶體的2×2像素陣列探測器研究[J].光電子激光,2010,21(5):639-643.