李曉晴，劉錫明，苗積臣，吳志芳

(清華大學(xué) 核能與新能源技術(shù)研究院 核檢測技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100084)

計(jì)算機(jī)斷層成像技術(shù)(CT)是一種先進(jìn)的非接觸無損檢測技術(shù)，能清晰、準(zhǔn)確地對物體成像，可反映物體的內(nèi)部結(jié)構(gòu)、密度、缺陷位置和形狀等特性[1]。經(jīng)過了近50年的發(fā)展，CT技術(shù)已廣泛應(yīng)用于醫(yī)學(xué)診斷、工業(yè)探傷、港口安檢等諸多領(lǐng)域。隨著平板探測器技術(shù)的不斷成熟，基于平板探測器的錐束CT成為了重要的工業(yè)無損檢測手段[2]。

影響錐束CT成像質(zhì)量的因素較多，其中CT系統(tǒng)的幾何參數(shù)精確獲取程度是影響CT圖像重建質(zhì)量的關(guān)鍵因素。按照是否需要特定模體，現(xiàn)有的參數(shù)校正方法可分為無模體法[3-6]與模體法[7-10]兩大類。用于獲取幾何參數(shù)的標(biāo)定模體的結(jié)構(gòu)材質(zhì)各異，典型的標(biāo)定模體有西北工業(yè)大學(xué)使用的細(xì)琴弦和帶豁口的鋼尺[11]、Ford與Cho等[4,12]使用的將多個(gè)小鋼球等角排列而成的雙鋼珠平行圓環(huán)。相較于通過反復(fù)迭代重建圖像以獲取幾何參數(shù)的無模體法，模體法具有參數(shù)獲取速度較快、結(jié)果不受初值影響等優(yōu)點(diǎn)，使用更為廣泛，但模體法存在模體不易制作或集成度差，使用時(shí)可能需要更換模體、操作繁瑣等問題。

為解決上述問題，本文設(shè)計(jì)一種易于制作、操作簡便、集成化的模體，提出一種參數(shù)快速獲取方法：平行雙絲法。本文首先介紹該方法的原理，利用幾何投影關(guān)系推導(dǎo)錐束CT幾何參數(shù)的計(jì)算公式，然后通過仿真實(shí)驗(yàn)評估方法的準(zhǔn)確性與魯棒性，最后使用真實(shí)模體，在高精度錐束CT系統(tǒng)平臺上進(jìn)行實(shí)驗(yàn)驗(yàn)證。

1 參數(shù)定義

對于錐束CT，本文定義射線源到探測器的距離為DSD、射線源到旋轉(zhuǎn)中心的距離為DSO，旋轉(zhuǎn)中心在探測器上的投影位置坐標(biāo)為(u0,v0)，探測器繞軸u=u0的旋轉(zhuǎn)角度記為σ，探測器繞軸v=v0的旋轉(zhuǎn)角度記為ζ，探測器在u-v平面內(nèi)繞點(diǎn)(u0,v0)旋轉(zhuǎn)的角度記為η。圖1為錐束CT幾何參數(shù)的定義。這些參數(shù)對重建圖像的質(zhì)量影響程度不同。根據(jù)文獻(xiàn)[13-14]，影響最大的兩個(gè)參數(shù)分別為u0和η，本文將重點(diǎn)關(guān)注這兩個(gè)參數(shù)。默認(rèn)探測器的旋轉(zhuǎn)角度σ和ζ為0°。

2 平行雙絲法原理

本文提出的方法使用的模體是兩根平行金屬絲，故將方法命名為平行雙絲法。兩根平行絲記為l1、l2，設(shè)計(jì)時(shí)l1、l2的距離d已知。

a——CT系統(tǒng)DSD、DSO參數(shù)定義；b——探測器繞軸u=u0的旋轉(zhuǎn)角度σ；c——探測器繞軸v=v0的旋轉(zhuǎn)角度ζ；d——探測器繞點(diǎn)(u0，v0)的旋轉(zhuǎn)角度η圖1 錐束CT幾何參數(shù)定義Fig.1 Geometric parameter of cone-beam CT

