劉躍進(jìn)，孔麗娜，白鴿，韓路長(zhǎng)，曹楊，羅和安

(湘潭大學(xué)化工學(xué)院化工過(guò)程模擬與優(yōu)化教育部工程研究中心，湖南湘潭411105)

CT(computed tomography)是利用輻射信號(hào)對(duì)被測(cè)對(duì)象進(jìn)行斷面掃描，采集通過(guò)被測(cè)對(duì)象衰減后的掃描信號(hào)，利用計(jì)算機(jī)手段對(duì)所采集信號(hào)進(jìn)行數(shù)學(xué)處理和圖像重構(gòu)，得到被測(cè)對(duì)象某種指標(biāo)的斷面圖像和數(shù)字分布，其發(fā)射信號(hào)一般為χ或γ射線，是一種無(wú)介入、非接觸、不受體系濃度和高溫高壓限制的探測(cè)方法，可以一次得到被測(cè)對(duì)象斷面二維測(cè)量結(jié)果，測(cè)量結(jié)果可以圖像和數(shù)字化同時(shí)顯示［1］.γ射線由于具有較強(qiáng)的穿透力在國(guó)外早已被應(yīng)用于多相流研究［2-4］，國(guó)內(nèi)雖然有關(guān)于 γ-CT 的研究報(bào)道［5-7］，但實(shí)際應(yīng)用測(cè)量成像的報(bào)道甚少［8-9］.在多相流領(lǐng)域，CT處理的對(duì)象是具有瞬時(shí)隨機(jī)不確定性分布的多相流體，這與處理具有靜態(tài)確定分布的多相物質(zhì)有一定區(qū)別，所涉及到的數(shù)學(xué)處理和圖像重構(gòu)方法也不盡相同.γ射線CT用于多相流測(cè)量，主要是測(cè)量多相流斷面分相分率分布，如氣液體系氣含率分布和氣固或者液固體系固含率分布［10］.

1 CT測(cè)量氣含率分布原理

當(dāng)γ射線輻射信號(hào)通過(guò)流體介質(zhì)時(shí)，符合Beer-Lambert定律:

式中:Φ為ρμ的乘積，是流體參數(shù);I0、I分別為放射源發(fā)射信號(hào)、探測(cè)器接受信號(hào)的強(qiáng)度.

對(duì)單相流體，分別有

下標(biāo)d、f分別代表流體分散相和連續(xù)相.

對(duì)于兩相流體，有

當(dāng)以ε表示兩相流體的分相分率(指體積分率);L為輻射所通過(guò)兩相流體距離，則有

由式(2)和(5)，可得

式(6)反映了利用輻射衰減定律來(lái)求兩相流體分相分率的基本原理，但這里求得的是輻射距離L上流體分相的平均分率，而分相分率ε存在著一定的分布.當(dāng)以f(l)表示流體參數(shù)Φ在任一輻射距離L上的分布，則有

式(7)是普遍意義上的Beer-Lambert輻射衰減定律，它將輻射信號(hào)的衰減變化與流體參數(shù)Φ的分布關(guān)聯(lián)起來(lái).

CT信號(hào)的發(fā)射與收集有多種方式，如圖1是一個(gè)輻射源和多個(gè)探測(cè)器旋轉(zhuǎn)的方式.按照CT裝置的空間分辨率，將兩相流體斷面劃分為很多像素小區(qū)間.圖1所示是一個(gè)輻射照的情況.

參照?qǐng)D2，容器斷面上任一割線L的方程可表示為

這樣，通過(guò)容器斷面任一輻射信號(hào)衰減方程可表示為

方程(9)實(shí)際上是一個(gè)圖像重構(gòu)問(wèn)題.數(shù)學(xué)上可以證明，如果P(r，θ)在無(wú)限平面上連續(xù)，且對(duì)r有連續(xù)的一階偏導(dǎo)數(shù)，通過(guò)Radon逆變換［11］，可以求得流體參數(shù)的唯一分布函數(shù)為

由于兩相流體積具有線性加和性，則

這樣可求得氣含率在容器斷面上的分布為

