李新龍，張昭，王嘉輝，張延勝，董芳

（1 河北大學(xué)質(zhì)量技術(shù)監(jiān)督學(xué)院，河北 保定 071002；2 計(jì)量?jī)x器與系統(tǒng)國(guó)家地方聯(lián)合工程研究中心，河北 保定071002；3 河北省能源計(jì)量與安全檢測(cè)技術(shù)重點(diǎn)實(shí)驗(yàn)室，河北 保定 071002）

彈狀流是一種典型的氣液兩相流，具有間歇性與不穩(wěn)定性的特點(diǎn)，廣泛存在于石油、天然氣、核電站等工業(yè)過(guò)程中[1]。彈狀流動(dòng)過(guò)程中，氣彈與液塞交替出現(xiàn)，會(huì)導(dǎo)致管道產(chǎn)生高強(qiáng)度振動(dòng)、高強(qiáng)度壓力、加快管道老化。精準(zhǔn)預(yù)測(cè)氣彈特性參數(shù)對(duì)于工業(yè)管道設(shè)計(jì)和安全操作至關(guān)重要。

眾多學(xué)者提出了不同測(cè)量方法對(duì)彈狀流進(jìn)行了實(shí)驗(yàn)研究，根據(jù)是否對(duì)流場(chǎng)產(chǎn)生擾動(dòng)可以分為侵入式和非侵入式[2]。侵入式包括差壓法[3-4]、電導(dǎo)探針法[5]，這些方法無(wú)疑會(huì)對(duì)流場(chǎng)造成干擾，尤其是微通道流動(dòng)測(cè)量中[6]。非侵入式主要包括高速攝像法[7-9]、紅外光學(xué)法[10]、電容法[11]和超聲法等[12]。與侵入傳感器相比，非侵入傳感器具有很大優(yōu)勢(shì)，但現(xiàn)有的非侵入傳感器也存在許多問題，如對(duì)管道材料有一定要求、操作難度大、響應(yīng)速度慢等。聲發(fā)射法作為非侵入法的一種，具有采樣頻率高、響應(yīng)速度快、操作簡(jiǎn)單、對(duì)流場(chǎng)無(wú)干擾等特點(diǎn)。

聲發(fā)射技術(shù)是工業(yè)生產(chǎn)中廣泛使用的無(wú)損檢測(cè)技術(shù)，用于發(fā)現(xiàn)和定位設(shè)備中的缺陷[13-14]，但越來(lái)越多的學(xué)者將聲發(fā)射技術(shù)應(yīng)用于多相流測(cè)量領(lǐng)域[15-16]。Husin等[17]利用聲發(fā)射技術(shù)采集單個(gè)氣泡開始到破裂時(shí)的信號(hào)，得到聲發(fā)射信號(hào)與氣泡大小和液體黏度相互關(guān)系。Zhao等[18]利用聲發(fā)射技術(shù)測(cè)量泡狀流和彈狀流的流動(dòng)噪聲，信號(hào)通過(guò)Hilbert Huang 變換、R/S 分析和流動(dòng)噪聲數(shù)學(xué)模型進(jìn)行了流型識(shí)別。聲發(fā)射技術(shù)可以實(shí)時(shí)采集氣液兩相流動(dòng)中由于氣液、氣固和液固之間相互作用產(chǎn)生的耦合噪聲，通過(guò)噪聲解耦可以獲取不同流型下氣液兩相流動(dòng)信息，實(shí)現(xiàn)流動(dòng)參數(shù)定量測(cè)量。該技術(shù)具有無(wú)侵入性、高采樣率的特點(diǎn)[19]。

本文設(shè)計(jì)了一種基于聲發(fā)射原理的彈狀流氣彈參數(shù)測(cè)量系統(tǒng)。利用鯨魚優(yōu)化算法變分模態(tài)分解對(duì)聲發(fā)射信號(hào)進(jìn)行分析，結(jié)合彈狀流流動(dòng)特性從能量和熵的角度，實(shí)現(xiàn)了水平管彈狀流噪聲信息的解耦、彈頻和氣彈長(zhǎng)度的測(cè)量。通過(guò)彈狀流彈頻與氣彈長(zhǎng)度的相關(guān)性分析，選取了8個(gè)特征值，并基于CatBoost 算法構(gòu)建了彈狀流彈頻、氣彈長(zhǎng)度的預(yù)測(cè)模型。

