基于改進卷積反投影模型的CT成像偽影消除方法*

韓寅奔， 董泉潤

(西北工業大學 a. 教務處， b. 動力與能源學院， 西安 710072)

隨著現代工業的發展，計算機斷層成像(computed tomography，CT)技術被越來越多地運用到無損檢測和逆向工程等行業中[1].CT重建是一種無損檢測技術，其利用精確準直的X射線在不破壞被測物體的前提下，根據被測物體不同組織結構吸收射線能力的差異性對被測物體進行斷層成像，能省時、無損地探測物體內部情況，獲取內部結構信息[2-4].

目前，常用的CT重建算法主要有迭代法和解析法兩種[5].迭代法能處理不同方式的采樣數據，且能獲得稀疏采樣數據的重建圖像，但該種方法需要較大的計算量，重建速度較慢[6-7]；解析法中最具代表性的是卷積反投影法，其通過處理完全等間距的投影數據來獲得更好的成像質量和更快的成像速度，但這種方法得到的重建圖像容易產生偽影[8-9].

相比于迭代法，卷積反投影法的應用范圍更廣.許多專家和學者提出了諸多方法來提高卷積反投影法的重建速度與成像質量，消除偽影.在提升重建速度方面，文獻[10]基于極坐標的特性，提出了一種PCBP重建算法，但其在坐標轉換過程中需要較大的計算量；文獻[11]基于正余弦函數的特性，同時重建16幅采樣圖像，明顯提升重建速度；文獻[12]則基于三角函數的對稱性，使用對稱優化算法同時重建64幅圖像.而在提升成像質量方面，文獻[13]通過分析不同偽影的種類和成因，分別消除各種醫學圖像的偽影；文獻[14]則基于周期性原理，使用正切函數來優化重建圖像，但這種方法得到的重建圖像仍存在部分偽影.

針對上述問題，本文提出了一種基于卷積反投影算法的CT成像偽影消除方法.首先分析并標定給定模板的測量參數；然后使用Radon和傅里葉變換得到反投影方程，并計算每一點吸收率與投影值的積分關系，從而得到較清晰的重建圖像.

1 參數標定

使用CT系統進行無損檢測時，保持射線源和探測器不動.CT系統探測模板與著色圖及其坐標系示意圖如圖1所示.為采集不同旋轉角度的投影信息，探測物包含橢圓和圓形兩個模板.然而，理論上旋轉中心與射線源焦點的連線是垂直于探測器的，但是在安裝CT系統時，難以實現這一要求，通常存在誤差，并導致重建圖像存在偽影.因此，需要對安裝好的CT系統進行參數標定，即借助模板標定CT系統的參數，并使用這些參數對未知結構的樣品進行成像.

本文首先繪制模板投影信號著色圖，得到模板旋轉軌跡.圖1b中用黃色表示投影信號強度為零，并用深淺不同的藍色表示不同強度的投影信號，強度越高，顏色越深；然后以著色圖左下角為原點，橫軸表示旋轉角度，縱軸表示其探測器端點的距離，建立直角坐標系；最后結合著色圖的幾何特征，建立其與實際照射情況的對應關系，并在此基礎上標定CT系統的旋轉中心、探測器單元間距離及探測方向等參數.

圖1 探測模板與著色圖及其坐標系Fig.1 Detection template，coloring diagram and coordinate system

1.1 探測器單元間距離

當X射線平行于橢圓長軸入射時，光在介質中走過的路徑最長，衰減幅度最大，接收到的信號最強.此時，橢圓模板在接收器上的投影長度最短且恰為短軸長，即圖1b中A、B兩點的縱坐標之差.測得AB之間共有L個單元格，即對應L-1個探測單元間距，故可得探測器單元間距離d=m/(L-1)，其中，m為橢圓模板的短軸長.

1.2 系統照射方向

模板每旋轉360°后將恢復為初始狀態，著色圖上的曲線則完成一個首尾相接的周期.由模板形狀可知，其旋轉0°～180°所得數據與旋轉180°～360°所得數據相同，因此，只需要將探測物轉過180°即可獲得實驗所需的所有數據.

