王 輝，康南生，楊 波

(1.中國(guó)科學(xué)院合肥智能機(jī)械研究所，安徽合肥230031;2.中國(guó)科學(xué)技術(shù)大學(xué) 自動(dòng)化系，安徽 合肥230027)

0 引言

計(jì)算機(jī)斷層掃描成像技術(shù)(computerized tomography，CT)是計(jì)算機(jī)技術(shù)與放射學(xué)相結(jié)合產(chǎn)生的一門(mén)新的成像技術(shù)，它是在無(wú)損狀態(tài)下獲得被檢測(cè)斷面的灰度信息，以圖像的形式清晰準(zhǔn)確直觀地展現(xiàn)被檢測(cè)物體內(nèi)部結(jié)構(gòu)特征情況、材料的密度、有無(wú)缺陷等，廣泛地應(yīng)用于醫(yī)學(xué)、工業(yè)等領(lǐng)域。

毒品、爆炸物對(duì)個(gè)人和社會(huì)都有極為嚴(yán)重的危害，如何在火車(chē)站，海港和飛機(jī)場(chǎng)等重要通道處快速準(zhǔn)確地識(shí)別出毒品、爆炸物非常重要。利用雙能量算法進(jìn)行毒品、爆炸物檢測(cè)獲得廣泛的應(yīng)用，根據(jù)不同材料，不同密度的物質(zhì)對(duì)不同能級(jí)的X射線衰減程度不同，對(duì)雙能量系統(tǒng)獲取的2種圖像數(shù)據(jù)進(jìn)行綜合考慮計(jì)算，可以得到與被識(shí)別物體密度相關(guān)的值，可以提高系統(tǒng)的物質(zhì)分辨能力，在毒品、爆炸物檢測(cè)領(lǐng)域中起到重要作用。

國(guó)內(nèi)外學(xué)者對(duì)雙能量重建算法做了相關(guān)的研究，但基本都是建立在單色譜雙能量的基礎(chǔ)上的，所得到的透射圖像與X射線穿過(guò)物質(zhì)的厚度無(wú)關(guān)，但是在實(shí)際應(yīng)用中，所用的系統(tǒng)基本都是連續(xù)譜X射線，這樣就必須考慮待測(cè)樣品的厚度變化對(duì)探測(cè)結(jié)果的影響。因此，常用的單色譜雙能量重建算法并不能在連續(xù)譜雙能量系統(tǒng)中得到應(yīng)用［1］。

Alvarez RE和Macovsky A于1976年首次提出經(jīng)典的雙能量X射線斷層掃描基本重建算法，它采用牛頓—拉夫遜迭代法去迭代2個(gè)非線性雙能量方程組來(lái)重建原子序數(shù)和物質(zhì)密度的分布，算法簡(jiǎn)單，但該方法耗時(shí)嚴(yán)重，而且，數(shù)字迭代過(guò)程中可能是發(fā)散的，因此，本文提出了一種新穎的基于雙能量曲線的重建算法，通過(guò)對(duì)探測(cè)器得到的高、低能圖像灰度值進(jìn)行曲線擬合，實(shí)現(xiàn)對(duì)物質(zhì)電子密度的仿真重建。實(shí)驗(yàn)結(jié)果表明:該算法效果很好，相對(duì)誤差較?。?］。

1 投影重建與FBP算法

投影重建的理論基礎(chǔ)是基于Lambert-Beers定律

式中 I0為入射X射線強(qiáng)度，I為出射X射線強(qiáng)度，Δx為厚度，μ為材料的線性衰減系數(shù)，通常μ隨著X射線的能量改變，并因所選材料的不同而不同。

濾波反投影(FBP)算法的基礎(chǔ)是傅立葉切片理論，物體在某個(gè)角度上得到的平行投影的傅立葉變換，等于該物體在同一角度下進(jìn)行的二維傅里葉變換的一條直線，如圖1所示。

圖1 傅立葉切片理論示意圖Fig 1 Diagram of Fourier slice theory

基于傅立葉切片理論［3］，只需要獲取各個(gè)方面的投影數(shù)據(jù)的一維傅立葉變換，并將它們拼在一起得到整個(gè)圖像的二維傅立葉變換，然后，做二維傅立葉反變換就可以得到重建圖像。但由于傅立葉空間中產(chǎn)生的采樣模式不是笛卡兒坐標(biāo)的，所以，投影采樣是在極坐標(biāo)柵格上，為了執(zhí)行二維傅立葉反變換，這些采樣不得不被柵格化到一個(gè)笛卡兒坐標(biāo)系中，然而，在頻率域空間的差值會(huì)給重建圖像造成漂移、偽影等干擾，所以，直接的差值傅立葉變換是不可取的。

物體的二維傅立葉變換是由許多一維傅立葉變換拼接起來(lái)的。由于通過(guò)投影采樣得到的一維傅立葉變換是長(zhǎng)條狀，所以，拼接的過(guò)程必須加權(quán)，這個(gè)加權(quán)的過(guò)程就是加權(quán)濾波，而整個(gè)重建過(guò)程，一般稱(chēng)為濾波反投影，如圖2所示［4］。

圖2 濾波反投影示意圖Fig 2 Diagram of filtered back-projection

2 雙能量重建算法

2.1 雙能量重建算法的基本原理

Alvarez R E和Macovsky A［5］在1976年發(fā)表的文章中指出給定物質(zhì)對(duì)光子能的線性衰減系數(shù)μ在30～100kV范圍內(nèi)，可以表達(dá)為康普頓散射和光電效應(yīng)衰減的線性組合，光電效應(yīng)衰減成分很大程度地取決于被探測(cè)物質(zhì)的原子序數(shù)，而康普頓散射給出的是物質(zhì)的電子密度分布，它正比于物質(zhì)的密度