1 氣彈參數(shù)測(cè)量傳感器

1.1 基于聲發(fā)射原理的氣彈參數(shù)測(cè)量傳感器設(shè)計(jì)

本文為了獲取彈狀流氣彈特性參數(shù)信息，設(shè)計(jì)了一種基于聲發(fā)射傳感器的彈狀流氣彈參數(shù)測(cè)量系統(tǒng)，如圖1 所示。4 個(gè)探頭C1、C2、C3、C4 被上下對(duì)稱固定在管徑為50mm的水平管道上，利用高真空油脂耦合劑使4 個(gè)探頭與管道壁面充分接觸。如圖1所示，SH-Ⅲ為信號(hào)調(diào)理裝置，對(duì)探頭測(cè)得的微弱電信號(hào)進(jìn)行放大、整形、濾波，并轉(zhuǎn)換成數(shù)字量傳輸給上位機(jī)。實(shí)驗(yàn)開始前設(shè)置聲發(fā)射采集閾值，作為要接收信號(hào)幅值的下限值。經(jīng)測(cè)量在空管時(shí)聲發(fā)射探頭所測(cè)周圍環(huán)境噪聲均在40dB 以下，因此為保持信號(hào)靈敏度的同時(shí)減少環(huán)境噪聲的影響，本次實(shí)驗(yàn)將采集閾值設(shè)置為40dB。

圖1 實(shí)驗(yàn)系統(tǒng)結(jié)構(gòu)

1.2 實(shí)驗(yàn)裝置

基于聲發(fā)射技術(shù)的彈狀流測(cè)量實(shí)驗(yàn)在河北大學(xué)高精度氣液兩相流實(shí)驗(yàn)裝置上進(jìn)行，實(shí)驗(yàn)測(cè)量裝置如圖2所示。

圖2 高精度氣液兩相流實(shí)驗(yàn)設(shè)備

氣相由空氣壓縮機(jī)產(chǎn)生，經(jīng)穩(wěn)壓罐干燥機(jī)提供穩(wěn)定氣源，流經(jīng)科里奧利流量計(jì)（DN8mm、DN40mm，U=0.1%，k=2）測(cè)量后進(jìn)入實(shí)驗(yàn)管段。液相通過(guò)水泵進(jìn)入管道，通過(guò)電磁流量計(jì)（DN10mm、DN32mm，U=0.2%，k=2） 和科里奧利流量計(jì)（DN40mm，U=0.1%，k=2）測(cè)量后通過(guò)引射器進(jìn)入實(shí)驗(yàn)管段后與氣相充分混合。氣液輸送管道上安裝調(diào)節(jié)開關(guān)實(shí)現(xiàn)對(duì)氣液流量的精準(zhǔn)控制。實(shí)驗(yàn)管道介質(zhì)溫度和壓力由溫度變送器（JWB/38Z/A，0～50℃，U=0.2%，k=2）和壓力變送器（JT-8016CRA，0～1MPa，U=0.5%，k=2）實(shí)時(shí)采集。聲發(fā)射探頭采樣頻率為100kHz，采樣時(shí)間為10s。本文的實(shí)驗(yàn)參數(shù)范圍見表1。通過(guò)水平管Baker 流型圖可以看出本文的實(shí)驗(yàn)流型均為彈狀流（圖3）。

表1 實(shí)驗(yàn)參數(shù)范圍

圖3 水平管Baker流型

2 實(shí)驗(yàn)數(shù)據(jù)分析

2.1 變分模態(tài)分解

