錐束計算機斷層成像系統快速定位方法

李明君，張定華，黃魁東，張順利，于清超

(西北工業大學 現代設計與集成制造技術教育部重點實驗室，陜西 西安710072)

計算機斷層成像(CT)通過對物體不同角度的射線投影測量而獲得被測物體內部斷層圖像信息，目前較實用的CT 系統主要是采用第3 代錐束CT(CBCT)掃描方式，其特點是射線源和平板探測器固定不動，被測物體繞旋轉軸旋轉一周，利用探測器采集的二維數字投影圖像序列，快速重構出掃描區域內的所有斷層圖像。它具有其它無損檢測方法無法替代的特點，在無損檢測領域具有廣闊的應用前景。

錐束CT 系統的硬件部分主要包括X 射線源、旋轉工作臺和平板探測器，無論采用哪種掃描方式，采用任何一種重建算法都要求知道定位參數的信息。在進行系統安裝定位時，應盡可能地保證系統的平面度和平行度，以減少定位參數測量過程的復雜度。然而，在錐束CT 硬件系統組裝時，由于機械定位誤差不可避免，探測器在安裝定位時，存在旋轉、偏轉和傾斜等現象;射線源和旋轉工作臺在安裝定位時，射線源焦點經旋轉中心難以投影到平板探測器的中心;射線源焦點到平板探測器的距離和射線源焦點到旋轉中心的距離(源心距)也會存在測量不夠精確的現象。錐束CT 硬件系統安裝定位時的空間幾何參數被稱作定位參數，主要包括射線源焦點到工作臺旋轉中心的距離Dso，射線源焦點到探測器的距離Dsd，射線源經旋轉工作臺轉軸到達平板探測器的投影中心(X0，Z0)以及探測器分別繞X、Y和Z 軸的角度α、β 和γ.若不能得到這些定位參數的精確值，將影響到重建算法中投影地址的計算，其誤差會引起CT 圖像中的偽影，導致圖像降質。因此在錐束CT 系統中的定位參數測定對重建圖像質量具有至關重要的意義。

對于以上問題，學者們進行了相關的研究。如針對投影中心的問題，文獻［1］研究了投影中心偏移對重建圖像質量的影響，提出了在圓柱表面軸向貼細鉛絲，在0°和180°投影，讀取鉛絲位置平均得到投影中心的方法。文獻［2］針對投影中心偏移的情況，提出了得到實際投影中心的算法，但并不完全適用于實際的復雜情況。針對源心距的問題，文獻［3］通過仿真實驗研究了射線源定位精度對CT 圖像質量的影響，只要射線源的位置精度誤差小于0.57°，重建圖像質量是可接受的。文獻［4］分析了射線源焦點到轉軸的距離的誤差對重建圖像質量的影響，當焦點距離遠大于被重建物體的尺寸時，其影響可忽略。文獻［5］利用自相關法對雙圓目標體二次成像，求得系統焦距值。針對探測器角度的問題，文獻［6］提出了一種基于細琴弦重建圖像評估的黃金分割迭代搜索法來確定定位參數，但搜索較耗時。在其它參數定位方面，文獻［7］研究了錐角對FDK［8］重建圖像質量的影響，得出FDK 算法只適用于錐角較小情況下(6°左右)的近似重建。文獻［9］分析了射線源到探測器的最小距離與檢測樣件的關系，在射線源到探測器距離選定的情況下，旋轉工作臺應該足夠靠近探測器。通過其它一些研究及仿真實驗，在經過按平行度和平面度要求定位安裝后，探測器的三個角度都很小，但其中繞Y 軸的角度β 對重建切片質量的影響較大，其余2 個角度的影響很小，可認為這2 個角度為0.

提出一種針對錐束CT 系統的快速定位方法，通過設計一個高精度的圓柱測具將無形的旋轉軸實體化，并通過此計算出定位參數中射線源焦點到工作臺旋轉軸的距離Dso，射線源焦點到探測器的距離Dsd，投影中心(X0，Z0)以及探測器繞Y 軸的角度β等主要定位參數。實驗表明，采用該方法測量校正后，使重建出的切片圖像偽影得到了明顯改善，證明了該方法的有效性。為錐束CT 的系統定位調試和圖像高質量重建提供了強有力的保障。

