基于GBDT的氣液兩相流相分布測量模型

曾思睿，孔明

（中國計量大學(xué)計量測試工程學(xué)院，浙江 杭州 310018）

氣液兩相流是兩相流中最常見的形式之一，涉及諸多工業(yè)領(lǐng)域[1-2]，例如化工行業(yè)中的吸收塔、蒸餾塔、冷凝器等[3-4]。氣液兩相在空間中的分布有不同形式，這些不同的分布形式會改變力學(xué)特性、傳熱特性以及傳質(zhì)特性等，因此在工業(yè)領(lǐng)域和科學(xué)研究領(lǐng)域，氣液兩相流的相分布測量技術(shù)具有重要意義。該技術(shù)能夠支持工業(yè)生產(chǎn)過程的流動狀態(tài)監(jiān)測和控制[5]，并為理論模型和仿真模型提供可靠數(shù)據(jù)。憑借這項技術(shù)，人們能夠?qū)崟r監(jiān)測和掌握工業(yè)流體的分布情況，從而實現(xiàn)精確控制和優(yōu)化工藝。

機器學(xué)習(xí)在解決復(fù)雜非線性回歸問題中表現(xiàn)出優(yōu)異特性，為兩相流流動參數(shù)測量模型的建立提供了一種新思路。國內(nèi)外學(xué)者常用的機器學(xué)習(xí)主要有隨機森林、BP神經(jīng)網(wǎng)絡(luò)、支持向量機（SVM）、梯度提升決策樹（GBDT）以及深度學(xué)習(xí)等算法，已應(yīng)用于氣液兩相流流型識別[6]、氣泡輪廓重建[7]、流量和相含率測量[8]之中，但都未對相分布參數(shù)進行研究。

綜合現(xiàn)有的氣液兩相流參數(shù)測量研究文獻，對管道界面內(nèi)氣泡的大小和位置的研究較少。戴振韜[9]提出提取光強分布特征的方法測量相分布，采用線性回歸來擬合相分布測量模型，需要指出的是，該方法使用的是單波長且僅限于水平管道。在此基礎(chǔ)上，曹鳴等[10]提出一種基于雙波長透射法對小通道豎直上升氣泡流進行研究，并利用BP 神經(jīng)網(wǎng)絡(luò)建立測量模型。由于BP 神經(jīng)網(wǎng)絡(luò)具有收斂速度慢、預(yù)測能力和訓(xùn)練能力矛盾等問題[11]，本文提出用泛化性能較好、預(yù)測精度更高的GBDT算法來訓(xùn)練模型，并預(yù)測實驗中氣泡的相分布。

1 基本原理

1.1 雙波長測量原理

光學(xué)檢測系統(tǒng)如圖1所示，光源發(fā)出包含兩種波長的激光光線，光線通過1mm 的狹縫變?yōu)槠瑺罟馐?，再通過分束鏡分成兩束互相垂直的片狀激光，與管道垂直相交，在相交位置形成截面，稱該截面為檢測截面。包含兩種波長的片狀激光在檢測截面內(nèi)發(fā)生折射，管道和管道內(nèi)兩相流對不同波長光的折射率存在較大差異，且不同大小的氣泡位于兩相流中的不同位置時也會導(dǎo)致激光呈現(xiàn)出不同的折射效果，因此影響傳感器上接收光強分布的因素有波長、氣泡大小和氣泡位置。用高速CMOS傳感器實時記錄氣泡運動到檢測截面的光強分布數(shù)據(jù)，并將這些數(shù)據(jù)傳送至計算機。通過分析和處理這些光強分布數(shù)據(jù)，能夠獲取光強分布與氣泡位置和大小的重要信息，提取這些特征，用大量仿真數(shù)據(jù)訓(xùn)練基于GBDT的相分布測量模型，將提取出的光強分布特征輸入模型中，對實驗中氣泡經(jīng)過檢測截面時的半徑和位置坐標進行預(yù)測，根據(jù)預(yù)測結(jié)果對氣泡中心位置運動軌跡進行追蹤。

圖1 光學(xué)檢測系統(tǒng)

