曹大泉，王雅霄，闕介民，孫翠麗，魏存峰，魏 龍

(1.中國科學(xué)院 高能物理研究所 核分析技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100049；

2.北京市射線成像技術(shù)與裝備工程技術(shù)研究中心，北京 100049)

X射線薄板層析成像(CL)系統(tǒng)是一種非同軸掃描的計(jì)算機(jī)斷層成像系統(tǒng)，專門針對板狀構(gòu)件進(jìn)行結(jié)構(gòu)成像和缺陷檢測。CL與常規(guī)錐束CT的差異在于轉(zhuǎn)臺的轉(zhuǎn)軸不再垂直于X射線源焦點(diǎn)與探測器中心的連線(稱為光軸)，而是成一定的傾斜角度。CL的掃描方式解決了CT掃描板狀構(gòu)件時(shí)的兩個(gè)難題：1) 某些角度穿透困難，圖像信噪比低；2) 高空間分辨率下無法完成掃描。另一方面，這種特殊的掃描結(jié)構(gòu)使得CL的投影數(shù)據(jù)不滿足數(shù)據(jù)完備性條件，因此，傳統(tǒng)的解析重建方法容易出現(xiàn)不同斷層間圖像信息混疊。

實(shí)際數(shù)據(jù)獲取系統(tǒng)與圖像重建數(shù)學(xué)模型之間的差異導(dǎo)致產(chǎn)生各種圖像偽影。由于CL的射線源同樣采用X光管，而非重建模型中的單色光，因此導(dǎo)致產(chǎn)生射束硬化偽影。硬化偽影的校正方法很多，如濾片預(yù)硬化法、多項(xiàng)式擬合法、蒙特卡羅法、雙能校正法等。這些校正方法均對硬化偽影有一定的校正效果，不過都有一定的限制條件。CL檢測的目標(biāo)樣品形狀多不規(guī)則(如電路板的焊點(diǎn))，不易提取準(zhǔn)確邊界，不適合采用重投影多項(xiàng)式擬合法[1]；材質(zhì)多樣且多為復(fù)合材料，缺乏先驗(yàn)知識，不宜使用蒙特卡羅法模擬校正[2]；濾片預(yù)硬化法簡單實(shí)用，但校正效果不及其他方法理想[3]；雙能校正法成本高，在實(shí)驗(yàn)室不易實(shí)現(xiàn)。考慮到射束硬化偽影廣泛存在于CL重建圖像中，嚴(yán)重降低了圖像質(zhì)量，因此需要一種不同于上述射束硬化偽影校正的方法來校正CL圖像的偽影。文獻(xiàn)[4]提出一種校正常規(guī)CT圖像射束硬化偽影的方法，該方法在SART重建算法中加入射束硬化物理過程，通過準(zhǔn)確表達(dá)射線投影獲取模型來校正CT圖像的硬化偽影。本工作進(jìn)一步研究該校正方法，將該方法用于CL圖像的硬化偽影校正。

1 校正方法及數(shù)據(jù)獲取

1.1 CL硬化偽影校正方法

SART算法是一種對ART算法進(jìn)行改進(jìn)的代數(shù)迭代重建算法[5]。與ART對每條射線依次更新不同，SART對每個(gè)投影角度下所有的射線同時(shí)更新，因此重建速度更快。SART迭代重建首先將被重建區(qū)域賦初值0，對該區(qū)域做前投影，將前投影與真實(shí)投影的差值進(jìn)行權(quán)重的反投影。然后將反投影的結(jié)果迭加到初始被重建區(qū)域，其體素更新過程可表示為：

λ

(1)

考慮到SART在更新每個(gè)重建點(diǎn)像素時(shí)，將每條射線上的衰減值按權(quán)重進(jìn)行累加，權(quán)重值即為射線在每個(gè)重建網(wǎng)格中的交線長度，因此可將射束硬化校正的模型引入SART重建的過程，稱該種圖像重建算法為SART-BHC算法。SART-BHC使用射線在重建點(diǎn)的有效衰減系數(shù)累加代替衰減系數(shù)的直接累加，可更準(zhǔn)確地反映多色光譜投影數(shù)據(jù)獲取的過程，進(jìn)而在重建圖像時(shí)抑制射束硬化偽影。有效衰減系數(shù)與多色光譜、射線在物體中穿過的長度及物體本身材質(zhì)相關(guān)，可表示為[6]：