1 定位參數測量與計算方法

1.1 圓桿測具的設計

為實現定位參數的測量與計算，設計了一標準圓桿測具，它由均勻鋁質圓桿和圓盤底座構成，其模型如圖1所示。圓盤底座用于確保圓桿在旋轉工作臺上放置平衡且垂直于旋轉工作臺。圓桿測具尺寸可根據實驗要求按模型定制，本文所采用圓桿測具相關尺寸為:圓桿直徑2R =20.00 mm，在圓桿上開2 個寬度為4.00 mm 的退刀槽，其直徑2r=16.00 mm，它們相距L=60.00 mm 且分別位于中心射束的兩側。整個圓桿測具除要求較高的尺寸精度處，還要求有較小的形位公差，主要包括圓度、同軸度、平行度和垂直度，這些一般采用磨床精加工后可以滿足要求。

圖1 圓桿測具模型Fig.1 Round bar measurement tool

標準圓桿測具的設計中，兩側環帶用于測量計算Dso，Dsd和Z0，垂直圓桿用于測量計算X0和β，下面將分別介紹這5 個參數的測量計算方法。

1.2 Dso，Dsd和Z0 的測量與計算

錐束CT 系統各部件的位置關系及坐標系的定義，如圖2所示。當把圓桿測具放置在旋轉工作臺中心時，由射線源發出的1、2、3 號射線分別與圓桿測具相交的位置，如圖3所示。

圖2 錐束CT 硬件系統定位關系Fig.2 Localization relationship in CBCT

圖3 圓桿測具成像關系Fig.3 Imaging relationship of round bar measurement tool

因為圓桿與射線源的距離遠大于圓桿直徑，且圓桿直徑較小，因此可認為這3 條線與水平方向所成角度是近似相等的，設其所成角度為θ，則射線1穿過圓桿測具的長度為2R/cosθ，則射線2 穿過圓桿測具的長度為(R+r)/cosθ，則射線3 穿過圓桿測具的長度為2r/cosθ.根據Beer 定律可知

式中:I0為入射X 射線強度;I1，I2，I3分別為射線1，2，3 處對應的透射X 射線強度;l1，l2，l3分別為射線1，2，3 處穿過的材料厚度;μ 為材料線性衰減系數。

對(1)式進行對數運算處理后可得

式中:p1，p2，p3分別為射線1，2，3 處對應的投影圖像中的灰度值。

由此可得射線2 所得到的灰度值是射線1 與射線3 所得到灰度值的中值。另一方面，射線1 所得到的灰度值在其附近是最高的，射線3 所得到的灰度值在其附近是最低的。為了得到3 條射線在平板探測器上的準確位置，并計算射線源焦點到旋轉中心的距離Dso，射線源焦點到探測器的距離Dsd，以及中心射束在平板探測器上的位置Z0，采用以下處理步驟:

1)將圓桿測具放置在旋轉工作臺的中心處，為了減少噪聲的影響，采用多幅平均的方式，即采集多幅投影圖像，然后將這些圖像按對應像素進行灰度疊加并平均，得到一幅投影圖像G1，并按重建時的投影圖像尺寸進行裁剪。

2)保持圓桿測具不動，將旋轉工作臺往射線源或探測器方向平移一定距離d，按步驟1)方法采集得到一幅投影圖像G2.

3)為降低噪聲影響，將求對數后的投影圖像G1 按行進行投影像素灰度疊加，得到一列投影。

4)以該列投影的中心為起始點，逐像素向上查找第一個局部灰度最大值像素點(對應于射線1)，將此點標記為Pmax，判斷準則為該像素的相鄰5～9個像素的灰度均小于該像素的灰度;從Pmax逐像素向上查找第一個局部灰度最小值像素點(對應于射線3)，將此點標記為Pmin，判斷準則為該像素的相鄰5～9 個像素的灰度均大于該像素的灰度。