圖2是光線在檢測截面內(nèi)的傳播示意圖，一束包含波長為λ1、λ2的光線入射到1mm 厚的管道，Rin為管道的內(nèi)徑，Rout為管道的外徑，入射點到管道中心的距離為l，空氣的折射率為n，波長λ1、λ2在水中的折射率為n11、n12，在石英玻璃管道中的折射率分別為ng1、ng2。以管道中心為原點建立坐標系，原點到檢測平面的距離為D。由斯涅爾定律可以得出以下結(jié)論：一束包含兩種波長不同的光經(jīng)過管道和管道內(nèi)的氣液兩相流后，它們的出射角度是不同的，這意味著波長λ1、λ2將分別入射到光檢測平面上的不同點，設(shè)這兩點之間的距離為Δδ。

圖2 全液相管道光線傳播路徑

以波長為λ1的一束光線為例，圖中a、c、e、g表示入射角，b、d、f、h表示折射角，波長為λ1的出射光線與x方向夾角為δ1。分別計算出兩波長出射光線與x方向的夾角δ1、δ2，最后計算得出同一位置入射的兩束不同波長的光入射到光檢測平面不同兩點的距離Δδ，計算過程如下。

入射角a可用入射點到管道中心距離l和管道的外徑Rout表示如式(1)。

根據(jù)斯涅耳定律，光線經(jīng)管道外壁折射后的折射角b如式(2)。

結(jié)合式(1)～式(4)得入射角e=d，折射角f=c；入射角g=b，折射角h=a。

波長為λ1的光線通過管道后出射光線與x方向的夾角δ1與各角度的幾何關(guān)系如式(5)。

綜合式(5)～式(8)，可計算同一位置入射的兩束不同波長的光線入射到光檢測平面兩點的距離Δδ如式(9)。

以上為同一位置發(fā)出的一束包含兩個波長光線的計算過程，可以通過計算并記錄上百萬條入射光線的傳播路徑來獲得檢測平面上雙波長光強分布的結(jié)果。對于某個特定波長的光線而言，光檢測平面上接收的光強分布可以看作是反射和折射效應(yīng)的疊加。為了描述光在空間中的光強分布，可以使用復(fù)振幅函數(shù)。當(dāng)光線入射角度為θ時，相應(yīng)的光檢測平面的復(fù)振幅函數(shù)可以表示如式(12)。

式中，S(θ)為復(fù)振幅函數(shù)；Sreflect(θ)為反射部分的復(fù)振幅函數(shù)；Srefract,p(θ)為折射部分的復(fù)振幅函數(shù)。當(dāng)p=0時，表示光線在管道內(nèi)發(fā)生全發(fā)射，Srefract,0(θ)表示反射光的復(fù)振幅函數(shù)；當(dāng)p≥1時，Srefract,p(θ)表示入射光線在管道內(nèi)經(jīng)過(p-0)次反射后出射光線的復(fù)振幅函數(shù)。

根據(jù)式(9)～式(11)可得出：由于管道和管道內(nèi)液體對不同波長光的折射率存在差異，因此同一位置發(fā)出不同波長的光線出射到檢測平面上的位置也不同。對上百萬條光線進行追跡，通過式(12)，對光線進行疊加，最終光檢測平面上兩個波長的光形成不同的光強分布。當(dāng)有氣泡經(jīng)過管道的檢測段時，光線在氣液界面發(fā)生折射，導(dǎo)致大量光線被折射到光檢測平面外，檢測平面無法接收到這部分的光強分布，光強分布曲線出現(xiàn)缺失部分，氣液分界面對于不同波長光線的折射能力不同，因此在相分布相同的情況下，兩個波長光強分布也呈現(xiàn)不同特質(zhì)，研究兩個不同波長光強之間的關(guān)系和變化規(guī)律可獲得更多相分布信息，文獻[10]中對單波長方法和雙波長方法進行了仿真實驗，對比相應(yīng)數(shù)據(jù)后可得出使用雙波長方法對氣泡相分布參數(shù)的預(yù)測平均絕對誤差和均方誤差均小于單波長方法所得。本文通過研究兩個不同波長光強分布的關(guān)系和變化規(guī)律，以獲取光強分布與相分布的相關(guān)信息。結(jié)合上述公式，波長間隔越大折射率相差越大，到達光檢測平面的間距也越大，再結(jié)合對市場現(xiàn)有激光器的調(diào)研，選定波長分別為445nm和635nm的雙波長激光光源。

