基于图像处理的低速横流中液体射流轨迹提取方法研究

兰 天, 孔令真, 陈家庆, 王奎升,*

(1. 北京化工大学 机电工程学院, 北京 100029; 2. 北京石油化工学院 机械工程学院, 北京 102617; 3. 深水油气管线关键技术与装备北京市重点实验室, 北京 102617)

0 引 言

液体射流注入横向运动气流是一种常见的气液雾化混合方法，在超燃冲压发动机[1]、农用喷雾器[2]和新型管式气液接触吸收器[3]中得到了广泛应用。液体射流柱在横向运动气流作用下经历弯曲、断裂以及破碎等过程，最终雾化成粒径更小的液滴群分布于来流气相中。由于不同工况下射流弯曲轨迹不同，所形成的雾化液滴群在射流方向上的粒径大小和分布存在差异，这将在很大程度上影响气液两相的接触比表面积，进而影响气液两相的混合效果[4]。

针对横向运动气流中的射流弯曲轨迹特征及其影响因素，国内外学者采用脉冲阴影、高速摄像等测量方法开展了一系列可视化实验研究[5-6]，但这些测量方法在原始图像采集过程中都存在一些问题，例如部分射流区域“高亮”、图像背景嘈杂、射流柱边界拍摄不清晰等，进而对图像后处理时射流轨迹的提取带来很大挑战，因此，需要对采集到的射流破碎图像进行处理，以获得更为准确的液体射流轨迹。Song等[7]采用一种简单的阈值分割方法从高速摄像的原始图像中提取射流轨迹，由于阈值选择存在主观判断，提取结果不能很好地切合射流柱轨迹轮廓。刘静[8]采用Photoshop和Origin软件提取射流纹影图像边界线，在处理大量实验图像时具有很大难度。Amighi[9]通过脉冲阴影获取原始图像，研究了不同阈值算法对图像分割结果的影响，发现大多数阈值算法或稀释图像而低估边界，或受背景噪声影响而高估雾化区域，较为适合的阈值算法是均值、三角形和Huang算法等基于单个强度分布的算法，但这些算法并未针对每帧图像的特点来自适应选择阈值。柴进[10]针对超声速横流条件设计了一种基于外轮廓的激波线提取法，认为阈值选取的合理性很大程度上决定了边缘提取的准确性，通过设定合理阈值，从滤波结果中去掉梯度值较小的部分边缘，可以精确提取激光阴影图像中的射流激波线。

横向运动气流中液体射流柱在气液相互作用下经历复杂的运动和变形过程，不同工况下表现出的射流破碎形态不同，这使得射流区域中的液体密度分布也发生变化，所造成的光照强度在背景和射流柱等不同区域存在差异[11]。在图像处理过程中，需要根据这种差异对成像后的图像进行一系列处理与优化，以消除图像采集过程中存在的问题，适应不同工况下的射流破碎图像，并最终准确得到描述液体射流破碎的轨迹参数[12]。

轨迹轮廓提取是计算机图像视觉的重要研究领域。由于目标区域的多样化，轨迹轮廓提取很难具有统一的处理方法，众多的轮廓提取均是基于目标区域特点提出针对性的处理方式。对于大量的实验图像而言，Photoshop、ImageJ等软件的处理过程较为通用化，且很难做到批量选择和处理。为此，本文结合图像处理技术，采用高速摄像仪记录横向运动气流中液体射流破碎的历经过程，提出低速横流作用下的液体射流弯曲轨迹提取方法，并结合MATLAB予以编程实现，从而准确得到射流弯曲轨迹，为进一步分析不同工况下射流弯曲轨迹的影响因素奠定基础。

1 实验装置及轨迹提取方法

1.1 实验装置