5)計算Pmax和Pmin像素灰度的平均值，由于Pmax與Pmin之間的投影像素灰度值為遞減變化，因此在這2 個像素之間查找與此灰度平均值最接近的2個像素點。

6)對這2 個像素點的灰度進行線性插值，根據灰度平均值可得到與射線2 所對應的像素點位置，它是一亞像素值，即計算精度可達到亞像素級，將此點標記為P1點。

7)采用與上述步驟3)～步驟6)完全類似的方法，以該列投影的中心為起始點，逐像素向下查找并計算可得到P2點。

8)對求對數后的投影圖像G2，采用與上述步驟3)～步驟7)類似的方法可得到P3點和P4點。

9)根據三角形的相似性(圖4)，可以得

圖4 Dso、Dsd、Z0 測量的成像幾何關系Fig.4 Imaging relationship for Dso，Dsd，Z0

1.3 X0 和β 的測量與計算

在測量X0和β 時，需要將圓桿測具放置在略偏離旋轉工作臺中心處，且旋轉工作臺軸線須包含在圓桿之內，另外使圓桿中心線與旋轉工作臺軸線所組成的平面與探測器平面接近平行，并標記此位置為0°，采集多幅投影圖像，然后將這些圖像按對應像素進行灰度疊加并平均，得到一幅投影圖像G3，并按重建時的投影圖像尺寸進行裁剪;接下來保持圓桿測具不動，將旋轉工作臺旋轉180°，同樣采集得到一幅投影圖像G4.

圖5 投影圖像G3 和投影圖像G4 的成像幾何關系Fig.5 Imaging relationship of projection images G3 and G4

投影圖像G3 和G4 的成像幾何關系，如圖5所示，其中圖5(a)為XY 空間視圖，圖5(b)為YZ 空間視圖。對于中心射束，當O1O2與探測器平面平行時，A1B1與A2B2重合，則顯然有A1B1=A2B2;而大多數情況下，O1O2與探測器平面不能完全平行時，在0°和180°位置處，OO1=OO2，直角△OO1H1與△OO2H2全等，得到O1H1=O2H2，對于2 個位置相同直徑的圓，由于弦心距相等，則弦長A1B1=A2B2，即穿過旋轉中心的X 射線在此2 個位置處穿過的零件長度是相等的。對于錐束射線成像中通過旋轉軸的其它射線，也可推導出A'1B'1=A'2B'2，即所有通過旋轉軸的射線，在0°和180°位置處所穿過的零件長度相等。因此，將多幅平均后的投影圖像G3 和G4 按對應像素進行灰度值相減并取絕對值，得到一幅減影圖像G5，則圖像G5 中旋轉軸投影所在位置的像素灰度值為局部最小值(0 或接近于0)，在其兩側的像素灰度值關于旋轉軸投影呈對稱分布。

為了得到旋轉軸投影直線在平板探測器上的數學表達式，以計算中心射束在平板探測器上的位置X0和探測器旋轉角β，采用以下處理步驟:

1)在減影圖像G5 在高度方向取等間距的若干行像素，對每一行像素，查找圓桿投影區域內的最小灰度值像素。

2)以最小灰度值像素為中心，在相應行上分別取其左右相鄰的2～7 個像素，得到連續的5～15 個像素。

3)以像素位置為自變量x，采用最小二乘二次多項式對這些像素灰度進行擬合，得到與行數量相同的多個二次多項式。

4)分別計算各個二次多項式的最小值，得到與最小值對應的X 值(解析計算，可得亞像素級精度)，結合所在行的Z 坐標，即得到與行數量相同的多個平板探測器上坐標點。

5)采用最小二乘直線擬合這些坐標點，得到的直線方程就是旋轉軸投影直線在平板探測器上的數學表達式。

6)計算直線傾斜角，可得探測器旋轉角β.

7)將前面計算所得的中心射束縱坐標值Z0代入直線方程，即可計算出中心射束橫坐標值X0(解析計算得亞像素級精度)。

