王小鑫,張艷麗,白 萌,席 煒,黨 博,胡紅利

(1.西安石油大學(xué) 陜西省油氣井測(cè)控技術(shù)重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710065;2. 西安交通大學(xué) 電力設(shè)備電氣絕緣國(guó)家重點(diǎn)實(shí)驗(yàn)室, 陜西 西安 710049)

兩相流普遍存在于石油、化工及制藥等工業(yè)領(lǐng)域的開(kāi)發(fā)和生產(chǎn)過(guò)程中。其流場(chǎng)分布參數(shù)演變復(fù)雜且在時(shí)間、空間上相互耦合,因此,研發(fā)能夠直觀顯示兩相流流動(dòng)狀態(tài)的測(cè)量技術(shù)至關(guān)重要[1-2]。電學(xué)層析成像(electrical tomography, ET)技術(shù)因其無(wú)輻射、適用性強(qiáng)、圖像成像快、檢測(cè)成本低等優(yōu)點(diǎn),在兩相流可視化測(cè)量中得到了廣泛的應(yīng)用。經(jīng)過(guò)多年的科學(xué)研究和實(shí)踐應(yīng)用,電學(xué)層析成像技術(shù)在傳感器、數(shù)據(jù)采集和圖像重建等方面取了大量研究成果[3-4]。其中,傳統(tǒng)ERT/EIT技術(shù)在測(cè)量過(guò)程中通常需要接觸待測(cè)介質(zhì),容易引起傳感器電極與被測(cè)流體之間發(fā)生物理和化學(xué)反應(yīng)[5]。電容層析成像(electrical capacitance tomography, ECT)技術(shù)則通過(guò)檢測(cè)場(chǎng)域中不同物質(zhì)對(duì)應(yīng)的介電常數(shù),重建出被測(cè)場(chǎng)域截面的介電常數(shù)分布,該技術(shù)主要適用于非導(dǎo)電物體組成的兩相流或者多相流檢測(cè)[6]。近年來(lái),鑒于含導(dǎo)電相多相流的非接觸式檢測(cè)需求,Wang和Xu等人將電容耦合非接觸電導(dǎo)檢測(cè)(capacitively coupled contactless conductivity detection,簡(jiǎn)稱C4D) 技術(shù)與ERT技術(shù)相結(jié)合,組成電容耦合電阻層析成像(CCERT)技術(shù)以實(shí)現(xiàn)非接觸式等效電阻測(cè)量[7-8];Jiang等研究一種電容耦合相位式介電譜層析成像技術(shù),用于含不同電導(dǎo)率(0.065～0.4 S/m)的固液兩相流層析成像[9];Gunes等提出位移電流相位層析成像技術(shù),相較于ECT/ERT,該技術(shù)用于連續(xù)相為有損介質(zhì)(電導(dǎo)率為5～50 mS/m)的混合流體時(shí)具有更寬的線性范圍[10]。

傳統(tǒng)電學(xué)層析成像的重構(gòu)算法部分通常需要借助Matlab軟件平臺(tái)實(shí)現(xiàn),不利于工業(yè)現(xiàn)場(chǎng)的實(shí)時(shí)監(jiān)測(cè)及工作人員的操作。NI公司創(chuàng)立的LabVIEW是目前國(guó)際上廣泛應(yīng)用的虛擬儀器開(kāi)發(fā)環(huán)境之一,已成為測(cè)量與測(cè)試領(lǐng)域的工業(yè)標(biāo)準(zhǔn),尤其在儀器控制、數(shù)據(jù)采集、數(shù)據(jù)分析及顯示等領(lǐng)域發(fā)揮了巨大優(yōu)勢(shì)。周英鋼等人研發(fā)了一種基于LabVIEW的直接三維電容層析成像系統(tǒng),其中, LabVIEW主要實(shí)現(xiàn)數(shù)據(jù)的采集[11];顧鵬同樣設(shè)計(jì)一款基于LabVIEW上位機(jī)的三維電阻抗成像系統(tǒng),該上位機(jī)主要負(fù)責(zé)數(shù)據(jù)的采集和預(yù)處理,并實(shí)現(xiàn)Matlab圖像重建程序的調(diào)用[12]。為了提高程序的執(zhí)行效率,劉琳將Matlab成像算法模塊打包成COM組件,供LabVIEW調(diào)用,進(jìn)而降低了成像的速度[13]。