如图1所示，所采用的实验装置系专门搭建用于研究液体射流注入横向运动气流中的射流雾化机理的实验平台。气相由空气压缩机以及与之相连的高压储气罐提供，储气罐最高储气压力0.6 MPa，采用温压补偿一体式涡街气体流量计(合肥精大，VAFTP-050-DC1-213-1.0-P1.6)测量通过的气体流量。水相为储液罐提供的超纯水。实验时，将高压气瓶的压缩空气通入储液罐的气相空间中，以气压产生稳定流动的液体流，再从射流口注入运动气流。采用微小量液体流量计(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 液体射流轨迹提取方法

液体射流注入横向运动气流中，经历复杂的相互作用，射流柱从发生弯曲、变形的射流初次破碎，到破裂、雾化的射流二次破碎，最终形成较为稳定的雾化液滴群分布于来流气相中。射流初次破碎状态关系到雾化液滴群的液滴数量、粒径尺寸和空间分布[13]。图3直观展示了液体射流破碎的主要研究内容，即射流区域内外边界、射流变形与破碎、雾化结构等。一般的，众多学者将射流区域外轮廓定义为液体射流轨迹，并以此进一步分析影响射流轨迹特征的因素[14]，这是因为外轮廓决定了沿射流方向上射流区域的最大分布区间。

因此，在本文研究中将射流区域外轮廓定义为射流弯曲轨迹，并通过高速摄像机拍摄射流破碎图像，提出基于MATLAB图像处理技术[15]的低速横流情况下的液体射流轨迹提取方法，以定量获取射流轨迹数据。

图3 液体射流破碎结构示意图

液体射流轨迹提取方法如图4所示。首先进行图像采集，采用高速摄像机拍摄不同工况下多帧连续射流破碎图像，镜头光圈调至固定值，调整镜头焦距并固定，确保不同工况下足够的进光量以及相同景深，根据实际拍摄的情况，灵活调整光源角度和光照强度；其次，对原始图像进行增强处理，即将原始图像灰度化处理后进行直方图均衡化，目的是进一步锐化增强所关注的射流破碎区域；然后，采用最佳直方图熵法(KSW熵法)及传统遗传算法进行图像阈值分割，该阈值方法可以自适应快速选择各帧图像的2个阈值，将射流区域同背景分割；最后，使用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]。

灰度图像可以通过灰度直方图来反映图像中的不同灰度等级分布，其中，横坐标为不同的灰度等级，纵坐标为不同灰度等级下像素点的个数。本文采用直方图均衡化方法进行图像增强，可使图像的灰度间距拉开或使图像的灰度分布趋于均匀，增大图像与背景的反差，使图像细节清晰，从而达到增强图像的目的[17]。

图5为工况Wej=496.46，q=117.64条件下图像增强处理前后的对比。图5(a)为实验观测的原始图像、均衡化后的图像以及去背景后的图像(从上至下)，图5(b)为图5(a)各图像对应的灰度直方图。从图5可以看出：经过图像增强处理，本文关注的液体射流区域的亮度得到很大程度的增强；通过MATLAB去背景函数Imsubtract，原始彩色图像中大小不一的壁面黏附液滴反射的光斑等噪点也得到了很大程度的去除。从灰度直方图还可以看出：原始彩色图像灰度分布于较窄区间，这会导致图像细节不够清晰、背景与目标区域混叠严重不易区分；在经过直方图均衡化后的灰度直方图中，灰度区间均衡分布，背景作用减弱，射流区域增强。

图5 图像增强处理结果

2.2 图像阈值分割

图像阈值分割是一种用于灰度图像的图像分割方法。对于本文实验拍摄的射流破碎图像而言，就是通过某种途径寻找射流破碎图像中适合的某个阈值，若图像中某处像素点的灰度低于该阈值，则视为背景区域，否则视为液体射流区域的一部分[10,17]。常用的全局阈值选取方法有峰谷法、最小误差法、最大类间方差法(Otsu)、最佳直方图熵法(KSW熵法)等。其中，Pun等[18]提出的最佳直方图熵法是一种以图像灰度直方图最佳熵为判断准则进行图像阈值分割的算法。从横向运动气流中的射流破碎高速摄像图像可以看出，在不同气液工况条件下，液体射流柱的破碎剧烈程度不同，这就使得在光源条件不发生变化的情况下，高速摄像机所拍摄的不同工况下的射流破碎图像完全不同，因此，需要在图像处理中依据不同射流破碎图像自身的灰度特征自适应地选取阈值分割的最佳阈值[19]。对于KSW最佳直方图熵法二阈值分割情况(即存在2个分割阈值S1和S2)，则有：

(2)

(3)

图6 图像像素直方图和KSW二阈值结果

图7 阈值分割后图像

2.3 边缘轮廓检测