設X射線與旋轉中心和模板上小圓圓心連線之間所夾銳角為φ，則小圓圓心在探測器上投影點距起始點的距離l可表示為

l=l0+d0sinφ

(1)

式中：l0為旋轉中心投影點到探測器起始點的距離；d0sinφ為旋轉中心與圓心連線在探測器上的投影長度.

由式(1)可知，模板與探測器始端的距離可以表示為三角函數形式.而式(1)則為著色圖中等寬曲線的數學描述，故可使用傅里葉級數擬合該曲線并求得其跨度q.由CT系統的工作原理可知，系統每次轉過的角度相同，故探測器每次轉過180°/q.

1.3 旋轉中心標定

如圖1b所示，C點對應等寬曲線的最高點，設θ1、θ2分別為A、C兩點在坐標圖上對應的橫坐標，A、C兩點的橫軸間距θm為

θm=θ1-θ2

(2)

在實際情況下，θm與光線轉過夾角相等，對應狀態幾何關系如圖2所示.

圖2 A、C兩點對應狀態的幾何關系Fig.2 Geometric relation of corresponding state between A and C points

由于探測器每次轉過180°/q，故通過觀察A、C兩點之間對應列的列數I，即可求得

θm=I180°/q

(3)

當等寬曲線位于C點時，模板圓心與旋轉中心的間距被完全投影到探測器上，表明此時圓心投影點距探測器始端最遠；而當位于其對稱位置時，距探測器始端最近.故只需計算出圖1b中等寬曲線兩極值點豎直方向上相差的單元格數目T，即可求出圓心與旋轉中心之間的距離D，即

l(T-1)=2D

(4)

若以模板在圖2的幾何狀態為基準，以左下角為原點建立直角坐標系，則旋轉中心的坐標值表示為

(5)

2 模型CT成像偽影消除

CT圖像斷層平面中某一點的密度值可以看作平面內所有經過該點的射線投影值之和的均值.然而觀測樣品各點的形貌與吸收率各不相同，有些點比較突出，經卷積反投影處理后，這些點依舊突出，并可據此還原斷層；但有些原本不為零的地方現在成為了零，即產生偽跡，因此，需要在反投影之前對投影進行濾波處理.

2.1 Radon變換及傅里葉變換

Pφ(ρ)=P(ρ，φ)

(6)

式中：w1、w2為傅里葉變換[15]的特征參量；F為傅里葉變換映射；P為投影函數p的傅里葉變換.可由傅里葉反變換求得待重建圖像，即

a(x，y)=F-1[E(w1，w2)]=

(7)

式中：α為實際某點投影角度；r為弧長.對于式(7)進行積分，可以整理為

g[rcos(α-φ)，φ]

(8)

故待建圖像轉化為

(9)

由前述推導過程可得

(10)

反投影得到托盤上各點的吸收率為

(11)

考慮到實際分析處理問題時需要將介質劃分為256×256的微元格點研究，故需對上述連續模型進行離散化處理，則有

(12)

濾波因子采樣得到

(13)

其取值為

(14)

可將離散條件下的反投影積分轉化為級數求和，即

(15)

式中，N為一個微元滿足的最小旋轉角的最大倍數.將式(15)離散化則有

(16)

式中，Δx、Δy為橫縱坐標最小變化量.

2.2 偽影消除

使用2.1節變換操作可確定絕大部分點的介質吸收率，然而，實際操作中發現有少部分點并沒有X射線穿過而無法確定其吸收率[16].若使用近鄰插值方法求取這些點的吸收率，雖仍能保持樣品的形貌特征，但重建圖像部分邊緣較模糊，并形成偽跡.因此為抑制偽跡出現，一方面要對源數據進行降噪處理；另一方面需要對變換后的數據進行插值賦值，減小數據的偏離幅度.

本文使用3階樣條法對式(16)得到的重建圖像進行插值處理.3階樣條法除了要求在兩端與給定數據值相等外，還要求樣條系統的光滑性，即內部鄰接數據點的1階及2階導數必須匹配.在外側端點1階導數沒有變化，2階導數為零.樣條函數定義為

(17)

其導數為

(18)

為滿足數據匹配條件，則有

(19)

