馬 敏 趙 亮 姬晶晶

（中國民航大學(xué)航空自動(dòng)化學(xué)院 天津 300300）

計(jì)算機(jī)層析成像(CT)是通過對物體進(jìn)行不同角度的射線投影測量而獲取物體橫截面信息的成像技術(shù)。CT技術(shù)的核心是投影重建圖像，其實(shí)質(zhì)是由掃描所得的投影數(shù)據(jù)反求出成像平面上每個(gè)點(diǎn)的衰減系數(shù)值。圖像重建方法大致分為[1]：變換法和級數(shù)展開法。濾波反投影算法是目前應(yīng)用最廣泛的基于變換法的圖像重建算法[2,3]，具有重建速度快、空間和密度分辨率高等優(yōu)點(diǎn)，缺點(diǎn)是對投影數(shù)據(jù)的完備性要求較高，否則會產(chǎn)生偽影。本文分析了濾波反投影重建算法產(chǎn)生偽影的原因，提出一種新的圖像重建方法，將 傳 統(tǒng)濾波反投影算法與Wiener濾波相結(jié)合。實(shí)驗(yàn)結(jié)果表明，該方法較好地克服了圖像偽影現(xiàn)象，圖像質(zhì)量更好。

1 重建圖像偽影

1.1 偽影概述[4]

CT圖像偽影又稱偽像，是指圖像重建過程中，不同類型的圖像互相干擾和其他各種非隨機(jī)干擾在圖像上的表現(xiàn)，使部分影像根本不反映被測物體的對應(yīng)區(qū)域。CT投影是整個(gè)CT過程的輸入，投影的好壞直接影響重建圖像的質(zhì)量。若投影數(shù)據(jù)先天不足，則無論何種算法都不能得到完美的重建圖像，得到的只是帶有偽影的圖像，且投影數(shù)據(jù)截?cái)嘣酱?，偽影越多?/p>

造成投影質(zhì)量差的原因主要為[5]幾何誤差、探測器缺陷、射線源不理想以及各種噪聲。幾何誤差因掃描方式而異，對二維CT圖像中心旋轉(zhuǎn)中心不重合，會造成“拖尾”偽影；探測器缺陷和探測器陣列元素響應(yīng)不一致，會造成環(huán)形偽影；射線源的多色性，會導(dǎo)致射線穿透試件時(shí)的硬化效應(yīng)，從而造成杯狀偽影；射線源能量歧離也會產(chǎn)生條紋狀偽影；噪聲主要包括隨機(jī)噪聲、散射噪聲等。在射線穿過物體時(shí)，由于像素對射線采樣的不足，通道信號會產(chǎn)生混疊失真，出現(xiàn)偽影。而濾波反投影重建的本質(zhì)是把取自有限物體空間的射線投影均勻地反投影到射線所及的無限空間的各點(diǎn)之上，因此產(chǎn)生的是星狀偽影。

1.2 偽影原理

本文以反投影重建算法為例來闡述偽影產(chǎn)生機(jī)理[6,7]，該算法可等效為以原圖像為輸入，重建后圖像為輸出的成像系統(tǒng)(圖1)。

圖1 反投影重建的等效成像系統(tǒng)Fig.1 Equivalent imaging system of back projection reconstruction.

該系統(tǒng)的點(diǎn)擴(kuò)展函數(shù)(PSF)可如下求得。設(shè)位于坐標(biāo)原點(diǎn)的點(diǎn)源δ(x,y)為x-y斷面中唯一像點(diǎn)，并設(shè)掃描方式為平移/旋轉(zhuǎn)。待重建圖像為f(x,y)，射線方向由旋轉(zhuǎn)坐標(biāo)軸xr與固定坐標(biāo)軸x的夾角確定，射線透射物體后的投影為pφ(xr)。

