陳日健，閆紅杰，劉柳，肖俊兵，宋彥坡

(中南大學 能源科學與工程學院，湖南 長沙，410083)

底吹熔池熔煉是指利用底吹氣體的噴吹作用，實現(xiàn)熔池內熔體的劇烈攪拌，強化熔池內氣液之間的傳質傳熱過程，從而達到提高熔煉反應速率的目的[1]。底吹熔池熔煉技術作為常用的熔池熔煉方法之一，已廣泛應用于鉛、銅冶煉等冶金工業(yè)生產過程。在底吹過程中，氣液兩相相界面積直接影響氣液的傳質及反應速率，而獲得氣泡群在液相中的直徑分布能得到氣液相界面積，因此，研究底吹過程中氣泡群的直徑分布規(guī)律具有十分重要的意義。目前，研究者針對單孔底吹模型中的實驗研究大多集中在單氣泡形狀[2?3]、上升速度[4?5]、上升軌跡[6]等方面，而對單孔底吹過程中氣泡群直徑分布的研究較少。大量研究表明[7?10]，圖像處理技術具有非接觸、操作簡單、實用性強、應用廣泛等優(yōu)點，是提取氣泡參數的重要手段。對于氣泡直徑分布的研究，研究者大都關注密集鼓泡過程中氣泡群的直徑分布情況。LAU等[11]運用圖像處理技術研究了“偽二維”模型下多孔鼓泡過程中氣泡群直徑分布狀況，探究了不同氣速、不同液面高度下的氣泡群直徑概率密度分布規(guī)律。MAJUMDER等[12]對密集鼓泡塔反應器內的氣泡群運動進行了研究，通過圖像處理得到氣泡群直徑分布、氣液相界面積與軸向位置、氣相速度之間的函數關系。GADALLAH等[13]使用高速攝影儀研究了蜂窩形分布器不同氣體流量下的氣泡聚并破碎過程，得到了不同的概率密度分布規(guī)律，同時分析了Sauter平均直徑、最大直徑與氣體流量的關系。王樂樂等[14]通過圖像處理技術提取曝氣池模型內的氣泡直徑分布情況，發(fā)現(xiàn)氣泡直徑分布具有良好的正態(tài)分布規(guī)律。GAILLARD等[15?17]運用圖像處理技術對氣泡群直徑分布進行了提取與分析。以上研究對象往往是直徑均勻的氣泡群，與之不同的是單孔底吹過程的氣泡直徑不一，而對于這類氣泡群直徑的提取仍存在一定難度。為此，本文通過搭建單孔底吹水模型實驗平臺，以底吹過程中直徑不一的氣泡群為對象，通過高速攝影儀拍攝氣泡群圖像，提出合適的圖像處理算法，解決圖像中氣泡的粘連現(xiàn)象，分類提取氣泡直徑，得到氣泡直徑分布情況，并分析氣泡Sauter平均直徑、最大直徑與氣體流量的關系。

1 實驗裝置

圖1 實驗裝置示意圖Fig. 1 Schematic diagram of experiment equipment

實驗裝置示意圖如圖1所示，包括供氣裝置、實驗主臺架、攝像系統(tǒng)三大系統(tǒng)。實驗以無色無毒、難溶于水的N2作為氣相，水為液相。供氣裝置中的N2以液態(tài)形式存儲在高壓儲氣罐中，通過減壓閥的減壓氣化使氣態(tài) N2通過連接管道和質量流量控制系統(tǒng)噴入實驗容器。其中，氣體流量分別通過型號為SEVEN-STAR?D07?19B 質 量 流 量 控 制 器 和D07?19BM 質量流量計進行控制和檢測。實驗主臺架采用1個長為0.15 m、寬為0.15 m、高為0.50 m的長方體有機玻璃容器，其中，底部壁面厚度為10 mm，側面壁面厚度均為 8 mm。在該玻璃容器的底部中心處開1個直徑為30 mm圓孔，將橡膠塞塞入，并在橡膠塞的中心處插入內徑為0.6 mm的不銹鋼注射針，針頭處均為平口。攝像系統(tǒng)采用美國 RedlakeTMMotionProX?3 高速攝影儀，實驗中設定拍攝頻率為1 kHz，圖像分辨率可達1 280×1 024像素，可將底吹過程的氣泡群圖片實時保存至計算機中。實驗中光源采用型號為M-TGD-01的LED投光燈，功率為50 W。為了增強拍攝效果，實驗中采用硫酸紙作為濾光裝置，保證拍攝區(qū)域光照均勻分布。實驗中獲得的氣泡圖像通過 Matlab軟件進行圖像處理，可提取氣泡直徑分布。由于圖像處理中得到的氣泡直徑等幾何參數都是以像素為單位，故需對其像素進行標定，將其轉換為實際直徑。實驗中，將1個標有刻度的直尺放在拍攝區(qū)域，確定拍攝直尺圖像的像素范圍，得到10 cm標定長度所對應的像素為1 142個。