變分模態(tài)分解（variational modal decomposition，VMD）是一種完全非遞歸、且具有完備數(shù)學(xué)理論支撐的信號(hào)處理方法。相比于小波變換，VMD克服了小波閾值的不同選取對(duì)信號(hào)分解、降噪時(shí)的影響；相比于經(jīng)驗(yàn)?zāi)B(tài)分解（empirical mode decomposition，EMD）分解，VMD 避免了EMD 分解時(shí)帶來(lái)的模態(tài)混疊和端點(diǎn)效應(yīng)，因此VMD 在聲發(fā)射數(shù)據(jù)的處理上具有明顯優(yōu)勢(shì)[20]。VMD 的核心思想是構(gòu)建和求解變分問題，假設(shè)原始信號(hào)F被分解為K個(gè)分量，保證分解序列為具有中心頻率有限帶寬的模態(tài)分量，同時(shí)各模態(tài)的估計(jì)帶寬之和最小，約束條件為所有模態(tài)之和與原始信號(hào)相等。相應(yīng)約束變分表達(dá)如式(1)。

式中，K為需要分解的模態(tài)個(gè)數(shù)（正整數(shù)）；{uk}、{ωk}分別對(duì)應(yīng)分解后第k個(gè)模態(tài)分量和中心頻率；δ(t)為狄拉克函數(shù)；*為卷積運(yùn)算符。

在該變分模型中，對(duì)uk(t)進(jìn)行Hilbert變換得到單邊譜，然后通過(guò)乘指數(shù)項(xiàng)exp(-jωkt)調(diào)節(jié)預(yù)估的中心頻率，并將其頻譜調(diào)整到相應(yīng)的基帶上。為了將原變分約束模型轉(zhuǎn)換為非變分約束模型，引入增廣Lagrange函數(shù)如式（2）。

式中，α為二次懲罰因子；λ為L(zhǎng)agrange因子。利用交替方向乘子算法不斷迭代求解出如下自適應(yīng)中心頻率及各IMF分量表達(dá)式(3)～式(5)。

從變分模態(tài)分解原理可以看出，VMD 的參數(shù)設(shè)置對(duì)分解結(jié)果有本質(zhì)的影響，不同的參數(shù)組合會(huì)產(chǎn)生不同的結(jié)果，并且手動(dòng)調(diào)整參數(shù)非常耗時(shí)。因此，本文利用鯨魚優(yōu)化算法來(lái)優(yōu)化VMD參數(shù)。

鯨魚優(yōu)化算法（whale optimization algorithm，WOA）是一種基于座頭鯨種群迭代進(jìn)化搜索的元啟發(fā)式優(yōu)化算法，算法通過(guò)模擬座頭鯨獨(dú)特的狩獵行為，實(shí)現(xiàn)參數(shù)尋優(yōu)。利用排列熵計(jì)算初始種群中每個(gè)個(gè)體的適應(yīng)度，衡量參數(shù)組合的分解效果。熵越小，則時(shí)間序列分布越有規(guī)律，表示VMD 處理得到的IMF包含更多的有效信息；反之，時(shí)間序列越接近隨機(jī)分布，IMF中噪聲成分更多。因此，當(dāng)排列熵最小時(shí)，對(duì)應(yīng)的參數(shù)最優(yōu)，得出實(shí)驗(yàn)最優(yōu)K值為17，α值為29872。鯨魚優(yōu)化算法變分模態(tài)分解流程如圖4所示。

圖4 鯨魚優(yōu)化算法變分模態(tài)分解流程

2.2 噪聲信號(hào)解耦分析

彈狀流在管道流動(dòng)過(guò)程中會(huì)產(chǎn)生流動(dòng)噪聲，流動(dòng)噪聲主要包括氣液相互作用噪聲、氣固耦合噪聲、液固耦合噪聲和流體內(nèi)部產(chǎn)生的輕微噪聲。采集的聲發(fā)射信號(hào)經(jīng)過(guò)VMD 分解為固定數(shù)量的本征模態(tài)信號(hào)后，彈狀流的特征信息寄存于各IMF分量中。在VMD 對(duì)聲發(fā)射信號(hào)f(t)進(jìn)行分解時(shí)，可由式(6)計(jì)算出不同IMF 分量信號(hào)的能量E1、E2、…、Ek。

式中，ui為第i個(gè)IMF分量；t為采樣時(shí)間。