本文設(shè)計(jì)了一種基于LabVIEW的非接觸電阻抗層析成像系統(tǒng),圖像重構(gòu)算法部分在LabVIEW上實(shí)現(xiàn),不需要額外調(diào)用Matlab程序,實(shí)現(xiàn)了數(shù)據(jù)的實(shí)時(shí)采集及圖像在線重構(gòu)。該系統(tǒng)通過(guò)阻抗測(cè)量電路獲取輸出電信號(hào)的實(shí)部、虛部和相位進(jìn)行分別成像,并結(jié)合主成分分析融合算法以提高重構(gòu)精度。

1 非接觸式電阻抗層析成像系統(tǒng)構(gòu)架

非接觸式電阻抗層析成像系統(tǒng)主要由12電極傳感器、測(cè)量采集電路以及用于圖像重建的PC機(jī)組成。其系統(tǒng)結(jié)構(gòu)如圖1所示。

圖1 電阻抗層析成像系統(tǒng)結(jié)構(gòu)

圖1中測(cè)量及采集電路主要由激勵(lì)信號(hào)發(fā)生電路、I/U(阻抗/電壓)測(cè)量電路、開(kāi)關(guān)控制電路、基于模擬乘法器的相敏解調(diào)電路、低通濾波器和模數(shù)轉(zhuǎn)換電路組成。測(cè)量及采集電路是整個(gè)系統(tǒng)的核心,連接陣列電極,控制每個(gè)電極的激勵(lì)、待測(cè)及接地狀態(tài),同時(shí),將得到阻抗測(cè)量值由串口傳輸?shù)絇C機(jī)完成可視化測(cè)量,在測(cè)量過(guò)程中通過(guò)PC端的上位機(jī)來(lái)進(jìn)行控制。

1.1 傳感器結(jié)構(gòu)及靈敏度矩陣計(jì)算

本系統(tǒng)采用12電極非接觸式陣列傳感器,其徑向截面結(jié)構(gòu)及單電極激勵(lì)的電勢(shì)分布如圖2所示。采用柔性電極板加工技術(shù)以保證其陣列電極的對(duì)稱性和均一性,電極片貼附在PVC管道周圍,并在外部包裹上接地金屬屏蔽罩,圖3為傳感器實(shí)物圖。其中,管道內(nèi)徑為46 mm,外徑為2 mm,柔性電極的張角為26°,電極的覆蓋率為86.6%,管道與屏蔽罩外殼之間的距離為8 mm。

圖2 電極結(jié)構(gòu)及電勢(shì)分布

圖3 12電極傳感器

電學(xué)層析成像系統(tǒng)的靈敏度矩陣一般由Geselowitz靈敏度理論計(jì)算得出,基于電勢(shì)分布法的靈敏度矩陣可表示為[14]

(1)

式中:Si,j(k)為傳感器電極i和電極j加激勵(lì)時(shí)場(chǎng)域內(nèi)第k個(gè)單元處的靈敏度;φi為電極i加激勵(lì)信號(hào)時(shí)的場(chǎng)內(nèi)電勢(shì);φj為電極j加激勵(lì)信號(hào)時(shí)的場(chǎng)內(nèi)電勢(shì);vk為第k個(gè)單元的體積;Ue是激勵(lì)電極的電壓。第k個(gè)單元格的電勢(shì)梯度φ(k)為

φ(k)=Exx+Eyy+Ezz

(2)

式中:x、y和z分別為x軸方向、y軸方向和z軸方向的單位矢量;Ex、Ey和Ez分別表示x軸方向、y軸方向和z軸方向的電場(chǎng)強(qiáng)度,本文借助Comsol有限元分析軟件獲得。進(jìn)一步在激勵(lì)電極和檢測(cè)電極的作用下,第k個(gè)單元的靈敏度矩陣可以表示為

(3)

本文所使用的靈敏度矩陣在空?qǐng)?管道內(nèi)部為空氣)狀態(tài)下獲取,靈敏度矩陣作為解決層析成像逆問(wèn)題數(shù)學(xué)建模中的權(quán)重?cái)?shù)矩陣,反映某一單元的介質(zhì)發(fā)生變化時(shí)引起的阻抗值的變化量。圖4為相鄰電極和相對(duì)電極加激勵(lì)時(shí)的敏感場(chǎng)分布圖。

圖4 靈敏場(chǎng)分布

1.2 測(cè)量電路