平行雙絲法的基本思想為：將模體放在工件轉(zhuǎn)臺上，保持雙絲與轉(zhuǎn)臺平面垂直，轉(zhuǎn)臺旋轉(zhuǎn)1周，獲取模體360°投影圖像，對投影圖像進(jìn)行數(shù)據(jù)處理，計(jì)算出DSD、DSO、u0、η4個(gè)參數(shù)。關(guān)鍵步驟推導(dǎo)如下。

2.1 u0、DSD計(jì)算方法

操作時(shí)將模體置于轉(zhuǎn)臺上，兩絲垂直于轉(zhuǎn)臺表面，使絲與旋轉(zhuǎn)軸平行。由幾何知識可知，在旋轉(zhuǎn)1周的過程中，有且僅有兩個(gè)位置使得兩平行絲的投影重合，且這兩個(gè)位置的投影射線與中心線夾角相等，記為α。圖2為平行絲投影重合位置示意圖，圖2中選取過射線中心線且垂直于旋轉(zhuǎn)軸的平面作為截面。A1、A2分別為第1根細(xì)絲l1的兩個(gè)重合位置，B1、B2分別為第2根細(xì)絲l2的兩個(gè)重合位置，兩個(gè)虛線圓分別為l1、l2旋轉(zhuǎn)1周的軌跡，O為旋轉(zhuǎn)中心，D為射線源與旋轉(zhuǎn)中心連線在探測器上的投影位置，S為射線源的位置。用A0代表細(xì)絲l1的初始位置，轉(zhuǎn)臺從初始位置旋轉(zhuǎn)到雙絲重合的第1個(gè)位置旋轉(zhuǎn)角度為θ1、旋轉(zhuǎn)到第2個(gè)重合位置的角度為θ2。

圖2 平行絲投影重合位置示意圖Fig.2 Coincidence position of projection for parallel wires

在理想位置下，SD⊥DP1、SD⊥DP2、DP1=DP2。因而射線源中心線在探測器上的投影橫坐標(biāo)u0等于細(xì)絲l1在探測器的投影橫坐標(biāo)uP1與細(xì)絲l2在探測器的投影橫坐標(biāo)uP2的平均值：

(1)

射線源到探測器的距離為：

(2)

2.2 DSO計(jì)算方法

找到雙絲投影重合位置后，尋找平行絲所在平面與中心線垂直時(shí)的投影位置。由幾何知識可知，旋轉(zhuǎn)1周的過程中，這樣的位置有且僅有兩個(gè)，如圖3所示。A3、A4分別為雙絲所在平面與中心線垂直時(shí)細(xì)絲l1的兩個(gè)位置，B3、B4分別為垂直時(shí)細(xì)絲l2的兩個(gè)位置。兩虛線圓分別為l1、l2旋轉(zhuǎn)1周的軌跡。用A0代表細(xì)絲l1的初始位置，記轉(zhuǎn)臺從初始位置旋轉(zhuǎn)到平行絲所在平面與中心線垂直的第1個(gè)位置時(shí)角度為θ3，旋轉(zhuǎn)到第2個(gè)位置時(shí)角度為θ4。

圖3 平行絲連線垂直于中心射線時(shí)投影示意圖Fig.3 Projection diagram of perpendicular to central ray position for parallel wires

在旋轉(zhuǎn)1周的過程中，兩平行絲距離是不變的，A1B1=A2B2=A3B3=A4B4=d，旋轉(zhuǎn)中心到兩絲連線的距離r不變。由等比定理可得：

(3)

2.3 不同投影位置的角度關(guān)系

通過對圖像中兩根絲投影重合度的判斷可得到θ1、θ2，需要通過θ1、θ2計(jì)算出α、θ3、θ4，下面給出角度關(guān)系的推導(dǎo)公式。為了使推導(dǎo)最為簡易，選取旋轉(zhuǎn)中心到雙絲連線的垂足位置進(jìn)行計(jì)算。圖4為角度關(guān)系示意圖，其中圓為旋轉(zhuǎn)中心到雙絲連線垂足的旋轉(zhuǎn)軌跡，C0為初始時(shí)垂足的位置，C1、C2為雙絲投影重合時(shí)垂足的位置，對應(yīng)的轉(zhuǎn)臺旋轉(zhuǎn)角度為θ1、θ2，C3、C4為雙絲平面與中心線垂直時(shí)垂足的位置，對應(yīng)的轉(zhuǎn)臺旋轉(zhuǎn)角度為θ3、θ4；圓的兩條切線表示雙絲投影重合時(shí)穿過雙絲的射線，與中心線的夾角為α。