(2)

其中：α為校正強(qiáng)度系數(shù)；β為能量相關(guān)系數(shù)；Eeff為有效衰減系數(shù)。

將式(2)代入SART迭代公式可得SART-BHC算法更新公式[4]：

(3)

其中，N為重建體素個(gè)數(shù)。

(4)

(5)

其中：SO為光源到物體旋轉(zhuǎn)中心的距離；OD為物體旋轉(zhuǎn)中心到探測器中心的距離；φ為中心射線與物體旋轉(zhuǎn)軸的夾角；θ為物體繞z″的旋轉(zhuǎn)角；di為第i條射線與探測器的交點(diǎn)。

在物體坐標(biāo)系下，得到每個(gè)投影角度下射線源及每個(gè)探測器像素的位置坐標(biāo)，就可利用射線驅(qū)動(dòng)的Siddon算法[7]計(jì)算系統(tǒng)矩陣ω。

與解析重建算法(如濾波反投影(FBP)重建算法)需要完備的投影數(shù)據(jù)相比，SART-BHC在投影數(shù)據(jù)不完備的情況下也能得到較好的重建結(jié)果，因此更適用于CL這種無法獲取完備投影數(shù)據(jù)的掃描系統(tǒng)。SART-BHC的射束硬化偽影校正方法通過將有效衰減系數(shù)引入迭代過程，不需被測物體和光譜的信息，不需對被測物體進(jìn)行刻度實(shí)驗(yàn)以建立硬化校正模型。與標(biāo)準(zhǔn)SART重建相比，不用增加額外的迭代次數(shù)。此外，由于以往錐形束CT射束硬化偽影校正方法無法適用于CL掃描結(jié)構(gòu)[8]，該校正方法可很好地應(yīng)用于CL這種特殊的掃描結(jié)構(gòu)。由于很多工業(yè)部件無法獲取準(zhǔn)確的成分信息，或無法進(jìn)行刻度實(shí)驗(yàn)以建立硬化偽影校正模型，本文提出的方法可很好滿足工業(yè)部件硬化偽影校正的需要。

圖1 CL系統(tǒng)幾何示意圖

1.2 模擬數(shù)據(jù)獲取

模擬數(shù)據(jù)不受實(shí)驗(yàn)硬件條件不足的限制，又可排除實(shí)驗(yàn)中各種因素的干擾，可方便準(zhǔn)確地驗(yàn)證校正算法。本文采用蒙特卡羅方法在Geant4環(huán)境下搭建CL掃描系統(tǒng)模型，獲取投影數(shù)據(jù)以驗(yàn)證算法的有效性。

首先，在Geant4環(huán)境中利用電子打靶獲取X射線光譜，模擬參數(shù)列于表1。模擬空間中的幾何關(guān)系示于圖2。

表1 模擬電子打靶主要參數(shù)

圖2 模擬電子打靶的幾何關(guān)系

高速電子(160 keV)轟擊在鎢靶上，大部分電子與鎢靶作用后散射出去，極少部分電子與鎢靶作用后產(chǎn)生X射線，將探測器收集到的X射線繪制成能譜，結(jié)果示于圖3。

圖3 160 keV電子打鎢靶產(chǎn)生的X光譜

圖4 CL掃描模型

按圖3光譜的分布發(fā)射一定數(shù)量的X光子，就相當(dāng)于在模擬環(huán)境中建立一連續(xù)譜分布的X射線源。將該光譜作為X射線源，在Gate環(huán)境中建立CL掃描模型(圖4)，模型參數(shù)列于表2。

表2 CL投影模型主要參數(shù)

利用上述模型，樣品繞轉(zhuǎn)軸旋轉(zhuǎn)360°，獲取間隔均勻的60幅投影圖，圖5示出其中1幅。

圖5 模擬CL掃描的投影圖

1.3 實(shí)驗(yàn)數(shù)據(jù)獲取

