曾思睿,孔明
(中國計量大學計量測試工程學院,浙江 杭州 310018)
氣液兩相流是兩相流中最常見的形式之一,涉及諸多工業領域[1-2],例如化工行業中的吸收塔、蒸餾塔、冷凝器等[3-4]。氣液兩相在空間中的分布有不同形式,這些不同的分布形式會改變力學特性、傳熱特性以及傳質特性等,因此在工業領域和科學研究領域,氣液兩相流的相分布測量技術具有重要意義。該技術能夠支持工業生產過程的流動狀態監測和控制[5],并為理論模型和仿真模型提供可靠數據。憑借這項技術,人們能夠實時監測和掌握工業流體的分布情況,從而實現精確控制和優化工藝。
機器學習在解決復雜非線性回歸問題中表現出優異特性,為兩相流流動參數測量模型的建立提供了一種新思路。國內外學者常用的機器學習主要有隨機森林、BP神經網絡、支持向量機(SVM)、梯度提升決策樹(GBDT)以及深度學習等算法,已應用于氣液兩相流流型識別[6]、氣泡輪廓重建[7]、流量和相含率測量[8]之中,但都未對相分布參數進行研究。
綜合現有的氣液兩相流參數測量研究文獻,對管道界面內氣泡的大小和位置的研究較少。戴振韜[9]提出提取光強分布特征的方法測量相分布,采用線性回歸來擬合相分布測量模型,需要指出的是,該方法使用的是單波長且僅限于水平管道。在此基礎上,曹鳴等[10]提出一種基于雙波長透射法對小通道豎直上升氣泡流進行研究,并利用BP 神經網絡建立測量模型。由于BP 神經網絡具有收斂速度慢、預測能力和訓練能力矛盾等問題[11],本文提出用泛化性能較好、預測精度更高的GBDT算法來訓練模型,并預測實驗中氣泡的相分布。
光學檢測系統如圖1所示,光源發出包含兩種波長的激光光線,光線通過1mm 的狹縫變為片狀光束,再通過分束鏡分成兩束互相垂直的片狀激光,與管道垂直相交,在相交位置形成截面,稱該截面為檢測截面。包含兩種波長的片狀激光在檢測截面內發生折射,管道和管道內兩相流對不同波長光的折射率存在較大差異,且不同大小的氣泡位于兩相流中的不同位置時也會導致激光呈現出不同的折射效果,因此影響傳感器上接收光強分布的因素有波長、氣泡大小和氣泡位置。用高速CMOS傳感器實時記錄氣泡運動到檢測截面的光強分布數據,并將這些數據傳送至計算機。通過分析和處理這些光強分布數據,能夠獲取光強分布與氣泡位置和大小的重要信息,提取這些特征,用大量仿真數據訓練基于GBDT的相分布測量模型,將提取出的光強分布特征輸入模型中,對實驗中氣泡經過檢測截面時的半徑和位置坐標進行預測,根據預測結果對氣泡中心位置運動軌跡進行追蹤。

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

圖2 全液相管道光線傳播路徑
以波長為λ1的一束光線為例,圖中a、c、e、g表示入射角,b、d、f、h表示折射角,波長為λ1的出射光線與x方向夾角為δ1。分別計算出兩波長出射光線與x方向的夾角δ1、δ2,最后計算得出同一位置入射的兩束不同波長的光入射到光檢測平面不同兩點的距離Δδ,計算過程如下。
入射角a可用入射點到管道中心距離l和管道的外徑Rout表示如式(1)。
根據斯涅耳定律,光線經管道外壁折射后的折射角b如式(2)。

結合式(1)~式(4)得入射角e=d,折射角f=c;入射角g=b,折射角h=a。
波長為λ1的光線通過管道后出射光線與x方向的夾角δ1與各角度的幾何關系如式(5)。

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

以上為同一位置發出的一束包含兩個波長光線的計算過程,可以通過計算并記錄上百萬條入射光線的傳播路徑來獲得檢測平面上雙波長光強分布的結果。對于某個特定波長的光線而言,光檢測平面上接收的光強分布可以看作是反射和折射效應的疊加。為了描述光在空間中的光強分布,可以使用復振幅函數。當光線入射角度為θ時,相應的光檢測平面的復振幅函數可以表示如式(12)。
式中,S(θ)為復振幅函數;Sreflect(θ)為反射部分的復振幅函數;Srefract,p(θ)為折射部分的復振幅函數。當p=0時,表示光線在管道內發生全發射,Srefract,0(θ)表示反射光的復振幅函數;當p≥1時,Srefract,p(θ)表示入射光線在管道內經過(p-0)次反射后出射光線的復振幅函數。
根據式(9)~式(11)可得出:由于管道和管道內液體對不同波長光的折射率存在差異,因此同一位置發出不同波長的光線出射到檢測平面上的位置也不同。對上百萬條光線進行追跡,通過式(12),對光線進行疊加,最終光檢測平面上兩個波長的光形成不同的光強分布。當有氣泡經過管道的檢測段時,光線在氣液界面發生折射,導致大量光線被折射到光檢測平面外,檢測平面無法接收到這部分的光強分布,光強分布曲線出現缺失部分,氣液分界面對于不同波長光線的折射能力不同,因此在相分布相同的情況下,兩個波長光強分布也呈現不同特質,研究兩個不同波長光強之間的關系和變化規律可獲得更多相分布信息,文獻[10]中對單波長方法和雙波長方法進行了仿真實驗,對比相應數據后可得出使用雙波長方法對氣泡相分布參數的預測平均絕對誤差和均方誤差均小于單波長方法所得。本文通過研究兩個不同波長光強分布的關系和變化規律,以獲取光強分布與相分布的相關信息。結合上述公式,波長間隔越大折射率相差越大,到達光檢測平面的間距也越大,再結合對市場現有激光器的調研,選定波長分別為445nm和635nm的雙波長激光光源。
梯度提升決策樹(gradient boosting decision tree,GBDT)屬于機器學習中的監督學習算法,是一種迭代決策樹算法,是一個高性能的非線性回歸預測算法[12],也可用于實現分類任務。其主要理念在于通過持續的迭代過程逐漸降低殘差,并通過梯度優化來構建多個回歸決策樹。最終,將所有回歸樹的結論總結起來形成最終模型,其目的在于同時降低模型的方差和偏差。除了具備樹模型較強的解釋性和對混合模型有效的處理能力,其另外的優勢在于強大的預測性能和出色的穩定性。
GBDT 算法的本質是一種以決策樹作為基函數的提升方法,具體如式(13)。
式中,決策樹用T(x;Φm)表示,其中Φm為決策樹的參數,M為決策樹的個數。模型算法采用向前分布算法。
首先確定初始提升樹F0(x)=0。
其次,利用向前分布算法,得到第m步的模型如式(14)。
式中,Fm-1(x)為當前模型。
之后,采用經驗風險最小化方法確定下一棵決策樹的參數Φm如式(15)。

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