2 圖像處理算法

圖像處理流程如圖2所示。首先將拍攝的原始圖像進行預處理，將其轉換為二值圖像，再將處理后的二值圖像轉換為標注矩陣提取氣泡周長及面積，獲取氣泡圓形度；通過設定圓形度范圍進行氣泡分類，將其分為單氣泡、粘連小氣泡、粘連大氣泡3類，運用改進分水嶺分割算法對粘連氣泡進行分割處理，最后分別對單氣泡及分割后的粘連氣泡進行標注矩陣變換，獲取氣泡直徑分布。

圖2 圖像處理流程圖Fig. 2 Flowchart of image processing

2.1 圖像預處理

圖像預處理是將高速攝影儀拍攝的原始圖像經過一系列預處理，包括灰度轉換、圖像濾波、背景去除、二值化、孔洞填充、邊緣檢測等步驟。以拍攝的某一區(qū)域圖像為例，圖像預處理過程如圖3所示。

為方便圖像連續(xù)處理，一般將 RGB格式的原始圖像(見圖3(a))轉換為灰度圖像(見圖3(b))。由于高速攝影設備或數據傳輸過程等都可能引入各種各樣的隨機噪聲，需對灰度圖像進行濾波處理，此處采用中值濾波處理，如圖 3(c)所示。在拍攝過程中，圖像不僅包括目標氣泡，而且可能包含多余的背景噪聲，這種噪聲無法通過濾波處理去除，可通過氣泡圖像與背景圖像的減法運算有效消除背景噪聲的影響(如圖 3(d)所示)，其運算關系式為[18]

式中：I(x,y)為減法運算后圖像像素點灰度；O(x,y)為氣泡圖像像素點灰度；B(x,y)為背景圖像像素點灰度。為了保證減法運算后圖像像素點灰度大于 0，需取差值的絕對值作為最終輸出的圖像灰度。

二值圖像指的是矩陣中像素點灰度只能為0或1時的圖像。閾值分割的主要思路是利用氣泡區(qū)域與背景區(qū)域灰度特性的不同，將圖像分成2類不同灰度級別的區(qū)域組合，選取合適的閾值，確定像素點氣泡區(qū)域或背景區(qū)域。其中，閾值的選擇尤為重要，OTSU[19]提出的一種自適應閾值確定方法能有效實現(xiàn)圖像的二值化(如圖3(e)所示)，稱為最大類間方差法或大津法，具體表達式如下：

式中：F(x,y)為輸出圖像像素點的灰度；I(x,y)為輸入圖像像素點的灰度；TOtsu為大津法計算的閾值。

在拍攝過程中由于存在光線問題，氣泡中心區(qū)域與背景區(qū)域像素接近，二值化處理時會出現(xiàn)氣泡中心區(qū)域灰度小于閾值從而被識別成背景的現(xiàn)象，此時，氣泡中心處可能會出現(xiàn)“孔洞”。由于后續(xù)過程需統(tǒng)計氣泡面積等參數，需對二值圖像進行填充處理(見圖3(f))，以便能精確地提取所需信息。通過對填充后的氣泡圖像進行邊緣檢測，可以得到清晰的氣泡邊緣圖像，如圖3(g)所示。

