基于圖像處理的低速橫流中液體射流軌跡提取方法研究

蘭 天, 孔令真, 陳家慶, 王奎升,*

(1. 北京化工大學 機電工程學院, 北京 100029; 2. 北京石油化工學院 機械工程學院, 北京 102617; 3. 深水油氣管線關鍵技術與裝備北京市重點實驗室, 北京 102617)

0 引 言

液體射流注入橫向運動氣流是一種常見的氣液霧化混合方法，在超燃沖壓發(fā)動機[1]、農用噴霧器[2]和新型管式氣液接觸吸收器[3]中得到了廣泛應用。液體射流柱在橫向運動氣流作用下經歷彎曲、斷裂以及破碎等過程，最終霧化成粒徑更小的液滴群分布于來流氣相中。由于不同工況下射流彎曲軌跡不同，所形成的霧化液滴群在射流方向上的粒徑大小和分布存在差異，這將在很大程度上影響氣液兩相的接觸比表面積，進而影響氣液兩相的混合效果[4]。

針對橫向運動氣流中的射流彎曲軌跡特征及其影響因素，國內外學者采用脈沖陰影、高速攝像等測量方法開展了一系列可視化實驗研究[5-6]，但這些測量方法在原始圖像采集過程中都存在一些問題，例如部分射流區(qū)域“高亮”、圖像背景嘈雜、射流柱邊界拍攝不清晰等，進而對圖像后處理時射流軌跡的提取帶來很大挑戰(zhàn)，因此，需要對采集到的射流破碎圖像進行處理，以獲得更為準確的液體射流軌跡。Song等[7]采用一種簡單的閾值分割方法從高速攝像的原始圖像中提取射流軌跡，由于閾值選擇存在主觀判斷，提取結果不能很好地切合射流柱軌跡輪廓。劉靜[8]采用Photoshop和Origin軟件提取射流紋影圖像邊界線，在處理大量實驗圖像時具有很大難度。Amighi[9]通過脈沖陰影獲取原始圖像，研究了不同閾值算法對圖像分割結果的影響，發(fā)現大多數閾值算法或稀釋圖像而低估邊界，或受背景噪聲影響而高估霧化區(qū)域，較為適合的閾值算法是均值、三角形和Huang算法等基于單個強度分布的算法，但這些算法并未針對每幀圖像的特點來自適應選擇閾值。柴進[10]針對超聲速橫流條件設計了一種基于外輪廓的激波線提取法，認為閾值選取的合理性很大程度上決定了邊緣提取的準確性，通過設定合理閾值，從濾波結果中去掉梯度值較小的部分邊緣，可以精確提取激光陰影圖像中的射流激波線。

橫向運動氣流中液體射流柱在氣液相互作用下經歷復雜的運動和變形過程，不同工況下表現出的射流破碎形態(tài)不同，這使得射流區(qū)域中的液體密度分布也發(fā)生變化，所造成的光照強度在背景和射流柱等不同區(qū)域存在差異[11]。在圖像處理過程中，需要根據這種差異對成像后的圖像進行一系列處理與優(yōu)化，以消除圖像采集過程中存在的問題，適應不同工況下的射流破碎圖像，并最終準確得到描述液體射流破碎的軌跡參數[12]。

軌跡輪廓提取是計算機圖像視覺的重要研究領域。由于目標區(qū)域的多樣化，軌跡輪廓提取很難具有統(tǒng)一的處理方法，眾多的輪廓提取均是基于目標區(qū)域特點提出針對性的處理方式。對于大量的實驗圖像而言，Photoshop、ImageJ等軟件的處理過程較為通用化，且很難做到批量選擇和處理。為此，本文結合圖像處理技術，采用高速攝像儀記錄橫向運動氣流中液體射流破碎的歷經過程，提出低速橫流作用下的液體射流彎曲軌跡提取方法，并結合MATLAB予以編程實現，從而準確得到射流彎曲軌跡，為進一步分析不同工況下射流彎曲軌跡的影響因素奠定基礎。

1 實驗裝置及軌跡提取方法

1.1 實驗裝置

如圖1所示，所采用的實驗裝置系專門搭建用于研究液體射流注入橫向運動氣流中的射流霧化機理的實驗平臺。氣相由空氣壓縮機以及與之相連的高壓儲氣罐提供，儲氣罐最高儲氣壓力0.6 MPa，采用溫壓補償一體式渦街氣體流量計(合肥精大，VAFTP-050-DC1-213-1.0-P1.6)測量通過的氣體流量。水相為儲液罐提供的超純水。實驗時，將高壓氣瓶的壓縮空氣通入儲液罐的氣相空間中，以氣壓產生穩(wěn)定流動的液體流，再從射流口注入運動氣流。采用微小量液體流量計(GICAR112)測量通過的液體流量。高速攝像機(IDT Motion Studio)幀率固定為6000 幀/s，配備Nikon微距鏡頭(AF-S Micro NIKKOR 60 mm f/2.8G ED)。為滿足高速攝像需求，實驗觀測采用LED光源，并調整合適的光源角度和光照強度，以獲得較為理想的原始圖像。

