王化祥，薛 倩

(天津大學(xué)電氣與自動(dòng)化工程學(xué)院，天津 300072)

工業(yè)過(guò)程中存在大量多相流測(cè)量問(wèn)題，高精度、實(shí)時(shí)在線測(cè)量對(duì)生產(chǎn)過(guò)程的控制與管理有重要意義．圖像重建算法是實(shí)現(xiàn)層析成像(computed tomography，CT)系統(tǒng)可視化測(cè)量的關(guān)鍵．隨著 CT技術(shù)的廣泛應(yīng)用，一些新算法如小波變換、神經(jīng)網(wǎng)絡(luò)模型、遺傳算法等被引入到圖像重建中，取得了較好效果．然而對(duì)于工業(yè)多相流 CT系統(tǒng)，投影數(shù)據(jù)少導(dǎo)致的圖像重構(gòu)病態(tài)問(wèn)題并未很好地解決．筆者針對(duì)設(shè)計(jì)的多相流 CT系統(tǒng)，改進(jìn)了投影矩陣計(jì)算和迭代重建算法，并通過(guò)實(shí)驗(yàn)證明該方法的有效性．

1 多相流CT系統(tǒng)簡(jiǎn)介

多相流 CT系統(tǒng)由γ射線源、探測(cè)器陣列、數(shù)據(jù)采集處理及圖像重建單元組成[1]．射線穿過(guò)多相流體后衰減，其衰減率與介質(zhì)吸收系數(shù)有關(guān)，測(cè)量射線的衰減程度即射線投影，經(jīng)數(shù)據(jù)處理及圖像重構(gòu)，可再現(xiàn)管內(nèi)多相流體分布．

如圖 1(a)所示，重建區(qū)域分成J nn= × 像素，依次編號(hào)為 1～J．在第 j像素中，射線吸收系數(shù)可視為常數(shù)，以像素值xj表示．j像素對(duì)i射線投影的貢獻(xiàn)權(quán)值以加權(quán)因子rij表示，則射線i總的射線投影為

式中：I為射線數(shù)．為保證成像質(zhì)量，I、J一般選取較大，所以投影矩陣計(jì)算復(fù)雜．但對(duì)每條射線，它僅與不超過(guò)2,n像素相交，因此R為大型稀疏矩陣．

圖1 多相流CT系統(tǒng)Fig.1 Multiphase flow CT system

2 擴(kuò)充投影矩陣

忽略射線寬度，并設(shè)其間距為τ，加權(quán)因子 rij以2種方法定義．

(2) rij＝i號(hào)射線在j號(hào)像素內(nèi)的長(zhǎng)度．

針對(duì)以上2種定義，提出投影數(shù)較少的CT系統(tǒng)扇束投影矩陣的2種算法，并比較二者的成像效果．

2.1 逐點(diǎn)插值算法

(1) 賦R初值為I×J的零矩陣；

(2) 將成像區(qū)域內(nèi)的像素E的直角坐標(biāo)(x, y)轉(zhuǎn)化成極坐標(biāo)(r,θ)；圖 2中 S0代表射線源， C1C2代表探測(cè)器陣列，由S0發(fā)出的射線經(jīng)過(guò)E，與C1′C2′交于A點(diǎn)，與C1C2交于B點(diǎn)；Os為C1C2的中點(diǎn)；C1′C2′//C1C2//EF．由△ A OS0～△ E FS0可得O與A間距離為

圖2 等距射線扇形束參數(shù)關(guān)系Fig.2 Relations among parameters of fan beam

若i是 1～17間的整數(shù)，則對(duì) rij賦值為 1，若i是1～17間的小數(shù)，則令 a = f loor(i)，(floor(i)取i的整數(shù)部分)，對(duì) ra+1,j賦值為t = i - a ，而 ra,j＝1- t．

(4) 依次計(jì)算全部像素點(diǎn)，循環(huán) J 次．

該方法可以圖3解析．由于投影稀疏，在射線1、2間存在大量無(wú)射線穿過(guò)的像素，若令其加權(quán)因子為零，這些像素值將得不到修正，當(dāng)該像素過(guò)多時(shí)將無(wú)法成像．現(xiàn)假設(shè)存在第1.4條射線穿過(guò)像素e，則e對(duì)射線1投影的貢獻(xiàn)權(quán)值為 0.6，對(duì)射線 2投影的貢獻(xiàn)權(quán)值為 0.4．這種賦值方法相當(dāng)于卷積反投影重建算法中的線性插值，因此稱其為插值算法．

(3) 由s算出經(jīng)過(guò)像素jx的射線號(hào)

2.2 逐線循環(huán)算法