a——FXE-225.99型微焦點(diǎn)X光機(jī)；b——PaxScan? 2520V型探測器；c——Huber轉(zhuǎn)臺；d——偏轉(zhuǎn)轉(zhuǎn)臺

本實(shí)驗(yàn)在1臺225 kV高分辨的CL系統(tǒng)上進(jìn)行，實(shí)驗(yàn)平臺示于圖6。CL系統(tǒng)搭建在光學(xué)平臺上，主要包括1臺FXE-225.99型微焦點(diǎn)X光機(jī)，1臺PaxScan? 2520V型探測器，1臺加裝1個(gè)很薄的有機(jī)玻璃圓筒的Huber轉(zhuǎn)臺，以及位于Huber轉(zhuǎn)臺下方用于控制偏轉(zhuǎn)角的偏轉(zhuǎn)轉(zhuǎn)臺。為了與模擬數(shù)據(jù)更好符合，加工一鋁質(zhì)圓盤，圓盤直徑30 mm，高10 mm。調(diào)整CL系統(tǒng)的參數(shù)(表3)，對鋁盤進(jìn)行CL掃描，得到圖7所示的投影圖。

表3 實(shí)驗(yàn)參數(shù)

圖7 實(shí)驗(yàn)投影圖

2 結(jié)果分析

2.1 模擬數(shù)據(jù)校正結(jié)果

為證明SART-BHC校正方法的可行性，作為對比，將蒙特卡羅模擬獲得的投影利用FBP算法重建成斷層圖。FBP是一種最普通的重建算法，由于其算法簡單、重建速度快，廣泛應(yīng)用于CT圖像重建領(lǐng)域。FBP的重建參數(shù)列于表4，重建的CL結(jié)果示于圖8。由圖8可看出，CL掃描重建出的斷層圖像有明顯的杯狀偽影。灰度曲線的底部較平，分析可知，這個(gè)杯狀偽影不僅來自于硬化，還有一部分來自于層間混疊效應(yīng)，即CL的傾斜掃描方式使得圖像重建時(shí)不同層之間的信息相互干擾。層間混疊效應(yīng)在本實(shí)驗(yàn)中體現(xiàn)在把本屬于鋁盤內(nèi)部體素對投影的貢獻(xiàn)，重建圖像時(shí)一部分丟失在鋁盤外部，這樣產(chǎn)生層間混疊偽影，同時(shí)降低了內(nèi)部的灰度，加重了杯狀偽影。

表4 模擬數(shù)據(jù)的FBP算法重建參數(shù)

a——鋁盤的斷層圖像；b——斷層圖像灰度曲線

利用SART-BHC算法重建蒙特卡羅模擬的投影數(shù)據(jù)，重建參數(shù)列于表5，SART-BHC算法重建結(jié)果示于圖9。由圖9可見，模體的斷層圖像不再存在杯狀偽影，灰度曲線的趨勢變?yōu)槠街?，與實(shí)際情況相符。因此，應(yīng)用SART-BHC算法對CL掃描方式所產(chǎn)生的杯狀偽影有很好的校正效果，證明對該算法改進(jìn)的正確性及該算法對CL掃描方式的杯狀偽影校正的有效性。

表5 模擬數(shù)據(jù)的SART-BHC算法重建參數(shù)

a——鋁盤的斷層圖像；b——斷層圖像灰度曲線

2.2 實(shí)驗(yàn)數(shù)據(jù)校正結(jié)果

為進(jìn)一步證明SART-BHC校正方法在CL中的實(shí)用性，對CL實(shí)驗(yàn)平臺上獲取的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行校正研究。首先將在CL實(shí)驗(yàn)平臺上獲取的實(shí)驗(yàn)數(shù)據(jù)利用FBP算法重建斷層圖像，重建參數(shù)與表4所列參數(shù)相同，重建結(jié)果示于圖10。為與模擬數(shù)據(jù)直觀對照，從整幅斷層圖像的正中選取一部分，并壓縮至與模擬圖像同等大小。從圖10可看出，實(shí)驗(yàn)數(shù)據(jù)的FBP重建結(jié)果與圖8所示的模擬數(shù)據(jù)結(jié)果一致，同樣有嚴(yán)重的杯狀偽影，嚴(yán)重降低了圖像質(zhì)量，應(yīng)予以校正。