圖1 橫流液體射流實驗平臺示意圖

如圖2(a)所示，高速攝像測量段為剖面尺寸50 mm×50 mm、長800 mm的有機玻璃方管，噴嘴安裝于圖2(a)中紅色矩形所示位置，噴嘴射流口端面與方管內壁平齊。實驗采用的噴嘴結構如圖2(b)所示，噴嘴射流口孔徑dj=1 mm，長徑比L/dj=4。實驗中，參考坐標系以射流口中心為原點，氣流方向為x軸方向，液體射流方向為y軸方向。

圖2 測量段和噴嘴結構示意圖

1.2 測量工況

實驗在常溫常壓下進行，實驗工況如表1所示。其中，橫向運動氣流速度為平均速度ug=Qg/A，Qg為測得的氣體流量，A為氣體管路的橫截面積；液體射流速度也為平均速度uj=Qj/a，Qj為測得的液體流量，a為液體管路的橫截面積。

表1 實驗工況Table 1 Test conditions

1.3 液體射流軌跡提取方法

液體射流注入橫向運動氣流中，經歷復雜的相互作用，射流柱從發(fā)生彎曲、變形的射流初次破碎，到破裂、霧化的射流二次破碎，最終形成較為穩(wěn)定的霧化液滴群分布于來流氣相中。射流初次破碎狀態(tài)關系到霧化液滴群的液滴數量、粒徑尺寸和空間分布[13]。圖3直觀展示了液體射流破碎的主要研究內容，即射流區(qū)域內外邊界、射流變形與破碎、霧化結構等。一般的，眾多學者將射流區(qū)域外輪廓定義為液體射流軌跡，并以此進一步分析影響射流軌跡特征的因素[14]，這是因為外輪廓決定了沿射流方向上射流區(qū)域的最大分布區(qū)間。

因此，在本文研究中將射流區(qū)域外輪廓定義為射流彎曲軌跡，并通過高速攝像機拍攝射流破碎圖像，提出基于MATLAB圖像處理技術[15]的低速橫流情況下的液體射流軌跡提取方法，以定量獲取射流軌跡數據。

圖3 液體射流破碎結構示意圖

液體射流軌跡提取方法如圖4所示。首先進行圖像采集，采用高速攝像機拍攝不同工況下多幀連續(xù)射流破碎圖像，鏡頭光圈調至固定值，調整鏡頭焦距并固定，確保不同工況下足夠的進光量以及相同景深，根據實際拍攝的情況，靈活調整光源角度和光照強度；其次，對原始圖像進行增強處理，即將原始圖像灰度化處理后進行直方圖均衡化，目的是進一步銳化增強所關注的射流破碎區(qū)域；然后，采用最佳直方圖熵法(KSW熵法)及傳統(tǒng)遺傳算法進行圖像閾值分割，該閾值方法可以自適應快速選擇各幀圖像的2個閾值，將射流區(qū)域同背景分割；最后，使用Sobel算子對閾值分割后的圖像進行射流邊緣檢測，結合凸包算法，對每個工況下處理后的100張圖像采取“先平均、后取點”的方式提取邊緣輪廓數據點。上述軌跡提取方法能夠對高速攝像機拍攝的多幀原始圖像進行一系列批量處理，準確得到最終的射流軌跡數據點，并以之擬合經驗公式。

圖4 液體射流軌跡提取方法框圖

2 液體射流軌跡提取方法的實現

2.1 圖像增強處理

圖像增強處理是對原始彩色圖像進行加工，使圖像的某些特征更加突出，從而獲得更為理想的目標圖像[16]。圖像處理一般基于灰度圖像，因為灰度化后的圖像矩陣維數下降，可以大幅提高運算速度。圖像灰度化處理有分量法、最大法、平均法、加權平均法等多種方法。本文采用MATLAB中的rgb2gray函數對原始圖像進行灰度化處理，即：

f(x,y)=0.30R(x,y)+0.59G(x,y)+

0.11B(x,y)

(1)

式中：(x,y)為原始彩色圖像中某一點的像素坐標；f(x,y)為該點彩色圖像灰度化后的灰度值；R、G、B分別為該點原始彩色圖像3個通道的紅色、綠色和藍色分量[17]。