如圖 3所示，將重建區(qū)域劃分為 n×n方格，依次編號(hào)為1，2，…，n2．若射線斜率的絕對(duì)值不大于1，設(shè)距離參數(shù) d表示射線與網(wǎng)格的交點(diǎn)到網(wǎng)格底邊的距離，d1為當(dāng)前的d，d2表示在x方向移動(dòng)1個(gè)網(wǎng)格后的d．若射線斜率的絕對(duì)值大于 1，則 d為橫向距離，每步在y方向移動(dòng)1個(gè)網(wǎng)格，在圖中相應(yīng)表示為 d1′、．5個(gè)源與對(duì)應(yīng)探測(cè)器陣列的組合依次編為 1～5組，第3組探測(cè)器1～17的縱坐標(biāo)為-8×16～8×16，如圖4所示．

圖3 成像區(qū)域坐標(biāo)Fig.3 Imaging region on xy-plane

圖4 5源模型投影示意Fig.4 Skech map of projections in 5 sources model

步驟1在-128～128之間每隔τ＝2插入1條虛擬射線，總射線數(shù)由 5×17＝85增至 5×129＝645，實(shí)現(xiàn)了信息擴(kuò)充，同時(shí)以射線穿過(guò)像素的長(zhǎng)度為投影系數(shù)，較式(3)更為精確；

步驟2根據(jù)幾何模型，計(jì)算每條射線方程：第3組的各條射線可由其起點(diǎn)與終點(diǎn)確定其方程，而 1、2、4、5組的各條射線可由第 3組的相應(yīng)射線旋轉(zhuǎn)相應(yīng)角度獲得，易得其點(diǎn)斜式方程．事實(shí)上，由于模型關(guān)于x軸對(duì)稱，只需計(jì)算第1、2組以及第3組的一半射線穿過(guò)的網(wǎng)格及長(zhǎng)度，其他射線穿過(guò)的網(wǎng)格號(hào)及長(zhǎng)度可由其對(duì)稱射線的計(jì)算結(jié)果直接獲得．例如關(guān)于x軸對(duì)稱的網(wǎng)格Ea、Eb，其關(guān)系式為

式中 re m(Ea- 1 ,n)表示取 ( Ea- 1 )/n 的余數(shù)．

步驟3根據(jù)射線斜率與源位置判斷射線i與重建區(qū)域相交邊界，算出射線i與邊界的交點(diǎn)坐標(biāo)，由此得到射線i穿過(guò)的第 1個(gè)網(wǎng)格(記為0E)及最后 1個(gè)網(wǎng)格(記為1E)；

步驟4射線斜率k分為k＞1，0≤k≤1，-1 ≤ k ＜0，k ＜ -1，分別通過(guò)4種不同的循環(huán)，由 d1、k、δ求出 d2，根據(jù) d2求出射線穿過(guò)的網(wǎng)格號(hào)及長(zhǎng)度，網(wǎng)格號(hào)存入矩陣 G，長(zhǎng)度存入矩陣 L．以第 2組k＞ 1 的情況為例，計(jì)算流程如圖 5所示，其中s =，j初值為1．其他情況可參照編程．

圖5 射線i穿過(guò)網(wǎng)格號(hào)及長(zhǎng)度計(jì)算流程Fig.5 Flow chart of loop for ray i

步驟 5若當(dāng)前網(wǎng)格號(hào) E = E1，射線i計(jì)算完成，否則計(jì)算新的 d2，求下一個(gè)網(wǎng)格號(hào)和長(zhǎng)度；

步驟6依次計(jì)算各條射線，循環(huán) I 次；

步驟 7以 G中元素為列標(biāo)，L 中對(duì)應(yīng)元素為值，將 G重排成 85行的矩陣 R，或用插值法對(duì)測(cè)量矢量P進(jìn)行擴(kuò)充，將射線投影值加權(quán)后賦予相鄰虛擬射線．

對(duì)每組虛擬112條射線，參考平行射束的投影矩陣算法[2-3]計(jì)算射線穿過(guò)的網(wǎng)格號(hào)及被網(wǎng)格所截長(zhǎng)度，較之于逐點(diǎn)循環(huán)，由于I遠(yuǎn)少于J，并利用實(shí)際系統(tǒng)的對(duì)稱性簡(jiǎn)化計(jì)算(以 5源系統(tǒng)為例，計(jì)算量減少1/2；若系統(tǒng)結(jié)構(gòu)關(guān)于x和y軸均對(duì)稱，則可減少3/4)，極大提高了計(jì)算速度，且像素越多，速度提高越明顯．

由于平均每條射線穿過(guò)n個(gè)像素，而上述 2種算法使投影矩陣具有遠(yuǎn)大于85n個(gè)非零值，故可視為投影矩陣的擴(kuò)充式算法．

