基于灰度图差分梯度的雾化角测量方法

魏永杰，李明凯，王浩然，张树日

（河北工业大学机械工程学院，天津300130）

引言

喷雾技术在燃油、消防、冶金、环保、航空航天、化工等领域有了越来越广泛的应用。由喷嘴喷出的液滴在周围空气流的作用下形成复杂的气液两相流，雾场特性分析是提高雾化效果的重要途径，其中雾化角是一项重要指标。刘祺等分析了在不同环境条件下对航空涡轮风扇发动机中所使用的离心式喷嘴喷雾特性，并开展了实验研究[1]。梁博健等以喷射角和射流流量作为评价指标，考虑多项参考因素的影响，研究高压水除鳞喷嘴的射流性能[2]。蒋仲安等采用一次雾化和二次多级雾化理论得出影响喷嘴雾化性能的主要因素，并通过实验分析各因素对喷嘴雾化性能的影响[3]。杨国华等基于高速摄影技术和图像处理方法，以不同螺旋升角作为变量，对氨法脱硫用螺旋形实心锥喷嘴的雾化特性进行了分析和预测[4]。郭鹏宇等分析了变压器灭火领域喷嘴的安装角度、喷嘴直径和喷雾压力等对喷嘴雾化特性的影响规律，得到了最佳的喷嘴安装角度[5]。贾卫东等使用激光相位多普勒粒子测速仪直接测定喷雾雾化角，精度较高，直观地显示了雾滴数分布状况[6]。

根据流体力学原理，液体以一定压力经喷嘴喷射后，与空气流相互作用，形成一定范围内的雾场。由于喷射出的是微米级的小颗粒，而喷雾流的速度较高，会扰动周围空气形成差压气流，尤其是多相流形成的细小雾滴随空气湍流飘洒，造成喷雾场的边缘不清晰，测量比较困难。

喷嘴雾化角即喷嘴形成的雾场边界所形成的夹角。除目视测量外，测量喷嘴雾化角主要是采用图像法[7]。除普通可见光图像法外，Wei Yijie 等采用高速相机对燃油喷雾使用粒子追踪测速法进行图像分析的同时，获得了喷雾的瞬态轮廓和雾化角[8]；Chen Run等根据激光散射吸收法(laser absorption-scattering,LAS)，采用紫外、可见光双光束成像技术分析了燃油喷雾的特性和形貌特征[9]；Gröger Karsten 等采用X 射线图像方法对喷雾的近场形貌进行了分析[10]。与上述方法相比，采用普通可见光图像法分析雾场边界和雾化角是最经济快速的方法，但由于喷雾场是由液滴颗粒组成的不连续流体场，采用通常的图像边缘算子，如Sobel算子、Canny 算子、Laplace 算子等直接在灰度图上提取到的雾场边界并不是很理想。

基于上述原因，一般采用的处理过程是通过对喷雾图像降噪滤波后进行二值化处理，然后根据图像形态学针对二值化的图像进行边缘提取和计算[11]。由于喷雾图像对比度较差且雾滴粒度、浓度分布不同，造成灰度图像不均匀，图像二值化过程比较困难。目前有多种算法进行复杂灰度图的二值化方法，如全局阈值二值化、自适应二值化、OTSU 算法等。但这些方法主要是针对图像进行后期处理，计算喷雾场边界缺少客观的评价依据。

根据图像形态学梯度原理，将灰度图像分别膨胀、腐蚀后与原图像相减得到梯度图像，理论上也可以得到图像边界，但由于雾场图像灰度值分布复杂，难以得到有效边界。

图像的差分梯度算法即针对有某种关联特性的两幅图像或经过处理的两幅图像，将对应像素或相邻某个区域的图像做差分，从而得到图像的特征[12]。

本文提出根据喷雾灰度图得到雾场范围扩大的二值化掩模板，并将此掩模板作用于灰度图像，针对原始灰度图进行迭代，并对喷雾范围边缘的平均灰度值进行梯度计算，认为梯度最大时的图像为实际喷雾图像，从而得到雾场边界，并可以进行喷雾角度计算。

1 算法原理

1.1 雾场边界特征分析

图1 是拍摄的喷雾图像和采用边缘算子直接处理灰度图的结果。图1（a）是拍摄的喷雾原图。该喷嘴是一个小型的细水雾喷嘴，喷雾由左到右水平喷出。由于喷雾浓度较小、气流影响较大，因此在喷射一定距离后雾场和背景的灰度值相差不大。为了更清晰地看清楚喷雾图像，采用直方图均衡化对图1（a）进行了处理，如图1（b）所示，图中可以看到上、下2 个灰度值较大的边缘。由于射流卷吸问题[13]，受气液两相流和空气剪应力的共同作用[14]，造成空气流在2 个边缘附近形成负压且细水雾中的较小颗粒受阻力大，飘散慢，停留时间长，故在2 个灰度值较大的边缘外各有一条灰度值较小且分布不连续的细水雾飘散带[15]。因此，认为灰度值较大的边缘为实际喷雾的有效边界范围。