1.2 梯度提升決策樹原理

梯度提升決策樹（gradient boosting decision tree，GBDT）屬于機器學(xué)習(xí)中的監(jiān)督學(xué)習(xí)算法，是一種迭代決策樹算法，是一個高性能的非線性回歸預(yù)測算法[12]，也可用于實現(xiàn)分類任務(wù)。其主要理念在于通過持續(xù)的迭代過程逐漸降低殘差，并通過梯度優(yōu)化來構(gòu)建多個回歸決策樹。最終，將所有回歸樹的結(jié)論總結(jié)起來形成最終模型，其目的在于同時降低模型的方差和偏差。除了具備樹模型較強的解釋性和對混合模型有效的處理能力，其另外的優(yōu)勢在于強大的預(yù)測性能和出色的穩(wěn)定性。

GBDT 算法的本質(zhì)是一種以決策樹作為基函數(shù)的提升方法，具體如式(13)。

式中，決策樹用T(x;Φm)表示，其中Φm為決策樹的參數(shù)，M為決策樹的個數(shù)。模型算法采用向前分布算法。

首先確定初始提升樹F0(x)=0。

其次，利用向前分布算法，得到第m步的模型如式(14)。

式中，F(xiàn)m-1(x)為當(dāng)前模型。

之后，采用經(jīng)驗風(fēng)險最小化方法確定下一棵決策樹的參數(shù)Φm如式(15)。

2 基于GBDT的相分布測量模型

2.1 雙波長特征提取

為了討論光強分布隨氣泡半徑和位置變化的規(guī)律，使用Trace Pro 光學(xué)仿真軟件建立檢測系統(tǒng)模型，進行氣泡流相分布測量仿真實驗。光源發(fā)出的光線經(jīng)過管道與不同相分布下兩相流折射，通過Trace Pro軟件自帶的光線追跡功能，得到了光檢測平面上的光強分布曲線。為了更好地模擬氣液兩相流的情況，引入均質(zhì)模型，并假設(shè)氣泡在任何位置的尺寸保持不變。將小通道截面內(nèi)的氣相分布情況近似描述為中心位置和尺寸可變的圓形，中心位置坐標為x和y，氣泡半徑為r，記相分布參數(shù)為(x,y,r)。圖3為小通道氣液兩相流模型部分光線的追跡。

圖3 小通道氣液兩相流模型光線追跡

圖4為氣泡不同相分布參數(shù)情況下光強分布對比，對比圖4(a)～(d)可以發(fā)現(xiàn)，在管道中有氣泡時，光強分布曲線出現(xiàn)缺失部分。由圖4(b)可看出，兩波長光強分布曲線的缺失部分兩側(cè)雙波長的間隔存在差異。對比圖4(b)、(c)可以發(fā)現(xiàn)當(dāng)氣泡半徑相同，中心位置不同時，缺失部分位置不同。對比圖4(c)、(d)可以發(fā)現(xiàn)當(dāng)氣泡中心位置相同，氣泡半徑不同時，缺失部分長度不同。綜上可以得出結(jié)論：氣泡的中心位置和氣泡半徑會影響雙波長光強分布曲線的缺失部分長度、缺失部分偏移量和雙波長間隔寬度，因此可以選擇以上三類特征參數(shù)作為特征值，對氣泡相分布參數(shù)進行預(yù)測。

圖4 氣泡不同相分布參數(shù)情況下光強分布對比