2.2 氣泡分類

圖3 圖像預處理過程Fig. 3 Image preprocessing procedures

為保證單氣泡直徑提取的準確性，同時考慮粘連氣泡直徑不一，需對二值化圖像中的氣泡進行分類處理。在氣泡分類前，需將二值化圖像轉化成標注矩陣，把圖像中的不同氣泡以不同的整數進行標記，從而確定二值圖像中的每個對象，如圖 4(a)所示。之后，通過函數計算每個氣泡的面積、周長，計算氣泡的圓形度。其中，氣泡面積為二值化填充圖像中每個氣泡連通域內像素的總和，周長為邊緣檢測后每個氣泡邊緣輪廓線上像素的總和。最后，依據圓形度的不同將其分為單氣泡、粘連小氣泡、粘連大氣泡3類。通過多次嘗試后，將圖像與圓形度進行對比分析后發(fā)現(xiàn)：當圓形度大于0.85時，單氣泡圖像為圖4(b)；當圓形度大于0.50小于0.85時，粘連小氣泡圖像為圖4(c)；當圓形度小于0.50時，粘連大氣泡圖像為圖4(d)。圓形度定義如下：

式中：R為圓形度；P為氣泡周長，mm；A為氣泡面積，mm2。當R=1時，氣泡為圓形；R越接近于1，則氣泡形狀越趨近于圓形。在分類過程中，單氣泡圓形度較大，更接近于 1，容易區(qū)分；粘連氣泡圓形度較小，且比較接近，需經過不斷嘗試后，依據處理過程中的經驗進行判斷，設定適當的圓形度閾值，區(qū)分大、小粘連氣泡。

圖4 氣泡分類Fig. 4 Bubble classifications

2.3 粘連氣泡分割

分水嶺算法[20]是一種經典的圖像分割算法，在粘連目標分割方面具有很好的適用性。本文采用改進后的分水嶺算法進行氣泡分割。首先對目標圖像和背景圖像進行標記，然后，運用極小值標定技術[21]修改梯度圖像，最后對修改后的梯度圖像進行分水嶺分割得到分割圖像。

分水嶺算法一般不對原圖像進行分割，而是對梯度圖像進行處理，這是因為提取梯度能夠得到圖像的邊緣信息，更易于構建分水嶺。因此，需要求取灰度圖像的梯度，其中，粘連小氣泡梯度提取如圖 5(a)所示。梯度定義如下：

式中：I?為梯度圖像；grad(.)為梯度運算；f(x,y)為原始圖像。

前景標記是指對二值圖像中目標氣泡區(qū)域進行標記處理，根據氣泡的形狀，使用合適直徑的圓形結構元素進行形態(tài)學腐蝕操作，使得目標區(qū)域的氣泡變小，解除氣泡的粘連狀態(tài)，達到前景氣泡圖像的標記效果，此處粘連小氣泡采取像素直徑為4的圓形結構進行腐蝕，如圖5(b)所示。二值圖像腐蝕操作運算定義為

式中：g為前景目標標記圖像；O為目標區(qū)域；S為結構元素；Θ為腐蝕算符。得到的集合是當結構元素S完全包括在目標區(qū)域O中時S中心像素位置的集合。

為得到背景標記圖像，需將二值圖像進行歐式距離變換[22]。將計算得到距離圖像進行分水嶺變換，然后求取分割邊界，尋找分水嶺界限，能夠得到背景標記圖像。借助極小值標定技術，可修改梯度圖像使其在前景標記和背景標記處為局部最小值，修改后的梯度為

圖5 粘連氣泡分割Fig. 5 Overlapping bubbles segmentation

式中： ?I′ 為修改后的梯度圖像；IMMIN(.)為形態(tài)學極小值標定操作；m為背景標記圖像；g為前景標記圖像。

最后對修改后的梯度圖像進行分水嶺變換，完成圖像的分割。大、小粘連氣泡通過改進分水嶺分割算法分割后的圖像分別如圖5(c)和圖5(d)所示。

2.4 氣泡直徑提取