2 定位參數測量及實驗驗證

將錐束CT 先進行預定位，硬件系統需保證一定的平面度及平行度等，射線源焦點與探測器之間的距離粗測約為1 350 mm，探測器與工作臺旋轉中心的距離約為250 mm.先將圓桿測具放置在錐束CT 系統的旋轉工作臺中心處，采集得到投影圖像G1;保持圓桿測具不動，將旋轉工作臺往射線源方向平移d =50.00 mm，采集得到投影圖像G2;然后將圓桿測具沿與探測器成像平面平行的方向平移3 mm左右，在0°位置采集得到投影圖像G3;保持圓桿測具不動，將旋轉工作臺旋轉180°，采集得到投影圖像G4.所得圖像按1 024 像素×1 024 像素尺寸進行裁剪，如圖6所示，下側黑條處為鉛板屏蔽。整個過程僅測量時需3～5 min，參數計算可自動完成，時間可忽略。最后利用上述采集的投影圖像，計算系統主要定位參數，包括射線源焦點到旋轉中心的距離Dso，射線源焦點到探測器的距離Dsd，中心射束在平板探測器上的照射位置(X0，Z0)以及探測器旋轉角β.

圖6 投影圖像G1、G2、G3、G4Fig.6 Projection images G1，G2，G3 and G4

按本文1.2 節中方法計算Dso，Dsd和Z0.后文中未標識單位處單位皆為像素，為保證計算結果準確性，系統對精度要求較高的計算結果都采用小數點后6 位有效數字。在求對數后的投影圖像G1 中計算得到P1=233.245 101，P2=817.660 577;在求對數后的投影圖像 G2 中計算得到P3=203.343 123，P4=846.659 914;根據2 次投影的相關 性，計 算 可 得 Dso=1 092.19 mm，Dsd=1 348.81 mm，Z0=529.931.然后，按本文1.3 節中方法計算X0和β.將投影圖像G3 和G4 按對應像素灰度值相減并取絕對值得到G5，如圖7所示;通過在G5 中取等間距的9 行像素，對每一行分別進行最小二乘二次多項式擬合，如其中第150 行擬合結果為g=1.906 8x2-1 932.771x+489 792，其最小值對應的x=506.81，則由該行得到的坐標為(506.81，150).將這些坐標點采用最小二乘直線擬合，擬合后所得到的直線方程Z=-342.699 3X +173 850.672 2，即為旋轉工作臺旋轉軸的投影直線在平板探測器上的數學表達式。則平板探測器的旋轉角為β°=-0.167 2°，將前面計算得到的中心射束縱坐標值Z0=529.931 代入直線方程，即可計算出中心射束橫坐標值X0=505.751 664.

圖7 減影圖像G5Fig.7 Subtraction image G5

為驗證上述所得的錐束CT 系統主要定位參數的準確性，對由該系統掃描的一條鋼絲采用這些參數進行重建，并與文獻［1］中方法計算所得投影中心X0=505 對比，其它參數按本文方法得到的結果進行計算。取重建切片的第512 層觀察其鋼絲重建灰度，如圖8所示，由灰度的收斂性可知，本文方法所得的鋼絲重建灰度的集中性較好;然后又分別采用2 種方法計算得到的定位參數對某型號陶芯膜模進行重建，如圖9所示，可見本文方法得到的切片圖像邊緣模糊度較好。以上結論也印證了文獻［1］的結論——測量的投影中心與實際的投影中心即使存在較小的誤差，都會使重建圖像變形并產生偽像。通過以上實驗結果，表明采用本文方法獲取的錐束CT 系統主要定位參數是準確可靠的。

圖8 鋼絲重建灰度分布圖Fig.8 Gray-scale distribution of wire reconstruction

圖9 陶芯膜模局部區域重建切片圖Fig.9 Local slice images of ceramic core wax mold

3 結論