式中:Φd、Φf應(yīng)由式(5)求出以消除容器壁對(duì)γ射線信號(hào)的衰減影響.

圖1 一個(gè)輻射照中射線通過(guò)容器像素示意圖Fig.1 The schematic of CT projections through pixel area of reactor in a view

圖2 幾何構(gòu)型的層析投影Fig.2 The geometrical configuration of tomographic

以上說(shuō)明了γ射線CT測(cè)量?jī)上嗔鞣稚⑾喾植嫉幕驹砗土黧w參數(shù)分布Φ的唯一可解性，式(10)、(12)即為氣含率分布求取式.但由于γ射線CT輻射信號(hào)的發(fā)射、衰減、接受和容器中兩相流體的瞬時(shí)分布均是隨機(jī)過(guò)程，因此式(10)、(12)很難直接求解，而必須采用數(shù)理統(tǒng)計(jì)方法來(lái)建立氣含率分布的數(shù)字成像模型.

2 CT測(cè)量氣含率分布數(shù)字成像模型

2.1 最大似然估計(jì)方法

對(duì)于方程(9)這樣的圖像重構(gòu)問(wèn)題，可以通過(guò)大量的實(shí)驗(yàn)測(cè)量數(shù)據(jù)P(r，θ)對(duì)f(x，y)進(jìn)行估計(jì)來(lái)實(shí)現(xiàn) CT 彩色數(shù)字圖像［12］，如過(guò)濾背投影法［13］、代數(shù)重構(gòu)法［14］和最大似然估計(jì)法［15］等.當(dāng)把 CT發(fā)射和接受信號(hào)過(guò)程看成Poisson分布、流體吸收衰減信號(hào)過(guò)程為二項(xiàng)分布時(shí)，適用最大似然估計(jì)法來(lái)進(jìn)行參數(shù)估計(jì).

設(shè)被測(cè)變量樣本空間為Y，其概率密度分布函數(shù)為g(Y，Φ)，Φ為待估計(jì)參數(shù).要求一個(gè)完整母空間X，其概率密度分布函數(shù)為f(X，Φ)，且與Y具有對(duì)應(yīng)關(guān)系h(X)=Y，Y是X的一個(gè)子空間，具體步驟為:

1)在已知實(shí)驗(yàn)樣本空間Y和參數(shù)Φ(開(kāi)始時(shí)需假定Φ的初值)條件下，基于發(fā)射和接受輻射信號(hào)的Poisson分布和輻射衰減信號(hào)的二項(xiàng)分布隨機(jī)模型，對(duì)完整母空間X進(jìn)行估計(jì):

2)在所估計(jì)X空間下，求其聯(lián)合概率密度分布函數(shù)，并對(duì)其聯(lián)合概率密度分布函數(shù)采用線性最大似然估計(jì)法求參數(shù)Φ的迭代式:

先求表達(dá)式:

令

求得

3)以Y、X和Φ代入式(14)判斷是否達(dá)到最大，如果是，則Φ即為所求;如果不是，重復(fù)步驟1)、2)，直至達(dá)到最大為止.

2.2 基本假定

假定1 在一輻射照內(nèi)以橫排像素為單元進(jìn)行參數(shù)估值.雖然這個(gè)假設(shè)與實(shí)際情況出入很大，但這樣做的好處一是使估值過(guò)程大大簡(jiǎn)化并得以進(jìn)行，二是這個(gè)誤差權(quán)重只占1%左右，因?yàn)樾D(zhuǎn)一周有近百個(gè)輻射照，而最后結(jié)果是對(duì)所有輻射照情況下所估值再求平均.

假定2 放射源發(fā)射信號(hào)數(shù)即為進(jìn)入像素橫排Gi的信號(hào)數(shù)，放射源發(fā)射信號(hào)數(shù)和探測(cè)器接受輻射衰減信號(hào)數(shù)yi均分別用平均值為d和λi為的Poisson概率分布來(lái)描述，而它們之間的信號(hào)傳遞則用二項(xiàng)概率分布來(lái)描述.