3 盲去卷積復(fù)原

設(shè) f為原始圖像，g為退化圖像，h為點(diǎn)擴(kuò)展函數(shù)(point spread function，PSF)，“*”表示卷積，則圖像退化模型為

圖像復(fù)原即由g、h對(duì)f進(jìn)行最佳估計(jì)．將 CT系統(tǒng)視為一個(gè)圖像退化系統(tǒng)，重建圖像為退化圖像，則可將其代入圖像復(fù)原算法．由于CT系統(tǒng)的點(diǎn)擴(kuò)展函數(shù)(PSF)未知，故選擇盲去卷積算法．盲去卷積僅利用退化圖像，即可同時(shí)估算PSF和清晰圖像[4]．

設(shè)f、g均為離散概率頻率函數(shù)，那么f、g在像素i上的值if，ig可以認(rèn)為是事件(假定收集到的單位光子為一個(gè)事件)在該像素點(diǎn)上發(fā)生的頻率數(shù)．由貝葉斯條件概率公式得

4 改進(jìn)重建算法

4.1 代數(shù)重建算法

Herman提出的代數(shù)重建算法(algebrai reconstruction technique，ART)迭代公式為式中：ik= k ( modI )+1，k為迭代次數(shù)，ik為射線號(hào)；pi為測(cè)得的射線i的投影，以探測(cè)的粒子數(shù)表示[5-6]．

基于式(11)重建，經(jīng)3I ～ 8I次迭代，結(jié)果較好[5].為進(jìn)一步提高成像質(zhì)量與實(shí)時(shí)性，將盲去卷積算法融入代數(shù)重建算法的迭代過(guò)程(記為 BD-ART)，流程如圖6所示．其中，ε為根據(jù)精度要求指定的正數(shù)，KI為總迭代次數(shù)．實(shí)驗(yàn)證明該方法經(jīng)1I ～ 2I次迭代即可取得優(yōu)于ART的效果．

圖6 加入盲去卷積的ART算法流程Fig.6 Flow chart of BD-ART

4.2 同乘代數(shù)重建算法

Byrne提出的同乘代數(shù)重建算法(SMART)迭代公式為

最小交叉熵算法適用于投影數(shù)據(jù)較少的場(chǎng)合，重建質(zhì)量?jī)?yōu)，但需要非均勻分布的先驗(yàn)信息．在缺乏先驗(yàn)信息情況下將其應(yīng)用于多相流 CT系統(tǒng)，需進(jìn)行適當(dāng)改進(jìn)．將盲去卷積算法融入 SMART的迭代過(guò)程(記為 BD-SMART)，以復(fù)原圖像為先驗(yàn)信息，進(jìn)行迭代計(jì)算，流程如圖7所示．其中，K為迭代次數(shù)．

圖7 加入盲去卷積的SMART算法流程Fig.7 Flow chart of BD-SMART

4.3 圖像處理

圖像處理步驟如下．

步驟 1將管道外部的像素賦值成任意背景值，以消除物場(chǎng)外的偽影．

步驟2利用盲去卷積估算的 PSF，對(duì)圖像進(jìn)行Richardson-Lucy 復(fù)原[9]．

步驟 3根據(jù)成像結(jié)果，選擇適當(dāng)閾值分割圖像[10]．

5 實(shí)驗(yàn)結(jié)果與分析

計(jì)算機(jī)配置為Pentium(R)4 2.93GHz CPU、504 MB內(nèi)存，模擬 4種油水兩相流流型(見(jiàn)圖 8)并測(cè)得投影數(shù)據(jù)，以MATLAB7.1為平臺(tái)驗(yàn)證本文方法．

圖8 油水兩相流流型Fig.8 Flow regimes of oil-water two phase flow

未擴(kuò)充的投影矩陣記為 R0，逐點(diǎn)插值法計(jì)算的擴(kuò)充投影矩陣記為 R1，逐線循環(huán)法計(jì)算的擴(kuò)充投影矩陣記為R2．

圖9為采用R0基于ART重建的層流圖像，幾乎無(wú)法由該圖像識(shí)別流型．對(duì)比表 1可見(jiàn)，擴(kuò)充投影矩陣后成像效果顯著改善．

圖9 采用R0重建層流圖像Fig.9 Laminar flow image reconstructed with R0

像素足夠多時(shí)，逐線循環(huán)法才比逐點(diǎn)插值法的速度明顯快，故表 1中，重建像素為 200×200的圖像．對(duì)比發(fā)現(xiàn)，采用 R2明顯提高了實(shí)時(shí)性；加入盲去卷積運(yùn)算后，成像質(zhì)量明顯提高；雖然盲去卷積運(yùn)算使迭代耗時(shí)增加，但由于所需迭代次數(shù)顯著減少，成像速度得到提高．與采用R1的ART相比，采用R2的BD-ART在實(shí)時(shí)性與成像質(zhì)量方面均有明顯優(yōu)勢(shì)．

