符號動(dòng)態(tài)濾波的氣液兩相流流型識別

宋情洋,劉鐵軍,陳宜倩

(中國計(jì)量大學(xué) 計(jì)量測試工程學(xué)院,浙江 杭州 310018)

近些年,許多國內(nèi)外專家在氣液兩相流流動(dòng)規(guī)律研究領(lǐng)域取得豐碩的成果.特別是在流態(tài)圖轉(zhuǎn)變和流型識別方面.隨著計(jì)算機(jī)控制理論的發(fā)展,神經(jīng)網(wǎng)絡(luò)因其操作簡單,能很好地處理相關(guān)問題,成為流型檢測的一個(gè)非常有用的方法.王妍芃[1]等將BP神經(jīng)網(wǎng)絡(luò)應(yīng)用于氣液兩相流流型的識別中,BP神經(jīng)網(wǎng)絡(luò)的缺點(diǎn)在于其對樣本的處理過程非常繁瑣并且十分不靈活.陳露陽[2]等人以渦街流量計(jì)為元件誘發(fā)鈍體繞流,通過管壁差壓法獲取氣液兩相流鈍體繞流的尾跡波動(dòng)信號,再采用集總經(jīng)驗(yàn)?zāi)B(tài)分解對信號進(jìn)行分解,利用Hilbert變換得到Hilbert邊際譜,并利用最大互相關(guān)系數(shù)法對固有模態(tài)函數(shù)進(jìn)行篩選,分別構(gòu)造出了與體積含氣率以及兩相雷諾數(shù)相關(guān)的兩種流型圖,得到了較高流型識別率.D-Markov 模型是基于馬爾科夫過程,用于復(fù)雜系統(tǒng)的模式識別的一種快速時(shí)間序列方法,最早由Ray在2004年提出.Chin、Ray等將該模型應(yīng)用于異常檢測領(lǐng)域中,并與傳統(tǒng)的統(tǒng)計(jì)方法和神經(jīng)網(wǎng)絡(luò)方法進(jìn)行對比分析,發(fā)現(xiàn)D-Markov模型的性能比上述兩種方法要好得多[3].Khatkhate、Gupta、Ray等將D-Markov模型應(yīng)用于聯(lián)軸器的異常檢測并取得了較好的結(jié)果[4].周迅、向馗等用D-Markov模型識別發(fā)動(dòng)機(jī)零部件的疲勞裂紋的發(fā)展模式,并且發(fā)現(xiàn)這個(gè)模型對模式轉(zhuǎn)變具有預(yù)測功能[5].

本文提出了一種基于符號動(dòng)態(tài)濾波的氣液兩相流流型識別方法,這是一種改進(jìn)數(shù)據(jù)序列系統(tǒng)處理信號的新方式.在對實(shí)數(shù)序列進(jìn)行符號化的過程中,需要除去一些無關(guān)的特征,保留能表征流型信息的一些特征,基于這些有用特征建立符號序列模型.通過符號序列分析方法將數(shù)據(jù)序列轉(zhuǎn)變?yōu)榉栃蛄泻?計(jì)算效率有了明顯的提高.對符號數(shù)據(jù)分析來說,噪聲的影響并不是很大,因此其魯棒性也得到了加強(qiáng).該方法是基于超聲波氣液兩相流參數(shù)檢測原理,采用一收一發(fā)式超聲換能器,兩相流體在管道流動(dòng)過程中不斷采集換能器接收端輸出信號強(qiáng)度大??；不同流型情況下?lián)Q能器輸出信號強(qiáng)度不同,對不同流型下實(shí)時(shí)的信號強(qiáng)度進(jìn)行統(tǒng)計(jì),分析得出不同流型下信號強(qiáng)度的差異.將差異數(shù)據(jù)運(yùn)用符號動(dòng)態(tài)濾波方法處理得出流型與信號差異信號的對應(yīng)關(guān)系；進(jìn)而為現(xiàn)有的氣液兩相流流行識別方法提供參考[6-7].