2.3 數(shù)學(xué)推導(dǎo)

放射源在一個(gè)輻射照時(shí)間Δt內(nèi)所發(fā)射信號(hào)數(shù)的平均值為

根據(jù)輻射信號(hào)衰減式(7)，任一探測(cè)器在Δt時(shí)間內(nèi)所接受信號(hào)數(shù)的平均值為

按照假設(shè)2，有

式中:M、N分別為一個(gè)輻射照內(nèi)的輻射線總數(shù)和橫排像素總數(shù);Gi、Ni分別表示與輻射線i初次相遇和最后離開(kāi)的橫排像素序號(hào)數(shù).

當(dāng)為xi，j的輻射信號(hào)進(jìn)入橫排像素j時(shí)，一部分被流體介質(zhì)所吸收，一部分又繼續(xù)輻射出去.按照假定2，離開(kāi)的輻射信號(hào)數(shù)為xi，j+1的概率是:

各橫排像素之間是相互獨(dú)立的，故所有像素區(qū)間輻射信號(hào)數(shù)的聯(lián)合概率密度為

取對(duì)數(shù)，得

式中:R與Φj為無(wú)關(guān)項(xiàng).按照式(7)，進(jìn)入橫排像素j的輻射信號(hào)數(shù)xi，j的期望值是:

又

在給定進(jìn)入橫派像素j的輻射信號(hào)數(shù)xi，j條件下，探測(cè)器接受輻射衰減信號(hào)數(shù)yi的概率是:

xi，j服從平均值為 γi，j的 Poisson 分布，而在確定

xi，j條件下的 yi則服從成功概率為 γyi/γi，j的二項(xiàng)分布，因此這兩者的聯(lián)合概率密度是:

給定yi條件下xi，j的概率密度是:

由式(30)可知，(xi，j-yi)在給定 yi條件下服從期望值為(γi，j-γyi)的 Poisson 分布，故有

式(32)用來(lái)完成式(13)即對(duì)X的估值.對(duì)式(23)中的Φj求偏導(dǎo)，并令其等于0，得

當(dāng)容器斷面S＜1時(shí)，有如下近似式

常見(jiàn)流體水的 Φ為0.086 cm，空間分辨率2 mm，ΔliΦj=0.017 2，滿足式(34)，這樣可解得

式(35)即為φj的迭代式.

通過(guò)式(25)～(27)、(32)、(35)可對(duì)參數(shù)Φ進(jìn)行反復(fù)估值，直至式(23)達(dá)到最大值為止.

改變輻射照的位置，又可對(duì)另外一組流體參數(shù)Φ矢量進(jìn)行類(lèi)似的估值.然后按照所有像素小區(qū)間的幾何位置，加和求平均，即得各像素小區(qū)間的流體參數(shù)Φ.對(duì)每個(gè)像素小區(qū)間，有

得

這樣就求出了每個(gè)小像素區(qū)間分相分率，其計(jì)算過(guò)程用Fortrun語(yǔ)言編程計(jì)算，并且通過(guò)Matlab實(shí)現(xiàn)分相分率的彩色圖像顯示和數(shù)字變化曲線.

3 實(shí)驗(yàn)條件

3.1 測(cè)量參數(shù)

實(shí)驗(yàn)測(cè)量參數(shù)如下:

1)100m Ci的137Cs為放射源，7個(gè)NaI晶體探測(cè)器扇形排列，每個(gè)探測(cè)器前面都裝有一個(gè)長(zhǎng)7 cm、直徑2.5 cm的準(zhǔn)直器，準(zhǔn)直器中間有一個(gè)寬2 mm、高5 mm、長(zhǎng)70 mm的矩形孔讓放射線通過(guò)到達(dá)探測(cè)器，以減少Compton漫反射效應(yīng).CT裝置的空間分辨率和密度分辨率分別是2 mm和40 kg/m3.

2)采用20 cm內(nèi)直徑的標(biāo)準(zhǔn)Rushton氣液攪拌釜，壁厚6.5 mm，擋板厚度1.5 mm，4 塊，攪拌漿軸直徑8 mm，在15 cm液高處掃描.