本系統(tǒng)測(cè)量電路的激勵(lì)源采用MAX038模擬信號(hào)發(fā)生器產(chǎn)生固定頻率為2 MHz的正弦交流信號(hào)。為了實(shí)現(xiàn)電阻抗實(shí)部與虛部的測(cè)量,引入了模擬乘法器相敏解調(diào)電路。模擬乘法器采用AD734構(gòu)成相敏解調(diào)電路,2個(gè)乘法器分別實(shí)現(xiàn)同向解調(diào)和正交解調(diào),如圖5所示。

圖5 實(shí)部、虛部測(cè)量結(jié)構(gòu)圖

正弦信號(hào)一路經(jīng)過(guò)I/U轉(zhuǎn)換電路調(diào)制得到包含阻抗信息的電壓信號(hào),一路通過(guò)移相器將正弦信號(hào)的相位超前或滯后90°,得到正交信號(hào),并將此信號(hào)作為乘法器的參考信號(hào),用于正交解調(diào);另一路正弦信號(hào)直接與乘法器相連,作為正交解調(diào)的參考信號(hào),用于同向解調(diào);最后經(jīng)過(guò)乘法解調(diào)后的信號(hào)經(jīng)過(guò)低通濾波器,分別得到包含阻抗的實(shí)部信息和虛部信息的直流電壓。此外,相位信息可由實(shí)部和虛部值進(jìn)一步計(jì)算得出。

設(shè)激勵(lì)信號(hào)Vi為

Vi=Asin(ωt+α)

(4)

式中:A為振幅;ω為角頻率;α為初始相角。

則I/U轉(zhuǎn)換電路輸出的時(shí)諧信號(hào)Vo為

AωRfCxcos(ωt+α)

(5)

式中:Rx為電極間有損介質(zhì)的等效電阻;Rf為放大電阻的反饋電阻;Cx為電極間有損介質(zhì)的等效電容。

Vo和Vi經(jīng)過(guò)第一個(gè)乘法器相乘可以得到

A2ωRfCxsin(ωt+α)cos(ωt+α)

(6)

通過(guò)低通濾波器后輸出的直流信號(hào)為

(7)

同理,將Vo和Vi經(jīng)過(guò)移相90°得到的正交信號(hào)相乘,并且經(jīng)過(guò)濾波器可以得到

(8)

得到相位的計(jì)算式為

(9)

2 基于LabVIEW的成像算法實(shí)現(xiàn)

2.1 層析成像逆問(wèn)題

圖像重建屬于逆問(wèn)題,即在已知結(jié)果的情況下,反推產(chǎn)生它的原因。阻抗層析的逆問(wèn)題是根據(jù)測(cè)得的阻抗信息來(lái)反演成像區(qū)域內(nèi)損耗因子的分布,求解逆問(wèn)題的實(shí)質(zhì)是求解靈敏度矩陣S的廣義逆S-1。類似EIT的歸一化線性模型可以表示為[15]

G=S-1Z

(10)

式中:G為介質(zhì)分布像素列向量(N×1);S-1為歸一化的靈敏度矩陣的廣義逆矩陣(N×M);Z為歸一化的阻抗向量(M×1)。通常情況下N?M,因此該問(wèn)題的解不唯一,存在欠定性,圖像重構(gòu)技術(shù)可以用于解決EIT的逆問(wèn)題。

本系統(tǒng)圖像重建部分在LabVIEW軟件上實(shí)現(xiàn),該軟件由NI公司研究開(kāi)發(fā),其優(yōu)勢(shì)在于LabVIEW使用的是圖形化編輯語(yǔ)言G編寫程序,以框圖的樣式完成功能實(shí)現(xiàn)[16]。本文在LabVIEW上編輯3種經(jīng)典圖像重構(gòu)算法,分別是LBP(linear back projection)算法、Landweber迭代算法、ART(algebraic reconstruction technique)代數(shù)重建算法。

2.2 成像算法的實(shí)現(xiàn)

2.2.1 LBP算法的實(shí)現(xiàn)

目前常用的LBP算法近似認(rèn)為STS=I(I為單位矩陣),則可用ST近似代替S的廣義逆S-1,LBP的表達(dá)式為[17]

G=STZ

(11)

式中:G為歸一化灰度向量;ST為靈敏度矩陣的轉(zhuǎn)置矩陣;Z為測(cè)量的阻抗向量。