由于VMD 分解的IMF 分量具有正交性，因此所有IMF 的能量之和恒等于原始信號(hào)f(t)中的總能量，各IMF 分量包含原始信號(hào)中的不同頻率成分，包含不同的信息值。以實(shí)驗(yàn)流量條件（usl=0.990m/s，usg=0.849m/s）為例，VMD分解后的流動(dòng)噪聲如圖5所示，將氣水兩相流噪聲分解為16 個(gè)本征模態(tài)分量和1個(gè)殘余模態(tài)分量，本征模態(tài)分量中心頻率由高到低排列。

在流體流過(guò)管道時(shí)，氣相、液相、管道相互作用產(chǎn)生的噪聲會(huì)以彈性波的形式被附著在管壁上的聲發(fā)射探頭捕捉到，噪聲內(nèi)有大量的有用信息。通過(guò)分析VMD 分解后本征模態(tài)分量的能量變化，從微觀角度對(duì)兩相流內(nèi)在機(jī)理進(jìn)行研究。圖6顯示了在彈狀流分析中能量，隨著流量的增加，氣液管道相互作用碰撞產(chǎn)生的流動(dòng)噪聲和卷繞效應(yīng)增強(qiáng)，產(chǎn)生更強(qiáng)的流動(dòng)噪聲，能量也更大。

為進(jìn)一步對(duì)彈狀流噪聲信號(hào)進(jìn)行分析引入信息熵的概念。信息熵是對(duì)信號(hào)的復(fù)雜性與不確定性的一種度量指標(biāo)，與微狀態(tài)的對(duì)數(shù)成正比。微狀態(tài)的數(shù)量越多，則越混亂，信息量越多，反之信息熵的減少則會(huì)提高組織水平，信號(hào)混亂程度越低。彈狀流流動(dòng)過(guò)程中氣相、液相、管道相互作用產(chǎn)生不同混亂程度噪聲信號(hào)。信息的獲取意味著概率分布在可能性之間的集中，信息熵的大小能間接反映出彈狀流流動(dòng)過(guò)程中系統(tǒng)復(fù)雜程度，系統(tǒng)越復(fù)雜，混亂程度越高，信息熵越大。彈狀流氣彈與液塞具有不同的混亂程度，利用信息熵對(duì)彈狀流進(jìn)行進(jìn)一步研究。信息熵定義如式(7)。

式中，p(xi)為隨機(jī)事件xi的概率。

聲發(fā)射探頭C1、C2 在不同流速下能量如圖7所示，不同流速下信息熵如圖8所示。IMF的能量與熵主要集中在前5個(gè)分量中，且前5個(gè)分量?jī)蓚€(gè)探頭能量與熵之間有較大差距。結(jié)合彈狀流流動(dòng)特性分析，在彈狀流流動(dòng)過(guò)程中，會(huì)出現(xiàn)前部液塞脫落，后部液塞不斷拾起，并且拾起過(guò)程中會(huì)在液塞前方產(chǎn)生一個(gè)混合區(qū)域，液塞脫落的過(guò)程中會(huì)不斷與管壁產(chǎn)生相互作用，此過(guò)程具有最高的混亂度和熵。IMF6～I(xiàn)MF12 間能量與熵?cái)?shù)值差距逐漸變小且存在近似相等的情況，由于上下探頭對(duì)彈狀流內(nèi)氣液相互作用噪聲采集條件相同，所以存在很小的差異。IMF12～I(xiàn)MF14 為液固噪聲，液固噪聲主要由液相與管道之間摩擦產(chǎn)生，與氣固噪聲和氣液相互作用噪聲相比系統(tǒng)混亂程度較低，具有最小的能量與熵且比較穩(wěn)定。通過(guò)對(duì)聲發(fā)射測(cè)量裝置設(shè)置噪聲門檻，避免了絕大部分噪聲干擾，但仍有小部分噪聲具有低頻性、不確定性，與IMF15、IMF16 具有相同特性，因此本文選擇前IMF14個(gè)分量進(jìn)行重構(gòu)濾除周圍噪聲。

圖7 不同流速下C1、C2探頭能量

圖8 不同流速下C1、C2探頭信息熵

3 氣彈特征參數(shù)測(cè)量

3.1 彈頻及氣彈長(zhǎng)度