3)探測(cè)器間距5°，整組探測(cè)器又以0.2°間隔旋轉(zhuǎn)4.8°形成一個(gè)輻射照.在一個(gè)輻射照內(nèi)一個(gè)探測(cè)器停留25個(gè)測(cè)量位置，每個(gè)測(cè)量位置停留5.7 s，采樣頻率20 Hz.這樣在一個(gè)測(cè)量位置采樣114個(gè)，最后在360°范圍內(nèi)均勻旋轉(zhuǎn)測(cè)量100個(gè)輻射照，共采集1 995 000個(gè)數(shù)據(jù)，形成17 500條輻射線信號(hào)，以滿足γ射線隨機(jī)信號(hào)采樣和CT測(cè)量裝置空間分辨率的要求.旋轉(zhuǎn)一周共花4 h之久，又是依據(jù)大量隨機(jī)信號(hào)進(jìn)行統(tǒng)計(jì)處理，因此測(cè)得的是氣含率分布的時(shí)均值.

3.2 掃描位置

在攪拌時(shí)，攪拌槳透平上、下方空間會(huì)形成一個(gè)負(fù)壓循環(huán)區(qū)，分布器上升的氣泡被攪拌槳葉片及液面剪切力打碎分散，形成大量小氣泡，一部分快速進(jìn)入此內(nèi)循環(huán)負(fù)壓空間，另一部分被旋轉(zhuǎn)液體甩向外圍.通過(guò)實(shí)驗(yàn)觀察，掃描位置被選定在形成相對(duì)均勻的氣含率斷面分布的3/4靜液高處.

3.3 CT掃描實(shí)驗(yàn)設(shè)計(jì)

氣液Rushton攪拌釜?dú)夂史植糃T掃描實(shí)驗(yàn)設(shè)計(jì)見(jiàn)表1.

表1 氣含率分布CT掃描實(shí)驗(yàn)Table 1 CT scan experiment of gas holdup distribution r/m

4 結(jié)果與討論

氣含率分布受攪拌釜尺寸及操作條件等因素的影響，為研究所構(gòu)建的CT測(cè)量氣含率分布彩色數(shù)字成像模型可行性及準(zhǔn)確性，本文對(duì)比研究了不同操作條件下氣含率斷面CT彩色數(shù)字圖像.

4.1 攪拌轉(zhuǎn)速的影響

圖3 18.88 L/min通氣量下不同攪拌轉(zhuǎn)速CT圖像Fig.3 CT images of gas holdup distribution with gas flow rate 18.88 L/min at different impeller rotating speeds

圖3為通氣量18.88 L/min，攪拌轉(zhuǎn)速不同氣含率斷面CT彩色數(shù)字圖像.其清楚表明CT圖像與釜實(shí)際尺寸完全一致，這說(shuō)明了數(shù)字成像算法的可行性;同時(shí)，CT彩色數(shù)字圖像也清楚形象地表現(xiàn)出3種不同攪拌轉(zhuǎn)速下不同攪拌區(qū)氣含率分布的差別.從CT彩色成像圖還可看出，隨著轉(zhuǎn)速的增加，氣含率主要在攪拌槳正上方空間增加較多.為了驗(yàn)證CT彩色數(shù)字圖像的準(zhǔn)確性，圖4列出了上述3種不同攪拌轉(zhuǎn)速下平均氣含率沿?zé)o因次半徑的分布曲線.

圖4 通氣量18.88 L/min下不同攪拌轉(zhuǎn)速氣含率沿?zé)o因次半徑分布曲線Fig.4 The effect of impeller rotating speed on gas holdup distribution at gas rate 18.88 L/min

由圖4可知，在同一通氣量下，不同攪拌轉(zhuǎn)速的徑向氣含率分布趨勢(shì)幾乎相同，從攪拌器中心位置到器壁均呈先增大后減小趨勢(shì)，當(dāng)無(wú)因次半徑在0.25～0.6之間時(shí)，達(dá)到最大，并且此區(qū)間氣含率分布隨攪拌轉(zhuǎn)速的增加而明顯增加，由氣含率斷面CT彩色數(shù)字圖像圖角度也可得出同樣結(jié)論，證明了實(shí)驗(yàn)測(cè)量及數(shù)字成像算法的準(zhǔn)確性.

