張立峰, 王 智

(華北電力大學(xué) 自動化系，河北 保定 071003)

1 引 言

氣液兩相流系統(tǒng)作為一個高度隨機(jī)、非線性的瞬時動態(tài)系統(tǒng),廣泛存在于自然環(huán)境與工業(yè)生產(chǎn)中,如火力發(fā)電、核電、油氣運輸、航空航天等領(lǐng)域[1,2]。兩相界面不穩(wěn)定的熱量、動量及質(zhì)量傳遞過程使得流動機(jī)理復(fù)雜,難以用精確的數(shù)學(xué)模型描述。流型是兩相流的流態(tài)特征之一,準(zhǔn)確實現(xiàn)流型識別對定量塑造兩相流流動、深入理解流動機(jī)理及工業(yè)設(shè)施的安全生產(chǎn)具有重要意義,然而目前尚未見合適的流動模型,對工藝參數(shù)的處理方法也不完善,因此,急需探索一種魯棒性強(qiáng)、準(zhǔn)確率高的流型識別框架。 流型測量方法大致可分為射線法[3]、光學(xué)法[4]、超聲波法[5]、差壓法[6]、電學(xué)法等,其中電學(xué)法可以獲得管道截面的相分布狀態(tài),且具有成本低、靈活性高等優(yōu)點,被認(rèn)為是一種有前途的多相流檢測技術(shù)。電學(xué)法主要包括電阻層析成像[7,27](electrical resistance tomography, ERT)、電容層析成像[8](electrical capacitance tomography, ECT)及電磁層析成像[9](electromagnetic tomography, EMT),其中ERT在氣液兩相流檢測中得到了廣泛應(yīng)用,為研究流體行為提供了一種非侵入、非放射性的檢測方法,有助于研究不同尺度與極端條件下的動態(tài)流動行為。ERT系統(tǒng)通過多點測量提供了較為全面的流動信息,使用特定的算法可以實現(xiàn)截面電導(dǎo)率分布的實時重建[10]。若能利用ERT系統(tǒng)獲取流動數(shù)據(jù),進(jìn)而即可通過數(shù)據(jù)處理、特征提取與分類模型實現(xiàn)流型識別。

對于傳感器測量值序列的處理方法有功率譜分析[11]、概率密度分析[12]、多尺度熵[13]、復(fù)雜網(wǎng)絡(luò)[14]、時頻分析[15]、統(tǒng)計分析[16]等,但其中大部分方法僅僅從時域或頻域單空間中認(rèn)識測量數(shù)據(jù),盡管時頻分析法兼顧了時域與頻域空間,但其忽視了ERT系統(tǒng)可以重建截面相分布的優(yōu)點,因此,為充分利用ERT系統(tǒng)測量信息,本文將重建圖像作為測量數(shù)據(jù)的補(bǔ)充,在多域空間中處理ERT系統(tǒng)采集的流動信息。ERT系統(tǒng)測量值序列蘊(yùn)含著介質(zhì)分布變化的時域信息,重建圖像提供了介質(zhì)分布在空間域的直觀表達(dá),而變換域有助于提取更加全面的圖像信息。離散傅里葉變換(discrete fourier transform, DFT)常被用于實現(xiàn)圖像由空域到頻域的變換,但DFT在計算中涉及到了復(fù)數(shù)運算,作為頻域的擴(kuò)展,Walsh變換[17]為在列率域中研究頻率提供了理論依據(jù),且具有計算過程簡單、能量集中的優(yōu)點。然而,提取的高維特征往往包含著冗余信息,采用合適的方法實現(xiàn)特征簡約對提高計算效率、提升識別準(zhǔn)確率有重要作用。

流形學(xué)習(xí)[18]方法已成功用于尋找高維空間中數(shù)據(jù)流形結(jié)構(gòu)的低維嵌入映射,從而實現(xiàn)維數(shù)簡約,廣泛用于生物醫(yī)學(xué)[19]、故障診斷[20]、圖像處理[21]等領(lǐng)域。流形學(xué)習(xí)方法可分為線性與非線性兩類,線性方法如獨立成分分析 (independent component analysis, ICA)、主成分分析 (principal components analysis, PCA)等,非線性方法有局部線性嵌入[22](locally Linear Embedding, LLE)、等度量映射 (isometric Mapping, Isomap)、統(tǒng)一流形逼近與投影[23](uniform manifold approximation and projection, UMAP)等。其中,UMAP基于k-近鄰思想,通過計算成對點在高維空間與低維映射空間中的距離度量并使用隨機(jī)梯度下降法(stochastic gradient descent, SGD)優(yōu)化結(jié)果實現(xiàn)降維,具有保留全局結(jié)構(gòu)、計算快速及非線性降維等優(yōu)點。

