基于等徑外切圓的工業(yè)CT系統(tǒng)空間分辨率測試

齊子誠，倪培君，姜 偉，郭智敏，高紅俐

（1.浙江工業(yè)大學(xué)機(jī)械工程學(xué)院，浙江杭州310014；2.中國兵器科學(xué)研究院寧波分院，浙江寧波315103）

1 引 言

工業(yè)計算機(jī)斷層成像技術(shù)（Industrial Com?puterized Tomography，ICT）具有限制條件少、成像直觀、分辨率高的優(yōu)點，廣泛應(yīng)用于航空航天、兵器及汽車等行業(yè)的產(chǎn)品內(nèi)部質(zhì)量檢測［1-3］。目前，國內(nèi)工業(yè)CT系統(tǒng)已成系列，包括微焦點（納米）工業(yè)CT、15MeV大型工業(yè)CT系統(tǒng)等，系統(tǒng)的空間分辨率覆蓋40 lp/mm（微焦點CT）到0.5 lp/mm（大型加速器CT），密度分辨率最小可達(dá)0.1%?？臻g分辨率作為表征工業(yè)CT系統(tǒng)性能的一個重要參數(shù)，與缺陷檢測能力密切相關(guān)。因此，實現(xiàn)工業(yè)CT系統(tǒng)空間分辨率的有效測試對CT系統(tǒng)設(shè)計、研制、驗收和應(yīng)用等方面有著重要意義。

目前，工業(yè)CT系統(tǒng)空間分辨率測試方法按原理［4］可分為周期性結(jié)構(gòu)測試法和調(diào)制傳遞函數(shù)（Modulation Transfer Function，MTF）測試法。周期性結(jié)構(gòu)測試法［5-6］是在CT圖像中直接觀察線對、圓（方）孔陣列等周期性結(jié)構(gòu)，利用可識別的最小結(jié)構(gòu)尺寸來確定系統(tǒng)的空間分辨率，具有結(jié)果直觀、檢測效率高等優(yōu)點，但是測試能力和精度完全依賴于模體中周期性結(jié)構(gòu)的空間頻率范圍、間隔及加工精度。隨著系統(tǒng)空間分辨率的提高，測試模體的加工要求、難度顯著提升，影響了該方法的實用性。MTF法具有模體結(jié)構(gòu)簡單、加工要求低、檢測分辨率高等優(yōu)點。根據(jù)模體形狀不同，該方法可分為圓盤法［7-9］和刃邊法［10-15］，這兩種方法均利用模體在CT圖像上的邊緣灰度信息，利用計算機(jī)提取邊緣圖像灰度分布進(jìn)行微分運(yùn)算形成邊緣擴(kuò)散函數(shù)（Edge Spread Function，ESF），再經(jīng)過傅里葉變換計算獲得MTF曲線，用于表征工業(yè)CT系統(tǒng)的空間分辨率。在實際應(yīng)用中邊緣圖像灰度提取、平均降噪、插值擬合等過程較為復(fù)雜，影響了方法的實用性?；谏鲜鲅芯楷F(xiàn)狀，本文使用一對直徑相同的圓柱體形成連續(xù)變間距測試模型，建立了一種可等價于周期性結(jié)構(gòu)（線對卡模體）的工業(yè)CT系統(tǒng)空間分辨率測試方法。

2 空間分辨率測試模型

2.1 周期性結(jié)構(gòu)測試模型

工業(yè)CT掃描過程如圖1所示。將周期性結(jié)構(gòu)（線對卡模體）固定于轉(zhuǎn)臺上，射線源產(chǎn)生X射線穿透模體并被探測器接收，模體隨轉(zhuǎn)臺旋轉(zhuǎn)一周（360°），與此同時探測器實時采集射線穿透模體后的衰減信號，并傳輸至上位機(jī)實現(xiàn)CT圖像重建。

圖1 條形模體在工業(yè)CT掃描示意圖Fig.1 Scanning of strip phantom（line-to-card）in indus?trial CT

圖2 CT成像過程Fig.2 Process of CT imaging