1 符號動(dòng)態(tài)濾波方法原理

聲波在介質(zhì)傳播過程中會伴隨著能量的傳遞,當(dāng)超聲波在氣液兩相流傳播中遇到氣相和液相分相界面會產(chǎn)生聲波的反射、散射和衍射，使換能器接收到超聲波幅值、相位、頻率產(chǎn)生變化.本文通過符號動(dòng)態(tài)濾波算法分析流經(jīng)氣液兩相流后被超聲波換能器接收信號幅值衰減情況，從而得到氣液兩相流的流型信息.對超聲信號處理采用符號動(dòng)態(tài)濾波法,其計(jì)算過程可以分為解析信號相空間分割和D-Markov模型兩大部分[8].

1.1 解析信號相空間分割

解析信號相空間分割是符號動(dòng)態(tài)濾波算法的一個(gè)非常重要的環(huán)節(jié),它主要分為3個(gè)部分,分別為Hilbert變換[9-10]、映射和分割,最終將一維數(shù)據(jù)序列轉(zhuǎn)換為符號序列.整個(gè)流程如圖1.

圖1 解析信號相空間分割流程圖Figure 1 Flow chart of phase separation of the analytical signal

設(shè)符號集為S={s1,s2,…,sk},那么該符號集的符號熵定義為

(1)

式(1)中:Pi表示輸出符號Si的概率.

為了表述符號熵隨符號集大小k的變化,引入一個(gè)變量h(k),定義為

h(k)=H(k+1)-H(k).

(2)

取閾值ε,當(dāng)h(k)<ε時(shí),可得到符號集大小k的取值.

1.2 概率有限狀態(tài)機(jī)

經(jīng)過1.1描述的步驟后,實(shí)現(xiàn)了從原始的時(shí)間序列到符號序列的轉(zhuǎn)換,最后獲得一連串的符號序列.對于得到的符號序列,需要將其輸入到D-Markov模型中進(jìn)行計(jì)算輸出狀態(tài)轉(zhuǎn)移概率矩陣,提取其左特征向量作為狀態(tài)轉(zhuǎn)移概率向量,從而對矩陣狀態(tài)轉(zhuǎn)移概率矩陣的相似性進(jìn)行度量.符號處理的流程圖如圖2.

圖2 符號序列處理流程圖Figure 2 Symbol sequence processing flow chart

D-Markov模型最早運(yùn)用于工程領(lǐng)域的一些檢測中.對于給定的一串符號序列,構(gòu)造一個(gè)長度為D的窗口,從最初的符號開始,將長度為D的窗口向左滑動(dòng),每次滑動(dòng)窗口產(chǎn)生一個(gè)新的符號,將這個(gè)新產(chǎn)生的符號加上前一個(gè)狀態(tài)的D-1個(gè)符號,構(gòu)成了一個(gè)新的狀態(tài).將原狀態(tài)的符號串的個(gè)數(shù)表示為N(sn-1sn-2…sn-D),將原狀態(tài)加上新產(chǎn)生的符號后形成的新的符號串的個(gè)數(shù)表示為N(snsn-1sn-2…sn-D).

2 氣液兩相流流型識別實(shí)驗(yàn)裝置設(shè)計(jì)

2.1 氣液兩相流流型識別裝置設(shè)計(jì)

超聲氣液兩相流流型識別實(shí)驗(yàn)系統(tǒng)流程圖如圖3,主要由信號發(fā)生模塊、氣/水源模塊、測量模塊和數(shù)據(jù)采集模塊四個(gè)模塊組成.其中,信號發(fā)生模塊包括信號發(fā)生器、功率放大器,氣/水源模塊包括氣泵、水箱、調(diào)節(jié)閥和流量計(jì),測量模塊包括透明管道和一發(fā)一收式的超聲換能器,數(shù)據(jù)采集模塊包括示波器、數(shù)據(jù)采集卡和PC上位機(jī).信號發(fā)生器產(chǎn)生一個(gè)激勵(lì)信號,經(jīng)過功率放大器放大后傳到超聲換能器的發(fā)射端,由此生成超聲信號注入到垂直管道中.同時(shí),打開調(diào)節(jié)閥,根據(jù)調(diào)節(jié)閥的開度調(diào)整氣泵和水箱分別注入到垂直管道的氣液兩相流量中.在垂直管道的另一端,由超聲換能器的接收端接收超聲波信號,在連接示波器觀測波形的同時(shí),連接數(shù)據(jù)采集卡并將接收到的數(shù)據(jù)上傳至PC上位機(jī).