為實現(xiàn)垂直上升管道氣液兩相流流型識別,本文提出了一種多域特征的流型識別方案:1)在內(nèi)徑50 mm的垂直上升管道上進(jìn)行動態(tài)試驗操作,并利用ERT系統(tǒng)監(jiān)測兩相流流動以獲取氣液兩相介質(zhì)分布變化信息;2)為融合ERT傳感器測量數(shù)據(jù)與ERT重建圖像蘊(yùn)含的流動分布信息,提出了一種多域特征處理方案,從時域、空域及列率域提取并量化流型特征。使用Walsh變換代替時頻變換將圖像由空域轉(zhuǎn)換至列率域簡化計算過程;3)為實現(xiàn)特征簡約,使用UMAP算法對特征進(jìn)行降維處理,以獲得低維空間中具有代表性的特征向量,并與ICA、LLE及Isomap降維方法進(jìn)行對比,驗證UMAP的有效性;4)搭建支持向量機(jī)(support vector machine, SVM)實現(xiàn)流型辨識,并進(jìn)一步研究低維映射空間維度對辨識性能的影響。對上述方案進(jìn)行實驗驗證,結(jié)果表明,該流型識別框架具有較高的穩(wěn)定性與識別效果。

2 ERT系統(tǒng)原理與兩相流實驗

2.1 電阻層析成像

電阻層析成像基于不同介質(zhì)(相)電導(dǎo)率不同的原理。通過1對或多對電極向被測區(qū)域施加激勵以形成敏感場,當(dāng)場內(nèi)介質(zhì)分布改變時,電導(dǎo)率分布發(fā)生變化從而引起電勢分布變化,這種改變可以通過測量邊界電極輸出所反應(yīng)。本文所用16電極ERT傳感器結(jié)構(gòu)如圖1所示,管道內(nèi)徑50 mm,圓形電極直徑5 mm,均勻排列在管壁上。

圖1 ERT傳感器結(jié)構(gòu)Fig.1 Structure of ERT sensor

ERT系統(tǒng)工作模式如圖2所示,采用電壓激勵電流檢測的單電極激勵、測量傳感策略。激勵信號為頻率256 kHz、峰-峰值4 V的正弦交流電壓。ERT系統(tǒng)工作時,首先向電極E1施加激勵,依次測量其余電極輸出,然后將激勵施加到電極E2,仍依次測量其余電極輸出,重復(fù)上述過程直到完成所有測量。1次完整掃描可以獲得16×15=240個測量數(shù)據(jù),ERT系統(tǒng)數(shù)據(jù)采集速度為120 幀/s。

ERT圖像重建是由測量值到電導(dǎo)率分布的反演過程,可以描述為式(1),其中I表示測量值,W為靈敏度矩陣,反應(yīng)電導(dǎo)率分布g與邊界測量值的映射關(guān)系。

I=Wg

(1)

線性反投影[24](linear back projection, LBP) 算法由于其快速性已廣泛用于電阻層析成像的在線實時重建,可以表達(dá)為式(2),

(2)

式中uI中所有元素為1。

2.2 兩相流實驗

在華北電力大學(xué)成像與智能感測實驗室的氣液兩相流實驗裝置上進(jìn)行實驗,實驗裝置示意圖如圖3所示。氣、水輸送管道由不銹鋼管與透明有機(jī)玻璃管組成,其中玻璃管安裝在透明直管段,用于觀測流型,記錄流動圖像。測試段安裝16電極ERT傳感器并連接ERT系統(tǒng),通過計算機(jī)保存數(shù)據(jù)。

圖3 實驗裝置Fig.3 Experimental device

實驗介質(zhì)為空氣與水?？諝馔ㄟ^空氣壓縮機(jī)壓縮后存儲在儲氣罐內(nèi),并依靠穩(wěn)壓閥將氣體出口壓強(qiáng)固定在0.5 MPa。存儲在水箱內(nèi)的水由水泵抽出,與氣體混合后繼續(xù)向上運動,在測試段完成數(shù)據(jù)采集后返回水箱構(gòu)成水循環(huán)回路,同時氣體排放至環(huán)境中。分別使用電磁流量計與浮子流量計監(jiān)測氣液兩相流量,準(zhǔn)確度為±0.5%與±1.5%。