圖4 轉(zhuǎn)臺旋轉(zhuǎn)角度關(guān)系Fig.4 Rotation angle relation of turntable diagram

依據(jù)初始位置的不同，角度關(guān)系有4種情況：

(4)

2.4 η的計(jì)算方法

探測器繞點(diǎn)(u0,v0)旋轉(zhuǎn)角度為η時(shí)，投影情況如圖5所示，圖5中虛線表示理想探測器位置，實(shí)線表示旋轉(zhuǎn)后探測器位置?？芍D(zhuǎn)后兩絲的投影平行，與探測器v′軸的夾角等于探測器繞點(diǎn)(u0,v0)的旋轉(zhuǎn)角度η。

圖5 探測器繞點(diǎn)(u0,v0)旋轉(zhuǎn)前后投影對比示意圖Fig.5 Projection before and after detector rotate about point (u0,v0)

實(shí)驗(yàn)獲得雙絲投影后，計(jì)算投影圖中兩絲與探測器v軸的夾角，即得到探測器繞點(diǎn)(u0,v0)的旋轉(zhuǎn)角度η。按照式(5)進(jìn)行坐標(biāo)變換，變換后的幾何參數(shù)和理想情況下是一致的。

(5)

考慮實(shí)際情況下絲的放置位置與轉(zhuǎn)臺旋轉(zhuǎn)軸可能不完全平行，為了消除這種影響，計(jì)算雙絲旋轉(zhuǎn)1周獲得的所有投影圖像的η值，并計(jì)算這些η值的平均值，平均值即為系統(tǒng)的η值。

在實(shí)際測量校正中，由于旋轉(zhuǎn)中心投影橫坐標(biāo)v0對重建圖像質(zhì)量影響很小，所以不對其進(jìn)行精確計(jì)算，簡單采用探測器坐標(biāo)v軸中值作為v0值。

3 基于仿真方法的誤差分析

為驗(yàn)證在探測器繞軸u=u0旋轉(zhuǎn)，或繞軸v=v0旋轉(zhuǎn)情況時(shí)方法的魯棒性，利用Matlab進(jìn)行仿真模擬計(jì)算。模擬計(jì)算的參數(shù)為：DSD=1 150 mm，DSO=900 mm，d=150 mm，旋轉(zhuǎn)中心到兩絲連線的距離為100 mm。

3.1 σ誤差分析

圖6a為探測器繞軸u=u0旋轉(zhuǎn)時(shí)細(xì)絲的投影，虛線表示轉(zhuǎn)前探測器的位置，實(shí)線表示轉(zhuǎn)后探測器的位置；圖6b為Matlab模擬出的探測器無旋轉(zhuǎn)和繞軸u=u0旋轉(zhuǎn)角度σ為5°時(shí)的雙絲投影?？煽闯觯瑑山z投影依然平行，但投影在探測器上的長度、兩投影間距離均有改變。

圖6 探測器繞軸u=u0旋轉(zhuǎn)前后平行絲投影對比Fig.6 Projection before and after detector rotate around axis u=u0

本文計(jì)算了探測器旋轉(zhuǎn)角度σ在0.5°～10°范圍內(nèi)變化時(shí)，對幾何參數(shù)u0、DSD、DSO的影響，表1列出幾何參數(shù)的計(jì)算結(jié)果與精確值的相對誤差。

3.2 ζ誤差分析

圖7a為探測器繞v=v0軸旋轉(zhuǎn)時(shí)細(xì)絲的投影，虛線表示旋轉(zhuǎn)前探測器位置，實(shí)線表示旋轉(zhuǎn)后探測器位置；圖7b為計(jì)算機(jī)模擬出的探測器無旋轉(zhuǎn)和繞軸v=v0旋轉(zhuǎn)角度ζ為5°時(shí)雙絲的投影對比。可看出，兩絲的投影不再平行，會因探測器的傾斜而產(chǎn)生一定的夾角。

表1 不同σ下幾何參數(shù)相對誤差Table 1 Relative error of geometric parameter in different σ

圖7 探測器繞v=v0軸旋轉(zhuǎn)前后平行絲投影對比Fig.7 Projection before and after detector rotate around axis v=v0