a——鋁盤的斷層圖像；b——斷層圖像灰度曲線

利用SART-BHC算法對CL的實(shí)驗(yàn)投影數(shù)據(jù)進(jìn)行重建，主要重建參數(shù)列于表6，校正結(jié)果示于圖11。從圖11的重建結(jié)果可看出，鋁質(zhì)圓盤的斷層圖像不再存在杯狀偽影，灰度曲線較圖10所示的結(jié)果更加均勻，與實(shí)際相符?？梢姡琒ART-BHC應(yīng)用在CL中可有效地校正杯狀偽影，得到更好的圖像質(zhì)量。

表6 實(shí)驗(yàn)數(shù)據(jù)的SART-BHC算法重建參數(shù)

a——鋁盤的斷層圖像；b——斷層圖像灰度曲線

3 結(jié)論

本文利用蒙特卡羅模擬及CL實(shí)驗(yàn)平臺分別獲取了CL掃描方式的模擬和實(shí)驗(yàn)的投影數(shù)據(jù)，并利用基于SART迭代重建的硬化偽影校正算法(SART-BHC算法)對兩組投影數(shù)據(jù)分別進(jìn)行重建，均得到了硬化偽影消除的斷層圖像。具體結(jié)論如下。

1) CL成像系統(tǒng)同樣存在射束硬化現(xiàn)象，重建圖像存在杯狀偽影。杯狀偽影部分來自射束硬化，也有一部分來自層間混疊效應(yīng)。

2) 通過對SART-BHC的研究及對CL系統(tǒng)矩陣的推導(dǎo)，實(shí)現(xiàn)了應(yīng)用于CL的改進(jìn)SART-BHC算法。該算法通過將有效衰減系數(shù)引入迭代過程，在不增加重建迭代次數(shù)的情況下，可以很好消除CL圖像的射束硬化偽影。

3) 與以往的硬化偽影校正算法不同，SART-BHC在校正過程中不需被測物體和光譜信息，不需對被測物體進(jìn)行刻度實(shí)驗(yàn)以建立硬化校正模型。由于很多工業(yè)部件無法獲取準(zhǔn)確的成分信息，或無法進(jìn)行刻度實(shí)驗(yàn)，SART-BHC可很好地滿足工業(yè)部件硬化偽影校正的需要。

4) SART-BHC在校正射束硬化偽影的同時(shí)，能很好地適用于CL這種投影數(shù)據(jù)采樣不完備的掃描幾何，因此可較好地抑制CL層間混疊。

參考文獻(xiàn)：

[1] ZHANG Q H, LU H N, YANG M. Beam hardening correction using polynomial fitting based on reprojection[J]. Optical Technique, 2005, 31(4): 633-638.

[2] 曾剛，郁忠強(qiáng)，閻永廉. 基于蒙特卡羅模擬的射束硬化校正方法[J]. 高能物理與核物理，2006，30(2)：178-182.

ZENG Gang, YU Zhongqiang, YAN Yonglian. Beam hardening correction based on Monte Carlo simulation[J]. High Energy Physics and Nuclear Physics, 2006, 30(2): 178-182(in Chinese).

[3] 趙維. 錐束CT射束硬化偽影校正[D]. 北京：中國科學(xué)院高能物理研究所，2012.

[4] BRABANT L, PAUWELS E, DIERICK M, et al. A novel beam hardening correction method requiring no prior knowledge, incorporated in an iterative reconstruction algorithm[J]. NDT & E International, 2012, 51: 68-73.

[5] KAK A C, SLANEY M. Principles of computerized tomographic imaging[M]. New York: IEEE Press, 1988.

[6] QUE J M, CAO D Q, ZHAO W, et al. Computed laminography and the reconstruction algorithm[J]. Chinese Phys C, 2012, 36(8): 777-783.

[7] SIDDON R. Fast calculation of the exact radiological path length for a three-dimensional CT array[J]. Med Phys, 1985, 12: 252-255.

[8] ZHAO W, FU G T, SUN C L, et al. Beam hardening correction for a cone-beam CT system and its effect on spatial resolution[J]. Chinese Phys C, 2011, 35(10): 978-985.