參考Taitel[25]流型圖改變氣液流量,同時通過透明直管段觀測流動圖像,利用流動結(jié)構(gòu)的形態(tài)學(xué)特征確定流型。氣體以小氣泡形式離散的分布在液相中為泡狀流;隨著氣相流速的增加,氣泡間的碰撞行為加劇,小氣泡聚合成直徑稍大的氣泡,偶爾出現(xiàn)直徑比管徑略小的“Taylor”氣泡,系統(tǒng)由泡狀流向段塞流過渡;進(jìn)一步增大氣相流速,氣泡以近于管道內(nèi)徑大小的大型彈狀氣塞形式出現(xiàn),氣塞與液塞交替發(fā)生并具有一定的偽周期性趨勢,流動系統(tǒng)發(fā)展為段塞流;當(dāng)氣塞表面液膜變得不穩(wěn)定時,長氣塞逐漸破碎,流動系統(tǒng)展現(xiàn)出了一定的震蕩特性,發(fā)展為嚴(yán)重段塞流。

使用ERT系統(tǒng)采集123種不同工況點流動數(shù)據(jù),其中泡狀流13組,泡狀-段塞過渡流型20組,段塞流42組,嚴(yán)重段塞流48組,每種工況點采集10次,每次采集10 s,流動圖像分別如圖4(a)～圖4(d)所示。

圖4 氣液兩相流流型Fig.4 Flow pattern of gas-liquid two-phase flow

為獲得流動中兩相分布的實時變化情況,依據(jù)式(2)使用LBP算法反演測量數(shù)據(jù),圖5展示了試驗中基于該算法的ERT系統(tǒng)重建圖像,圖中像素灰度值較高部分代表氣相,灰度值較低部分表示液相。

圖5 LBP算法重建圖像Fig.5 Image reconstruction based on LBP algorithm

3 流型識別框架

基于ERT系統(tǒng)采集的流動數(shù)據(jù),采用多域特征提取的流型識別方案,如圖6所示。該方案從測量值序列與LBP重建圖像兩個角度分析,由時域、空域、列率域提取量化特征,將數(shù)據(jù)集劃分為訓(xùn)練集與測試集兩類,為實現(xiàn)特征簡約,使用UMAP算法對訓(xùn)練集高維特征做降維處理以獲取低維流形嵌入空間,并使用訓(xùn)練集降維數(shù)據(jù)訓(xùn)練SVM分類器,最后使用測試集數(shù)據(jù)驗證分類效果。

圖6 流型識別框架Fig.6 Structure of flow pattern recognition

3.1 Walsh-Hadamard變換

Walsh變換[17]可以將函數(shù)映射到正交、歸一化的Walsh空間,實現(xiàn)了圖像由空域到列率域的轉(zhuǎn)換。不同于Fourier變換中頻率的概念,列率是一種廣義的頻率,其代表信號單位時間內(nèi)跨越零點次數(shù)的1/2。Walsh函數(shù)系的取值僅有-1與+1,由于其運算只涉及實數(shù)的加減而沒有復(fù)數(shù)的乘法,因此相對于Fourier變換更為快速簡單。

對于二維離散圖像信號G(x,y),x∈[0,M-1],y∈[0,N-1],其Walsh變換為:

(3)

式中:X為圖像G的列率譜;u∈[0,M-1];v∈[0,N-1];φ(x,y,u,v)為Walsh基函數(shù)。

當(dāng)M=2m,N=2n時,可采用具有簡單遞推關(guān)系的Hadamard排列作為基函數(shù)。即

(4)

式中bi(x) 為x倒序自然二進(jìn)制碼的第i位,僅可取0或1。

對于2k階Hadamard排列可以表示為:

(5)

式中?表示克羅內(nèi)克積。

3.2 特征量化

對ERT系統(tǒng)采集的高維測量值序列降維處理,將每1幀數(shù)據(jù)壓縮至一維以簡化分析:

(6)

式中:Iij為氣液兩相流時第i次激勵下第j次測量的電流值;Iij0為Iij在管道內(nèi)充滿水時的測量數(shù)據(jù);則I可表示為1次完整掃描下測量電極對氣液兩相介質(zhì)分布變化的平均響應(yīng)。

進(jìn)一步提取時間序列時域特征并量化。使用均值tμ表征數(shù)據(jù)集中趨勢,使用標(biāo)準(zhǔn)差tσ表征測量數(shù)據(jù)離散程度,使用偏度tskew表征數(shù)據(jù)分布的對稱性,使用峰度tkurt表征數(shù)據(jù)分布的平滑、陡峭程度。

