底吹过程中基于图像处理技术的气泡直径分布特性

陈日健，闫红杰，刘柳，肖俊兵，宋彦坡

(中南大学 能源科学与工程学院，湖南 长沙，410083)

底吹熔池熔炼是指利用底吹气体的喷吹作用，实现熔池内熔体的剧烈搅拌，强化熔池内气液之间的传质传热过程，从而达到提高熔炼反应速率的目的[1]。底吹熔池熔炼技术作为常用的熔池熔炼方法之一，已广泛应用于铅、铜冶炼等冶金工业生产过程。在底吹过程中，气液两相相界面积直接影响气液的传质及反应速率，而获得气泡群在液相中的直径分布能得到气液相界面积，因此，研究底吹过程中气泡群的直径分布规律具有十分重要的意义。目前，研究者针对单孔底吹模型中的实验研究大多集中在单气泡形状[2−3]、上升速度[4−5]、上升轨迹[6]等方面，而对单孔底吹过程中气泡群直径分布的研究较少。大量研究表明[7−10]，图像处理技术具有非接触、操作简单、实用性强、应用广泛等优点，是提取气泡参数的重要手段。对于气泡直径分布的研究，研究者大都关注密集鼓泡过程中气泡群的直径分布情况。LAU等[11]运用图像处理技术研究了“伪二维”模型下多孔鼓泡过程中气泡群直径分布状况，探究了不同气速、不同液面高度下的气泡群直径概率密度分布规律。MAJUMDER等[12]对密集鼓泡塔反应器内的气泡群运动进行了研究，通过图像处理得到气泡群直径分布、气液相界面积与轴向位置、气相速度之间的函数关系。GADALLAH等[13]使用高速摄影仪研究了蜂窝形分布器不同气体流量下的气泡聚并破碎过程，得到了不同的概率密度分布规律，同时分析了Sauter平均直径、最大直径与气体流量的关系。王乐乐等[14]通过图像处理技术提取曝气池模型内的气泡直径分布情况，发现气泡直径分布具有良好的正态分布规律。GAILLARD等[15−17]运用图像处理技术对气泡群直径分布进行了提取与分析。以上研究对象往往是直径均匀的气泡群，与之不同的是单孔底吹过程的气泡直径不一，而对于这类气泡群直径的提取仍存在一定难度。为此，本文通过搭建单孔底吹水模型实验平台，以底吹过程中直径不一的气泡群为对象，通过高速摄影仪拍摄气泡群图像，提出合适的图像处理算法，解决图像中气泡的粘连现象，分类提取气泡直径，得到气泡直径分布情况，并分析气泡Sauter平均直径、最大直径与气体流量的关系。

1 实验装置

图1 实验装置示意图Fig. 1 Schematic diagram of experiment equipment

实验装置示意图如图1所示，包括供气装置、实验主台架、摄像系统三大系统。实验以无色无毒、难溶于水的N2作为气相，水为液相。供气装置中的N2以液态形式存储在高压储气罐中，通过减压阀的减压气化使气态 N2通过连接管道和质量流量控制系统喷入实验容器。其中，气体流量分别通过型号为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的不锈钢注射针，针头处均为平口。摄像系统采用美国 RedlakeTMMotionProX−3 高速摄影仪，实验中设定拍摄频率为1 kHz，图像分辨率可达1 280×1 024像素，可将底吹过程的气泡群图片实时保存至计算机中。实验中光源采用型号为M-TGD-01的LED投光灯，功率为50 W。为了增强拍摄效果，实验中采用硫酸纸作为滤光装置，保证拍摄区域光照均匀分布。实验中获得的气泡图像通过 Matlab软件进行图像处理，可提取气泡直径分布。由于图像处理中得到的气泡直径等几何参数都是以像素为单位，故需对其像素进行标定，将其转换为实际直径。实验中，将1个标有刻度的直尺放在拍摄区域，确定拍摄直尺图像的像素范围，得到10 cm标定长度所对应的像素为1 142个。

2 图像处理算法

图像处理流程如图2所示。首先将拍摄的原始图像进行预处理，将其转换为二值图像，再将处理后的二值图像转换为标注矩阵提取气泡周长及面积，获取气泡圆形度；通过设定圆形度范围进行气泡分类，将其分为单气泡、粘连小气泡、粘连大气泡3类，运用改进分水岭分割算法对粘连气泡进行分割处理，最后分别对单气泡及分割后的粘连气泡进行标注矩阵变换，获取气泡直径分布。

图2 图像处理流程图Fig. 2 Flowchart of image processing

2.1 图像预处理

图像预处理是将高速摄影仪拍摄的原始图像经过一系列预处理，包括灰度转换、图像滤波、背景去除、二值化、孔洞填充、边缘检测等步骤。以拍摄的某一区域图像为例，图像预处理过程如图3所示。