工業(yè)CT掃描形成的檢測圖像是空間域里物體經(jīng)過系統(tǒng)調(diào)制后得到的結(jié)果。因此，可以利用調(diào)制度反映物體細(xì)節(jié)在背景中的對比度，進(jìn)而表征CT系統(tǒng)的空間分辨率。線對卡模體中不同空間頻率結(jié)構(gòu)的CT成像前后對比如圖2所示?？梢娊?jīng)過CT成像后結(jié)構(gòu)的空間頻率不變，對比度降低了，形成了衰減的近似正弦波。在各個頻段的周期性結(jié)構(gòu)中，分別將圖像最亮處（灰度最大值）記為Imax，圖像最暗處（灰度最小值）記為Imin，則調(diào)制度MT可表示為［16-17］：

為探求線對卡模體的CT成像過程，先建立周期性條形結(jié)構(gòu)的數(shù)學(xué)模型。在理想情況下，條形結(jié)構(gòu)（線對卡）可表示為矩形脈沖函數(shù)。為研究方便，分別對波峰、波谷進(jìn)行討論，如圖3所示。條形結(jié)構(gòu)中的波峰p(x)、波谷q(x)分別表示為［18］：

圖3 條形模體表現(xiàn)形式Fig.3 CT manifestation of strip phantom

其中：n∈Z，2ai為波峰（波谷）矩形的脈沖寬度，ε(x)為 單 位 階 躍 信 號 ，定 義 為由于工業(yè)CT成像系統(tǒng)是一種線性非時變系統(tǒng)［12］，因此其成像導(dǎo)致的對比度下降可近似為高斯退化函數(shù)模型g(x)。成像過程可表示為輸入函數(shù)與高斯退化函數(shù)的卷積，則條形結(jié)構(gòu)中波峰的工業(yè)CT成像過程表示為：

同理，波谷成像過程表示為：

在不考慮噪聲的情況下，卷積退化后形成的波峰（波谷）極值必然在x=0處，則波峰Imax（波谷Imin）表示為：

分析各個條形結(jié)構(gòu)之間的影響關(guān)系，簡化波峰Imax（波谷Imin）的表達(dá)式。假設(shè)工業(yè)CT成像過程的高斯退化函數(shù)模型為：

根據(jù)標(biāo)準(zhǔn)正態(tài)分布規(guī)律可知，在3倍標(biāo)準(zhǔn)偏差的區(qū)間內(nèi)，g(x)的積分函數(shù)面積約占整體面積的99.73%。當(dāng)波峰（波谷）寬度2ai等于標(biāo)準(zhǔn)差σ時，在3倍標(biāo)準(zhǔn)偏差的區(qū)間內(nèi)包含了3組線對，此時的MT值約為1.1%，可見對波峰（波谷）幅值影響較為顯著的是相鄰的兩組線對。因此，當(dāng)1.1%≤MT≤100%時，波峰Imax可近似為：

將式（6）、式（8）代入式（1），即可獲得線對卡CT成像后調(diào)制度MT的表達(dá)式：

2.2 等徑外切圓結(jié)構(gòu)模型分析

雙球法［19］是Carmignato等人提出的一種三維CT系統(tǒng)空間分辨率評估方法，主要是利用一個簡單結(jié)構(gòu)的參考試塊（兩個具有相同公稱直徑的點接觸校準(zhǔn)球體），對該參考試塊進(jìn)行CT重建，觀察接觸點的失真情況。在指定誤差范圍內(nèi)，利用可測量的最小球體直徑來確定CT系統(tǒng)的空間分辨率。該方法的優(yōu)點在于參考試塊結(jié)構(gòu)簡單，但是對于具有微米分辨率的工業(yè)CT系統(tǒng)則需要直徑在微米級的球體，球體試塊制造、校準(zhǔn)及使用難度顯著提升。Zanini等［20］提出了一種基于雙球的改進(jìn)空間分辨率測試方法，采用較大直徑的球體，利用局部自適應(yīng)算法提取球體邊緣，計算接觸點扭曲區(qū)域的高度來表征CT系統(tǒng)的空間分辨率，但是在水平方向上由于射束硬化、噪聲等原因，接觸點附近灰度分布特別復(fù)雜，影響扭曲區(qū)域高度的測量精度和穩(wěn)定性，測量結(jié)果發(fā)生振蕩。本文構(gòu)造了一種垂吊式圓柱體線接觸模體，用于改善CT檢測效果，測量結(jié)果等價于線對卡法，更適用于線陣工業(yè)CT系統(tǒng)的空間分辨率測試。

如圖4所示，兩個具有相同公稱直徑(D)的校準(zhǔn)柱體，使它們在一條線上相互接觸。由于CT重建過程造成接觸區(qū)域圖像畸變，該畸變區(qū)域的面積隨著CT系統(tǒng)空間分辨率的提高而減小。因此，兩個相互接觸的校準(zhǔn)柱體可作為一種狹縫寬度連續(xù)變化的矩形波測試卡［21］。

