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

(1. 北京化工大學(xué) 機(jī)電工程學(xué)院, 北京 100029; 2. 北京石油化工學(xué)院 機(jī)械工程學(xué)院, 北京 102617; 3. 深水油氣管線關(guān)鍵技術(shù)與裝備北京市重點(diǎn)實(shí)驗(yàn)室, 北京 102617)

0 引 言

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

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

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

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

1 實(shí)驗(yàn)裝置及軌跡提取方法

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

如圖1所示，所采用的實(shí)驗(yàn)裝置系專(zhuān)門(mén)搭建用于研究液體射流注入橫向運(yùn)動(dòng)氣流中的射流霧化機(jī)理的實(shí)驗(yàn)平臺(tái)。氣相由空氣壓縮機(jī)以及與之相連的高壓儲(chǔ)氣罐提供，儲(chǔ)氣罐最高儲(chǔ)氣壓力0.6 MPa，采用溫壓補(bǔ)償一體式渦街氣體流量計(jì)(合肥精大，VAFTP-050-DC1-213-1.0-P1.6)測(cè)量通過(guò)的氣體流量。水相為儲(chǔ)液罐提供的超純水。實(shí)驗(yàn)時(shí)，將高壓氣瓶的壓縮空氣通入儲(chǔ)液罐的氣相空間中，以氣壓產(chǎn)生穩(wěn)定流動(dòng)的液體流，再?gòu)纳淞骺谧⑷脒\(yùn)動(dòng)氣流。采用微小量液體流量計(jì)(GICAR112)測(cè)量通過(guò)的液體流量。高速攝像機(jī)(IDT Motion Studio)幀率固定為6000 幀/s，配備N(xiāo)ikon微距鏡頭(AF-S Micro NIKKOR 60 mm f/2.8G ED)。為滿足高速攝像需求，實(shí)驗(yàn)觀測(cè)采用LED光源，并調(diào)整合適的光源角度和光照強(qiáng)度，以獲得較為理想的原始圖像。

圖1 橫流液體射流實(shí)驗(yàn)平臺(tái)示意圖

如圖2(a)所示，高速攝像測(cè)量段為剖面尺寸50 mm×50 mm、長(zhǎng)800 mm的有機(jī)玻璃方管，噴嘴安裝于圖2(a)中紅色矩形所示位置，噴嘴射流口端面與方管內(nèi)壁平齊。實(shí)驗(yàn)采用的噴嘴結(jié)構(gòu)如圖2(b)所示，噴嘴射流口孔徑dj=1 mm，長(zhǎng)徑比L/dj=4。實(shí)驗(yàn)中，參考坐標(biāo)系以射流口中心為原點(diǎn)，氣流方向?yàn)閤軸方向，液體射流方向?yàn)閥軸方向。

圖2 測(cè)量段和噴嘴結(jié)構(gòu)示意圖

1.2 測(cè)量工況

實(shí)驗(yàn)在常溫常壓下進(jìn)行，實(shí)驗(yàn)工況如表1所示。其中，橫向運(yùn)動(dòng)氣流速度為平均速度ug=Qg/A，Qg為測(cè)得的氣體流量，A為氣體管路的橫截面積；液體射流速度也為平均速度uj=Qj/a，Qj為測(cè)得的液體流量，a為液體管路的橫截面積。

表1 實(shí)驗(yàn)工況Table 1 Test conditions

1.3 液體射流軌跡提取方法

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

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

圖3 液體射流破碎結(jié)構(gòu)示意圖

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

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

2 液體射流軌跡提取方法的實(shí)現(xiàn)

2.1 圖像增強(qiáng)處理

圖像增強(qiáng)處理是對(duì)原始彩色圖像進(jìn)行加工，使圖像的某些特征更加突出，從而獲得更為理想的目標(biāo)圖像[16]。圖像處理一般基于灰度圖像，因?yàn)榛叶然蟮膱D像矩陣維數(shù)下降，可以大幅提高運(yùn)算速度。圖像灰度化處理有分量法、最大法、平均法、加權(quán)平均法等多種方法。本文采用MATLAB中的rgb2gray函數(shù)對(duì)原始圖像進(jìn)行灰度化處理，即：

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