在LabVIEW中利用數(shù)學(xué)線性代數(shù)矩陣A×B控件搭建LBP算法,以阻抗信息的實(shí)部計(jì)算結(jié)果為例,如圖6所示。

圖6 LBP算法模塊

2.2.2 Landweber算法的實(shí)現(xiàn)

Landweber算法是電學(xué)層析成像中較為經(jīng)典的圖像重建算法之一[18],Landweber算法最終求解的目標(biāo)為

(12)

以灰度值G為自變量,則目標(biāo)函數(shù)為

(13)

f(G)的梯度為

f(G)=ST(SG-Z)

(14)

由最優(yōu)化理論可知,最速下降法選擇負(fù)梯度方向?yàn)榈较騕19],因此Landweber算法的表達(dá)式為

Gk+1=Gk+αST(Z-SGk)

(15)

式中:α為增益因子,α=2/λmax,λmax為STS的最大特征值。

Landweber迭代算法的實(shí)現(xiàn),以實(shí)部計(jì)算為例,如圖7所示。該算法主要通過(guò)for循環(huán)和數(shù)學(xué)線性代數(shù)矩陣A×B控件實(shí)現(xiàn),為了實(shí)現(xiàn)迭代計(jì)算,需要通過(guò)移位寄存器將上一次的計(jì)算結(jié)果轉(zhuǎn)變?yōu)橄乱淮蔚某跏贾?。迭代次?shù)不同,則圖像重構(gòu)的顯示結(jié)果也不同。同時(shí),在Landweber圖像重構(gòu)算法中,計(jì)算迭代步長(zhǎng)α,計(jì)算過(guò)程見(jiàn)圖8。

圖7 Landweber迭代算法模塊

圖8 α計(jì)算值模塊

2.2.3 ART算法的實(shí)現(xiàn)

ART算法主要是在重構(gòu)范圍中給定初值[20],將數(shù)據(jù)殘差按照射線方向均勻地反投影,直到符合要求。

首先建立重建數(shù)據(jù)和投影數(shù)據(jù)之間的線性方程,如式(16)所示

(16)

式中:M為投影射線總數(shù);wij為權(quán)重因子,反映了像元fi對(duì)投影pi的貢獻(xiàn)。ART迭代算法中圖像重建的任務(wù)即由w和p求取f。由于N?M,采用迭代的方法計(jì)算,先給定一初值f(0),再按照式(17)逐次迭代,

(17)

結(jié)合ART數(shù)學(xué)模型公式與電學(xué)層析成像系統(tǒng)模型可以得到ART算法迭代公式,

(18)

式中:Sk表示靈敏度矩陣的第k行,ART算法迭代以第個(gè)k方程為基礎(chǔ),即在第k步迭代中,只有第k個(gè)歸一化值用于更新歸一化介電常數(shù)的分布[21]。

在LabVIEW中ART算法主要通過(guò)2個(gè)for循環(huán)實(shí)現(xiàn),如圖9所示。

圖9 ART算法模塊

2.3 基于PCA的圖像融合

為了進(jìn)一步提高圖像精度,本文基于主成分分析方法(PCA)對(duì)實(shí)部、虛部和相位圖像像素灰度值進(jìn)行融合。主成分分析法是一種使用最廣泛的數(shù)據(jù)降維算法[22],主要步驟如下。

1)Gr、Gi、Gp分別為歸一化后實(shí)部、虛部和相位的圖像灰度值向量,構(gòu)建歸一化后的圖像數(shù)據(jù)矩陣G=[GrGiGp]。

2)計(jì)算標(biāo)準(zhǔn)化后矩陣G的協(xié)方差矩陣。

(19)

式中:cov(GrGi) 為向量Gr和Gi的協(xié)方差。

3) 計(jì)算C的特征值和對(duì)應(yīng)的特征向量。選擇主特征向量對(duì)應(yīng)于最大特征值,確定主特征向量η。

4) 根據(jù)主特征向量計(jì)算加權(quán)系數(shù),完成圖像融合,得到最終灰度值向量。

(20)

在LabVIEW中PAC融合算法如圖10所示。

圖10 PCA圖像融合模塊

2.4 圖像顯示

本設(shè)計(jì)利用圖形-三維圖形-散點(diǎn)控件來(lái)顯示圖像,該控件輸入的數(shù)據(jù)格式主要為矩陣,通過(guò)連接橫、縱坐標(biāo)值和算法的結(jié)果進(jìn)行圖像顯示,如圖11所示。結(jié)合成像系統(tǒng)功能需求,該系統(tǒng)的操作面板設(shè)計(jì)如圖12所示。