本文計(jì)算了探測器旋轉(zhuǎn)角度ζ在0.5°～10°范圍內(nèi)變化時(shí)，對幾何參數(shù)u0、DSD、DSO的影響，表2列出了幾何參數(shù)計(jì)算結(jié)果與精確值的相對誤差。

表2 不同ζ下幾何參數(shù)的相對誤差Table 2 Relative error of geometric parameter in different ζ

由上述模擬結(jié)果可知，本方法受探測器繞u=u0軸旋轉(zhuǎn)角度的影響略高于受繞v=v0軸旋轉(zhuǎn)角度的影響，但影響均非常小，在肉眼可見的角度偏差內(nèi)對參數(shù)計(jì)算的影響可忽略不計(jì)，證明該方法魯棒性較好。

4 實(shí)驗(yàn)驗(yàn)證

4.1 實(shí)驗(yàn)方案

為檢驗(yàn)方法的實(shí)際效果，實(shí)際制作了模體并在高精度工業(yè)CT系統(tǒng)上進(jìn)行實(shí)驗(yàn)驗(yàn)證。選取0.3 mm鉭絲作為模體細(xì)絲，外殼采用有機(jī)玻璃，兩絲距離d為100 mm。系統(tǒng)參數(shù)及實(shí)驗(yàn)參數(shù)為：X射線源能量，320 keV；探測器像素尺寸，0.2 mm；探測器像素?cái)?shù)，2 048×2 048；掃描投影數(shù)，2 000。

4.2 操作步驟

實(shí)驗(yàn)操作步驟為：1) 將模體垂直放置于轉(zhuǎn)臺平面上，轉(zhuǎn)臺旋轉(zhuǎn)360°進(jìn)行投影，得到各角度的投影值；2) 根據(jù)投影圖中雙絲與v軸的夾角，計(jì)算探測器繞點(diǎn)(u0,v0)的旋轉(zhuǎn)角度η；3) 找到雙絲投影重合時(shí)的投影圖，記錄轉(zhuǎn)臺旋轉(zhuǎn)到兩個(gè)重合位置的旋轉(zhuǎn)角度θ1、θ2，記錄投影圖中雙絲投影的橫坐標(biāo)uP1、uP2；4) 根據(jù)式(4)計(jì)算兩平行絲所在平面垂直于中心線時(shí)的轉(zhuǎn)臺旋轉(zhuǎn)角度α、θ3、θ4；5) 找出轉(zhuǎn)臺旋轉(zhuǎn)角度為θ3的投影圖，記錄投影圖中兩絲的距離P3P4；6) 找出轉(zhuǎn)臺旋轉(zhuǎn)角度為θ4的投影圖，記錄投影圖中兩絲的距離P4P5；7) 根據(jù)式(1)～(3)，計(jì)算出幾何參數(shù)DSD、DSO和旋轉(zhuǎn)中心投影橫坐標(biāo)u0。

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

由于探測器v軸像素為2 048，故將v0設(shè)為1 024。表3列出幾何參數(shù)設(shè)計(jì)值與實(shí)驗(yàn)值的對比，可看出，實(shí)驗(yàn)值與設(shè)計(jì)值吻合較好。

表3 幾何參數(shù)設(shè)計(jì)值與實(shí)驗(yàn)值對比Table 3 Comparison between designed and experimental geometric parameters

圖8為實(shí)驗(yàn)值重建圖形，可看出，圖像質(zhì)量清晰度較好，說明平行雙絲法獲取的參數(shù)滿足實(shí)際需求。

圖8 實(shí)驗(yàn)值重建圖形Fig.8 Reconstruction image of experiment result

5 結(jié)論

針對現(xiàn)有錐束CT參數(shù)獲取模體法存在的模體不集成、制作要求較高、參數(shù)獲取過程操作較為繁瑣等問題，本文設(shè)計(jì)了一種操作簡便的集成化參數(shù)獲取方法：平行雙絲法。該方法將模體旋轉(zhuǎn)1周，獲取360°的投影圖像，即可準(zhǔn)確獲取DSD、DSO、u0、η。研究結(jié)果表明，該方法所使用的雙絲模體結(jié)構(gòu)簡單、易于制作、操作簡單、實(shí)用性強(qiáng)，參數(shù)獲取結(jié)果精確度高、魯棒性好，能較好地滿足高精度錐束工業(yè)CT實(shí)際應(yīng)用需求。