通過圖像預處理、氣泡分類、粘連氣泡分割步驟可得到單氣泡、粘連大氣泡、粘連小氣泡分割圖像，分別將其轉換為標注矩陣后提取氣泡的面積，從而得到氣泡的等效直徑。假設氣泡的二維圖像均為圓形，則具有相同面積氣泡的等效直徑定義為

式中：de為氣泡的等效直徑，mm；A為氣泡面積，mm2。

3 結果與分析

3.1 氣泡直徑分布

圖6 使用分割算法前、后氣泡數目對比Fig. 6 Comparisom of number of bubbles with and without segmentation algorithm

使用分割算法前、后氣泡數目的對比如圖6所示。從圖6可以看出：使用分割算法前、后，氣泡數目均隨氣體流量的增大而增大，但使用分割算法后的氣泡數目較使用前明顯增多，且這種增長趨勢隨著氣體流量的增大更加明顯。這是由于使用分割算法之前，氣泡之間的粘連使得氣泡被識別成單氣泡來處理，導致氣泡數目無法真實體現(xiàn)；當氣體流量為 134 mL/min時，氣相呈單氣泡狀態(tài)，無需對氣泡進行分割處理。但隨著氣體流量增大，氣泡之間的聚并粘連現(xiàn)象增多，未使用分割算法的氣泡數目與實際值之間的誤差也會增大，從而影響氣泡直徑分布的真實性。因此，本文采用改進分水嶺分割算法對氣泡群圖像進行處理，能夠減少粘連氣泡的影響，使獲得的氣泡直徑分布更接近實際情況。

通過提取氣泡直徑，得到不同氣體流量下以氣泡數量表示的氣泡直徑分布，如圖7所示。從圖7可以看出：水中小直徑氣泡居多，中等直徑氣泡和大直徑氣泡較少，氣泡數量與氣泡直徑的關系可以近似地看成呈雙曲線分布。結合氣泡流型變化進行分析可知：當氣體流量為267 mL/min時，氣泡不再是以單氣泡流型進入液相，而是逐漸向無規(guī)則狀態(tài)轉變，此時的氣泡直徑不再均一，開始出現(xiàn)氣泡聚并破碎現(xiàn)象，水中出現(xiàn)了少數小氣泡；隨著氣體流量增大，氣泡的聚并破碎現(xiàn)象加強，氣泡直徑依然不均一，氣泡聚并產生更多的大氣泡，同時，大氣泡的破碎也產生更多的小氣泡，此時氣泡的周期性已經喪失，處于無規(guī)則的混沌狀態(tài)；當氣體流量達到800 mL/min時，氣泡高度聚并產生更大的氣泡，氣泡直徑可達到13 mm左右，同時破碎程度也更加劇烈，導致小氣泡數目急劇增多；隨著氣體流量繼續(xù)增加，氣相狀態(tài)將逐漸由鼓泡狀態(tài)向氣體噴射狀態(tài)轉變。

圖7 不同氣體流量下氣泡直徑分布Fig. 7 Bubble size distribution under different gas flow rates

3.2 氣泡平均直徑與最大直徑

據已得到的氣泡直徑可以統(tǒng)計不同氣體流量下的氣泡直徑特性，如表1所示。其中：Sauter平均直徑d32是表征氣液界面?zhèn)髻|過程的重要參數[23]，其值為具有相同總體積與總表面積比值的氣泡的平均直徑；最大直徑dmax為氣泡群中最大氣泡的直徑；C為Sauter平均直徑d32與最大直徑dmax的比值。Sauter平均直徑定義式如下：

式中：d32為Sauter平均直徑，mm；di為單個氣泡的直徑，mm；N為總氣泡數量。

表1 不同氣體流量下的氣泡直徑特性Table 1 Bubble size characteristics under different gas flow rates

圖8 dmax和d32與氣體流量的關系Fig. 8 Relationship among gas flow rate, dmax and d32