圖11 圖像顯示控件

3 結(jié)果分析

電阻抗層析成像系統(tǒng)的實(shí)物如圖13所示。

圖13 電阻抗層析成像系統(tǒng)實(shí)物圖

實(shí)驗(yàn)中所使用的介質(zhì)分別為氣和水,共3組不同的介質(zhì)分布情況。將內(nèi)徑10 mm、管壁厚度1 mm、長(zhǎng)度300 mm的亞克力管裝滿水(上下端口封閉)放入傳感器中,模擬氣水兩相流。第1組:放置1根亞克力管,背景介質(zhì)為空氣;第2組:放置2根亞克力管,背景介質(zhì)為空氣;第3組:放置3根亞克力管,背景介質(zhì)為空氣。不同介質(zhì)分布如圖14所示。

圖14 不同介質(zhì)分布

對(duì)3組測(cè)量對(duì)象分別使用LabVIEW編輯的LBP、Landweber和ART算法實(shí)現(xiàn)圖像重構(gòu),并采用重構(gòu)圖像的相關(guān)系數(shù)Rcc對(duì)重構(gòu)精度進(jìn)行評(píng)價(jià)[23],結(jié)果如表1～6所示。表1、3、5為圖像重構(gòu)結(jié)果,表2、4、6為重構(gòu)圖像相關(guān)系數(shù)。

表1 單根棒圖像重建

表2 單根棒圖像重建相關(guān)系數(shù)

表3 2根棒圖像重建

表4 2根棒圖像重建相關(guān)系數(shù)

表6 3根棒圖像重建相關(guān)系數(shù)

從表1、3、5可以發(fā)現(xiàn),基于LabVIEW上位機(jī)軟件的NEIT系統(tǒng)可以有效實(shí)現(xiàn)數(shù)據(jù)的實(shí)時(shí)采集及介質(zhì)分布的圖像重構(gòu),并結(jié)合相關(guān)系數(shù)綜合分析不同算法及不同介質(zhì)分布下的重構(gòu)效果。根據(jù)重構(gòu)圖像及相關(guān)系數(shù)結(jié)果分析,在本文所搭建的實(shí)驗(yàn)條件下, ART算法在不同流型下,成像效果均不錯(cuò);Landweber算法比較適合兩根柱狀流,在其他兩種流型下,出現(xiàn)的偽影較多;LBP算法雖然計(jì)算方法簡(jiǎn)單,但誤差最大,成像效果不佳,不能正確地反映管道內(nèi)的介質(zhì)分布。此外,結(jié)合表2、4、6可以看出,無(wú)論是實(shí)部、虛部還是相位的單獨(dú)成像,其成像效果均受介質(zhì)分布或重構(gòu)算法影響較大。經(jīng)PCA算法融合后的圖像重構(gòu)效果更穩(wěn)定,且其相關(guān)系數(shù)普遍高于3個(gè)分量單獨(dú)成像相關(guān)系數(shù)的均值,因此,該融合算法的加入可以降低成像效果受介質(zhì)分布的影響,提高成像精度。

4 結(jié)語(yǔ)

本文提出一種基于LabVIEW的非接觸式電阻抗層析成像系統(tǒng),該系統(tǒng)的圖像重構(gòu)算法(LBP、Landweber、ART)全部由LabVIEW編譯,無(wú)需額外調(diào)用Matlab。此外,引入主成分分析融合算法,綜合利用輸出電信號(hào)的實(shí)部、虛部和相位的成像信息,提高成像精度。通過(guò)搭建的氣-水兩相流靜態(tài)實(shí)驗(yàn)平臺(tái)驗(yàn)證,該系統(tǒng)可以有效實(shí)現(xiàn)數(shù)據(jù)的實(shí)時(shí)采集及介質(zhì)分布的圖像重構(gòu)。進(jìn)一步結(jié)合相關(guān)系數(shù)分析,在成像過(guò)程中加入主成分分析圖像融合算法,有助于提高介質(zhì)分布重建的穩(wěn)定性及精度。該上位機(jī)軟件同樣可以為其他電學(xué)法多相流層析成像系統(tǒng)的數(shù)據(jù)實(shí)時(shí)采集及在線重構(gòu)提供借鑒。