表1 基于ART與BD-ART的成像結(jié)果Tab.1 Results obtained using ART and BD-ART

考慮本文中成像區(qū)域尺寸與分辨率要求，取 n＝58即可，對(duì)于58×58像素的圖像，R1與R2耗時(shí)相當(dāng)，而 R1計(jì)算簡(jiǎn)單，故表 2中采用 R1．對(duì)比表 2第 2、4行，可見(jiàn) BD-SMART明顯提高了成像質(zhì)量．第 3行為使用文獻(xiàn)[1]提出的3倍信息擴(kuò)充SMART算法(記為IE-SMART)重建的結(jié)果，對(duì)比BD-SMART的重建圖像，后者重建的偽影更少，速度更快，而前者使圖像邊緣趨于平滑．將2種方法結(jié)合(記為IE-BDSMART)，對(duì)較復(fù)雜流型重建的圖像質(zhì)量更佳，但實(shí)時(shí)性有所下降．

BD-SMART的重建結(jié)果經(jīng)圖像處理后，提取二值圖像面積，計(jì)算截面的分相含率，如表3所示．

表2 基于SMART與改進(jìn)SMART的成像結(jié)果Tab.2 Results obtained using SMART and improved SMART

表3 由重建圖像計(jì)算分相含率Tab.3 Phase volume fraction obtained with reconstructed images%

6 結(jié) 語(yǔ)

本文提出針對(duì)5源17探測(cè)器的工業(yè)多相流CT系統(tǒng)的擴(kuò)充投影矩陣計(jì)算，并將盲去卷積算法融入代數(shù)重建算法與同乘代數(shù)重建算法的迭代過(guò)程．利用改進(jìn)的算法對(duì)多相流 CT系統(tǒng)進(jìn)行圖像重建，將采用相同投影矩陣不同迭代過(guò)程、不同投影矩陣相同迭代過(guò)程的重建圖像精度和收斂速度進(jìn)行比較．實(shí)驗(yàn)結(jié)果表明，本文方法利用少量投影數(shù)據(jù)即可取得滿意的成像效果，能夠準(zhǔn)確識(shí)別流型，重建圖像的分相含率誤差大致在±1%以內(nèi)．

［1］王化祥，王 琦，郝魁紅. 最小交叉熵圖像重建算法[J]. 儀器儀表學(xué)報(bào)，2009，30(12)：2574-2579.

Wang Huaxiang，Wang Qi，Hao Kuihong. Reconstruction images based on minimum cross-entropy[J].Chinese Journal of Scientific Instrument，2009，30(12)：2574-2579(in Chinese).

［2］Herman G，Meyer L. Algebraic reconstruction can be made computationally efficient[J]. IEEE Transactions onMedical Imaging，1993，12(3)：600-609.

［3］張順利，張定華，李 山，等. ART算法快速圖像重建研究[J]. 計(jì)算機(jī)工程與應(yīng)用，2006，42(24)：1-3.

Zhang Shunli，Zhang Dinghua，Li Shan，et al. Research of fast image reconstruction on ART algorithm[J].Computer Engineering and Applications，2006，42(24)：1-3(in Chinese).

［4］Levente Kovács，Tamás Szirányi. Relative focus map estimation using blind deconvolution[J].Optics Letters，2005，30(22)：3021-3023.

［5］莊天戈. CT原理與算法[M]. 上海：上海交通大學(xué)出版社，1992.

Zhuang Tiange.Principles and Algorithms of CT[M].Shanghai：Shanghai Jiao Tong University Press，1992(in Chinese).

［6］Herman G T，Lent A. Iterative reconstruction algorithm[J].Computers in Biology and Medicine，1976，6(4)：273-294.

［7］Byrne C L. Iterative image reconstruction algorithms based on cross-entropy minimization[J].IEEE Transactions on Image Processing，1993，2(1)：96-103.

［8］Boer P T，Kroese D P，Mannor S，et al. A tutorial on the cross-entropy method[J].Annals of Operations Research，2005，134(1)：19-67.

［9］Masschaele B，Dierick M，Van Hoorebeke L，et al. Neutron CT enhancement by iterative de-blurring of neutron transmission images[J].Nuclear Instruments and Methods in Physics Research A，2005，542(1/2/3)：361-366.

［10］Golosio B，Masala G L，Piccioli A，et al. A novel multithreshold method for nodule detection in lung CT[J].Medical Physics，2009，36(8)：3607-3618.