圖3 實(shí)驗(yàn)系統(tǒng)流程圖Figure 3 Experimental system flow chart

2.2 實(shí)驗(yàn)測量設(shè)備及條件

本文使用的是型號為DG1022Z的信號發(fā)生器,正弦波范圍1 μHz～25 MHz,正弦波幅值范圍1.0 mVpp～10 Vpp.采用的功率放大器是TYH-200M中頻功率放大器,可放大增益為0～10.超聲換能器的型號為DYU-200K-01F,其中心頻率為(200±30) kHz,量程為0.02～1 m,工作溫度為-40～80 ℃,耐壓值不大于0.5 MPa.示波器是tekronix生產(chǎn)的DPO 2024B示波器.由于本文發(fā)生的超聲信號的頻率為200 kHz,按照一般的選型規(guī)則,數(shù)據(jù)采集卡的采樣頻率應(yīng)為所需采集信號頻率的5～10倍,那么數(shù)據(jù)采集卡的采樣頻率應(yīng)該要達(dá)到1 M以上,因此本文采用的數(shù)據(jù)采集卡為美國NI公司生產(chǎn)的USB6363 DAQ數(shù)據(jù)采集卡,16位AD分辨率,2 MS/s采樣速率,10 ns時(shí)間分辨率.

2.3 實(shí)驗(yàn)方案設(shè)計(jì)及數(shù)據(jù)預(yù)處理

氣液兩相流流型識別實(shí)驗(yàn)主要包括產(chǎn)生流型、發(fā)送超聲波信號以及接收采集超聲信號三個(gè)部分,其中泡狀流、彈狀流和環(huán)狀流分別進(jìn)行多組實(shí)驗(yàn),每組實(shí)驗(yàn)時(shí)長為2 min.由于純水狀態(tài)下信號并未發(fā)生衰減,所以做1組實(shí)驗(yàn).

本文在垂直透明管道中設(shè)置了四種氣液兩相流流型,分別為純水狀態(tài)、泡狀流、彈狀流和環(huán)狀流.純水狀態(tài)、泡狀流、彈狀流和環(huán)狀流的衰減信號分別如圖4-7.

圖4 純水狀態(tài)衰減信號圖Figure 4 Pure water state attenuation signal diagram

圖5 泡狀流衰減信號圖Figure 5 Bubble flow attenuation signal diagram

圖6 彈狀流衰減信號圖Figure 6 Slug flow attenuation signal diagram

圖7 環(huán)狀流衰減信號圖Figure 7 Annular flow attenuation signal diagram

由圖4-7可見,在不同流型下,超聲換能器接收端接收到的信號強(qiáng)度差異比較明顯.純水狀態(tài)下,超聲信號并沒有發(fā)生衰減;泡狀流狀態(tài)下,當(dāng)超聲信號遇到小氣泡時(shí),會產(chǎn)生較小程度上的衰減;彈狀流狀態(tài)下,當(dāng)遇到彈狀大氣泡時(shí),超聲信號在一段時(shí)間內(nèi)會產(chǎn)生較大程度上的衰減,遇到小氣泡時(shí),產(chǎn)生較小程度上的衰減;而環(huán)狀流狀態(tài)下,超聲信號則呈現(xiàn)出長時(shí)間內(nèi)很大程度的衰減.

3 基于符號動(dòng)態(tài)濾波的識別方法

3.1 流型的特征值提取

