基于等徑外切圓的工業CT系統空間分辨率測試

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

（1.浙江工業大學機械工程學院，浙江杭州310014；2.中國兵器科學研究院寧波分院，浙江寧波315103）

1 引 言

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

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

2 空間分辨率測試模型

2.1 周期性結構測試模型

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

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

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

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

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

圖3 條形模體表現形式Fig.3 CT manifestation of strip phantom

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

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

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

分析各個條形結構之間的影響關系，簡化波峰Imax（波谷Imin）的表達式。假設工業CT成像過程的高斯退化函數模型為：

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

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

2.2 等徑外切圓結構模型分析

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

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

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

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

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

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

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

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

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

其中W≥1。

2.3 等徑外切圓與周期性結構等效分析

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

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

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

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

其中k為修正系數，約為1.1～1.2。

式（17）轉換可得到：

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

在工程應用中，通常采用10%調制度下的線對寬度來比較工業CT系統的空間分辨率，因此將MT=10%代入式（16）并進行簡化，得到：

對式（19）進行多項式求根，取大于1的實數解，可得W≈2.810 49。

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

因此，10%的調制度下CT系統的空間分辨率為：

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

為了快速獲取雙圓模體在中軸線上的灰度分布，設計了一種圖像處理及測量算法實現中軸線灰度分布的提取及MTF曲線計算。基于等徑外切圓的工業CT系統MTF自動測試方法主要包含幾個重要的步驟：CT圖像二值化；模體質心計算；雙圓模體圓心計算；圓心中軸線上圖像灰度提取；計算CT系統的MTF曲線。數據處理流程圖如圖7所示。

圖7 CT圖像數據處理流程Fig.7 Data processing of CT images

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

其中num為滿足條件的i或j的數量。

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

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

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

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

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

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

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

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

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

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

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

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

3 實 驗

3.1 試塊的制作

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

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

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

3.2 有效性分析

為驗證本文方法的有效性，將它與線對卡法［16-17］進行比較，研究本文方法的測試結果受空間分辨率、噪聲水平、穩定性和雙圓模體直徑變化的影響，以及與線對卡法測試結果的差異。如圖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

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

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

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

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

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

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

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

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

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

4 結 論

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