圖4 雙圓模體與空間分辨率的關(guān)系Fig.4 Relationship between double circular phantom and spatial resolution

圖5 雙圓模體平面Fig.5 Double circular phantom plane

如圖5所示，兩個相鄰的圓柱形結(jié)構(gòu)平行于圓心線段的間距可表示為：

與圓心線段距離為h的灰度分布函數(shù)可表示為：

其中：ah為與圓心線段距離為h的雙圓間距半寬。需要注意，f(x，ah)表示歸一化的灰度值，且0≤f(x，ah)≤1。

將x=0代入式（11）可得不同間距下中軸線上的CT圖像灰度分布，有：

由于正態(tài)分布函數(shù)g(x)的積分無解析解，引入g(x)積分近似計算公式［22］：

其中W≥1。

2.3 等徑外切圓與周期性結(jié)構(gòu)等效分析

將線對卡CT成像后的調(diào)制度MT引入g(x)積分近似計算公式G(x)，即將式（13）代入式（9），則MT可表示為：

如圖6所示，當(dāng)線對卡的線對波峰（波谷）寬度2ai與雙圓模體的間距2ah相等時，將式（14）代入式（15）可得：

則線對波峰（波谷）寬度為2ah（2ai）對應(yīng)的空間頻率f（單位：線對數(shù)/毫米）可表示為：

圖6 雙圓模體與線對卡寬度相等Fig.6 Width of double-circle phantom equals to the one of line-to-card

其中k為修正系數(shù)，約為1.1～1.2。

式（17）轉(zhuǎn)換可得到：

將式（18）與式（16）聯(lián)立即可獲得利用雙圓模體測得的等效線對卡MTF曲線。

在工程應(yīng)用中，通常采用10%調(diào)制度下的線對寬度來比較工業(yè)CT系統(tǒng)的空間分辨率，因此將MT=10%代入式（16）并進(jìn)行簡化，得到：

對式（19）進(jìn)行多項式求根，取大于1的實數(shù)解，可得W≈2.810 49。

將W≈2.810 49代入式（14），求得：

因此，10%的調(diào)制度下CT系統(tǒng)的空間分辨率為：

2.4 等徑外切圓MTF自動測量算法流程

為了快速獲取雙圓模體在中軸線上的灰度分布，設(shè)計了一種圖像處理及測量算法實現(xiàn)中軸線灰度分布的提取及MTF曲線計算?；诘葟酵馇袌A的工業(yè)CT系統(tǒng)MTF自動測試方法主要包含幾個重要的步驟：CT圖像二值化；模體質(zhì)心計算；雙圓模體圓心計算；圓心中軸線上圖像灰度提取；計算CT系統(tǒng)的MTF曲線。數(shù)據(jù)處理流程圖如圖7所示。

圖7 CT圖像數(shù)據(jù)處理流程Fig.7 Data processing of CT images

假設(shè)雙圓模體所在CT圖像f(x，y)的尺寸為m×n，其中x為圖像數(shù)組橫坐標(biāo)，y為圖像數(shù)組縱坐標(biāo)。如圖8（a）和8（b）所示，利用自動閾值分割方法（clustering）轉(zhuǎn)換為二值化圖像b(x，y)，則b(x，y)的質(zhì)心C位置坐標(biāo)(xc，yc)可表示為：

其中num為滿足條件的i或j的數(shù)量。

為了能夠準(zhǔn)確識別雙圓心位置，需要對兩圓進(jìn)行有效分離。如圖8（c）所示，對b(x，y)進(jìn)行二值形態(tài)學(xué)處理，獲得雙圓邊緣圖像e(x，y)，表示為：

其中：Θ表示腐蝕，即b(x，y)被結(jié)構(gòu)B所腐蝕采用K-means聚類方法計算雙圓圓心坐標(biāo)，其主要步驟如下：

（1）首先，搜索e(x，y)中值為1且與質(zhì)心C距離最遠(yuǎn)的點，設(shè)為N。初始化類A={C}，并將質(zhì)心C的位置作為類A的初始中心位置；初始化類B={N}，并將點B的位置作為類B的初始中心位置；