4.2 通氣量的影響

圖5為攪拌轉(zhuǎn)速800 r/min，不同通氣量的氣含率斷面分布CT彩色成像圖.如圖所示，在3種通氣量下，氣含率分布均隨通氣量的增大而增大，通氣量對(duì)氣含率的影響明顯大于攪拌轉(zhuǎn)速對(duì)其的影響.圖6列出了上述3種不同通氣量下平均氣含率沿?zé)o因次半徑的分布曲線.

如上圖所示，在整個(gè)無(wú)因次半徑區(qū)間，氣含率分布均隨通氣量的增大而增大，當(dāng)無(wú)因次半徑在0.25到0.6之間時(shí)，增加趨勢(shì)更為顯著.此結(jié)果與圖5所呈現(xiàn)的結(jié)果基本一至，進(jìn)一步證明了實(shí)驗(yàn)測(cè)量及數(shù)字成像算法的準(zhǔn)確性.

圖5 攪拌轉(zhuǎn)速800 r/min下不同通氣量下CT圖像Fig.5 CT images of gas holdup distribution with impeller rotating speed 800 r/min at different gas flow rates

圖6 攪拌轉(zhuǎn)速800 r/min下不同通氣量平均氣含率沿?zé)o因次半徑變化曲線Fig.6 The effect of gas rate on gas holdup distribution at impeller rotating speed 800 r/min

4.3 分布器的影響

圖7為通氣量9.44 L/min，攪拌轉(zhuǎn)速500 r/min下2種分布器的氣含率斷面CT彩色數(shù)字圖像.觀察CT彩色數(shù)字圖像可知，16孔0.5 mm孔徑情況下，從圓環(huán)氣體分布器上升至攪拌槳間的氣泡較8孔1.5 mm孔徑下略密略小，但兩者差別不大.為進(jìn)一步定量比較，圖8用Excel列出了上述2種不同分布器下平均氣含率沿?zé)o因次半徑的分布曲線.

圖7 通氣量9.44 L/min，攪拌轉(zhuǎn)速500 r/min下不同分布Fig.7 CT images of gas holdup distribution at gas rate 9.44 L/min，impeller rotating speed 500 r/min with different spargers

圖8 通氣量9.44 L/min，攪拌轉(zhuǎn)速500 r/min下分布器平均氣含率沿?zé)o因次半徑變化曲線Fig.8 The effect of gas sparger on gas holdup distribution at gas rate 9.44 L/min，impeller rotating speed 500 r/min

如圖8所示，2種分布器下平均氣含率沿?zé)o因次半徑分布幾乎一致，只有攪拌槳區(qū)有略微差別，與圖7所觀察結(jié)果一致.

5 結(jié)論

1)構(gòu)建了γ射線CT測(cè)量Rushton攪拌釜?dú)夂史植嫉牟噬珨?shù)字成像模型.根據(jù)137Cs γ射線發(fā)射信號(hào)和接受信號(hào)的Poisson分布模型和兩相流體吸收衰減輻射信號(hào)的二項(xiàng)分布模型，在給出接受信號(hào)實(shí)驗(yàn)樣本空間條件下，采用最大似然參數(shù)估計(jì)方法，對(duì)容器斷面氣液兩相流體吸收衰減輻射γ射線信號(hào)的二維空間進(jìn)行了參數(shù)估計(jì)，推導(dǎo)和建立了其參數(shù)估計(jì)的數(shù)學(xué)模型.