掃描坐標(biāo)系統(tǒng)見圖2，虛線代表射線投影方向，其方向始終與yr平行，當(dāng)φ角從0°旋轉(zhuǎn)到90°時(shí)，得到一系列投影pφi(xr)xr=rcos(θ–φ)。

圖2 平移/旋轉(zhuǎn)掃描方式的坐標(biāo)系統(tǒng)Fig.2 Coordinate system of tran s lation/rotation scanning mode.

設(shè)輸入圖像為點(diǎn)源，利用 反 投 影重建算法，重建圖像即為該系統(tǒng)的點(diǎn)擴(kuò)展函數(shù)[5]。

式中，φ0為α(φ)=0的唯一解，令α(φ)=rcos(θ–φ)，則α(φ0)=rcos(θ–φ0)=0。

因r≠0，故sin(θ–φ0)=1。所以

可見，相應(yīng)于反投影重建算法的系統(tǒng)，它的點(diǎn)擴(kuò)展函數(shù)不是δ-函數(shù)，雖在r=0處能反映原圖是點(diǎn)源的情況，但在r≠0處，像素值不為 0 ，這是此種算法存在偽影的本質(zhì)。

2 Wiener濾波概述

Wiener濾波屬于反卷積算法，是從噪聲中提取信號的一種濾波方法，由Wiener[8]提出，在一、二維信號處理和圖像復(fù)原領(lǐng)域得到廣泛應(yīng)用。在濾波反投影算法(Filtered Back Projection, FBP)對圖像進(jìn)行重建時(shí)，濾波函數(shù)的選擇對消除反投影所產(chǎn)生的星形偽影及圖像重建質(zhì)量影響很大，目前常用R-L濾波函數(shù)和S-L濾波函數(shù)。R-L濾波函數(shù)優(yōu)點(diǎn)是形式簡單、實(shí)用，重建圖像輪廓清晰，缺點(diǎn)為重建圖像有明顯的振蕩現(xiàn)象，且當(dāng)投影數(shù)據(jù)含有噪音時(shí)，重建質(zhì)量比較差。S-L濾波函數(shù)優(yōu)點(diǎn)是對投影數(shù)據(jù)的高頻成分有抑制作用，重建圖像的震蕩響應(yīng)小，缺點(diǎn)是圖像分辨率不高。相比而言，Wiener濾波法頻域形式簡單，計(jì)算效率高，復(fù)原效果良好，且抗噪性能優(yōu)良，在重建圖像的精度和速度方面都有明顯的優(yōu)勢。

以一線性系統(tǒng)為例來闡述Wiener濾波原理[9]，若它的單位樣本響應(yīng)為h(n)，Wiener濾波器的輸入-輸出關(guān)系如圖3所示。

圖3 Wiener濾波原理示意圖Fig.3 Schematics of the principle of multistage Wiener filter.

由圖3，x(n)為隨機(jī)輸入信號，且

其中，s(n)表示信號，v(n)表示噪聲，則輸出y(n)為

y(n)為s(n)的估計(jì)值，用 ?(n)表示，即

用e(n)表示真實(shí)值與估計(jì)值之間的誤差，即

顯然，e(n)是一隨機(jī)變量，可能為正或負(fù)。因此，選擇用它的均方誤差來表達(dá)誤差，所謂均方誤差最小即它的平方的統(tǒng)計(jì)期望最?。?/p>

3 Wiener?Hopf方程的求解[10]

按最小均方誤差準(zhǔn)則確定 Wiener濾波器的沖激響應(yīng)h(n)，令ξ(n)對h(j)的導(dǎo)數(shù)等于零，即得式中，Rxs(m)=E(x(n)s(n+m))是s(n)與x(n)的互相關(guān)函數(shù)，Rxx(m) =E(x(n)x(n+m))是x(n)的自相關(guān)函數(shù)。式(8)所列為 W iener濾波器的標(biāo)準(zhǔn)方程或Wiener-Hopf方程。若已知Rxs(m)和Rxx(m)，則解此方程即可求出Wiener濾波器的沖激響應(yīng)。