圖8所示為氣泡Sauter平均直徑d32、最大直徑dmax與氣體流量變化的關系。由圖8可知：氣泡Sauter平均直徑與最大直徑均隨著氣體流量的增大而增大；當氣體流量為 134 mL/min時，氣相以單氣泡形式存在，最大直徑略大于Sauter平均直徑；隨著氣體流量增加，氣泡出現(xiàn)聚并、破碎現(xiàn)象，氣泡最大直徑增加幅度明顯，但最大直徑與Sauter平均直徑變化的趨勢較一致。圖9所示為氣泡C與氣體流量變化的關系。從圖9可以看出：當氣體流量為134 mL/min時，C為0.89，此時，氣泡以單氣泡形式上升，氣泡最大直徑主要取決于氣泡的形變；隨著氣體流量增大，氣泡發(fā)生聚并破碎，C穩(wěn)定在0.72左右，此時，氣泡最大直徑取決于氣泡的聚并融合過程。雖然此時氣泡破碎導致小氣泡的數量增加，但由于其數量量級小，Sauter平均直徑主要仍取決于氣泡的聚并。因此，當氣泡發(fā)生聚并破碎時，在一定流量范圍內也可通過提取最大直徑來估算Sauter平均直徑。

圖9 C與氣體流量的關系Fig. 9 Relationship between gas flow rate and C

4 結論

1)采用本文的圖像處理算法能對水中不同直徑范圍的氣泡進行分類提取，同時減少了氣泡粘連的影響，處理結果更接近實際情況。

2)當氣泡出現(xiàn)聚并破碎現(xiàn)象時，水中小氣泡數量較多，中等直徑和大直徑氣泡較少，氣泡數量與氣泡直徑的關系可看成呈雙曲線分布；隨著氣體流量增大，氣泡聚并破碎程度更加劇烈，小氣泡急劇增多，同時出現(xiàn)更大的氣泡。

3)不同氣體流量下氣泡的 Sauter平均直徑與最大直徑隨著氣體流量的增大而增大，且變化趨勢基本相同，兩者的比值C隨著氣體流量的增大穩(wěn)定在0.72左右。

[1]劉柳, 閆紅杰, 周孑民, 等. 氧氣底吹銅熔池熔煉過程的機理及產物的微觀分析[J]. 中國有色金屬學報, 2012, 22(7):2116?2124.LIU Liu, YAN Hongjie, ZHOU Jiemin, et al. Mechanism of copper smelting process by oxygen bottom blowing and microanalysis of smelting products[J]. The Chinese Journal of Nonferrous Metals, 2012, 22(7): 2116?2124.

[2]SURESH M, MANI A. Experimental studies on bubble characteristics for R134a-DMF bubble absorber[J]. Experimental Thermal and Fluid Science, 2012, 39: 79?89.

[3]LIU L, YAN H, ZHAO G. Experimental studies on the shape and motion of air bubbles in viscous liquids[J]. Experimental Thermal and Fluid Science, 2015, 62: 109?121.

[4]閆紅杰, 趙國建, 劉柳, 等. 靜止水中單氣泡形狀及上升規(guī)律的實驗研究[J]. 中南大學學報(自然科學版), 2016, 47(7):2513?2520.YAN Hongjie, ZHAO Guojian, LIU Liu, et al. Experimental study on shape and rising behavior of single bubble in stagnant water[J]. Journal of Central South University(Science and Technology), 2016, 47(7): 2513?2520.

[5]薛婷, 孟欣東, 張濤. 氣液兩相流中氣泡形態(tài)及運動特征參數提取[J]. 光電子.激光, 2010, 21(8): 1218?1221.XUE Ting, MENG Xindong, ZHANG Tao. Extraction of bubble shape and motion feature parameters in the gas-liquid two-phase flow[J]. Journal of Optoelectronics. Laser, 2010, 21(8):1218?1221.

[6]AMIRNIA S, BRUYN J R D, BERGOUGNOU M A, et al.Continuous rise velocity of air bubbles in non-Newtonian biopolymer solutions[J]. Chemical Engineering Science, 2013,94: 60?68.

[7]李娜. 氣泡運動的數字圖像測量方法研究[D]. 北京: 首都經濟貿易大學安全與環(huán)境工程學院, 2011: 42?49.LI Na. Research on digital image measurement method of bubble motion[D]. Beijing: Capital University of Economics and Business. Safety and Environmental Engineering College, 2011:42?49.