圖3 小通道氣液兩相流模型光線追跡
圖4為氣泡不同相分布參數情況下光強分布對比,對比圖4(a)~(d)可以發現,在管道中有氣泡時,光強分布曲線出現缺失部分。由圖4(b)可看出,兩波長光強分布曲線的缺失部分兩側雙波長的間隔存在差異。對比圖4(b)、(c)可以發現當氣泡半徑相同,中心位置不同時,缺失部分位置不同。對比圖4(c)、(d)可以發現當氣泡中心位置相同,氣泡半徑不同時,缺失部分長度不同。綜上可以得出結論:氣泡的中心位置和氣泡半徑會影響雙波長光強分布曲線的缺失部分長度、缺失部分偏移量和雙波長間隔寬度,因此可以選擇以上三類特征參數作為特征值,對氣泡相分布參數進行預測。

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

表1 仿真時相分布參數取值范圍

圖5 雙波長特征提取

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

表2 超參數及其選擇范圍
GirdSearch 的交叉驗證次數設定為5,超參數的選擇基于驗證數據集的擬合度最大化。當某一組超參數的5 次驗證數據的平均值最大時,認為這組超參數是最優的。在GBDT 模型中,整個模型的擬合能力一定程度上由決策樹的數量決定,該數量還影響著模型的復雜度。在選擇了超參數的基礎上,對擁有不同決策樹數量的模型進行了性能分析。最終,通過對比不同參數下的模型表現,設置了以下超參數:學習率為0.2,樹的最大深度為2,劃分所需最小樣本個數為2,最大迭代次數為3000。
隨機選取上述方法得到的22778組雙波長特征值數據集的90%作為訓練集對GBDT 模型進行訓練,利用訓練好的模型對剩余10%的測試集進行相分布預測。以橫坐標為相分布參數的真實值,縱坐標為預測值,虛線為±5%的誤差限,繪制誤差帶圖,當點在兩虛線范圍內時,說明預測結果較好,使用GBDT算法對測試集相分布的預測結果如圖6所示,預測值大部分都在±5%極限誤差內。為了驗證測量模型的優越性,與已有文獻中使用的BP 神經網絡進行對比,用相分布參數的均方誤差(mean square error,MSE)作為預測模型的評價指標。兩種預測模型的均方誤差對比見表3。通過對比可以發現,使用GBDT對氣泡相分布預測時,相較于BP 神經網絡均方誤差減小了33.33%,驗證GBDT算法可實現更高精度測量。

表3 BP神經網絡與GBDT預測均方誤差對比

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

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

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

圖9 氣泡中心位置運動軌跡
針對已有文獻中提出的雙波長透射法測量氣液兩相流相分布參數的技術,提出了用GBDT算法建立測量模型,用仿真實驗得到的光強分布與相分布對應的數據集,提取三類光強分布特征值作為輸入,對模型進行訓練后,用訓練好的模型對實驗中氣泡流的相分布進行預測,具體結論如下。
(1)用Trace Pro 建立仿真模型,研究雙波長光強分布曲線的三類特征值與相分布參數的關系,用GBDT建立測量模型,用仿真數據集對模型進行訓練,隨機抽取數據集的10%模擬實驗測量研究對相分布參數進行預測,預測值大部分都在±5%極限誤差內。
(2)搭建了實驗平臺,通過向液柱中注入空氣的方法得到氣泡流,用訓練好的模型對實驗采集的氣泡流進行相分布預測,實現了對氣泡中心位置運動軌跡的追蹤。
(3)本文使用的GBDT算法對氣泡相分布預測的定位精度優于已有文獻,但對于氣泡粒徑的預測精度還有待提高,未來發展可以對粒徑大小預測算法進行優化,提高粒徑預測精度。此外該方法可以對工業生產過程中兩相流流動狀態進行在線監測。該技術的測量原理是不同介質對不同波長光的折射率不同,出射點也不同,因此還可以將該技術用于顆粒流和折射率不同的液液兩相流等。