（2）在e(x，y)中尋找值為1的點，即雙圓的邊緣位置，計算該點與類A和類B的中心點的距離。當(dāng)該點與類A的中心距離小于類B的中心距離時，將該點計入類A，重新計算類A中所有點的中心位置；反之亦然。

圖8 雙圓模體圓心自動提取結(jié)果Fig.8 Automatic extraction of center of double circle phantom

（3）迭代計算步驟（2）直至e(x，y)中所有值為1的點處理完畢，即可獲得圓心坐標(biāo)分別為L=(xL，yL)和R=(xR，yR)。

如圖8（d）所示，利用圓心坐標(biāo)(xL，yL)和(xR，yR)計算中軸線l上的點坐標(biāo)(xl，yl)，可表示為：

計算中軸線l上每個點到圓心連線段LR的距離hl，則有：

其 中：A=xR?xL，B=(yL?yR)/(xL?xR)，

將hl代入式（17）即可獲得相應(yīng)的線對空間頻率。

在雙圓模體CT圖像f(x，y)中，采用數(shù)值微分法提取點坐標(biāo)(xl，yl)對應(yīng)位置的像素灰度值f(xl，yl)。將灰度值f(xl，yl)歸一化并進(jìn)行分段線性擬合處理，代入式（14）、式（16）即可獲得調(diào)制度MT。圖9分別表示中軸線的灰度分布和MTF曲線的計算結(jié)果。

圖9 中軸線灰度分布及MTF曲線Fig.9 Gray scale distribution of central axis and MTF curve

3 實 驗

3.1 試塊的制作

選用牌號為GCr15的軸承鋼作為雙圓模體材料，該材質(zhì)具有淬透性好、硬度高、耐磨性強(qiáng)和尺寸穩(wěn)定等優(yōu)點。首先在原料條上進(jìn)行斷切，得到兩個外形尺寸均為Ф9 mm×16 mm的坯料；進(jìn)行研磨加工及表面拋光，得到兩個Ф8 mm×15 mm的圓柱模體，接著用SGM15/14型高溫爐進(jìn)行淬火及低溫回火，使其硬度達(dá)到45～50；對其表面進(jìn)行防銹處理，最后采用三坐標(biāo)儀對模體尺寸進(jìn)行計量，驗證其尺寸公差是否滿足工業(yè)CT空間分辨率測試的要求。采用垂吊方式保證模體圓柱側(cè)面緊密貼合，如圖10所示。

圖10 雙圓模體示意圖Fig.10 Double circular phantom

實驗前，將雙圓模體置于20℃條件下恒溫2 h；實驗時，將模體置于CT設(shè)備轉(zhuǎn)臺上，選取模體的中部區(qū)域進(jìn)行工業(yè)CT掃描成像。實驗所用的設(shè)備為高能加速器工業(yè)CT成像系統(tǒng)，由6 MeV加速器射線源、608通道線陣探測器、雙立柱運(yùn)動平臺和Ф1 m轉(zhuǎn)臺組成。設(shè)備探測器單元水平方向尺寸為0.3 mm（固定），相鄰?fù)ǖ浪介g隔為1.3 mm，探測器單元垂直方向尺寸（切片厚度）可利用準(zhǔn)直器進(jìn)行調(diào)節(jié)，范圍為0.25～5 mm，最大劑量率約為800 cGy/min×1 m，焦點尺寸為2 mm，采用三代CT掃描方式。

3.2 有效性分析

為驗證本文方法的有效性，將它與線對卡法［16-17］進(jìn)行比較，研究本文方法的測試結(jié)果受空間分辨率、噪聲水平、穩(wěn)定性和雙圓模體直徑變化的影響，以及與線對卡法測試結(jié)果的差異。如圖11所示，線對卡寬度分別為0.1/0.15/0.2/0.25/0.3/0.4/0.5 mm。

圖11 線對卡模體及CT圖像Fig.11 Wire pair card phantom and CT image

微動次數(shù)是指探測器單元通過小范圍機(jī)械運(yùn)動形成過采集插值實現(xiàn)超分辨率重建，是影響高能工業(yè)CT系統(tǒng)空間分辨率的主要因素。微動裝置采用絲桿螺母副傳動、交叉滾子滑動模塊導(dǎo)向?qū)崿F(xiàn)高精度小范圍平移，將探測器相鄰單元間距按微動次數(shù)進(jìn)行等分，利用光柵尺測量探測器單元的實際位置用于軟件算法的插值計算。因此，通過改變工業(yè)CT系統(tǒng)檢測時的微動次數(shù)可形成不同空間分辨率的CT圖像。圖12（a）～12（c）分別給出了2次、5次和10次微動獲得的線對卡和雙圓模體CT圖像，實驗采用的切片厚度為1 mm、圖像分辨率為4 096×4 096（像元×像元）、成像范圍為200 mm×200 mm。采用兩種方法測得10%調(diào)制度下的空間分辨率如圖13所示，可見兩種方法獲得的10%MTF的空間分辨率基本吻合，最大誤差在10%以內(nèi)。