灰度圖像可以通過灰度直方圖來反映圖像中的不同灰度等級分布，其中，橫坐標為不同的灰度等級，縱坐標為不同灰度等級下像素點的個數。本文采用直方圖均衡化方法進行圖像增強，可使圖像的灰度間距拉開或使圖像的灰度分布趨于均勻，增大圖像與背景的反差，使圖像細節(jié)清晰，從而達到增強圖像的目的[17]。

圖5為工況Wej=496.46，q=117.64條件下圖像增強處理前后的對比。圖5(a)為實驗觀測的原始圖像、均衡化后的圖像以及去背景后的圖像(從上至下)，圖5(b)為圖5(a)各圖像對應的灰度直方圖。從圖5可以看出：經過圖像增強處理，本文關注的液體射流區(qū)域的亮度得到很大程度的增強；通過MATLAB去背景函數Imsubtract，原始彩色圖像中大小不一的壁面黏附液滴反射的光斑等噪點也得到了很大程度的去除。從灰度直方圖還可以看出：原始彩色圖像灰度分布于較窄區(qū)間，這會導致圖像細節(jié)不夠清晰、背景與目標區(qū)域混疊嚴重不易區(qū)分；在經過直方圖均衡化后的灰度直方圖中，灰度區(qū)間均衡分布，背景作用減弱，射流區(qū)域增強。

圖5 圖像增強處理結果

2.2 圖像閾值分割

圖像閾值分割是一種用于灰度圖像的圖像分割方法。對于本文實驗拍攝的射流破碎圖像而言，就是通過某種途徑尋找射流破碎圖像中適合的某個閾值，若圖像中某處像素點的灰度低于該閾值，則視為背景區(qū)域，否則視為液體射流區(qū)域的一部分[10,17]。常用的全局閾值選取方法有峰谷法、最小誤差法、最大類間方差法(Otsu)、最佳直方圖熵法(KSW熵法)等。其中，Pun等[18]提出的最佳直方圖熵法是一種以圖像灰度直方圖最佳熵為判斷準則進行圖像閾值分割的算法。從橫向運動氣流中的射流破碎高速攝像圖像可以看出，在不同氣液工況條件下，液體射流柱的破碎劇烈程度不同，這就使得在光源條件不發(fā)生變化的情況下，高速攝像機所拍攝的不同工況下的射流破碎圖像完全不同，因此，需要在圖像處理中依據不同射流破碎圖像自身的灰度特征自適應地選取閾值分割的最佳閾值[19]。對于KSW最佳直方圖熵法二閾值分割情況(即存在2個分割閾值S1和S2)，則有：

(2)

(3)

圖6 圖像像素直方圖和KSW二閾值結果

圖7 閾值分割后圖像

2.3 邊緣輪廓檢測

圖像處理中的邊緣定義為圖像局部特性的不連續(xù)處，如灰度的突變、顏色的突變、紋理結構的突變等。對于灰度圖像，邊緣是圖像中灰度變化較為劇烈的位置，亦即通常而言的信號發(fā)生奇異變化的位置[21]。橫向運動氣流中的液體射流柱存在明顯的軌跡特征[22]，從2.2節(jié)的閾值分割結果可以明顯看出，在射流破碎圖像中，射流柱迎風面沿邊緣切向的像素變化平緩，而法向(垂直于邊緣)的像素變化劇烈。射流柱迎風面、背風面與背景圖像之間具有明顯的灰度梯度，依據這種灰度梯度特性，可以對圖像的邊緣進行檢測和提取，并最終定量得到射流穿透軌跡特征。

常用于邊緣輪廓檢測的方法有Prewitt算子、Roberts算子、Canny算子、Sobel算子等，不同算子是基于圖像具體的突變特征或不同的計算方法得到。Sobel算子是一階微分算子(離散差分算子)，利用臨近(上下、左右)的像素值(灰度值)來計算當前像素值(灰度值)。由于圖像邊緣的灰度變化較快，呈現出階躍現象，Sobel算子邊緣檢測正是利用這一特點，對圖像中的每個像素計算其周圍相鄰像素的差分值，來判斷該像素點是否屬于邊緣[23]。由于本文關注的射流彎曲軌跡為射流柱外輪廓信息特征，其灰度值與背景存在很大的梯度，因此采用Sobel算子對閾值分割后的圖像提取邊緣輪廓。

圖8(a)是利用MATLAB中的Sobel算子對閾值分割后的圖像進行處理的結果；圖8(b)是提取的射流柱外輪廓與原始圖像的對比，外輪廓以紅色顯示，可以發(fā)現射流柱外輪廓能夠最大限度地保留，與原始圖像存在較好的切合性。