为方便图像连续处理，一般将 RGB格式的原始图像(见图3(a))转换为灰度图像(见图3(b))。由于高速摄影设备或数据传输过程等都可能引入各种各样的随机噪声，需对灰度图像进行滤波处理，此处采用中值滤波处理，如图 3(c)所示。在拍摄过程中，图像不仅包括目标气泡，而且可能包含多余的背景噪声，这种噪声无法通过滤波处理去除，可通过气泡图像与背景图像的减法运算有效消除背景噪声的影响(如图 3(d)所示)，其运算关系式为[18]

式中：I(x,y)为减法运算后图像像素点灰度；O(x,y)为气泡图像像素点灰度；B(x,y)为背景图像像素点灰度。为了保证减法运算后图像像素点灰度大于 0，需取差值的绝对值作为最终输出的图像灰度。

二值图像指的是矩阵中像素点灰度只能为0或1时的图像。阈值分割的主要思路是利用气泡区域与背景区域灰度特性的不同，将图像分成2类不同灰度级别的区域组合，选取合适的阈值，确定像素点气泡区域或背景区域。其中，阈值的选择尤为重要，OTSU[19]提出的一种自适应阈值确定方法能有效实现图像的二值化(如图3(e)所示)，称为最大类间方差法或大津法，具体表达式如下：

式中：F(x,y)为输出图像像素点的灰度；I(x,y)为输入图像像素点的灰度；TOtsu为大津法计算的阈值。

在拍摄过程中由于存在光线问题，气泡中心区域与背景区域像素接近，二值化处理时会出现气泡中心区域灰度小于阈值从而被识别成背景的现象，此时，气泡中心处可能会出现“孔洞”。由于后续过程需统计气泡面积等参数，需对二值图像进行填充处理(见图3(f))，以便能精确地提取所需信息。通过对填充后的气泡图像进行边缘检测，可以得到清晰的气泡边缘图像，如图3(g)所示。

2.2 气泡分类

图3 图像预处理过程Fig. 3 Image preprocessing procedures

为保证单气泡直径提取的准确性，同时考虑粘连气泡直径不一，需对二值化图像中的气泡进行分类处理。在气泡分类前，需将二值化图像转化成标注矩阵，把图像中的不同气泡以不同的整数进行标记，从而确定二值图像中的每个对象，如图 4(a)所示。之后，通过函数计算每个气泡的面积、周长，计算气泡的圆形度。其中，气泡面积为二值化填充图像中每个气泡连通域内像素的总和，周长为边缘检测后每个气泡边缘轮廓线上像素的总和。最后，依据圆形度的不同将其分为单气泡、粘连小气泡、粘连大气泡3类。通过多次尝试后，将图像与圆形度进行对比分析后发现：当圆形度大于0.85时，单气泡图像为图4(b)；当圆形度大于0.50小于0.85时，粘连小气泡图像为图4(c)；当圆形度小于0.50时，粘连大气泡图像为图4(d)。圆形度定义如下：

式中：R为圆形度；P为气泡周长，mm；A为气泡面积，mm2。当R=1时，气泡为圆形；R越接近于1，则气泡形状越趋近于圆形。在分类过程中，单气泡圆形度较大，更接近于 1，容易区分；粘连气泡圆形度较小，且比较接近，需经过不断尝试后，依据处理过程中的经验进行判断，设定适当的圆形度阈值，区分大、小粘连气泡。

图4 气泡分类Fig. 4 Bubble classifications

2.3 粘连气泡分割

分水岭算法[20]是一种经典的图像分割算法，在粘连目标分割方面具有很好的适用性。本文采用改进后的分水岭算法进行气泡分割。首先对目标图像和背景图像进行标记，然后，运用极小值标定技术[21]修改梯度图像，最后对修改后的梯度图像进行分水岭分割得到分割图像。

分水岭算法一般不对原图像进行分割，而是对梯度图像进行处理，这是因为提取梯度能够得到图像的边缘信息，更易于构建分水岭。因此，需要求取灰度图像的梯度，其中，粘连小气泡梯度提取如图 5(a)所示。梯度定义如下：

式中：I∇为梯度图像；grad(.)为梯度运算；f(x,y)为原始图像。

前景标记是指对二值图像中目标气泡区域进行标记处理，根据气泡的形状，使用合适直径的圆形结构元素进行形态学腐蚀操作，使得目标区域的气泡变小，解除气泡的粘连状态，达到前景气泡图像的标记效果，此处粘连小气泡采取像素直径为4的圆形结构进行腐蚀，如图5(b)所示。二值图像腐蚀操作运算定义为