彈狀流頻率為在單位時(shí)間內(nèi)通過(guò)管道內(nèi)固定點(diǎn)彈單元的平均數(shù)量，彈頻公式如式(8)所示。彈狀流是典型的間歇性兩相流，實(shí)驗(yàn)過(guò)程中，氣彈流過(guò)時(shí)，噪聲信號(hào)主要由氣固噪聲和氣液噪聲構(gòu)成，液塞流過(guò)時(shí)的噪聲信號(hào)主要由液固噪聲構(gòu)成。如圖8分析結(jié)果可知，氣彈處噪聲信號(hào)能量大于液塞處，即聲發(fā)射時(shí)域信號(hào)幅值更大。以u(píng)sl=0.283m/s、usg=0.849m/s實(shí)驗(yàn)工況點(diǎn)為例，其聲發(fā)射信號(hào)強(qiáng)度隨時(shí)間變化幅值如圖9所示，氣彈、液塞交替時(shí)具有較為明顯的突變信號(hào)，以此為特征判斷氣彈頭、氣彈尾、液塞尾。已知實(shí)驗(yàn)頻率為100kHz，采樣時(shí)間為10s，氣彈流過(guò)時(shí)通過(guò)圖像得出彈頭彈尾橫坐標(biāo)，得出彈狀流彈頻公式。不同流動(dòng)條件下的氣彈頻率如圖10 所示，從圖中可以看出彈狀流彈頻隨氣相表觀流速增大而減小，隨液相表觀流速增大而增大。

圖9 工況點(diǎn)（usl=0.283m/s，usg=0.849m/s）時(shí)聲發(fā)射時(shí)域重構(gòu)信號(hào)

圖10 不同流速下彈狀流彈頻

式中，f為彈頻；X為選擇點(diǎn)的橫坐標(biāo)值；S為彈單元個(gè)數(shù)。

氣彈長(zhǎng)度對(duì)氣液兩相孔隙率研究有較大幫助，氣彈長(zhǎng)度為彈狀流平移速度與氣彈通過(guò)聲發(fā)射探頭時(shí)間的乘積，彈狀流氣彈平移速度已有大量學(xué)者進(jìn)行了研究。Wang等[21]在50mm管徑下對(duì)彈狀流進(jìn)行了測(cè)量，將測(cè)量到的氣彈平移速度與Bendiksen 提出的預(yù)測(cè)模型[22]進(jìn)行了比較，結(jié)果表明在較低的混合物速度Vm范圍內(nèi)，測(cè)量值與Bendiksen 的預(yù)測(cè)值之間的一致性非常好。本文所選工況點(diǎn)弗勞德數(shù)Fr最大值為3.4，利用式(9)求得彈速。從聲發(fā)射時(shí)域圖中得出每個(gè)工況下氣彈和液塞平均時(shí)間，進(jìn)而求出氣彈長(zhǎng)度。不同流動(dòng)條件下的氣彈長(zhǎng)度值如圖11所示。從圖11中可以看出，氣彈長(zhǎng)度隨氣相表觀流速增大而增大，隨液相表觀流速增大而減小。

圖11 不同流速下彈狀流氣彈長(zhǎng)度

3.2 基于CatBoost算法的預(yù)測(cè)模型

從C1、C2 聲發(fā)射時(shí)域信號(hào)中提取以下6 個(gè)特征值作為信號(hào)特征進(jìn)行進(jìn)一步分析，對(duì)于信號(hào)x（n個(gè)元素的）可定義以下參數(shù)，見表2。

表2 聲發(fā)射信號(hào)特征公式