圖8 邊緣輪廓檢測圖像

檢測到射流柱外輪廓后，為定量分析不同工況下的射流軌跡特征并探究其影響因素，需要提取射流柱外輪廓的數據點。本文實驗中高速攝像圖像的像素尺寸為1024 pixel×1280 pixel，且射流柱外輪廓由一系列具有灰度值的像素點構成，這些點構成一個點集，因此，可以結合凸包概念(Convex Hull，將給定二維平面上的點集最外層的點連接為凸多邊型，以包含點集中所有的點[24])得到包圍這些點集的最外層數據點，其連線即為射流軌跡曲線。液體射流柱在橫向運動氣流中具有明顯的軌跡特征，經上述處理后，可以得到更為清晰的射流破碎輪廓。本文結合MATLAB中的凸包算法Convhull及邊緣點優(yōu)化和稠密化，得到射流柱外輪廓數據點集。圖9中的虛線為射流柱外輪廓數據點集的提取結果，圖10為多幀圖像的輪廓點提取對比，可以看出，提取出的數據點能夠很好地切合射流柱外輪廓。

3 軌跡提取方法的適應性驗證

液體射流注入橫向運動氣流中的射流破碎是一個氣液相互作用的復雜過程。Vich[25]、Birouk[26]等研究了低速橫流作用下的液體射流破碎模式特征，認為存在弱橫流作用破碎、拱形破碎和袋式破碎等3種典型射流破碎模式。這是由于在橫流氣速變化的情況下，運動氣流氣動力不同，使得射流柱破碎狀態(tài)存在差異。圖5和圖7～9是典型的袋式破碎模式，圖11(a)是弱橫流作用破碎模式，圖11(b)是拱形破碎模式(上為原始圖像，下為輪廓檢測提取后的圖像)。從圖中可以看出，射流破碎軌跡和破碎狀態(tài)在不同工況下存在差異，邊緣輪廓提取后的圖像與原始圖像能夠很好地吻合。圖11(a)的弱橫流作用破碎模式中，橫向氣流作用影響較弱，射流柱彎曲程度較好，在通過閾值處理和邊緣檢測后對射流柱外輪廓提取數據點時，數據點集能夠很好地切合射流柱外輪廓。在圖11(b)的拱形破碎模式中，射流柱在遠離射流孔的位置存在多處褶皺，其提取的外輪廓為射流區(qū)域最外層的邊界點集。

圖9 輪廓點提取后圖像

圖10 多幀圖像的輪廓點提取對比

圖11 不同破碎模式下圖像處理結果

根據本文提出的低速橫流條件下射流破碎軌跡提取方法，對得到的外輪廓數據點進行nlinfit非線性擬合，并結合液氣動量通量比q、液體韋伯數Wej和雷諾數Rej等無量綱參數，得出射流穿透軌跡經驗公式，其中擬合函數選擇：

(4)

式中，a=[a(1)a(2)a(3)a(4)a(5)]為擬合系數，dj=1 mm。

圖12為q=41.33、Wej=60.82工況下采用上述圖像處理方法得到的數據點非線性擬合結果。圖12(a)為擬合曲線與提取數據點的對比，圖12(b)為擬合得到的y值與軌跡提取的y值的對比，從圖中可以看出，通過該擬合函數得到的預測值較為準確，其中擬合系數a=[2.6404 -1.6543 0.4655 -0.0797 -0.0657]。

圖12 q=41.33、Wej=60.82工況下射流軌跡擬合結果

另外，利用本文提出的圖像處理方法，對不同工況下的射流破碎圖像進行后期處理，并對得到的外輪廓數據點進行nlinfit非線性擬合，得到預測低速橫流作用下的射流破碎軌跡經驗公式：

(5)

式中，擬合系數a系取平均后圓整得到。

圖13為不同工況下擬合得到的y值和軌跡提取的y值對比。可以看出，利用上述方法得到的結果和擬合得到的結果具有較好的一致性。這可以為研究不同工況條件對射流柱彎曲程度的影響因素以及霧化混合結構設計提供參考。

圖13 不同工況下射流軌跡擬合結果

4 結 論

本文提出了基于圖像處理技術的射流軌跡提取方法，該方法結合直方圖均衡化、最佳直方圖熵法、Sobel算子和凸包算法等，能夠有效適應低速橫流作用下不同射流破碎模式的射流軌跡提取，對提取后的彎曲軌跡數據點進行非線性擬合，進而獲得準確預測射流彎曲軌跡的經驗公式。該方法可以對高速攝像機拍攝的原始彩色圖像進行一系列批量處理，從而更加方便地討論液體射流軌跡的影響因素，為氣液霧化混合結構的設計提供參考。