借助Trace Pro 軟件計算兩個不同波長光線在不同相分布下，管道截面內(nèi)光線的傳播路徑，得到光檢測平面上不同波長光強分布曲線，用Matlab軟件與Trace Pro 建立通信，實現(xiàn)批量仿真，仿真時設(shè)置的相分布參數(shù)取值見表1，得到了22778 組光強分布與相分布參數(shù)對應(yīng)的數(shù)據(jù)集。對雙波長光強分布數(shù)據(jù)集進行歸一化處理，以全液相時的光強值為最大值，以全氣相時的光強值為最小值，進行最大最小歸一化，圖5 以0.9 為閾值處雙波長，間隔長度計算方式，橫坐標為像素位置，長度為1024。分別計算出445nm和635nm在固定閾值處歸一化光強分布曲線的缺失部分長度ΔL1,445nm、ΔL1,635nm，缺失部分偏移量ΔL2,445nm、ΔL2,445nm，計算得雙波長間隔寬度ΔL3,left、ΔL3,right，最終建立了共22778 組不同相分布與三類雙波長特征值相對應(yīng)的特征值數(shù)據(jù)集。

表1 仿真時相分布參數(shù)取值范圍

圖5 雙波長特征提取

2.2 測量模型的建立

本文借助GBDT來研究相分布每個參數(shù)和雙波長特征值的關(guān)系，選擇適當(dāng)?shù)臎Q策樹模型是影響GBDT模型泛化能力的關(guān)鍵因素之一。為此采用了GirdSearchCV 的方法調(diào)整GBDT 的超參數(shù)，通過在幾個主要超參數(shù)的取值范圍內(nèi)搜索最佳參數(shù)組合，研究了多種參數(shù)組合，使用建立的相分布特征值數(shù)據(jù)集進行測試，將所有驗證結(jié)果取平均即為該組參數(shù)組合的得分。建立模型的首要工作是確定模型的決策樹數(shù)目，用網(wǎng)格搜索和交叉驗證相結(jié)合的方法對其他超參數(shù)進行調(diào)優(yōu)，本文的GBDT超參數(shù)調(diào)優(yōu)設(shè)計見表2。

表2 超參數(shù)及其選擇范圍

GirdSearch 的交叉驗證次數(shù)設(shè)定為5，超參數(shù)的選擇基于驗證數(shù)據(jù)集的擬合度最大化。當(dāng)某一組超參數(shù)的5 次驗證數(shù)據(jù)的平均值最大時，認為這組超參數(shù)是最優(yōu)的。在GBDT 模型中，整個模型的擬合能力一定程度上由決策樹的數(shù)量決定，該數(shù)量還影響著模型的復(fù)雜度。在選擇了超參數(shù)的基礎(chǔ)上，對擁有不同決策樹數(shù)量的模型進行了性能分析。最終，通過對比不同參數(shù)下的模型表現(xiàn)，設(shè)置了以下超參數(shù)：學(xué)習(xí)率為0.2，樹的最大深度為2，劃分所需最小樣本個數(shù)為2，最大迭代次數(shù)為3000。

2.3 仿真數(shù)據(jù)分析

隨機選取上述方法得到的22778組雙波長特征值數(shù)據(jù)集的90%作為訓(xùn)練集對GBDT 模型進行訓(xùn)練，利用訓(xùn)練好的模型對剩余10%的測試集進行相分布預(yù)測。以橫坐標為相分布參數(shù)的真實值，縱坐標為預(yù)測值，虛線為±5%的誤差限，繪制誤差帶圖，當(dāng)點在兩虛線范圍內(nèi)時，說明預(yù)測結(jié)果較好，使用GBDT算法對測試集相分布的預(yù)測結(jié)果如圖6所示，預(yù)測值大部分都在±5%極限誤差內(nèi)。為了驗證測量模型的優(yōu)越性，與已有文獻中使用的BP 神經(jīng)網(wǎng)絡(luò)進行對比，用相分布參數(shù)的均方誤差（mean square error，MSE）作為預(yù)測模型的評價指標。兩種預(yù)測模型的均方誤差對比見表3。通過對比可以發(fā)現(xiàn)，使用GBDT對氣泡相分布預(yù)測時，相較于BP 神經(jīng)網(wǎng)絡(luò)均方誤差減小了33.33%，驗證GBDT算法可實現(xiàn)更高精度測量。

表3 BP神經(jīng)網(wǎng)絡(luò)與GBDT預(yù)測均方誤差對比

圖6 GBDT預(yù)測結(jié)果

3 實驗與結(jié)果

