气液两相流含气率图像检测方法

2013-08-24王秀英

计算机工程与应用 2013年5期

王秀英

上海新侨职业技术学院计算机信息系，上海 200237

在化工生产中，气液两相流的分析研究对改善生产效率和生产工艺都有很大的意义。比如化工设备中有许多都存在着气液、液液两相流，其中气泡、液滴大小的分布对研究设备中的能量传递、反应过程，以及工艺和设备的最优化设计与节能，具有十分重要的意义。在气液分散操作过程中，气泡的大小分布也是衡量过程结果的一个重要参数。以往的方法，如，电导探针法[1-2]、毛细管光电法[3-4]、光散射法[5]、摄像法[6-7]等，存在复杂、误判率大、支出成本高等缺点，因此本文开发了一套基于机器视觉技术的气液两相流在线检测系统，实现了气液两相流在线检测过程的自动化。

1 实验装置及方法

气液两相流在线检测系统的硬件使用CCD高速摄像仪摄取通过毛细管抽吸的动态气泡图像，然后由软件对采集的图像进行处理和分析，实现气液两相流中含气率的测量。实验系统如图1所示，主要分为4个部分：（1）毛细管抽吸部分。实验采用混有空气的水作为研究对象，采用直径为5 mm透明有机玻璃管固定到水槽中，毛细管的另一端使用抽吸泵，以1 m/s的速度对水槽里的水进行抽吸。（2）光源部分。经过在自然光、白炽光、荧光灯光三种光源下拍摄动态图像效果进行对比，发现采用定制的荧光灯及40 kHz以上的高频电子镇流器能提供稳定的光源，可以得到较好质量的图像。（3）图像采集部分。选用了台湾MINTRON公司的MTV-1802BC型1/2'黑白摄像机，单场技术，快门选在1/2 000 s以上，采集的图像清晰，无拖尾现象。（4）背景选择部分。经过比对墙壁和复印纸作为背景拍摄的图像，最后选用能够提供较均匀反射的墙壁作为背景。所得到的灰度图像如图2（a）所示。

图1 实验装置示意图

2 根据最小二乘法拟合毛细管的内外边界

2.1 图像处理

通过CCD拍摄的图像为512像素×512像素的灰度图像，为了能够得到毛细管和气泡的左右边界函数表达式，需要处理为二值图像[8-10]。（1）首先采用Laplacian算子与原图像进行叠加的方法进行边缘加强，见图2（b）；（2）然后利用Sobel算子对图2（b）中的图像进行边缘检测，见图2（c）；（3）最后利用熵法对图2（c）中的图像二值化处理，见图2（d）。经过上述处理的二值化图像，有较清晰的边界，接下来根据最小二乘法拟合毛细管的内外垂直边界，然后根据计算出的毛细管的左右边界，进行图像的校正。下面重点描述最小二乘法拟合毛细管的内外边界。

图2 图像处理

2.2 最小二乘拟合

在科学实验分析中，常常需要从一组实验的观测数据来确定变量之间的一个近似解析表达式。一般采用拟合的办法，即在确定的函数类Φ=S pan{φ0(x),φ1(x),…,φn(x)}中，求一个函数φ(x)∈Φ“最佳”地拟合已知的数据(xi,yi)(i=0,1,…,N)。所谓“最佳”的标准通常是要求φ(xi)与 yi之间差的平方和，即

取得最小值。按即采用最小二乘拟合[11]函数。

更一般地，最小二乘拟合是求φ(x)∈Φ，使

最小，式（2）中 wi＞0(i=0,1,…,N)是给定的常数，称为权，它表示数(xi,yi)的比重。

若φ0,φ1,…,φn是线性无关的，那么对于任意的φ∈Φ，都可由{φ0,φ1,…,φn}线性表出，即

于是问题化为求出这样的a0,a1,…,an，使δ2取得最小值。

若用向量、矩阵记号，令

那么，式（1）可改写为：

而式（2）可写成：

其中W=diag{w0,w1,...,wN}是N+1阶对角矩阵。

这样就把数据的最小二乘拟合归结为矩阵论中的最小二乘问题。不难验证，a是最小二乘问题公式（3）的解的充分必要条件是：它满足正规方程组

求数据的最小二乘拟合函数的关键，是建立并求解正规方程组，而这又取决于基函数。本文所感兴趣的是毛细管的内外边界线，其数学解析表达式为直线，其中φ0(x)=1，φ1(x)=x。

2.3 拟合毛细管边界