[8]HU Y, MENG X, SUN W. A novel quantificational assessment method of frothers effect on bubble characteristics[J]. Journal of Central South University of Technology, 2011, 18(6):1910?1916.

[9]AVACHAT S, ZHOU M. High-speed digital imaging and computational modeling of dynamic failure in composite structures subjected to underwater impulsive loads[J].International Journal of Impact Engineering, 2015, 77: 147?165.

[10]唐朝暉, 劉金平, 桂衛(wèi)華, 等. 基于數字圖像處理的浮選泡沫速度特征提取及分析[J]. 中南大學學報(自然科學版), 2009,40(6): 1616?1622.TANG Chaohui, LIU Jinping, GUI Weihua, et al. Froth bubbles speed characteristic extraction and analysis based on digital image processing[J]. Journal of Central South University(Science and Technology), 2009, 40(6): 1616?1621.

[11]LAU Y M, DEEN N G, KUIPERS J A M. Development of an image measurement technique for size distribution in dense bubbly flows[J]. Chemical Engineering Science, 2013, 94:20?29.

[12]MAJUMDER S K, KUNDU G, MUKHERJEE D. Bubble size distribution and gas-liquid interfacial area in a modified downflow bubble column[J]. Chemical Engineering Journal,2006, 122(1/2): 1?10.

[13]GADALLAH A H, SIDDIQUI K. Bubble breakup in co-current upward flowing liquid using honeycomb monolith breaker[J].Chemical Engineering Science, 2015, 131: 22?40.

[14]王樂樂, 雍曉東, 李然, 等. 圖像處理技術在氣泡特性研究中的應用[J]. 四川大學學報(工程科學版), 2012(S2): 188?192.WANG Lele, YONG Xiaodong, LI Ran, et al. Image processing technology in the application of bubble characteristics research[J]. Journal of Sichuan University (Engineering Science Edition), 2012(S2): 188?192.

[15]GAILLARD T, HONOREZ C, JUMEAU M, et al. A simple technique for the automation of bubble size measurements[J].Colloids and Surfaces A: Physicochemical and Engineering Aspects, 2015, 473: 68?74.

[16]BESAGNI G, INZOLI F. Bubble size distributions and shapes in annular gap bubble column[J]. Experimental Thermal and Fluid Science, 2016, 74: 27?48.

[17]KARN A, SHAO S, ARNDT R E A, et al. Bubble coalescence and breakup in turbulent bubbly wake of a ventilated hydrofoil[J]. Experimental Thermal and Fluid Science, 2016, 70:397?407.

[18]LIU L, YAN H, ZHAO G, et al. Experimental studies on the terminal velocity of air bubbles in water and glycerol aqueous solution[J]. Experimental Thermal and Fluid Science, 2016, 78:254?265.

[19]OTSU N. A threshold selection method from gray-level histograms[J]. IEEE Trans SMC, 1979, 9(1): 62?66.

[20]VINCENT L, SOILLE P. Watersheds in digital spaces: an efficient algorithm based on immersion simulations[J]. IEEE Trans on Pattern Analysis and Machine Intelligence, 1991, 13(6):583?598.

[21]SOILLE P. 形態(tài)學圖像分析: 原理與應用[M]. 2版. 北京: 清華大學出版社, 2008: 208?211.SOILLE P. Morphological image analysis: principles and applications[M]. 2nd ed. Beijing: Tsinghua University Press,2008: 208?211.

[22]戴丹. 基于改進分水嶺算法的粘連顆粒圖像分割[J]. 計算機技術與發(fā)展, 2013, 23(3): 19?22.DAI Dan. Image segmentation of adhering particle based on improved watershed algorithm[J]. Computer Technology and Development, 2013, 23(3): 19?22.

[23]FERREIRA A, PEREIRA G, TEIXEIRA J A, et al. Statistical tool combined with image analysis to characterize hydrodynamics and mass transfer in a bubble column[J].Chemical Engineering Journal, 2012, 180: 216?228.