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

（1 河北大學質量技術監(jiān)督學院，河北 保定 071002；2 計量儀器與系統(tǒng)國家地方聯(lián)合工程研究中心，河北 保定071002；3 河北省能源計量與安全檢測技術重點實驗室，河北 保定 071002）

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

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

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

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

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

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

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

圖1 實驗系統(tǒng)結構

1.2 實驗裝置

基于聲發(fā)射技術的彈狀流測量實驗在河北大學高精度氣液兩相流實驗裝置上進行，實驗測量裝置如圖2所示。

圖2 高精度氣液兩相流實驗設備

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

表1 實驗參數(shù)范圍

圖3 水平管Baker流型

2 實驗數(shù)據(jù)分析

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

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

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

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

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

從變分模態(tài)分解原理可以看出，VMD 的參數(shù)設置對分解結果有本質的影響，不同的參數(shù)組合會產生不同的結果，并且手動調整參數(shù)非常耗時。因此，本文利用鯨魚優(yōu)化算法來優(yōu)化VMD參數(shù)。

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

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

2.2 噪聲信號解耦分析

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

式中，ui為第i個IMF分量；t為采樣時間。

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

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

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

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

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

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

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

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

3.1 彈頻及氣彈長度

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

圖9 工況點（usl=0.283m/s，usg=0.849m/s）時聲發(fā)射時域重構信號

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

式中，f為彈頻；X為選擇點的橫坐標值；S為彈單元個數(shù)。

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

圖11 不同流速下彈狀流氣彈長度

3.2 基于CatBoost算法的預測模型

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

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

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

圖12 聲發(fā)射信號特征值相關性分析

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

表3 CatBoost參數(shù)設置

圖13 預測模型

為了評估預測結果，采用平均絕對百分比誤差（MAPE）作為統(tǒng)計指標，其定義如式(10)所示。誤差分析如圖14 所示，彈狀流頻率預測模型MAPE為5.12%，95.25%實驗點的相對偏差都在±15%范圍內；氣彈長度預測模型MAPE 為7.77%，90.48%實驗點的相對偏差都在±15%范圍內。

圖14 誤差分析

式中，yi為經分析所得實驗值；?i為預測值。

4 結論

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