式中：g为前景目标标记图像；O为目标区域；S为结构元素；Θ为腐蚀算符。得到的集合是当结构元素S完全包括在目标区域O中时S中心像素位置的集合。

为得到背景标记图像，需将二值图像进行欧式距离变换[22]。将计算得到距离图像进行分水岭变换，然后求取分割边界，寻找分水岭界限，能够得到背景标记图像。借助极小值标定技术，可修改梯度图像使其在前景标记和背景标记处为局部最小值，修改后的梯度为

图5 粘连气泡分割Fig. 5 Overlapping bubbles segmentation

式中： ∇I′ 为修改后的梯度图像；IMMIN(.)为形态学极小值标定操作；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可以看出：使用分割算法前、后，气泡数目均随气体流量的增大而增大，但使用分割算法后的气泡数目较使用前明显增多，且这种增长趋势随着气体流量的增大更加明显。这是由于使用分割算法之前，气泡之间的粘连使得气泡被识别成单气泡来处理，导致气泡数目无法真实体现；当气体流量为 134 mL/min时，气相呈单气泡状态，无需对气泡进行分割处理。但随着气体流量增大，气泡之间的聚并粘连现象增多，未使用分割算法的气泡数目与实际值之间的误差也会增大，从而影响气泡直径分布的真实性。因此，本文采用改进分水岭分割算法对气泡群图像进行处理，能够减少粘连气泡的影响，使获得的气泡直径分布更接近实际情况。

通过提取气泡直径，得到不同气体流量下以气泡数量表示的气泡直径分布，如图7所示。从图7可以看出：水中小直径气泡居多，中等直径气泡和大直径气泡较少，气泡数量与气泡直径的关系可以近似地看成呈双曲线分布。结合气泡流型变化进行分析可知：当气体流量为267 mL/min时，气泡不再是以单气泡流型进入液相，而是逐渐向无规则状态转变，此时的气泡直径不再均一，开始出现气泡聚并破碎现象，水中出现了少数小气泡；随着气体流量增大，气泡的聚并破碎现象加强，气泡直径依然不均一，气泡聚并产生更多的大气泡，同时，大气泡的破碎也产生更多的小气泡，此时气泡的周期性已经丧失，处于无规则的混沌状态；当气体流量达到800 mL/min时，气泡高度聚并产生更大的气泡，气泡直径可达到13 mm左右，同时破碎程度也更加剧烈，导致小气泡数目急剧增多；随着气体流量继续增加，气相状态将逐渐由鼓泡状态向气体喷射状态转变。

图7 不同气体流量下气泡直径分布Fig. 7 Bubble size distribution under different gas flow rates

3.2 气泡平均直径与最大直径

据已得到的气泡直径可以统计不同气体流量下的气泡直径特性，如表1所示。其中：Sauter平均直径d32是表征气液界面传质过程的重要参数[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平均直径；随着气体流量增加，气泡出现聚并、破碎现象，气泡最大直径增加幅度明显，但最大直径与Sauter平均直径变化的趋势较一致。图9所示为气泡C与气体流量变化的关系。从图9可以看出：当气体流量为134 mL/min时，C为0.89，此时，气泡以单气泡形式上升，气泡最大直径主要取决于气泡的形变；随着气体流量增大，气泡发生聚并破碎，C稳定在0.72左右，此时，气泡最大直径取决于气泡的聚并融合过程。虽然此时气泡破碎导致小气泡的数量增加，但由于其数量量级小，Sauter平均直径主要仍取决于气泡的聚并。因此，当气泡发生聚并破碎时，在一定流量范围内也可通过提取最大直径来估算Sauter平均直径。

图9 C与气体流量的关系Fig. 9 Relationship between gas flow rate and C

4 结论

1)采用本文的图像处理算法能对水中不同直径范围的气泡进行分类提取，同时减少了气泡粘连的影响，处理结果更接近实际情况。

2)当气泡出现聚并破碎现象时，水中小气泡数量较多，中等直径和大直径气泡较少，气泡数量与气泡直径的关系可看成呈双曲线分布；随着气体流量增大，气泡聚并破碎程度更加剧烈，小气泡急剧增多，同时出现更大的气泡。

3)不同气体流量下气泡的 Sauter平均直径与最大直径随着气体流量的增大而增大，且变化趋势基本相同，两者的比值C随着气体流量的增大稳定在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]闫红杰, 赵国建, 刘柳, 等. 静止水中单气泡形状及上升规律的实验研究[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]薛婷, 孟欣东, 张涛. 气液两相流中气泡形态及运动特征参数提取[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]. 北京: 首都经济贸易大学安全与环境工程学院, 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]唐朝晖, 刘金平, 桂卫华, 等. 基于数字图像处理的浮选泡沫速度特征提取及分析[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. 形态学图像分析: 原理与应用[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]. 计算机技术与发展, 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.