雙能DR物質識別演算法在CT成像系統中的應用

作者：湯昕燁， 張麗， 陳志強， 高河偉， 張國偉

【摘要】  本文介紹一種適用於計算機斷層成像（CT）的雙能透視（DR）物質識別演算法，此演算法基於單能CT重建映像，通過對CT映像的分割及對斷層幾何資訊的提取，分塊重建材料有效原子序數及電子密度的分布。結合各種掃描軌跡的CT成像系統，可以實現有效物質材料識別。同時，對比於傳統的雙能CT方法，本方法結合單能CT重建映像，改善了雙能DR識別效果，能夠實現較為精確的物質識別，在安全檢查等應用領域有著現實意義。

【關鍵詞】  雙能DR； CT成像； 物質識別

    Abstract: A dualenergy digital radiography (DR) method assisted by segmentation of single energy reconstruction image is proposed for material recognition in Xray 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 dualenergy detectors in the system of CT scanning, which shows its significance in application in the field of security inspection.

    Key words: dualenergy 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機是電子束掃描機（Electronbeam 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)是探測器的能量響應函數，μ是線性衰減係數。

    因此，材料的物質特性可以通過已知的雙能量投影資料以及對其中線性衰減係數的基本分解進行重建。分別基於重建映像和投影資料正弦圖的postreconstruction方法及prereconstruction方法是兩種經典的雙能量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/20，散射反應截面

    aσSC≈ZKN(E)，

    式中，α=1/137為精細結構常數，m0為電子品質，0=83π(e2/mc2)2，KN(E)表示康普頓散射的KleinNishina截面。

    由此可以定義物質的特徵係數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)

    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, Energyselective Reconstructions in Xray Computerized Tomography ［J］. Phys. Med. Biol., 1976, 21(5): 733-744.

　　［3］Guowei Zhang, Li Zhang, Zhiqiang Chen, An HL Curve Method for Material Discrimination of Dual Energy Xray 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 HL Curve Method ［A］.Conf. Record of 2006 IEEE Medical Imaging Conference［C］.2006， M14462.

　　［5］莊天戈，CT原理與演算法［M］.上海交通大學出版社，1992.

　　［6］Sidky E Y, Zou Y, and Pan X. Volume image reconstruction from a straightline source trajectory, IEEE Symposium on NSSMIC ［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. Dualenergy Xray CT with CdTe array and its extension ［J］. Nuclear Instruments and Methods in Physics Research, 2005, A (548): 72-77.