图像处理中的边缘定义为图像局部特性的不连续处，如灰度的突变、颜色的突变、纹理结构的突变等。对于灰度图像，边缘是图像中灰度变化较为剧烈的位置，亦即通常而言的信号发生奇异变化的位置[21]。横向运动气流中的液体射流柱存在明显的轨迹特征[22]，从2.2节的阈值分割结果可以明显看出，在射流破碎图像中，射流柱迎风面沿边缘切向的像素变化平缓，而法向(垂直于边缘)的像素变化剧烈。射流柱迎风面、背风面与背景图像之间具有明显的灰度梯度，依据这种灰度梯度特性，可以对图像的边缘进行检测和提取，并最终定量得到射流穿透轨迹特征。

常用于边缘轮廓检测的方法有Prewitt算子、Roberts算子、Canny算子、Sobel算子等，不同算子是基于图像具体的突变特征或不同的计算方法得到。Sobel算子是一阶微分算子(离散差分算子)，利用临近(上下、左右)的像素值(灰度值)来计算当前像素值(灰度值)。由于图像边缘的灰度变化较快，呈现出阶跃现象，Sobel算子边缘检测正是利用这一特点，对图像中的每个像素计算其周围相邻像素的差分值，来判断该像素点是否属于边缘[23]。由于本文关注的射流弯曲轨迹为射流柱外轮廓信息特征，其灰度值与背景存在很大的梯度，因此采用Sobel算子对阈值分割后的图像提取边缘轮廓。

图8(a)是利用MATLAB中的Sobel算子对阈值分割后的图像进行处理的结果；图8(b)是提取的射流柱外轮廓与原始图像的对比，外轮廓以红色显示，可以发现射流柱外轮廓能够最大限度地保留，与原始图像存在较好的切合性。

图8 边缘轮廓检测图像

检测到射流柱外轮廓后，为定量分析不同工况下的射流轨迹特征并探究其影响因素，需要提取射流柱外轮廓的数据点。本文实验中高速摄像图像的像素尺寸为1024 pixel×1280 pixel，且射流柱外轮廓由一系列具有灰度值的像素点构成，这些点构成一个点集，因此，可以结合凸包概念(Convex Hull，将给定二维平面上的点集最外层的点连接为凸多边型，以包含点集中所有的点[24])得到包围这些点集的最外层数据点，其连线即为射流轨迹曲线。液体射流柱在横向运动气流中具有明显的轨迹特征，经上述处理后，可以得到更为清晰的射流破碎轮廓。本文结合MATLAB中的凸包算法Convhull及边缘点优化和稠密化，得到射流柱外轮廓数据点集。图9中的虚线为射流柱外轮廓数据点集的提取结果，图10为多帧图像的轮廓点提取对比，可以看出，提取出的数据点能够很好地切合射流柱外轮廓。

3 轨迹提取方法的适应性验证

液体射流注入横向运动气流中的射流破碎是一个气液相互作用的复杂过程。Vich[25]、Birouk[26]等研究了低速横流作用下的液体射流破碎模式特征，认为存在弱横流作用破碎、拱形破碎和袋式破碎等3种典型射流破碎模式。这是由于在横流气速变化的情况下，运动气流气动力不同，使得射流柱破碎状态存在差异。图5和图7～9是典型的袋式破碎模式，图11(a)是弱横流作用破碎模式，图11(b)是拱形破碎模式(上为原始图像，下为轮廓检测提取后的图像)。从图中可以看出，射流破碎轨迹和破碎状态在不同工况下存在差异，边缘轮廓提取后的图像与原始图像能够很好地吻合。图11(a)的弱横流作用破碎模式中，横向气流作用影响较弱，射流柱弯曲程度较好，在通过阈值处理和边缘检测后对射流柱外轮廓提取数据点时，数据点集能够很好地切合射流柱外轮廓。在图11(b)的拱形破碎模式中，射流柱在远离射流孔的位置存在多处褶皱，其提取的外轮廓为射流区域最外层的边界点集。

图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算子和凸包算法等，能够有效适应低速横流作用下不同射流破碎模式的射流轨迹提取，对提取后的弯曲轨迹数据点进行非线性拟合，进而获得准确预测射流弯曲轨迹的经验公式。该方法可以对高速摄像机拍摄的原始彩色图像进行一系列批量处理，从而更加方便地讨论液体射流轨迹的影响因素，为气液雾化混合结构的设计提供参考。