0.11B(x,y)

(1)

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

灰度圖像可以通過(guò)灰度直方圖來(lái)反映圖像中的不同灰度等級(jí)分布，其中，橫坐標(biāo)為不同的灰度等級(jí)，縱坐標(biāo)為不同灰度等級(jí)下像素點(diǎn)的個(gè)數(shù)。本文采用直方圖均衡化方法進(jìn)行圖像增強(qiáng)，可使圖像的灰度間距拉開(kāi)或使圖像的灰度分布趨于均勻，增大圖像與背景的反差，使圖像細(xì)節(jié)清晰，從而達(dá)到增強(qiáng)圖像的目的[17]。

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

圖5 圖像增強(qiáng)處理結(jié)果

2.2 圖像閾值分割

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

(2)

(3)

圖6 圖像像素直方圖和KSW二閾值結(jié)果

圖7 閾值分割后圖像

2.3 邊緣輪廓檢測(cè)

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

常用于邊緣輪廓檢測(cè)的方法有Prewitt算子、Roberts算子、Canny算子、Sobel算子等，不同算子是基于圖像具體的突變特征或不同的計(jì)算方法得到。Sobel算子是一階微分算子(離散差分算子)，利用臨近(上下、左右)的像素值(灰度值)來(lái)計(jì)算當(dāng)前像素值(灰度值)。由于圖像邊緣的灰度變化較快，呈現(xiàn)出階躍現(xiàn)象，Sobel算子邊緣檢測(cè)正是利用這一特點(diǎn)，對(duì)圖像中的每個(gè)像素計(jì)算其周?chē)噜徬袼氐牟罘种?，?lái)判斷該像素點(diǎn)是否屬于邊緣[23]。由于本文關(guān)注的射流彎曲軌跡為射流柱外輪廓信息特征，其灰度值與背景存在很大的梯度，因此采用Sobel算子對(duì)閾值分割后的圖像提取邊緣輪廓。

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

圖8 邊緣輪廓檢測(cè)圖像

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

3 軌跡提取方法的適應(yīng)性驗(yàn)證

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

圖9 輪廓點(diǎn)提取后圖像

圖10 多幀圖像的輪廓點(diǎn)提取對(duì)比

圖11 不同破碎模式下圖像處理結(jié)果

根據(jù)本文提出的低速橫流條件下射流破碎軌跡提取方法，對(duì)得到的外輪廓數(shù)據(jù)點(diǎn)進(jìn)行nlinfit非線性擬合，并結(jié)合液氣動(dòng)量通量比q、液體韋伯?dāng)?shù)Wej和雷諾數(shù)Rej等無(wú)量綱參數(shù)，得出射流穿透軌跡經(jīng)驗(yàn)公式，其中擬合函數(shù)選擇：

(4)

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

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

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

另外，利用本文提出的圖像處理方法，對(duì)不同工況下的射流破碎圖像進(jìn)行后期處理，并對(duì)得到的外輪廓數(shù)據(jù)點(diǎn)進(jìn)行nlinfit非線性擬合，得到預(yù)測(cè)低速橫流作用下的射流破碎軌跡經(jīng)驗(yàn)公式：

(5)

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

圖13為不同工況下擬合得到的y值和軌跡提取的y值對(duì)比?？梢钥闯?，利用上述方法得到的結(jié)果和擬合得到的結(jié)果具有較好的一致性。這可以為研究不同工況條件對(duì)射流柱彎曲程度的影響因素以及霧化混合結(jié)構(gòu)設(shè)計(jì)提供參考。

圖13 不同工況下射流軌跡擬合結(jié)果

4 結(jié) 論

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