2)運(yùn)用所構(gòu)建的氣含率分布的彩色數(shù)字成像模型，對(duì)Rushton氣液攪拌釜在不同通氣量、不同攪拌速率及不同氣體分布器狀態(tài)下的氣含率進(jìn)行了圖像重構(gòu)，通過(guò)CT彩色數(shù)字圖像與平均氣含率沿?zé)o因次半徑變化實(shí)驗(yàn)值的比較研究，證明所構(gòu)建的CT測(cè)量氣含率分布彩色數(shù)字成像的模型有效可行，所建CT彩色數(shù)字圖像清晰準(zhǔn)確、對(duì)比度好，分辨率高，受噪聲影響小.

［1］CHAOUKI J，LARACHI F.Noninvasive tomographic and velocimetric monitoring of multiphase flows［J］.Ind ＆ Eng Chem Research，1997，36(11):4476-4503 .

［2］THIYAGARAJAN T K，SATYAMURTHY P，DIXIT N S，et al.Void fraction profile measurements in two-phase mercury-nitrogen flows using gamma-ray attenuation method［J］.Experimental Thermal and Fluid Science，1995，10(3):347-35.

［3］BOYER C，F(xiàn)ANGET B.Measurement of liquid flow distribution in trickle bed reactor of large diameter with a new gamma-ray tomographic system［J］.Chemical Engineering Science，2002，57(7):1079-1089.

［4］MARKUS S，GUNTHER H.Liquid flow texture analysis in trickle bed reactors using high-resolution gamma ray tomography［J］.Chemical Engineering Journal，2008，140(1/2/3):332-340.

［5］蘇顯素，楊學(xué)恒，趙炬.我國(guó)首臺(tái)層析γ射線機(jī)的性能及其應(yīng)用［J］.無(wú)損檢測(cè)，1995，17(9):245-247.

SU Xiansu，YANG Xueheng，ZHAO Ju.Properties of the first domestic γ-ray ICT and its applications［J］.Nondestructive Testing，1995，17(9):245-247.

［6］廖志勇，鄒宏.γ射線及其應(yīng)用［J］.宜賓科技，2006，4(14):32-34.

LIAO Zhiyong，ZOU Hong.γ-ray and its applications［J］.Yibinkeji，2006，4(14):32-34.

［7］代榮.γ射線CCD探測(cè)技術(shù)研究［D］.北京:中國(guó)工程物理研究院，2007:1-5.

DAI Rong.γ-ray CCD detection technology［D］.Beijing:China Academy of Engineering Physics，2007:1-5.

［8］王化祥，戴鵬，郝魁紅，等.基于Geant4的多相流CT系統(tǒng)優(yōu)化設(shè)計(jì)［J］.天津大學(xué)學(xué)報(bào)，2008，41(11):1293-1298.

WANG Huaxiang，DAI Peng，HAO Kuihong，et al.Optimum design of multiphase CT system based on Geant4 software toolkits［J］.Journal of Tianjin University，2008，41(11):1293-1298.

［9］劉躍進(jìn)，韓路長(zhǎng)，羅和安，等.Cs-137γ射線CT測(cè)量標(biāo)準(zhǔn)Rushton攪拌釜?dú)夂史植技皵?shù)字圖像重建［J］.高校化學(xué)工程學(xué)報(bào)，2006(4):648-652.

LIU Yuejin，HAN Luchang，LUO Hean，et al.Measurement of gas hold-up distribution and digital color image reconstruction for standard gas-liquid rushton stirred tank with Cs-137 γ-CT［J］.Journal of Chemical Engineering of Chinese Universities，2006(4):648-652.

［10］RADON J.über die bestimmung von funktionen durch ihre integralwerte l?ngsgewissermanningfaltigkeiten ［J］.Berichte Saechsische Akad Wiss，1917，69:262-277.

［11］莊天戈.CT原理與算法［M］.上海:上海交通大學(xué)出版社，1992:97-101.

［12］SHEPP L A，LOGAN B F.The Fourier reconstruction of the head section［J］.IEEE Trans on Nuclear Sci，1974，NS-21(3):21.

［13］HERMAN G T.Image reconstruction from projections:implementation and applications［M］. Springer-Verlag，1979.

［14］LANGE K，CARSON R.EM reconstruction algorithms for emission and transmission tomography［J］.Journal of Computer Assisted Tomography，1984，8(2):306-316.