設(shè)濾波器沖激響應(yīng)序列的長度為N[11]，沖激響應(yīng)矢量為

濾波器輸入數(shù)據(jù)矢量為

則濾波器的輸出為

Wiener-Hopf 方程也可寫成

其中

式(13)為s(n)與x(n)的互相關(guān)函數(shù)，是N維列矢量；R=E?x(n)xT(n)? 是x(n)的自相關(guān)函數(shù)，是N階方陣。利用求逆矩陣的方法直接求解式( 1 2 )，得

式中，opt表示“最佳”，即FIR Wiener濾波器的沖激響應(yīng)。

利用Wiener濾波器以平滑圖像邊緣，它對于圖像變化越明顯的區(qū)域所起的作用越強(qiáng)。實(shí)驗(yàn)結(jié)果顯示，在程序中加入Wiener濾波后得到圖像更清晰。

4 仿真及實(shí)驗(yàn)結(jié)果

4.1 數(shù)值仿真

為驗(yàn)證該法的有效性，進(jìn)行了數(shù)值仿真實(shí)驗(yàn)。本文重建圖像分辨率為 560×420。對三種氣液兩相流流型進(jìn)行仿真。程序在 Intel(R) Core(TM)2 Duo CPU、1 G內(nèi)存的PC上運(yùn)行。操作系統(tǒng)為Windows XP系統(tǒng)，程序運(yùn)行平臺為MAT LAB7.0。由投影數(shù)據(jù)分別對泡狀流、核心流和層流采用不同圖像重建算法進(jìn)行圖像重建，其中濾波反投影重建算法采用的是R-L濾波函數(shù)，兩種方法迭代次數(shù)均為20。仿真重建圖像如圖4所示。

圖4 FBP與改進(jìn)的FBP算法重建的圖像Fig.4 Images reconstructed using the FBP and improved FBP algorithm.

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

所用多相流CT實(shí)驗(yàn)裝置的射線源為1.11×1010Bq(300 mCi)241Am源，γ射線能量為59.5 keV，探測器陣列由17個(gè)CdZnTe探測器組成，單個(gè)探測器的尺寸為3 mm×7 mm×3 mm。將直徑為29、22.4、13.4 mm的有機(jī)玻璃棒插入充滿水的流體管道中，模擬三個(gè)泡狀體的泡狀流；將Φ60 mm有機(jī)玻璃管套在Φ82 mm流體管內(nèi)，將水充入內(nèi)管道模擬核心流；用Φ82 mm有機(jī)玻璃管模擬流體管道，在該管道中充入一半水模擬層流。從7個(gè)角度采集實(shí)驗(yàn)數(shù)據(jù)，在 M atlab7.0環(huán)境下分別用濾波反投影重建算法和加Wiener濾波的反投影重建算法作圖像重建。

為比較改進(jìn)前后算法的仿真結(jié)果，用空泡份額及相對誤差對重建圖像的質(zhì)量進(jìn)行評估?？张莘蓊~是指氣液兩相流中，氣相所占截面積與總流通截面積之比。以流體管道中含有三個(gè)氣泡為例計(jì)算空泡份額，并與實(shí)驗(yàn)?zāi)Ｐ椭械膶?shí)際比例進(jìn)行比較。根據(jù)已知實(shí)驗(yàn)?zāi)M的流體管道管徑和三個(gè)模擬泡狀氣體管徑數(shù)據(jù)，可算出三個(gè)模擬泡狀體管體的截面占流體管道截面的比例為Pa=0.1251、Pb=0.0746、Pc=0.0267。用圖像處理技術(shù)算出的各模擬泡狀管截面所成像所占流體管道截面成像的比例用P'a、P'b、P'c表示，空泡份額及相對誤差比較結(jié)果見表1。