图1（c）、图1（d）、图1（e）和图1（f）分别是采用4 种边缘算子对图1（b）灰度图进行处理的结果。但这4 种算子只得到雾场边缘的1 个范围或1 个局部，无法得到清晰完整的边缘，更无法判断实际雾场的有效边界。

图1 喷雾灰度图和边缘检测结果Fig.1 Grayscale image of spray and boundary detection results

图2 是将图1（b）膨胀后与原图差分得到的梯度图像。图2 中可以看到模糊的边界，与图1 中Canny 算子、Robert 算子得到的边缘相似，但仍不能判断实际边界。

图2 形态学梯度处理结果Fig.2 Gradient processing results by morphology

通过上述分析可知，这些方法是在一定灰度值变化范围内寻找图像灰度差异点，并认为这些点是图像的边界。因此，会找到多个边界点，从而得到边界是1 个范围，而不是某条确定的边界线。

1.2 灰度图差分梯度算法

为了得到有效的边界线，采用在灰度图上直接统计一定喷雾边缘范围内的灰度值，取灰度差的最大值作为判据提取边界的方法。

雾场和背景相比，雾场是高灰度值，而背景是低灰度值。为此，首先在原始灰度图像上，通过图像形态学的膨胀操作获得比实际雾场区域大的假设雾场区，然后采用全阈值二值化得到掩模板图像，其中假设的雾场区域为“白”，背景区域为“黑”，并将此二值化图作为迭代算法的初始掩模板。

直方图均衡化并不改变实际雾场的分布，为观察方便，从获得初始掩模板的步骤开始，均采用图1（b）直方图均衡化的图像作为原始图像。将该掩模板迭代进行腐蚀操作，依次得到边界范围逐渐减小的掩模板；采用相邻的2 个掩模板分别作用于图1（c）的图像，将得到的2 个图像做差分后，计算差分灰度图像的累计灰度值；将2 个相邻差分图像的平均灰度之差作为灰度梯度。

由于飘散的细水雾从有效雾场边界向外浓度逐渐降低，因此该灰度梯度从迭代开始逐渐增大；而在有效雾场范围内，灰度梯度依次减小，在飘散的细水雾颗粒和有效雾场边界处有最大值。因此，通过迭代可以得到该有效边界。

图3所示为迭代原理。由图3 可知，在图像处理整个过程中，只有模板是二值化的图像。将原始图像多次膨胀后取全阈值二值化，得到模板1，将此模板“罩”在原始灰度图上，得到喷雾区的灰度图1，同时将喷雾区外的灰度值置为0，并计算喷雾区的平均灰度值1；将二值化模板1 再腐蚀一次，作为二值化模板2，重复上述过程，得到喷雾区的灰度图2 和平均灰度值2。与灰度值1 相减，得到灰度图的梯度1；同理，得到第n次迭代的梯度n。由上面的分析可知，当梯度n取得最大值时停止迭代，得到喷雾有效边界。

图3 灰度图梯度迭代原理框图Fig.3 Block diagram of gradient iteration on grayscale images

2 实验及分析

采用上述提出的方法，针对图1（a）实际采集的雾场图像，经直方图均衡化后进行迭代处理，得到雾场有效边界，如图4所示。

图4（a）是迭代算法采用的掩膜板，是将图1（b）的灰度图通过膨胀操作后采用全局二值化获得。比较图4（a）和图1（b），喷雾范围扩大，尤其从左侧喷嘴位置观察更加明显。喷雾是由满足一定分布的液滴粒子构成，在拍摄的原始图像和直方图均衡化图像中，在喷雾区会有一定的孔洞，这些孔洞是由于该处不存在液滴粒子或液滴粒子数量较少形成的。图像形态学的膨胀操作恰好具有填补孔洞的作用，因此，图4（a）是连续性较好的图像。一般来说，对于雾场范围的图像测量方法需要采用平滑滤波，本文采用膨胀法提高了连通性，且后续步骤针对灰度图进行迭代，所以未采用平滑滤波方法。

图4 喷雾边界迭代结果Fig.4 Iterative results of spray boundary

图4（b）是针对灰度图进行迭代的结果。比较图4（b）和图1（b），处理后的结果与图1（b）中较亮的边界一致，因此认为此范围是雾场喷射的边界，而不是细水雾飘散形成的边界。比较图4（a）和图4（b），由于在获得初始掩模板时，在图4（a）右侧中部有1 小块未经膨胀去除的黑斑。因此，在迭代过程中，将掩膜板逐次腐蚀后在图4（b）右侧对应实际喷雾位置生成了1 块错误的小黑块。同时，在图4（a）掩膜板下端中间位置，代表扩大的假设喷雾区出现了较大的直线条，经多次腐蚀操作后，对应的图4（b）代表实际喷雾边界的下端中间位置出现了较大误差。迭代过程中的梯度值如表1所示。