綜上所述，定位參數的準確測量是保證錐束CT系統應用的必要前提和保障。本文首先分析了錐束CT 中定位參數對重建切片質量的影響，綜述了現有的定位參數確定方法的技術現狀及不足，然后針對錐束CT 系統中硬件定位參數快速測量的問題，設計了一個標準圓桿測具，并提出了基于此測具的定位測量計算方法。該方法首先通過前后移動圓桿測具，利用成像幾何關系測量計算Dso，Dsd和Z0;然后偏移圓桿測具后，通過0°和180°成像相減獲得減影圖像，利用中心射束穿越距離相等測量計算X0和β，實現了錐束CT 硬件系統的精確定位及準確測量;最后通過實驗結果表明，采用該裝置及其配套參數測量方法進行定位、測量及校正后，獲取的重建切片圖像質量有明顯的改善，本文方法實施簡單且提高了錐束CT 系統定位調試的精確度和效率。

References)

［1］王召巴.旋轉中心偏移對CT 重建圖像質量的影響［J］.兵工學報，2001，22(3):323 -326.WANG Zhao-ba.Effect of center deviation on CT reconstruction images［J］.Acta Armamentarii，2001，22(3):323 - 326.(in Chinese)

［2］蔡玉芳，王玨，楊德鴻.CT 三代掃描旋轉中心偏離的卷積反投影算法［J］.重慶大學學報，2002，25(9):55 -58.CAI Yu-fang，WANG Jue，YANG De-hong.Convolution backprojection reconstruction for ⅢCT scanning with a shift center of rotation［J］.Journal of Chongqing University，2002，25 (9):55 -58.(in Chinese)

［3］叢鵬，李志鵬，鄔海峰.射線源的定位精度對工業CT 圖像質量的影響［J］.原子能科學技術，2004，38(增刊):244 -247.CONG Peng，LI Zhi-peng，WU Hai-feng.Source position error influence on industry CT image quality［J］.Atomic Energy Science and Technology，2004，38(Suppl):244 - 247.(in Chinese)

［4］路宏年，楊民，張莉.射線源焦點距離誤差對CT 重建質量影響分析［J］.兵工學報，2003，24(1):65 -67.LU Hong-nian，YANG Min，ZHANG Li.A study on the reconstruction bias originating from error of focal distance of X-ray source［J］.Acta Armamentarii，2003，24(1):65 -67.(in Chinese)

［5］楊民，孫翠麗，李保磊.利用圖像相關法精確測量3D-CT 成像系統焦距［J］.北京航空航天大學學報，2007，33(10):1216 -1219.YANG Min，SUN Cui-li，LI Bao-lei.Method for measuring the focal distance of 3D-CT imaging system with high accuracy using image auto-correlation［J］.Journal of Beijing University of Aeronautics and Astronautics，2007，33(10):1216 -1219.(in Chinese)

［6］劉曉鵬，張定華，張豐收，等.錐束CT 系統安裝參數確定技術研究［J］.計算機工程與應用，2006，20:170 -173.LIU Xiao-peng，ZHANG Ding-hua，ZHANG Feng-shou，et al.Research on tracing technique of precise installation parameters in CBVCT system［J］.Computer Engineering and Applications，2006，20:170 -173.(in Chinese)

［7］孫宏宇，王召巴.三維CT 在錐角增大時對重建圖像質量的影響［J］.測試技術學報，2002，16(4):299 -301.SUN Hong-yu，WANG Zhao-ba.The effect of cone angle being increased on the 3D CT reconstruction capability［J］.Journal of Test and Measurement Technology，2002，16(4):299 -301.(in Chinese)

［8］Feldkamp L A，Davis L C，Kress J W.Practical cone-beam algorithm［J］.J Opt Soc Am A，1984，1 (6):612 -619.

［9］張朝宗.工業CT 技術參數對性能指標的影響——兼談如何選擇工業CT 產品［J］.無損檢測，2007，29(1):48 -52，29(2):84 -87，29(3):142 -145，172.ZHANG Chao-zong.Effect of technical data of an industrial CT on its performance:a guide for industrial CT system selection［J］.Nondestructive Testing，2007，29(1):48 -52，29(2):84 -87，29(3):142 -145，172.(in Chinese)