式中 ρ為物質(zhì)密度，g/cm3;A為原子質(zhì)量，n=4.5。

設(shè)S(E)為入射能量譜，則連續(xù)譜X射線CT投影公式［6］

設(shè)SH(E)，SL(E)為入射的高、低能量譜，則

雙能量算法實(shí)質(zhì)就是通過(guò)探測(cè)到的雙能量灰度值TH，TL值［7］，由以上兩式得到A1，A2的值，然后，利用 FBP重建算法重建得到 a1，a2，最后，由 a1=ρZn/A，a2=ρZ/A，可以得到被探測(cè)物體的原子序數(shù)和電子密度ρZ/A，從而可以得到物體的相對(duì)密度。

2.2 基于雙能量曲線的重建算法

在雙能量CT系統(tǒng)中，當(dāng)探測(cè)有效原子序數(shù)為Z的某樣品時(shí)，當(dāng)樣品的厚度tm變化時(shí)，可得到一系列點(diǎn)集(TH，TL)，每一組(Z，tm)都對(duì)應(yīng)一組灰度值(TH，TL)，形成一條TH-TL曲線，不同的Z值對(duì)應(yīng)不同的曲線。因此，雙能量曲線重建算法要首先標(biāo)定出各樣品的雙能量曲線圖TH-TL。

算法的具體實(shí)現(xiàn)步驟:

1)曲線標(biāo)定:標(biāo)定出若干違禁物品的雙能量曲線圖，在給定的高，低能量值下，對(duì)5種厚度的待探測(cè)樣品進(jìn)行分別透射，由得到的5種厚度下的高，低雙能量值，繪出待探測(cè)樣品的雙能量曲線圖TH-TL。

2)反向查找:當(dāng)對(duì)某樣品進(jìn)行探測(cè)時(shí)，由探測(cè)到的雙能量值(TH，TL)，根據(jù)已標(biāo)定出的雙能量曲線圖TH-TL，反向去查找得到該樣品的有效原子序數(shù)和厚度信息，即(Z，tm)。

3)計(jì)算投影數(shù)據(jù)集:由得到的Z，tm值計(jì)算投影數(shù)據(jù)集

4)FBP重建:由得到的投影數(shù)據(jù)集A1，A2，利用FBP重建算法重建出 a1，a2:a1=ρZn/A，a2=ρZ/A。

5)計(jì)算待測(cè)樣品密度:a1，a2是與待測(cè)樣品有效原子序數(shù)和電子密度相關(guān)的量，由此可以得到待測(cè)樣品的密度。

3 實(shí)驗(yàn)與結(jié)果

取樣品氯化鈉(Zeff=14.88)，糖(Zeff=8.876)，TNT(Zeff=7.104)并分別裝在自制的直徑為6，8，10，12，14 mm的圓柱形泡沫器皿中，分別做雙能量透射，實(shí)驗(yàn)設(shè)備為skyscan1076 CT機(jī)。

本文中雙能量系統(tǒng)選取40，100 keV，10 W，掃描角度范圍為0°～360°，旋轉(zhuǎn)角度為1°，空間分辨率為18μm。由得到的投影數(shù)據(jù)，繪制4種物質(zhì)的TH-TL曲線圖，如圖3所示。

圖3 雙能量曲線圖Fig 3 Dual-energy curves

標(biāo)定出物質(zhì)的TH-TL曲線圖后，就可以對(duì)該樣品的任意厚度進(jìn)行透射，通過(guò)曲線圖可以反向查找得到樣品的任意探測(cè)點(diǎn)的厚度，進(jìn)而可以計(jì)算出樣品的密度，得到的數(shù)據(jù)如表1。

由表1結(jié)果可知，本文提出的算法可以很好地得到被探測(cè)物質(zhì)的密度，精度較高，相對(duì)誤差小于3%。

4 結(jié)論

本文針對(duì)連續(xù)譜雙能量X射線CT系統(tǒng)，提出了一種新穎的基于雙能量曲線的重建算法，以得到被探測(cè)物質(zhì)密度，并結(jié)合實(shí)驗(yàn)證明了該算法具有較高的精確度，相對(duì)誤差小于3%，但現(xiàn)實(shí)系統(tǒng)中存在噪聲，因此，降低這些噪聲，可減小相對(duì)誤差。

表1 實(shí)驗(yàn)數(shù)據(jù)Tab 1 Experimental data

［1］孫麗娜.X射線安檢設(shè)備中物質(zhì)分類(lèi)識(shí)別方法的研究［D］.沈陽(yáng):東北大學(xué)，2003:2.

［2］Chuang K S，Huang H K.A fast dual-energy computational method using isotransmission lines and table look up［J］.Med Phys，1987，14(2):186－192.

［3］Jiang Hsieh.計(jì)算機(jī)斷層成像技術(shù)原理、設(shè)計(jì)、偽像和進(jìn)展［M］.張朝宗，譯.北京:科學(xué)出版社，2006.

［4］莊天戈.CT原理與算法［M］上海:上海交通大學(xué)出版社，1992:7.

［5］Alvarez E，Macovski A.Energy-selective reconstructionsin X-ray computerized tomogrphy［J］.Phys Med Biol，1976，21(5):733－744.

［6］Tsunoo T，Torikoshi M，ohno Y.Measurement of electron density and effective atomic number using dual-energy X-ray CT［J］.Medical Physics，2008，35(1):4924－4932.

［7］Zhang Guowei，Chen Zhiqiang，Zhang Li.Exact reconstruction for dual energy computed tomography using an H—L curve method［J］.IEEE Nuclear Science Symposium Conference Record，2006，6:3485－3488.