(7)

(8)

(9)

(10)

式中:T為時間序列長度;I(i)為第i次激勵的電流值。

(11)

式中Gi為LBP算法重建出的第i幅圖像。

提取紋理均值Sμ表示紋理分布的全局均勻程度,使用紋理方差Sσ表示紋理分布局部變化的劇烈程度,使用圖能量SE表示紋理分布的粗粒程度。

(12)

(13)

(14)

式中g(shù)為像素點個數(shù)。

為更全面地反應(yīng)圖像的紋理特征,進(jìn)一步將其轉(zhuǎn)化至列率域。由于LBP算法得到的圖像為僅有812個像素點的不規(guī)則圓形,因此在進(jìn)行Walsh-Hadamard變換之前需對其補(bǔ)0填充至32×32的形狀,記為G′。提取列率譜系數(shù)矩陣X′的一階矩Cμ、L1范數(shù)CL1與L2范數(shù)CL2描述紋理全局特征。

X′=Wal(G′)

(15)

(16)

(17)

(18)

式中:Wal表示W(wǎng)alsh-Hadamard變換;j=1,2,…,32;λ為X′TX′的最大特征值。綜上則可得到每個樣本的特征向量F=[tμ,tσ,tskew,tkurt,Sμ,Sσ,SE,Cμ,CL1,CL2]。

3.3 UMAP

UMAP算法[23]基于黎曼幾何與代數(shù)拓?fù)淅碚?常用于高維特征的降維與可視化,其能在降維過程中最大程度保留數(shù)據(jù)原始特征,且運算快速。下面為該方法的主要步驟描述。

給定RM維空間中含有n個樣本的向量集X={xi},i=1,2,…,n,尋找每個樣本點xi的k個近鄰點集合ηi。對于每個點xi,尋找其最近鄰點與距離:

ρi=min(d(xi,xj)|xj∈ηi,d(xi,xj)>0)

(19)

(20)

式中:ρi為xi與第1個最近鄰點的距離;d(xi,xj) 為點xi與xj之間的距離;σi為xi最近鄰點的半徑。進(jìn)一步便可得到無向加權(quán)圖G,其鄰接矩陣B定義為:

B=A+AT-A°AT

(21)

式中°表示哈達(dá)瑪積。A中元素由相應(yīng)的有向圖的權(quán)重給出:

(22)

式中w((xi,xj)) 為點xi與xj的連接權(quán)重。

UMAP在低維空間中使用力定向圖布局算法,通過在迭代中向邊或頂點上施加作用力與反作用力實現(xiàn),與模擬退火算法中所使用的方法相似。設(shè)yi與yj是高維空間中點xi與xj的低維表示,則yi與xj間的作用力定義為:

(23)

式中a、b為超參數(shù)。

排斥力定義為:

(1-w((xi-xj)))×(yi-yj)

(24)

式中ε為一個極小數(shù),默認(rèn)為0.001。

3.4 支持向量機(jī)

SVM[26]在解決小樣本、非線性問題方面具有優(yōu)勢,其思想為在數(shù)據(jù)空間中構(gòu)造最優(yōu)超平面實現(xiàn)分類。采用成對分類法構(gòu)建多分類模型,對于含有L個類別的數(shù)據(jù)集,該策略在每兩個類別之間構(gòu)建一個二分類SVM進(jìn)而實現(xiàn)多分類,則共需訓(xùn)練L(L-1)/2個二分類器。

給定含有n個樣本的訓(xùn)練集X={(x1,y1),(x2,y2),…,(xn,yn)},xi∈Rk,i=1,2,…,n,k為特征空間維度,yn={1,2,…,L},表示樣本類別。

對于一個二分類模型,使用高斯核函數(shù)時判決函數(shù)為:

(25)

(26)

(27)

式中:αi為拉格朗日乘數(shù);b為閾值;σ為可選參數(shù)。

4 結(jié)果分析

將本文所提流型識別模型應(yīng)用于辨識垂直上升管道的氣液兩相流流型,將樣本特征歸一化至[0,1]之間作為模型輸入,模型輸出為代表流型的標(biāo)簽值,即泡狀流、泡狀-段塞過渡流型、段塞流及嚴(yán)重段塞流的標(biāo)簽值分別為1、2、3、4。設(shè)置訓(xùn)練集與測試集比例為6∶4。使用分類準(zhǔn)確率定量評價模型識別效果:

(28)

式中:SACC,i為類別i的分類準(zhǔn)確率;TP,i為類別i分類正確的樣本個數(shù);Ti為類別i的測試集樣本總數(shù)。

為驗證UMAP算法有效性,使用LLE、Isomap、ICA作為對比,同時研究了降維維度對訓(xùn)練集分類準(zhǔn)確的影響,如表1所示。此時LLE算法鄰居數(shù)量設(shè)置為15,UMAP與Isomap算法的鄰居數(shù)量均為5。

表1 不同分類模型在訓(xùn)練集上的分類準(zhǔn)確率Tab.1 Classification accuracy of different classification models on the training set (%)

由表1可以發(fā)現(xiàn)UMAP-SVM的分類性能在訓(xùn)練集上展現(xiàn)出最優(yōu)效果,且降維空間的維度對分類結(jié)果影響較小。觀察其他分類模型可以發(fā)現(xiàn)保留的流形空間維度并不是越高越好,過高的維度可能會帶來無效特征,這不僅不利于降低模型計算量,甚至?xí)?dǎo)致分類準(zhǔn)確率略有下降。為確定UMAP算法的最佳近鄰數(shù)量,并保證其有一定的魯棒性,將流形空間設(shè)置為2維,通過改變近鄰點數(shù)量觀察UMAP-SVM分類器在測試集與訓(xùn)練集上的分類效果,結(jié)果如圖7所示。

由圖7可以發(fā)現(xiàn),鄰居數(shù)量對UMAP-SVM分類器在訓(xùn)練集上的分類性能影響較小,但會對模型泛化能力產(chǎn)生較大影響,在鄰居數(shù)量為7時模型魯棒性達(dá)到最佳效果,此時模型對測試集數(shù)據(jù)集的分類結(jié)果如表2所示,混淆矩陣如圖8所示。

表2 UMAP-SVM框架測試集分類準(zhǔn)確率Tab.2 Classification accuracy of UMAP-SVM framework on test set

圖8 UMAP-SVM辨識結(jié)果Fig.8 Identification results of UMAP-SVM

由表2可以發(fā)現(xiàn)UMAP-SVM流型識別模型在每種流型上的分類準(zhǔn)確率均可達(dá)到94.8%以上,平均分類準(zhǔn)確率可達(dá)95.7%,擁有較好的辨識性能,且具有較高的穩(wěn)定性。觀察圖8可以發(fā)現(xiàn),由于間歇流系統(tǒng)穩(wěn)定性較差,湍流動能相對泡狀流更強(qiáng),參數(shù)變化范圍更廣,易于其他流型產(chǎn)生交叉重疊,因此對間歇流的識別準(zhǔn)確率相對于泡狀流要低,特別在嚴(yán)重段塞流時,液膜不穩(wěn)定使得氣塞破碎現(xiàn)象相對明顯,系統(tǒng)具有一定的震蕩特性,使得測量參數(shù)與其他流型的重疊范圍更廣,因此其余流型樣本更易錯分為嚴(yán)重段塞流,這種現(xiàn)象可以通過增加訓(xùn)練集樣本數(shù)量,建立更加完備的樣本集來減弱,同時使用合適的去噪方法降低測量值序列噪聲也對提高模型辨識性能有利。

5 結(jié) 論

準(zhǔn)確識別氣液兩相流流型對研究流動機(jī)理及優(yōu)化產(chǎn)能設(shè)備具有重要意義。在垂直上升管道氣液兩相流設(shè)備上進(jìn)行動態(tài)試驗操作,使用ERT系統(tǒng)獲取流動中兩相介質(zhì)變化信息,進(jìn)而提出了一種基于多域特征的流型識別方案,該方案的主要貢獻(xiàn)如下:

1)為充分利用ERT系統(tǒng)采集的流型信息,提出了一種兼顧時域、空域與列率域的特征提取方案。使用列率域代替頻域簡化了運算過程,提高了計算效率。

2)使用UMAP算法對高維特征降維處理,得到了具有代表性的低維特征向量,為分類器高準(zhǔn)確率別提供了有效特征。

3)搭建SVM分類器實現(xiàn)流型識別,研究了UMAP-SVM模型的最佳參數(shù),結(jié)果表明,當(dāng)流形空間為2維,最近鄰點數(shù)量為7時,該模型性能最佳,對測試集的平均分類準(zhǔn)確率可達(dá)95.7%。