表1 灰度梯度迭代数据Table 1 Results of gradient iteration on grayscale images

由表1 可知，在迭代第8 次时，灰度梯度有最大值，此时即雾场有效边界。通过比较第8 次迭代以后的数据，尤其是第8 次和第9 次迭代的数据可知，由于雾场的连续性特性，在边界附近的值较大，而雾场内部液滴粒子分布较为一致，因此灰度梯度逐渐减小；在雾场外部，由于细水雾的飘洒，使雾场有效边界外的附近区域中，细水雾液滴逐渐减小，且细水雾液滴少、与喷射的水雾液滴不存在连续性，因此第7 次迭代数值为2.503，远小于第8 次迭代数值100.094，且由第7 次到第1 次依次减小。因此，将第8 次迭代的灰度图作为识别得到的雾场区域图，即图4（b）对应的灰度图。

表1 数据中，第2、4 次迭代数据略小于前一次迭代数据，第14 次迭代数据略小于后一次迭代数据，与雾场实际分布有波动有关。

为了与传统的针对二值化图像处理方法进行比较，采用二值化方法对图1（a）的灰度原始图先进行二值化边界识别，如图5所示。

图5 传统二值化雾场边界处理结果Fig.5 Results of traditional binary spray field boundary

图5（a）是直接观察和对比原始图设定阈值的结果；图5（b）是选取比图5（a）更大阈值的结果；图5（c）是采用OTSU 算法自动选取阈值的结果。从图5 的3 个结果可知，难以采用数值判据来分析雾场边界。

为验证上述方法的有效性，对另外3 种不同浓度的喷嘴雾场进行迭代运算。图6 是喷雾边界处理结果。其中图6（a）、图6（b）、图6（c）是拍摄的原始图，对应的雾场浓度依次减小。如图6（a）所示，在产生水喷雾场时，施加的压力越大，雾滴能完全破碎，雾场越浓；如图6（b）和图6（c）所示，施加的压力越小，雾滴不能破碎，且喷出的水量少、雾场稀。这两种雾场浓度依次减小，拍摄图片的灰度值也减小，同时雾滴因未完全破碎呈近细流状，雾场边界更模糊，直接用边界提取也更加困难。

图6（d）、图6（e）、图6（f）分别是采用本文方法的处理结果，可见随压力减小，雾场区域逐渐减小；由于边界附近雾滴少，因此图6（e）、图6（f）中得到的雾场边界更偏离直线。

图6 不同浓度雾场对比结果Fig.6 Comparison results of spray field with different density

为了拟合喷雾角，将图4（b）进行二值化处理，如图4（c）所示。由于处理过程中，已将图4（b）非雾场区域置0。因此，二值化后的图4（c）与图4（b）雾场区域边界是相同的。

提取图4（c）的边界点S，进行分段曲线拟合，再分析曲线斜率，分别找到上、下两条曲线斜率变化的最大值。考虑到喷雾角一般为锥形，因此取斜率变化最大值左边的上、下两段曲线，分别进行直线拟合，取两条拟合直线的夹角作为喷雾角。经计算，得到该喷雾角为82.65°，如图7所示。

图7 喷雾角计算结果Fig.7 Calculation results of spray angle

图7 中，取斜率变化最大值并截取由喷雾开始位置之后的平直段，得到两条由边界点构成的曲线。图中的两条直线分别是对上、下边界曲线拟合得到的直线，图中坐标代表的均是像素数。

由于雾场本身的动态特性且受周围气流扰动影响较大，因此对雾场边界和角度进行准确测量和标定还比较困难。本文提出用灰度梯度的最大值作为判据和图像处理后进行分界角度计算，提供了一种数值判断依据。

上述方法是针对深色背景下的水雾场图像进行处理，当背景较复杂时，用本方法与光谱图像的主成分分析方法等结合使用，也可以实现雾场边界定位。

3 结论

由于喷雾场液滴的非连续性以及两相流的作用，造成雾场空间分布动态波动特性明显，且雾场边界不清晰，灰度值差别较大。将灰度图像进行二值化再提取边界或计算喷雾角，缺少定量评判依据，得到的结果与实际雾场有偏离。

根据灰度图梯度算法，除掩模板处理的二值化图像外，其他处理步骤均是针对灰度图进行处理，并通过构造灰度图梯度，当灰度图梯度为最大值时，认为喷雾图像边界是实际喷雾场的有效边界，在此基础上进行喷雾角拟合计算。该方法提出了一种采用数值判断喷雾场有效边界的计算方法，通过灰度梯度最大值确定雾场边界，从而通过程序自动实现雾场边界提取和喷雾角计算，避免传统方法中需要通过直接观察缺乏评判依据的问题。