實驗裝置如圖7所示，實驗選用凌云光型號為BFS-U3-51S5C-C 高速CMOS 傳感器，像元大小為3.45μm×3.45μm，采集分辨率為2448×2048，最大幀頻可達200幀/秒，用SPINVIEW調(diào)整相機參數(shù)并采集光強分布數(shù)據(jù)。光源有635nm、520nm 和445nm 三個波長通道，實驗中使用的是波長為445nm 和635nm 的兩通道。采用Zolix 的精密電動直線滑臺PA050，搭配SC300-2B 型號的位移控制器，通過控制電動位移臺推動注射器的速度來控制進入管道的氣相流量。氣相浮子流量計型號分別為LZB-3WB，相應(yīng)量程范圍為10～100mL/min。

圖7 雙波長檢測實驗裝置

用最小二乘法分別提取出635nm和445nm兩單波長光強分布圖，單波長光強分布提取效果如圖8所示。分別提取出兩單波長的光強分布曲線后用上文提到的最大最小歸一化對三類特征值進行特征值提取，由于實驗中CMOS采集的信號具有噪聲，需要對歸一化后的圖像進行平滑處理后再計算出曲線的缺失部分長度、缺失部分偏移量和雙波長間隔寬度。計算出三類特征值將其輸入到由仿真數(shù)據(jù)訓(xùn)練好的模型中，對實驗數(shù)據(jù)的相分布進行預(yù)測。

圖8 單波長提取效果

設(shè)置電動位移臺的驅(qū)動速度為2mm/s，此時氣相浮子流量計測量的氣相速度為30～40mL/min。相機以200幀/秒的速度記錄一段時間內(nèi)氣泡的運動狀態(tài)，從實驗采集的雙波長光強分布中提取出三類特征值，計算出三類特征值后輸入到由22778 組仿真數(shù)據(jù)訓(xùn)練好的GBDT 模型中，得到相分布的預(yù)測值，圖9為根據(jù)相分布預(yù)測值建立的連續(xù)50幀氣泡中心位置運動軌跡，其中圖9(a)為檢測平面X，圖9(b)為檢測平面Y，橫軸為相對位置，單位為mm，縱軸為幀序數(shù)，實現(xiàn)了對氣泡中心位置運動軌跡的追蹤。

圖9 氣泡中心位置運動軌跡

4 結(jié)論

針對已有文獻中提出的雙波長透射法測量氣液兩相流相分布參數(shù)的技術(shù)，提出了用GBDT算法建立測量模型，用仿真實驗得到的光強分布與相分布對應(yīng)的數(shù)據(jù)集，提取三類光強分布特征值作為輸入，對模型進行訓(xùn)練后，用訓(xùn)練好的模型對實驗中氣泡流的相分布進行預(yù)測，具體結(jié)論如下。

（1）用Trace Pro 建立仿真模型，研究雙波長光強分布曲線的三類特征值與相分布參數(shù)的關(guān)系，用GBDT建立測量模型，用仿真數(shù)據(jù)集對模型進行訓(xùn)練，隨機抽取數(shù)據(jù)集的10%模擬實驗測量研究對相分布參數(shù)進行預(yù)測，預(yù)測值大部分都在±5%極限誤差內(nèi)。

（2）搭建了實驗平臺，通過向液柱中注入空氣的方法得到氣泡流，用訓(xùn)練好的模型對實驗采集的氣泡流進行相分布預(yù)測，實現(xiàn)了對氣泡中心位置運動軌跡的追蹤。

（3）本文使用的GBDT算法對氣泡相分布預(yù)測的定位精度優(yōu)于已有文獻，但對于氣泡粒徑的預(yù)測精度還有待提高，未來發(fā)展可以對粒徑大小預(yù)測算法進行優(yōu)化，提高粒徑預(yù)測精度。此外該方法可以對工業(yè)生產(chǎn)過程中兩相流流動狀態(tài)進行在線監(jiān)測。該技術(shù)的測量原理是不同介質(zhì)對不同波長光的折射率不同，出射點也不同，因此還可以將該技術(shù)用于顆粒流和折射率不同的液液兩相流等。