馬 敏 趙 亮 姬晶晶

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

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

1 重建圖像偽影

1.1 偽影概述[4]

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

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

1.2 偽影原理

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

以一線性系統(tǒng)為例來(lái)闡述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ī)輸入信號(hào)，且

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

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

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

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

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

按最小均方誤差準(zhǔn)則確定 Wiener濾波器的沖激響應(yīng)h(n)，令ξ(n)對(duì)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)序列的長(zhǎ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濾波器以平滑圖像邊緣，它對(duì)于圖像變化越明顯的區(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。對(duì)三種氣液兩相流流型進(jìn)行仿真。程序在 Intel(R) Core(TM)2 Duo CPU、1 G內(nèi)存的PC上運(yùn)行。操作系統(tǒng)為Windows XP系統(tǒng)，程序運(yùn)行平臺(tái)為MAT LAB7.0。由投影數(shù)據(jù)分別對(duì)泡狀流、核心流和層流采用不同圖像重建算法進(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，探測(cè)器陣列由17個(gè)CdZnTe探測(cè)器組成，單個(gè)探測(cè)器的尺寸為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é)果，用空泡份額及相對(duì)誤差對(duì)重建圖像的質(zhì)量進(jìn)行評(píng)估?？张莘蓊~是指氣液兩相流中，氣相所占截面積與總流通截面積之比。以流體管道中含有三個(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表示，空泡份額及相對(duì)誤差比較結(jié)果見表1。

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

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

5 結(jié)語(yǔ)

本文提出了一種應(yīng)用于CT成像系統(tǒng)圖像重建算法的改進(jìn)方案，針對(duì)濾波反投影算法產(chǎn)生偽影的問(wèn)題，通過(guò)采取將Wiener濾波技術(shù)融入傳統(tǒng)算法的方法提高重建圖像的質(zhì)量。此方法對(duì)有泡狀流、核心流和層流等典型流型的模擬管道進(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]. 中國(guó)醫(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