因?yàn)榉柤痥的大小無關(guān)乎信號測量時(shí)間段內(nèi)具體的信號波動(dòng)數(shù)值,所以可以在整個(gè)實(shí)驗(yàn)過程中截取任意時(shí)間段、任意流型的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行分析,以此來確定k的大小.而對于閾值ε的選擇既要考慮到計(jì)算準(zhǔn)確度同時(shí)也要兼顧超聲信號采集速度和噪音影響等因素.

取泡狀流的一個(gè)時(shí)間段內(nèi)的數(shù)據(jù)集來分析,在當(dāng)h(k)<ε時(shí),得到k的取值,本文選擇ε<0.25,則得到k=5.

一般來說,閾值較小,則求得的符號集就偏大,也就是劃分得更加細(xì)致,這可能會增加噪聲,導(dǎo)致錯(cuò)誤符號出現(xiàn)的概率也增大,而且計(jì)算量也會增大.如果閾值取得較大,求得的符號集就較小,也許不能較好地體現(xiàn)原始時(shí)間序列的所包含的關(guān)于流型的信息.所以,在選擇閾值的時(shí)候,需要兼顧計(jì)算速度和準(zhǔn)確度.

本文選用均勻分割的方式對解析信號進(jìn)行分割.其中,純水狀態(tài)、泡狀流、彈狀流和環(huán)狀流的部分?jǐn)?shù)據(jù)分組結(jié)果如圖8-11.

圖8 純水狀態(tài)分組結(jié)果Figure 8 Pure water status grouping results

圖9 泡狀流分組結(jié)果Figure 9 Bubble flow grouping results

圖10 彈狀流分組結(jié)果Figure 10 Slug flow grouping results

圖11 環(huán)狀流分組結(jié)果Figure 11 Ring stream grouping results

接下來,需要將分組完成后的數(shù)據(jù)進(jìn)行符號化.將0～0.2之間的數(shù)據(jù)命名為a,0.2～0.4之間的數(shù)據(jù)命名為b,0.4～0.6之間的數(shù)據(jù)命名為c,0.6～0.8之間的數(shù)據(jù)命名為d,0.8～1之間的數(shù)據(jù)命名為e.由此,得到各流型對應(yīng)的符號序列.

3.2 建立流型識別模型

Ray在應(yīng)用過程中發(fā)現(xiàn),取D=1和D=2時(shí),所產(chǎn)生的結(jié)果相近[3].所以,本文取D=1.那么狀態(tài)集Q={a,b,c,d,e}.然后根據(jù)3.1小節(jié)選擇的較優(yōu)的符號集大小,可以構(gòu)造出性能較好的D-Markov模型.選取部分實(shí)驗(yàn)數(shù)據(jù)作為訓(xùn)練樣本,將訓(xùn)練樣本輸入D-Markov模型,得出可以體現(xiàn)每種流型信號波動(dòng)情況的基準(zhǔn)狀態(tài)轉(zhuǎn)移矩陣.同理,得到能代表測試樣本的狀態(tài)轉(zhuǎn)移矩陣,與每種流型的基準(zhǔn)狀態(tài)轉(zhuǎn)移矩陣進(jìn)行相似性度量.相似度高說明其波動(dòng)信號更接近與之相對應(yīng)的流型,從而達(dá)到流型識別的目的.

3.3 相似性度量

根據(jù)狀態(tài)轉(zhuǎn)移矩陣,可以分別求得訓(xùn)練樣本的左特征向量.其中純水狀態(tài)訓(xùn)練樣本的左特征向量記為q1,將5組泡狀流訓(xùn)練樣本的左特征向量記為q21,q22,q23,q24,q25對這5個(gè)左特征向量求其歐拉距離,得到相似度,其統(tǒng)計(jì)如表1.

表1 泡狀流的狀態(tài)轉(zhuǎn)移矩陣的相似度