研究不同CT圖像分辨率對MTF測試結(jié)果的影響，實驗采用2次微動、切片厚度為1 mm、成像范圍為200 mm×200 mm，如圖14所示。增加圖像分辨率有助于提高檢測系統(tǒng)的空間分辨率，但是當(dāng)空間分辨率接近系統(tǒng)極限時，增加圖像分辨率對系統(tǒng)空間分辨率的提高效果不明顯且會大量消耗硬件內(nèi)存。從測量結(jié)果可見，兩種方法均體現(xiàn)出圖像分辨率對MTF測試結(jié)果的影響，在高圖像分辨率下（8 192×8 192），本文方法測量值（10%調(diào)制度）略高于線對卡測量值，原因在于線對卡法不同空間分辨率線對的間隔較大，10%調(diào)制度下的實際空間分辨率需要通過插值擬合計算獲得，而雙圓法的等效線寬間隔在10%調(diào)制度下遠(yuǎn)小于線對卡法。因此，雙圓法的測量精度高于線對卡法，同時在日常檢測時應(yīng)選擇合適的圖像分辨率以獲得最佳的空間分辨率及最優(yōu)的系統(tǒng)開銷。

圖12 不同微動條件下CT圖像Fig.12 CT images with different fretting frequencies

圖13 不同微動頻率下測試結(jié)果比較Fig.13 Comparison of test results with different fretting frequencies

圖14 不同圖像分辨率下測試結(jié)果比較Fig.14 Comparison of test results with different image resolutions

圖15 不同切片厚度下測試結(jié)果比較Fig.15 Comparison of test results with different slice thicknesses

研究不同噪聲水平對MTF測試結(jié)果的影響，實驗采用2次微動、圖像分辨率為4 096×4 096（像元×像元）、成像范圍為200 mm×200 mm，如圖15所示。切片厚度是影響CT圖像噪聲水平的主要因素，前期實驗表明切片厚度越大，圖像噪聲水平越低。可見，噪聲對本文方法具有一定的影響，噪聲越大與線對卡法的結(jié)果差異越大，原因在于本文方法通過對中軸線上灰度分布進(jìn)行分段線性擬合來降低圖像噪聲的干擾，降噪效果有限。通過分析，兩者誤差控制在10%以內(nèi)，基本滿足檢測要求。

圖16 不同直徑模體的測試結(jié)果比較Fig.16 Comparison of test results of phantom with differ?ent diameters

研究不同直徑模體對MTF測試結(jié)果的影響，實驗采用5次微動、圖像分辨率為4 096×4 096（像元×像元）、成像范圍為200 mm×200 mm、切片厚度為1 mm，雙圓模體直徑分別為3，6，8，20，26，32 mm，如圖16所示。可見，圓柱體直徑的變化不會顯著影響空間分辨率測量結(jié)果，噪聲、偽影是引起不同直徑模體測量結(jié)果波動的主要原因。

4 結(jié) 論

本文針對傳統(tǒng)線對卡法測量工業(yè)CT系統(tǒng)空間分辨率的局限性，推導(dǎo)了線對卡法測量模型及等徑外切圓變間距模型，建立了兩者的對應(yīng)等價關(guān)系，構(gòu)建了基于等徑外切圓的MTF測量方法，分析了測量結(jié)果受空間分辨率、圖像分辨率、噪聲水平及模體直徑的影響規(guī)律，并采用GCr15軸承鋼模體實驗驗證了本方法的有效性。實驗結(jié)果表明，本方法的測量結(jié)果接近線對卡法，誤差控制在10%以內(nèi)，具有模體加工要求低、測量過程簡單的優(yōu)點，可推廣用于CT系統(tǒng)的設(shè)計、研制、驗收和應(yīng)用等方面。下一步將針對等徑外切圓法的降噪方法及偽影對測量結(jié)果的影響開展研究，使本方法的適用性更廣。