杨志龙，田文斌，张珍，王志英，王一伟
（1 中国农业大学工学院，北京 100083；
2 石家庄铁道大学省部共建交通工程结构力学行为与系统安全国家重点实验室 石家庄 050043；
3 中国科学院力学研究所流固耦合系统力学重点实验室，北京 100190）
气液两相流是多相流领域中备受关注的研究方向，广泛存在于水利、化工、能源工程、航空航天等多种工业领域[1-3]。当气相以离散相形式分布于液体连续相中时，多以气泡的形态出现[4]。气泡-液体系统是气液两相体系中最为重要的一种流动形式。化工反应中利用气液两相流作为界面传质的载体有利于加快界面之间物质与热量的交换[5]，例如在鼓泡塔反应器中通入气体产生气泡增大相间接触面积以提高反应效率[6]。气泡的尺寸等参数影响着相间传质系数，决定着反应的传质传热效率[7-8]。而在水下航行体的研究中，可通过主动通气手段改变航体近壁面的流体状态获取预期的力学效应[9]。例如，在航行体表面开设气孔通入气体，产生气泡进而融合形成气腔可有效调节近壁面液体流动状态，从而实现减阻降载，提高航行稳定性[10]。因此，准确快速提取气泡轮廓和特征参数，是各个领域中气泡相关流动机制揭示的基础。
光学成像的方法作为非侵入式的测量方法，其优点是不会干扰流场，可以获得更加准确的数据。其中高速摄像技术结合数字图像处理的方法是目前气泡测量领域中的一种主流技术[11]。高速摄像能以大于1000 帧/s 的时间分辨率拍摄百万级像素的气泡图像，通过数字图像处理可以直观地测量出气泡的参数，如尺寸、位置分布和形状等信息。在拍摄的气泡图像中，由于拍摄条件的限制，图像中会存在多种问题影响测量的准确度，其中相互重叠粘连的气泡在分析过程中造成的误差影响最大。对于重叠粘连气泡的分割仍然是气泡图像处理领域中的一个难点，针对重叠粘连物体的分割，国内外学者提出了多种基于数字图像处理的处理方法。目前的主要方法按原理可以分为凹点检测方法和基于标记的分水岭算法。凹点检测方法通过气泡边缘曲线的曲率信息寻找凹点，并利用凹点进行分割。潘志成等[12]采用k余弦曲率算法寻找凹点并进行匹配连接，实现气泡进行分割。Zhong 等[13]采用凹点检测结合弗里曼链码的方法对气泡进行匹配与分割。Chen等[14]采用凹点分割圆弧，对圆弧信息进行聚类分割气泡。基于标记的分水岭算法通过设置不同区域的标记减少分水岭算法造成的过分割[15]。Lau 等[16]和Fu等[17]通过阈值分割提取气泡中心明亮的区域作为标记进行分水岭算法的分割，统计气泡的参数信息。Peng等[18]采用k-means 聚类提取浮选泡沫的反光作为前景标记并进行分水岭算法分割。张文飞等[19]结合距离变换图像与连通域分析分离粘连物体作为标记进行分水岭分割。葛翔等[20]运用形态学开运算的方法分离粘连雾滴作为前景标记进行分水岭分割，进行树叶雾滴数量统计。上述方法在取得一定效果的同时也存在以下问题：①基于凹点检测的分割方法一般适用于圆度较高的气泡，对于形态不规则的气泡检测与匹配较为困难；
②基于标记的分水岭算法中使用的气泡中心通常有规律地出现明亮区域或反光区域；
③通过传统的距离变换和形态学的方法只能在粘连程度低的检测物上获得较好的效果。
本文基于高速摄像技术，以单孔排气为例，获取了不同来流条件和通气量下的气泡图像。针对传统图像处理方法对形态不规则气泡及气泡重叠导致的气泡分割识别率低的问题，本文改进了基于标记的分水岭算法的标记提取方法，通过比较距离变换图像的局部极值点抑制了伪标记，对气泡进行了准确分割。通过最小二值化椭圆拟合方法对不规则轮廓进行重构，实现了气泡特征参数的准确提取，研究了不同流动参数下气泡的演化规律。
为了研究不同液体流速和通气量对气泡行为特征的影响，在一小型闭式循环水槽开展了相关实验，如图1所示。该装置实验段寸为400mm×150mm×150mm（长×宽×高）的有机透明玻璃。实验段上方敞开连通大气，流速调整范围为0.01～0.3m/s。采用气体控制装置实现了通气，该装置由流量控制阀、气体过滤阀、流量计、气管、过滤阀等组成。实验中来流速度采用0.1m/s和0.2m/s两种，在每种流速下设置有通气量100mL/min、200mL/min、300mL/min 和400mL/min 四组，共有8 种工况进行对照实验。
图1 实验装置示意图
实验中采用高速全流场显示系统获得了气泡形态特征，高速摄像机采用的是Phantom v2512，其分辨率最高为1280×800，最高拍摄速率可达100 万帧/s，完全满足实验要求，本次实验采用的采集频率为1000Hz。为了实现定量化，保证高速摄影机的位置与姿态以及焦距等不发生变化，使用不锈钢直尺作为已知长度参考物对摄像机像素进行标定，计算得到像素当量为0.1176mm/pixel。
2.1 气泡处理难点概述
图2 给出了高速摄影机拍摄的气泡图像样本，体现了图像处理过程中的问题与难点，包括：①复杂的气泡外形和组合；
②气泡中心处的复杂反光；
③图像背景的空间不均匀性；
④气泡的弱边缘。
图2 原气泡图像
具体来说，红色区域为气泡上升过程中发生变形，蓝色区域中气泡发生不同程度的融合，变为重叠粘连气泡，且圆度较低使得气泡粘连时的形态复杂。此时，通过边缘曲率分割的凹点检测等方法难以适用于这类气泡。绿色区域是气泡表面形状不平滑导致的复杂反光，通过阈值分割的方法来提取气泡中心区域有一定的困难。黄色区域内是不均匀背景，且随着气泡体积增加，光照到气泡上会产生阴影。紫色区域气泡弱边缘容易导致二值化过程中气泡的部分轮廓丢失。可见，图像中这些复杂的干扰因素，采用传统气泡处理方法有一定的局限性。图3为两种传统图像前景提取方法[21]的效果：一是在利用滤波算子对图像进行预处理后，直接使用大津法（OTSU）[22]进行二值化的方法，这造成如图3(a)红圈所示位置出现了轮廓丢失；
二是使用图像增强的方法加强边缘处的对比度再进行二值化，如图3(b)所示，边缘信息得到加强的同时背景方块处的阴影部分也被增强，对后续图像的二值化造成阻碍，且一部分阴影由气泡遮挡形成，其出现位置不固定，难以直接通过减去拍摄背景的差影法消除。由上可见，这两种传统方法都难以有效准确地提取复杂图像中的气泡特征。
图3 传统气泡提取方法
2.2 气泡图像处理方法
图4给出了本文图像处理流程，分为三个主要部分。
图4 图像处理流程
（1）前景提取。通过Canny算法提取图像纹理信息补充原图轮廓，转化为二值图，再用形态学方法提取出气泡的投影区域。
（2）图像分割。对获得的二值图像进行距离变换，利用极值点获取气泡的前景标记，制作掩模后进行分水岭分割粘连气泡。
（3）拟合重构。采用最小二乘法对分水岭算法得到的气泡边界曲线进行拟合重构，提取气泡的位置与尺寸等参数信息。
2.2.1 前景提取
根据图像背景中阴影边缘过渡平缓的特点，本文利用Canny边缘检测算法重新设计了前景提取方法。为了提高图像质量，首先采用中值滤波消除图像中的噪点，对图像进行预处理，再通过四个步骤实现前景提取：第一，使用Canny算法对图像进行检测[23]，提取气泡边缘轮廓以及内部纹理，并将其叠加在被检测图像上，图像的弱边缘处得到了加强，气泡整体轮廓得到了补充，但是由于背景的不均匀性，背景中的部分轮廓也被检测出，如图5(a)所示；
第二，使用大津法对叠加图像进行自适应阈值分割，将其转为二值图像，如图5(b)所示；
第三，针对阈值分割后气泡内部的孔洞，利用泛洪算法对其内部进行填充，如图5(c)所示，且从图中可以看出，背景中仍存在微小黑点，由于这些点为水中的一些细小杂质，所以通过连通域阈值的设定将其过滤掉；
第四，通过形态学开运算将气泡边缘的毛刺与凸起消除，气泡投影区域被完整地从背景中提取出来，如图5(d)所示。
图5 前景提取流程
2.2.2 图像分割
基于分水岭算法的图像分割，关键在于单个气泡的前景标记。基于前景提取的二值图像[图6(a)]，利用欧式距离变换可将二值图像转换为灰度级图像。欧式距离是图像中两像素之间的直线距离，距离边缘越远的像素灰度值越高[24-25]。图6(b)给出了气泡前景的距离变换图像，可见随着气泡中心距离边缘的欧式距离增加，其灰度值增加，中心处达到最大。
图6 距离变换图
基于距离变换图像，采用阈值分割提取大于阈值的区域可作为单个物体的前景标记，这种方法往往适合在大小相当、形状规整（圆形、方形）且只在边缘处有粘连的物体，如圆形细胞[26]、大豆种子[27]等。但对于这种大小形状都具有一定差异的气泡，其中心位置的灰度值之间也会有较大差异，难以通过设置全局阈值或采用连通域最大值的方法进行标记提取。本文将利用距离图像的局部极值点进行标记作为分水岭算法的注水点，查找像素局部极值点的判断式，如式(1)所示。
图7(a)给出了提取极值点的示意图，由于形状的不规则，气泡内部会出现多个局部极值点，出现伪标记，直接在距离变换图像上使用分水岭算法会导致严重的过分割。因此，需要在局部极值点的基础上增加判断条件来抑制部分局部极值点（伪标记）的出现。在查找完局部极值点再设定一个范围L，将所有的极值点与周围一定范围内的像素再次进行比较寻找峰值点，峰值点判断条件见式(2)～式(4)。
图7 极值点的抑制
式中，Pp为保留的峰值点；
Pn为剔除的局部极值点。
当局部极值点与设定范围内的某个像素灰度值相等时将局部极值点Pm置零，避免在设定范围内同时出现多个极值点或出现多个极值点相互抑制剔除的情况。图7(b)给出了多次判断后极值点示意图，此时大部分伪标记被抑制，实现了前景有效标记。
从图7(b)中右侧框内可以看出一个气泡中仍然可能存在多个峰值点，为了不造成过分割，一个气泡值内应只有一个标记。通过观察发现，属于同一个气泡的峰值点距离较近，将所有的峰值点进行适当膨胀，可将距离较近的点融合为一个前景标记，如图8(a)右侧框内所示。再使用形态学操作将二值图膨胀与所获取的前景标记取差，得到图8(b)所示的白色区域，即不确定区域，分水岭算法的分割线将在这个区域内产生。
图8 峰值点融合与不确定区域
将不确定区域标记为0，前景标记区域设置为1，使用分水岭算法，从标记区域开始漫出，在不确定区域中相遇形成分水岭完成气泡的分割，分割结果如图9所示。
图9 分水岭算法分割结果
同时，将该结果与基于阈值分割[27]和基于腐蚀运算[28]的标记获取方法对比，如图10 所示。可以观察到，腐蚀运算通过腐蚀粘连处达到分离气泡的目的，但腐蚀运算在腐蚀粘连处的同时也会腐蚀气泡本身，这会对气泡标记造成破坏。对于在距离变换图像上进行阈值分割的方法，不同尺寸气泡中心处的灰度值差异较大，而气泡粘连处与气泡中心处的灰度值差异较小，难以设定一个阈值将气泡中心区域提取，无法分离气泡造成欠分割。而本文方法通过对距离变换图像极值点之间的比较，抑制了伪标记的出现，寻找到气泡中心区域的最优点并进行了一定程度的融合，这使得本文方法在处理重叠、尺寸差异大以及形状不规整的气泡时比前两种方法具有更好的效果。
2.2.3 拟合重构
由于基于分水岭算法得到气泡形状并不规则，为提取其特征参数并进行统计分析，需要将其拟合为规则形状。在此采用基于边界的最小二乘椭圆拟合算法，即将待拟合的边缘曲线视为一个个相连的独立点，作为观测值，理论值为需要拟合的椭圆曲线函数，用一个损失函数去表征理论值与观测值之间的残差平方和，找到一条使得损失函数最小的椭圆曲线，进而以最小误差实现气泡的拟合重构[29]，如图11所示。
图11 最小二乘椭圆拟合原理
一般的二次曲线F可表达为式(5)。
点ti(xi,yi)是轮廓曲线上的N(N＞5)个检测点，根据最小二乘原理，椭圆二次曲线的拟合方程参数通过代数距离平方和的最小化来求解，拟合曲线代数平方和计算为式(6)。
求得符合条件的向量k使得C(k)取得最小值，即可得到椭圆方程参数。根据椭圆方程可以计算椭圆的位置与大小等参数。(Xc,Yc)为拟合椭圆在图像上的几何中心坐标，用于描述椭圆在图像上的位置，其计算分别为式(7)、式(8)。
A为椭圆的长轴直径，B为短轴直径，用于重构气泡的形状，计算为式(9)、式(10)。
θ是重构椭圆的长轴倾斜角度，为长轴与水平轴线的逆时针角度，计算为式(11)。
基于图9利用分水岭算法得到的气泡轮廓分割曲线，对其进行最小二乘椭圆拟合，可计算得到椭圆的中心坐标及长短轴等信息。图12 给出了绘制在图像上的重构气泡轮廓曲线，可以看出该拟合结果较好。
图12 气泡的最小二乘椭圆拟合重构
3.1 气泡数目统计分析
为了评估本文所提出算法的气泡分割效果，进行了气泡数目的统计实验验证。本文对实验中8种工况下的气泡图像（200张图片）分别采用了上述的3 种方法进行识别，并与人工计数结果进行对比。腐蚀运算难以分离粘连程度高的气泡，阈值分割则无法提取对于尺寸形状差异大且重叠粘连形式复杂的气泡标记，导致错误识别。本文方法能够较好地解决上述问题，相较于腐蚀运算方法识别准确率提高了22.7%，相较于阈值分割方法准确率提高了13.6%，如表1所示。
表1 传统算法与本文算法对比
3.2 气泡参数分析
基于气泡图像处理，提取气泡椭圆特征参数，标记数代表为气泡数，椭圆质心坐标表示气泡位置参数，坐标系的横、纵轴分别对应X轴与Y轴，并计算气泡的等效直径。图13 和图14 分别给出了两种来流速度0.1m/s和0.2m/s，不同通气量下气泡上升演化过程。发现在同一横向来流速度下，通气量为100mL/min时，表现为较为独立的气泡形式，其运动轨迹呈线性上升。随着通气量的增加，逐渐出现重叠与粘连的气泡，但气泡的整体位置分布并未随着进气量的增加发生明显变化。当液体流速从0.1m/s增加到0.2m/s时，相同通气量下气泡的整体位置发生改变，表明液体流速是气泡运动轨迹改变的主要因素[30]。
图13 0.1m/s液体流速下气泡位置分布
图14 0.2m/s液体流速下气泡位置分布
为了表明流速和通气量对气泡特征的影响，分别使用上述三种方法以及人工方法进行了测量对比。图15 和图16 分别给出了不同流动参数下平均气泡统计个数和气泡尺寸，从图中可以看出本文方法与人工识别获取的数据更为吻合，说明了本文方法可准确获取数据。从图中本文方法获取数据看出，当通气量从100mL/min增加到400mL/min，气泡数量略微增加，其中通气量200mL/min与300mL/min相比无明显增长。而气泡的平均尺寸随着通气量增加明显增加，尤其当通气量从200mL/min 增加到300mL/min时，增长幅度最大。表明气泡的生长过程中，通气量的增加促进了气泡的生长，使其尺寸增加，其次才是气泡的数量的增加。在相同通气量下，横向来流速度0.2m/s时气泡的平均尺寸相较于流速0.1m/s的情况更小，但是0.2m/s的流速下气泡平均数量更多。这是因为在气泡的生长过程中，随着气体的通入气泡产生的浮力随着体积增大而增大。当浮力大于气泡表面张力时，气泡脱离附着的气孔向上浮动，液体横向来流速度会给气泡施加一个横向的拖拽力加速气泡的脱离，流速越大则拖拽力越大，单个气泡生长的时间要比同一通气量下要短，造成气泡的平均尺寸更小，但是平均气泡数量更多[31]。
图15 不同液体流速（V）下气泡平均个数随通气量的变化
图16 不同液体流速（V）下气泡平均尺寸随通气量的变化
本文基于高速摄像获得了不同流动参数下的气泡特征，针对气泡图像中重叠粘连气泡造成的欠分割导致气泡特征参数提取误差问题，改进了分水岭算法的前景标记提取方法，实现了气泡的准确分割和拟合重构，提取了气泡特征参数并进行对比分析，得到以下结论。
（1）在对于大小形状都具有一定差异的重叠粘连气泡的处理上，通过抑制与融合距离变换图像极值点的方法可有效抑制伪标记，实现较为精准的气泡分割。本文方法相较于腐蚀运算分离气泡的方法识别准确率提高了22.7%，相较于对距离变换图像进行阈值分割的方法识别准确率提高了13.6%。
（2）通过统计气泡椭圆拟合重构获取的特征参数，发现随着通气量的增加，气泡的尺寸显著增大，但气泡的数量只是略微增加，通气量的增加主要促进了单个气泡的生长过程，其次是影响气泡的数量，而气泡整体分布并未随着通气量的改变有明显变化，通气量在100～400mL/min范围内变化时对于气泡的运动无明显影响。
（3）随着液相流速的改变气泡的整体位置分布发生改变，气泡的运动主要受液体流速影响。相同通气量下液体流速越高气泡数量越多，而平均尺寸更小，横向来流液体提供的拖拽力缩短了单个气泡的生长时间，使得气泡的数量更多而尺寸更小。
猜你喜欢 分水岭极值气泡 柠檬气泡水欣漾(2024年2期)2024-04-27极值点带你去“漂移”新世纪智能(数学备考)(2021年10期)2021-12-21SIAU诗杭便携式气泡水杯新潮电子(2021年7期)2021-08-14浮法玻璃气泡的预防和控制对策建材发展导向(2021年13期)2021-07-28极值点偏移拦路，三法可取河北理科教学研究(2020年3期)2021-01-04一类“极值点偏移”问题的解法与反思中学数学杂志(2019年1期)2019-04-032019，一定是个分水岭！兽医导刊(2019年1期)2019-02-21冰冻气泡儿童故事画报·发现号趣味百科(2019年9期)2019-02-02借助微分探求连续函数的极值点广东技术师范大学学报(2016年5期)2016-08-22“华北第一隧”——张涿高速分水岭隧道贯通城市道桥与防洪(2014年1期)2014-02-27