濾波反投影重建算法與基于 Wiener濾波的反投影重建算法均采用20次迭代，對泡狀流、核心流和層流進(jìn)行圖像重建，重建時(shí)間見表2。

表1 空泡份額及相對誤差比較Table 1 Comp a r ison of void fra c tio n and relative e rrors.

表2 FBP與改uter time of the FBP算法的重建 oved FBP and impr algo進(jìn)的 時(shí)間Table 2 Comp FBP ns).(20次迭代)rithm (20 iteratio

從實(shí)驗(yàn)結(jié)果對比可見，本文所用改進(jìn)的重建算法的重建效果明顯好于傳統(tǒng)的濾波反投影的重建效果。在近似相等的時(shí)間內(nèi)，濾波反投影重建圖像雖然也可以重建圖像，但表面有許多細(xì)條紋和亮紋，有明顯的 偽 影，而基于維納濾波的反投影重建算法重建的圖像則清晰得多，表層均勻，無雜質(zhì)。由表2和表3，加入維納濾波后，重建圖像明顯偽影減少，更接近原始圖像，因此有一定的可行性和應(yīng)用性。

5 結(jié)語

本文提出了一種應(yīng)用于CT成像系統(tǒng)圖像重建算法的改進(jìn)方案，針對濾波反投影算法產(chǎn)生偽影的問題，通過采取將Wiener濾波技術(shù)融入傳統(tǒng)算法的方法提高重建圖像的質(zhì)量。此方法對有泡狀流、核心流和層流等典型流型的模擬管道進(jìn)行圖像重建，實(shí)驗(yàn)結(jié)果表明，方法較好地克服了偽影，成像的邊界更加清晰，獲得較好的圖像重建效果。

1 Gilbert B K , Welch B, Holmes D R. Rapid execution of fan beam image reconstruction algorithms using efficient computational techniques and special-purpose[J]. IEEE Trans Biomedical Engineering, 1981, 28(2): 98–115

2 劉杰, 施寅, 阮秋琦. CT快速圖像重建算法研究[J]. 中國醫(yī)學(xué)物理學(xué)雜志, 2003, 3: 15–17 LIU Jie, SHI Yin, RUAN Qiuqi. The study of fast image reconstruction[J]. Chinese Journal of medical Physics,2003, 3: 15–17

3 Lucas G P, Albusaidi K. An axial scanning system for investigating vertical gas-liquid flows in the slug and bubbly-slug transition regimes[J]. Flow Measurement and Instrumentation, 1998, 9: 9–35

4 Luenberger D G. Linear and nonlinear programming[J].MA: Addison-Lesley, 1989, 5: 19–25

5 莊天戈. CT原理與算法[M]. 上海: 上海交通大學(xué)出版社, 1992: 34–35,77 ZHUANG Tiange. Principles andAlgorithms of CT[M].Shanghai: Shanghai Jiaotong UniversityPress, 1992:34–35,77

6 Fossa M. Design and performance of a conductance probe for measuring the liquid fractionin two-phase gas-liquid flows[J]. Flow Measurement and Instrumentation, 1998, 9:103–109

7 Thorn R, Hamme E A, Johanen G A. Recent developments inhere-phase flow measurement[J].Measurement Science Technology, 1997, 8(7): 691–701

8 Klaus Muelle. Anti-aliased three-dimentional cone-beam reconstruction of low-constrast objects with algebraic methods[J]. IEEE Transactiol Imagining, 1999, 5(2):116–134

9Alexander D P, Zayed M R. Adaptive filtering primer with MATLAB[J]. Chinese Journal of medical Physics,2006, 10: 50–75

10 Johansen G A, Froystein T, Hjertaker B T. A dual sensor flow imaging topographic system[J]. Measurement Sci Tech, 1996, 7(2): 297–307

11 Jiang Hsieh. Computer tomography principle, design,artifacts and recent advances[J]. SPIE of America, 2003,5(10): 10–12