王小鑫,張艷麗,白 萌,席 煒,黨 博,胡紅利
(1.西安石油大學(xué) 陜西省油氣井測控技術(shù)重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710065;2. 西安交通大學(xué) 電力設(shè)備電氣絕緣國家重點(diǎn)實(shí)驗(yàn)室, 陜西 西安 710049)
兩相流普遍存在于石油、化工及制藥等工業(yè)領(lǐng)域的開發(fā)和生產(chǎn)過程中。其流場分布參數(shù)演變復(fù)雜且在時(shí)間、空間上相互耦合,因此,研發(fā)能夠直觀顯示兩相流流動(dòng)狀態(tài)的測量技術(shù)至關(guān)重要[1-2]。電學(xué)層析成像(electrical tomography, ET)技術(shù)因其無輻射、適用性強(qiáng)、圖像成像快、檢測成本低等優(yōu)點(diǎn),在兩相流可視化測量中得到了廣泛的應(yīng)用。經(jīng)過多年的科學(xué)研究和實(shí)踐應(yīng)用,電學(xué)層析成像技術(shù)在傳感器、數(shù)據(jù)采集和圖像重建等方面取了大量研究成果[3-4]。其中,傳統(tǒng)ERT/EIT技術(shù)在測量過程中通常需要接觸待測介質(zhì),容易引起傳感器電極與被測流體之間發(fā)生物理和化學(xué)反應(yīng)[5]。電容層析成像(electrical capacitance tomography, ECT)技術(shù)則通過檢測場域中不同物質(zhì)對應(yīng)的介電常數(shù),重建出被測場域截面的介電常數(shù)分布,該技術(shù)主要適用于非導(dǎo)電物體組成的兩相流或者多相流檢測[6]。近年來,鑒于含導(dǎo)電相多相流的非接觸式檢測需求,Wang和Xu等人將電容耦合非接觸電導(dǎo)檢測(capacitively coupled contactless conductivity detection,簡稱C4D) 技術(shù)與ERT技術(shù)相結(jié)合,組成電容耦合電阻層析成像(CCERT)技術(shù)以實(shí)現(xiàn)非接觸式等效電阻測量[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)場的實(shí)時(shí)監(jiān)測及工作人員的操作。NI公司創(chuàng)立的LabVIEW是目前國際上廣泛應(yīng)用的虛擬儀器開發(fā)環(huán)境之一,已成為測量與測試領(lǐng)域的工業(yè)標(biāo)準(zhǔn),尤其在儀器控制、數(shù)據(jù)采集、數(shù)據(jù)分析及顯示等領(lǐng)域發(fā)揮了巨大優(yōu)勢。周英鋼等人研發(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)通過阻抗測量電路獲取輸出電信號的實(shí)部、虛部和相位進(jìn)行分別成像,并結(jié)合主成分分析融合算法以提高重構(gòu)精度。
非接觸式電阻抗層析成像系統(tǒng)主要由12電極傳感器、測量采集電路以及用于圖像重建的PC機(jī)組成。其系統(tǒng)結(jié)構(gòu)如圖1所示。
圖1 電阻抗層析成像系統(tǒng)結(jié)構(gòu)
圖1中測量及采集電路主要由激勵(lì)信號發(fā)生電路、I/U(阻抗/電壓)測量電路、開關(guān)控制電路、基于模擬乘法器的相敏解調(diào)電路、低通濾波器和模數(shù)轉(zhuǎn)換電路組成。測量及采集電路是整個(gè)系統(tǒng)的核心,連接陣列電極,控制每個(gè)電極的激勵(lì)、待測及接地狀態(tài),同時(shí),將得到阻抗測量值由串口傳輸?shù)絇C機(jī)完成可視化測量,在測量過程中通過PC端的上位機(jī)來進(jìn)行控制。
本系統(tǒng)采用12電極非接觸式陣列傳感器,其徑向截面結(jié)構(gòu)及單電極激勵(lì)的電勢分布如圖2所示。采用柔性電極板加工技術(shù)以保證其陣列電極的對稱性和均一性,電極片貼附在PVC管道周圍,并在外部包裹上接地金屬屏蔽罩,圖3為傳感器實(shí)物圖。其中,管道內(nèi)徑為46 mm,外徑為2 mm,柔性電極的張角為26°,電極的覆蓋率為86.6%,管道與屏蔽罩外殼之間的距離為8 mm。
圖2 電極結(jié)構(gòu)及電勢分布
圖3 12電極傳感器
電學(xué)層析成像系統(tǒng)的靈敏度矩陣一般由Geselowitz靈敏度理論計(jì)算得出,基于電勢分布法的靈敏度矩陣可表示為[14]
(1)
式中:Si,j(k)為傳感器電極i和電極j加激勵(lì)時(shí)場域內(nèi)第k個(gè)單元處的靈敏度;φi為電極i加激勵(lì)信號時(shí)的場內(nèi)電勢;φj為電極j加激勵(lì)信號時(shí)的場內(nèi)電勢;vk為第k個(gè)單元的體積;Ue是激勵(lì)電極的電壓。第k個(gè)單元格的電勢梯度φ(k)為
φ(k)=Exx+Eyy+Ezz
(2)
式中:x、y和z分別為x軸方向、y軸方向和z軸方向的單位矢量;Ex、Ey和Ez分別表示x軸方向、y軸方向和z軸方向的電場強(qiáng)度,本文借助Comsol有限元分析軟件獲得。進(jìn)一步在激勵(lì)電極和檢測電極的作用下,第k個(gè)單元的靈敏度矩陣可以表示為
(3)
本文所使用的靈敏度矩陣在空場(管道內(nèi)部為空氣)狀態(tài)下獲取,靈敏度矩陣作為解決層析成像逆問題數(shù)學(xué)建模中的權(quán)重?cái)?shù)矩陣,反映某一單元的介質(zhì)發(fā)生變化時(shí)引起的阻抗值的變化量。圖4為相鄰電極和相對電極加激勵(lì)時(shí)的敏感場分布圖。
圖4 靈敏場分布
本系統(tǒng)測量電路的激勵(lì)源采用MAX038模擬信號發(fā)生器產(chǎn)生固定頻率為2 MHz的正弦交流信號。為了實(shí)現(xiàn)電阻抗實(shí)部與虛部的測量,引入了模擬乘法器相敏解調(diào)電路。模擬乘法器采用AD734構(gòu)成相敏解調(diào)電路,2個(gè)乘法器分別實(shí)現(xiàn)同向解調(diào)和正交解調(diào),如圖5所示。
圖5 實(shí)部、虛部測量結(jié)構(gòu)圖
正弦信號一路經(jīng)過I/U轉(zhuǎn)換電路調(diào)制得到包含阻抗信息的電壓信號,一路通過移相器將正弦信號的相位超前或滯后90°,得到正交信號,并將此信號作為乘法器的參考信號,用于正交解調(diào);另一路正弦信號直接與乘法器相連,作為正交解調(diào)的參考信號,用于同向解調(diào);最后經(jīng)過乘法解調(diào)后的信號經(jīng)過低通濾波器,分別得到包含阻抗的實(shí)部信息和虛部信息的直流電壓。此外,相位信息可由實(shí)部和虛部值進(jìn)一步計(jì)算得出。
設(shè)激勵(lì)信號Vi為
Vi=Asin(ωt+α)
(4)
式中:A為振幅;ω為角頻率;α為初始相角。
則I/U轉(zhuǎn)換電路輸出的時(shí)諧信號Vo為
AωRfCxcos(ωt+α)
(5)
式中:Rx為電極間有損介質(zhì)的等效電阻;Rf為放大電阻的反饋電阻;Cx為電極間有損介質(zhì)的等效電容。
Vo和Vi經(jīng)過第一個(gè)乘法器相乘可以得到
A2ωRfCxsin(ωt+α)cos(ωt+α)
(6)
通過低通濾波器后輸出的直流信號為
(7)
同理,將Vo和Vi經(jīng)過移相90°得到的正交信號相乘,并且經(jīng)過濾波器可以得到
(8)
得到相位的計(jì)算式為
(9)
圖像重建屬于逆問題,即在已知結(jié)果的情況下,反推產(chǎn)生它的原因。阻抗層析的逆問題是根據(jù)測得的阻抗信息來反演成像區(qū)域內(nèi)損耗因子的分布,求解逆問題的實(shí)質(zhì)是求解靈敏度矩陣S的廣義逆S-1。類似EIT的歸一化線性模型可以表示為[15]
G=S-1Z
(10)
式中:G為介質(zhì)分布像素列向量(N×1);S-1為歸一化的靈敏度矩陣的廣義逆矩陣(N×M);Z為歸一化的阻抗向量(M×1)。通常情況下N?M,因此該問題的解不唯一,存在欠定性,圖像重構(gòu)技術(shù)可以用于解決EIT的逆問題。
本系統(tǒng)圖像重建部分在LabVIEW軟件上實(shí)現(xiàn),該軟件由NI公司研究開發(fā),其優(yōu)勢在于LabVIEW使用的是圖形化編輯語言G編寫程序,以框圖的樣式完成功能實(shí)現(xiàn)[16]。本文在LabVIEW上編輯3種經(jīng)典圖像重構(gòu)算法,分別是LBP(linear back projection)算法、Landweber迭代算法、ART(algebraic reconstruction technique)代數(shù)重建算法。
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為測量的阻抗向量。
在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所示。該算法主要通過for循環(huán)和數(shù)學(xué)線性代數(shù)矩陣A×B控件實(shí)現(xiàn),為了實(shí)現(xiàn)迭代計(jì)算,需要通過移位寄存器將上一次的計(jì)算結(jié)果轉(zhuǎn)變?yōu)橄乱淮蔚某跏贾?。迭代次?shù)不同,則圖像重構(gòu)的顯示結(jié)果也不同。同時(shí),在Landweber圖像重構(gòu)算法中,計(jì)算迭代步長α,計(jì)算過程見圖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對投影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算法主要通過2個(gè)for循環(huán)實(shí)現(xiàn),如圖9所示。
圖9 ART算法模塊
為了進(jìn)一步提高圖像精度,本文基于主成分分析方法(PCA)對實(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的特征值和對應(yīng)的特征向量。選擇主特征向量對應(yīng)于最大特征值,確定主特征向量η。
4) 根據(jù)主特征向量計(jì)算加權(quán)系數(shù),完成圖像融合,得到最終灰度值向量。
(20)
在LabVIEW中PAC融合算法如圖10所示。
圖10 PCA圖像融合模塊
本設(shè)計(jì)利用圖形-三維圖形-散點(diǎn)控件來顯示圖像,該控件輸入的數(shù)據(jù)格式主要為矩陣,通過連接橫、縱坐標(biāo)值和算法的結(jié)果進(jìn)行圖像顯示,如圖11所示。結(jié)合成像系統(tǒng)功能需求,該系統(tǒng)的操作面板設(shè)計(jì)如圖12所示。
圖11 圖像顯示控件
電阻抗層析成像系統(tǒng)的實(shí)物如圖13所示。
圖13 電阻抗層析成像系統(tǒng)實(shí)物圖
實(shí)驗(yàn)中所使用的介質(zhì)分別為氣和水,共3組不同的介質(zhì)分布情況。將內(nèi)徑10 mm、管壁厚度1 mm、長度300 mm的亞克力管裝滿水(上下端口封閉)放入傳感器中,模擬氣水兩相流。第1組:放置1根亞克力管,背景介質(zhì)為空氣;第2組:放置2根亞克力管,背景介質(zhì)為空氣;第3組:放置3根亞克力管,背景介質(zhì)為空氣。不同介質(zhì)分布如圖14所示。
圖14 不同介質(zhì)分布
對3組測量對象分別使用LabVIEW編輯的LBP、Landweber和ART算法實(shí)現(xiàn)圖像重構(gòu),并采用重構(gòu)圖像的相關(guān)系數(shù)Rcc對重構(gòu)精度進(jìn)行評價(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ì)算方法簡單,但誤差最大,成像效果不佳,不能正確地反映管道內(nèi)的介質(zhì)分布。此外,結(jié)合表2、4、6可以看出,無論是實(shí)部、虛部還是相位的單獨(dú)成像,其成像效果均受介質(zhì)分布或重構(gòu)算法影響較大。經(jīng)PCA算法融合后的圖像重構(gòu)效果更穩(wěn)定,且其相關(guān)系數(shù)普遍高于3個(gè)分量單獨(dú)成像相關(guān)系數(shù)的均值,因此,該融合算法的加入可以降低成像效果受介質(zhì)分布的影響,提高成像精度。
本文提出一種基于LabVIEW的非接觸式電阻抗層析成像系統(tǒng),該系統(tǒng)的圖像重構(gòu)算法(LBP、Landweber、ART)全部由LabVIEW編譯,無需額外調(diào)用Matlab。此外,引入主成分分析融合算法,綜合利用輸出電信號的實(shí)部、虛部和相位的成像信息,提高成像精度。通過搭建的氣-水兩相流靜態(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ù)分析,在成像過程中加入主成分分析圖像融合算法,有助于提高介質(zhì)分布重建的穩(wěn)定性及精度。該上位機(jī)軟件同樣可以為其他電學(xué)法多相流層析成像系統(tǒng)的數(shù)據(jù)實(shí)時(shí)采集及在線重構(gòu)提供借鑒。