相關(guān)性系數(shù)可以用于研究變量之間的線性相關(guān)程度，本文選擇皮爾遜相關(guān)系數(shù)研究聲發(fā)射信號(hào)特征量與彈狀流彈頻、氣彈長(zhǎng)度之間的相關(guān)性。如圖12所示為C1、C2探頭噪聲信號(hào)的最大值、最小值、平均值、標(biāo)準(zhǔn)差、偏斜度、峭度與彈狀流彈頻和氣彈長(zhǎng)度之間的相關(guān)性（C3、C4與C1、C2具有相似測(cè)量條件，只選擇C1、C2進(jìn)行分析）。相關(guān)系數(shù)絕對(duì)值越接近于1，相關(guān)性越強(qiáng)，越接近于0則相關(guān)性越弱。由圖12可以看出，C1探頭標(biāo)準(zhǔn)差和峭度、C2探頭標(biāo)準(zhǔn)差與彈頻之間具有最好的相關(guān)性；C1探頭標(biāo)準(zhǔn)差和峭度、C2探頭標(biāo)準(zhǔn)差和峭度與氣彈長(zhǎng)度之間具有最好的相關(guān)性。將這些相關(guān)性好的特征變量作為預(yù)測(cè)模型的特征輸入，建立彈狀流特征參數(shù)預(yù)測(cè)模型。

圖12 聲發(fā)射信號(hào)特征值相關(guān)性分析

CatBoost 是一種基于對(duì)稱決策樹（oblivious trees）算法，參數(shù)少、支持類別型變量和高準(zhǔn)確性的GBDT 框架，能高效合理地處理類別型特征。CatBoost 用了組合類別特征，利用特征之間的聯(lián)系，極大豐富了所用特征的維度。CatBoost 采用排序提升的方式替換傳統(tǒng)算法中梯度估計(jì)方法進(jìn)而減輕了梯度偏差和預(yù)測(cè)偏移的問題，減少了過(guò)擬合的發(fā)生，進(jìn)而提高算法的準(zhǔn)確性和泛化能力[23]。因此選用CatBoost來(lái)預(yù)測(cè)彈狀流頻率與氣彈長(zhǎng)度，訓(xùn)練集和測(cè)試集數(shù)據(jù)為4∶1，CatBoost 參數(shù)見表3，預(yù)測(cè)模型如圖13所示。

表3 CatBoost參數(shù)設(shè)置

圖13 預(yù)測(cè)模型

為了評(píng)估預(yù)測(cè)結(jié)果，采用平均絕對(duì)百分比誤差（MAPE）作為統(tǒng)計(jì)指標(biāo)，其定義如式(10)所示。誤差分析如圖14 所示，彈狀流頻率預(yù)測(cè)模型MAPE為5.12%，95.25%實(shí)驗(yàn)點(diǎn)的相對(duì)偏差都在±15%范圍內(nèi)；氣彈長(zhǎng)度預(yù)測(cè)模型MAPE 為7.77%，90.48%實(shí)驗(yàn)點(diǎn)的相對(duì)偏差都在±15%范圍內(nèi)。

圖14 誤差分析

式中，yi為經(jīng)分析所得實(shí)驗(yàn)值；?i為預(yù)測(cè)值。

4 結(jié)論

利用聲發(fā)射技術(shù)對(duì)彈狀流噪聲信號(hào)進(jìn)行測(cè)量，使用鯨魚優(yōu)化算法變分模態(tài)分解技術(shù)對(duì)彈狀流噪聲信號(hào)進(jìn)行了處理，從能量和熵的角度進(jìn)行分析，完成了對(duì)水平管彈狀流噪聲信號(hào)的解耦。分析得出彈狀流彈頻隨氣相表觀流速增大而減小，隨液相表觀流速增大而增大；氣彈長(zhǎng)度隨氣相表觀流速增大而增大，隨液相表觀流速增大而減小。對(duì)噪聲信號(hào)進(jìn)行重構(gòu)，提取重構(gòu)后信號(hào)的能量、信息熵、標(biāo)準(zhǔn)差、峭度等作為特征值進(jìn)行相關(guān)性分析，選擇相關(guān)性高的特征值作為特征輸入建立了彈狀流彈頻和氣彈長(zhǎng)度的預(yù)測(cè)模型，彈狀流頻率預(yù)測(cè)模型MAPE為5.12%，95.25%實(shí)驗(yàn)點(diǎn)的相對(duì)偏差都在±15%范圍內(nèi)；氣彈長(zhǎng)度預(yù)測(cè)模型MAPE為7.77%，90.48%實(shí)驗(yàn)點(diǎn)的相對(duì)偏差都在±15%范圍內(nèi)。