为了去除毛细管以左及以右的噪音，需要将找到毛细管的左右边界，本文通过最小二乘拟合法计算出毛细管的左右边界的函数表达式。

2.3.1 拟合毛细管左边界

从左到右对图3（a）二值化图像进行扫描拟合毛细管的左边界。

图3 毛细管内、外径边界

具体的方法是：

（1）采用间隔行取点法，减少噪声

由于噪声的分布是随机的，而毛细管的左边界的图像是按一定直线连续的，所以采用间隔行取点法可以提高拟合曲线的精度。

（2）从左向右扫描，找到这样第一个点，它满足在其±5个像素点范围内灰度值为零的个数多于某一个阈值（经验值）。

（3）为了减少墙壁上固有纹路对实验的影响，采用从不同的列开始按（2）取拟合点，在拟合出的多条直线中，将直线附近存在灰度值为零最多的那条确定为毛细管的左边界直线。由于从固定的列开始选取拟合点，所以在用CCD拍摄图像时毛细管要尽可能的垂直。对于这幅512像素×512像素的图像，从不同的列（具体的是第5列、第60列、第160列）开始按上述方法（1）（2）分别拟合出3条直线。然后计算在这3条直线±1范围内灰度值为零（即图像上的黑色）的像素点的个数，最多的就是所求的直线，从而提高拟合曲线的准确度。图4为分别从第5列、第60列、第160列拟合的左边界，其中圈出的部分是拟合到的毛细管的左边界。通过计算，图4（c）为所求直线。

图4 从不同的列开始拟合的效果图

2.3.2 拟合毛细管的右边界

以上述拟合出的毛细管左边界加10个像素点为出发点，以毛细管左边界加100个像素点为终点，在这个范围内从右向左扫描，按照拟合左边界方法的（1）和（2）得到毛细管右边界的拟合直线。

2.3.3 拟合气泡的左边界

以拟合出的毛细管左边界加10个像素点为出发点，以毛细管左边界加20个像素点为终点，在这个范围内从左向右扫描，按照拟合左边界方法的（1）和（2）得到毛细管右边界的拟合直线。

2.3.4 拟合气泡的右边界

方法同拟合气泡的左边界。只是以拟合出的毛细管右边界为基础向左加10个像素点为出发点，以毛细管右边界向左加20个像素点为终点，从右向左扫描，拟合出气泡的右边界。

2.4 气液两相流图像的校正

根据已经得到的毛细管和气泡的左右边界，对图像进行校正，剔除毛细管左右的噪音，为后续一阶导数法确定气泡的上下边界打下基础。

（1）在图3（b）中得到的四个边界方程（包括毛细管的左右边界，气泡的左右边界）选择一个较准确的斜率。理论上，这四个方程的斜率应该是一致的，但实际中，由于光线和拍摄角度的关系导致在拟合直线时会产生一些误差，可以在四个斜率中选择一个最接近的斜率。

方法是：分别用四个斜率作为毛细管左边界的斜率，后统计在这个斜率周围5个像素点范围内灰度值为零点像素点的个数，以个数最多的斜率定为最接近的斜率。

（2）将拟合出的毛细管左右边界内的图像保留，并校正为垂直线，其余部分置为255，即为白色，如图5所示。

图5 经校正后的图像

3 一种基于一阶导数的确定气泡和水柱上下边界方法

3.1 局部极值点

在用函数的方法解决理论和实际问题时，常要分析函数的增量Δy，在Δx→0的条件下，d y=y′d x是与之等价的量，在分析时可用以代替Δy，从而显示了微分的作用。但当Δx是个有限的改变量（不是无穷小）时，一般d y与Δy不相等，而只是近似。这样，在进行数学的分析时，由于没有等量关系就很难发挥作用。为了解决这个问题，在函数y=f(x)满足适当加强了的条件后，对自变量的有限增量Δx（为强调其有限而不是无穷小，有时就记Δx=b-a），对应的函数增量成立等式为：

因这里ξ是介于a，b之间的某个中间值，故常称这个结果为中值公式，也常称为拉格朗日公式或有限增量公式。

（1）局部极值点[12]

若函数 y=f(x)在点 x0处存在邻域 N(x0)，使 f(x)在每一点x处有定义，且对一切x∈̂(x0)成立：

则称 f(x)在 x0有一（局部）极小值（或极大值），称 x0为f(x)的（局部）极小点（或极大点）。

（2）局部极值点的必要条件

根据费马定理，若x0是 y=f(x)的（局部）极值点，则当f′(x0)存在时必有 f′(x0)=0，则局部极值点的必要条件是函数 f(x)的极值点必定是它的驻点或不可微点。