內部數據點1階、2階導數匹配，故在兩函數銜接點gn(i2)處有

(20)

在外側端點1階導數沒有變化，故2階導數為0，即

(21)

聯立式(18)～(21)可解得具體的插值函數，并可將[x1，x3]區間內任一點帶入相應坐標求得其值.

3 算例分析

本文使用圖1a所示的模板標定CT系統的參數，使用文中提出的方法求得探測器單元間距離l為0.280 4 mm，旋轉中心的坐標為(58.633 1 mm，42.395 3 mm)，著色圖中等寬曲線在整個數據點區間內跨度大約為180°，探測器每次轉過1°，傅里葉變換特征參量w1、w2分別取30和60，濾波因子由式(14)求得.

使用標定后的系統參數進行測試實驗，得到如圖3所示的反投影結果.從圖3可以看出，得到的重建圖較為模糊，且輪廓周圍泛有過渡白色.使用3階樣條法進行插值處理后得到的結果如圖4所示，處理后的圖像明顯清晰，內部分布有兩個彼此分離的橢圓空洞、兩個相交的更高吸收率的橢圓形小區域.

圖3 未經去噪、插值處理的圖像Fig.3 Images without denoising and interpolating processes

圖4 降噪、插值后樣品最終圖像Fig.4 Final image of sample after noise reduction and interpolation

本文也比較了消除偽影前后10個給定位置處的介質吸收率，如表1、2所示.根據表1結果可以看出，部分吸收率為負值，但其值均較小，表明實際吸收率應為0，存在一定的測量誤差.比較表1、2的吸收率可知，消除偽影后能測得每個位置的吸收率，并能得到更清晰的重建圖像.

表1 偽影消除前給定10個位置處吸收率Tab.1 Absorption rates at given 10 positions before artifact removal

表2 偽影消除后給定10個位置處吸收率Tab.2 Absorption rates at given 10 positions after artifact removal

3.1 算法精度分析

實驗中的誤差將降低測試精度，這些誤差主要包括截斷誤差、公式化誤差、舍入誤差和測量誤差.在參數標定時，本文討論不涉及多項式近似，且忽略了次要條件的影響，并不會大幅度影響最終的結果.因此，本文算法的主要誤差來源為給定實驗數據、計算軟件的舍入誤差以及數據收集過程中不精確性引入的測量誤差.下面以圓形遮光區域為例，推導吸光區尺寸與誤差的定量關系.

設k為圓遮擋的X射線數，設z為真實值與測量值之間的偏差，則有

kd+z=2L

(22)

式中，2L為測量值，故相對誤差為

(23)

由于一般情況下誤差會控制在較小范圍內，即kd≈2L，則有

(24)

由式(24)可知，隨著物體半徑L的增大，誤差η在減小，表明增大模板吸光區半徑，可有效減小誤差.如將半徑擴大為原來的二倍，則誤差為

可見半徑擴大為原來的二倍，誤差幾乎變為原來的1/2，因此，對于小圓半徑相對于整個平板較小的情況，參數的測量精度還有較大的提升空間.但由于偏差z值較小，即使誤差減半，也不會對最終結果影響過大，故本文測得的參數精度是可以接受的.

3.2 算法穩定性分析

本文假設X射線與橢圓中心、小圓圓心的連線恰好垂直，如圖2所示.但實際上取樣是非連續的，只能近似為垂直關系，由于探測器間距較近，這種近似導致的誤差雖會隨具體情況變化，但其值較小，在一定程度上是可接受的.

本文測試了以旋轉中心坐標為例的參數值對系統參數的敏感性.將式(5)中的旋轉中心橫縱坐標對夾角θm求導，得到

(25)

4 結 論

本文提出了一種基于卷積反投影算法的CT成像偽影消除方法.利用實物與投影的對應關系，根據投影方便、快捷地求出探測器間距等待系統參數；使用卷積反投影技術和三次樣條插值算法，精確地求出每一點的吸收率并得到較清晰的重建圖像.算例仿真測試的精確度、穩定性分析結果表明，提出的算法具有簡明易懂、操作方便和精確度高等諸多優點，且其相較于傳統方法具有更高的運算效率.