可見,5組泡狀流訓(xùn)練樣本的狀態(tài)轉(zhuǎn)移矩陣之間的距離很小,分布在0.020 5～0.074 1之間.這可以說明相對應(yīng)的5個(gè)狀態(tài)轉(zhuǎn)移矩陣之間的相似度較高.由于q21相對于其它4個(gè)左特征向量,與每個(gè)特征向量的距離最小,所以將q21作為泡狀流的基準(zhǔn)特征向量,與之對應(yīng)的基準(zhǔn)狀態(tài)轉(zhuǎn)移矩陣為Π21.

將5組彈狀流的左特征向量記為q31,q32,q33,q34,q35分別對任意兩個(gè)特征向量求歐拉距離,得到對應(yīng)相似度,如表2.

表2 彈狀流的狀態(tài)轉(zhuǎn)移矩陣的相似度

可見,5組彈狀流訓(xùn)練樣本的狀態(tài)轉(zhuǎn)移矩陣之間的距離較小,分布在0.054 0～0.139 2之間.這可以說明相對應(yīng)的5個(gè)狀態(tài)轉(zhuǎn)移矩陣之間的相似度較高.由于q34相對于其它4個(gè)左特征向量,與每個(gè)特征向量的距離最小,所以將q34作為彈狀流的基準(zhǔn)特征向量,與之對應(yīng)的基準(zhǔn)狀態(tài)轉(zhuǎn)移矩陣為Π34.

將5組環(huán)狀流訓(xùn)練樣本的左特征向量記為q41,q42,q43,q44,q45分別對任意兩個(gè)特征向量求歐拉距離,得到對應(yīng)相似度,如表3.

表3 環(huán)狀流的狀態(tài)轉(zhuǎn)移矩陣的相似度

可見,5組環(huán)狀流訓(xùn)練樣本的狀態(tài)轉(zhuǎn)移矩陣之間的距離較小,分布在0.037 1～0.080 1之間.這可以說明相對應(yīng)的5個(gè)狀態(tài)轉(zhuǎn)移矩陣之間的相似度較高.由于q41相對于其它4個(gè)左特征向量,與每個(gè)特征向量的距離最小,所以將q41作為泡狀流的基準(zhǔn)特征向量,與之對應(yīng)的基準(zhǔn)狀態(tài)轉(zhuǎn)移矩陣為Π41.

接下來,將純水左特征向量q1、泡狀流基準(zhǔn)特征向量q31、彈狀流基準(zhǔn)特征向量q31和環(huán)狀流基準(zhǔn)特征向量q41,進(jìn)行相似度測量,分別計(jì)算它們兩兩之間的歐拉距離,如表4.

表4 流型狀態(tài)轉(zhuǎn)移矩陣的相似度

可見,4種不同流型的基準(zhǔn)狀態(tài)轉(zhuǎn)移矩陣之間的歐氏距離相較于同種流型下的各組狀態(tài)轉(zhuǎn)移矩陣的歐氏距離要大得多.于是可以得出,在不同流型下,采集到的信號強(qiáng)度差異較大,而同種流型下采集到的信號強(qiáng)度較為相近.因此,相似性度量可以為流型識別提供一種新的思路.

4 結(jié)語

為了更加深入地認(rèn)識氣液兩相流流動(dòng)結(jié)構(gòu)及相應(yīng)流動(dòng)參數(shù)的測量,本文搭建了一套垂直上升的氣液兩相流流型識別實(shí)驗(yàn)系統(tǒng),對管徑為50 mm的透明管道進(jìn)行了氣液兩相流流型的識別實(shí)驗(yàn),其中超聲換能器用于獲取垂直上升氣液兩相流的波動(dòng)信號.然后運(yùn)用符號動(dòng)態(tài)濾波方法對實(shí)驗(yàn)過程采集到超聲波信號分析處理,得到能表征兩相流信號波動(dòng)情況的狀態(tài)轉(zhuǎn)移矩陣,對矩陣進(jìn)行相似性度量可判別不同流型的狀態(tài)轉(zhuǎn)移矩陣之間的差異程度,從而為以后氣液兩相流流型識別提供了一種可行性思路.