（3）单调性

对于在一个区间上连续、可微的函数 y=f(x)，利用其导数可简单地判明它在这个区间上的单调性。

设函数 y=f(x)在闭区间[a,b]连续，在开区间(a,b)可导，若

则 f(x)在[a,b]上严格单调上升（下降）。

（4）利用单调性判定法得出驻点为极值点的1阶充分条件

设 y=f(x)在 x0的某一邻域内可导，且 f′(x0)=0，又若存在 x0的左邻域 N-(x0,δ)(x0-δ,x0)及右邻域N+(x0,δ)(x0,x0+ δ)，能使

①当x∈ N-(x0,δ)时 f′(x)＜ 0，而x∈ N+(x0,δ)时 f′(x)＞0，则x0为 f(x)的局部极小点。

②当x∈ N-(x0,δ)时 f′(x)＞ 0，而x∈ N-(x0,δ)时 f′(x)＜0，则x0为 f(x)的局部极大点。

③当 x∈N(x0,δ)时，f′(x)保持符号不变，则 x0不是极值点。

3.2 利用灰度平均值的一阶导数确定气泡和水柱上下边界

再次分析经过校正的灰度图像。由于在摄像时光源是从上斜向下照射的，所以从图像上可以看出在气泡的上端较亮，而水柱的末端较暗（如图5所示）。根据这个特点，在有效的范围内分析每行的灰度平均值。图6是以图5的行为横坐标，以图5每一行的灰度平均值为纵坐标。从图6可知，尽管每行的灰度平均值都不一样（由于照明是自上而下的，所以使得整体图像表现为上面亮下面暗），但通过考察它的一阶导数（见图7）可知，存在着一定的规律。即，迎光部分表现为灰度值单调上升后单调下降，而背光部分表现为单调下降后又单调上升，见图7圈出的部分。这说明在这个范围内存在着极值点，根据费马定理，这个极值点应该是一阶导数为零点。为此，考察这个灰阶图像的一阶导数情况。

图6 每行取灰度平均值的曲线图

图7 灰度平均值的一阶导数图

尽管照明具有不均匀性，但根据光照的连续性，考察一阶导数的情况（由于为离散函数，所以考察的是后一行灰度平均值与前一行灰度平均值的差值）。图7即为图6的一阶导数图。由图7可知，大部分的差值是围绕着零点的（这体现了光照的连续性），只是有些点部分是围绕着零点震动（图中圈出部分）。分析这些过零点部分，发现其一阶导数先是大于零后小于零的部分对应气泡的亮的部分，反之为暗的部分。

确定气泡上下边界的流程图如图8所示。

图8 确定气泡上下边界流程图

4 测量结果

4.1 气泡上下边界的测量结果

经过上述步骤，测得每个气泡的上下边界，见图9。

在图9中“亮”表示气泡迎光部分的所在行号，“暗”表示气泡背光部分的所在行号，未编号的横线为修正的迎光部分或背光部分。其中图9（a）背景为复印纸，光源为荧光灯光，图9（b）背景为墙壁，光源为荧光灯光。通过比对图9（a）和（b），虽然背景使用复印纸也可以清晰地测量出每个气泡的上下边界，但对灯光的角度有特殊的要求，不易掌握，所以实验以反射均匀的墙壁作为背景。

图9 气泡上下边界示意图

4.2 测量结果

根据4.1节的结果，再根据标定的毛细管的内外径与在图像中测到的毛细管的内外径的比例关系，可以根据：

其中r为气泡的半径，r1为毛细管的内径，l为气泡的长度。计算出气泡的直径，进而得到气液比、含气率、气泡平均直径参数。图10为连续1 000幅图像气泡直径的概率密度分布曲线，其中横坐标d为气泡直径，单位为毫米，纵坐标 f(d)为气泡直径的概率。从图10看出，所测得的气泡直径小于5 mm，符合玻璃管内径为5 mm的限制，且气泡大小分布基本呈现正态分布，也和气泡随机产生的这一特性相符，说明系统测量方法测得结果在可信范围内。

图10 气泡直径的概率密度分布曲线

重复性是测试系统可靠性的重要指标，在相同的光照（高频荧光灯）、相同的背景（墙壁）、相同的气泡速率（1 m/s以上）条件下，表1给出了5次测量（每次200幅图像）结果。其中气泡的平均直径在3.2～3.4 mm之间，气泡平均直径的随机误差为5.47%，含气率在50%～54%之间，可见本系统具有良好的重复性及可靠性。