基于机器视觉的流体速度量化研究

陈金鹏 李贵山 李双科 马宏锋

(兰州工业学院电气工程系，甘肃 兰州 730050)

0 引言

不管是普通的流体视觉化，还是由计算机辅助的流体视觉化，它们在流体动力学中对于观察和监控流体结构都是非常有用的工具。流体视觉化通过在流体中引入染色物或烟气形式的污染物(即被动示踪体)，使流体的结构变得可视。

根据不同流体的密度，只有选择合适的污染物，才能确保很好地跟踪流体的流动，视觉化才能集中于被动示踪体以及它在时间和空间上的浓度的变化，从而识别和确认特定流体的特性。

对于用数码相机采集的流体/被动示踪体流动的图像，可以进一步对其进行处理，以用于不同目的的研究，如涡流分离分析［1］、回流区位置确认［2］、监控气蚀结构［3］、流体流动对沉积的影响［4］等。在这个阶段出现的问题是能否借助计算机视觉获取流体的定量特性，也就是说，能否采用被动示踪体将速度场这样的流体运动量从一组流体图像中抽取出来。

目前，流体速度场测量主要采用粒子图像速度测量仪。该测量仪通过在流体中引入固体微粒，采用机器视觉的方法进行测量。借助机器视觉对流体进行定性和定量的研究，至少应同时使用两种方法或系统［5－6］，但这样的组合系统往往很复杂。此外，大部分普通的测量方法都不能测量非平稳或周期性流体结构的流体速度。

本文仅通过采用视觉方法对加入示踪体的流体的速度场进行测量，并在实际使用中检验这种方法的可行性。该方法是利用一种已知的物理关系，把污染物浓度和流体速度结合起来的方法。

1 研究方法

在所描述的研究方法中，最重要的步骤是把速度向量场与浓度的标量场联系起来，后者用灰度图像表示。通过浓度N(即单位体积流体中污染体的分子数)表征污染体引入流体的情况。如果采用如下速度供给污染体分子，则可以使x方向上N的梯度保持不变。

式中:Φ为X方向的扩散流量;D为分子扩散系数。

式(1)为菲克扩散定律。如果考虑非稳态扩散过程，基本体积中浓度N的变化率可以用流出和流入该体积流体边界的污染体的微分表示，则式(1)可以改写为［7］:

当有流体运动(即对流)时，用总的导数代替时间偏导数，则式(2)可进一步展开为:

式(3)即为著名的对流扩散方程，它表示示踪体/污染体与流体动力学的基本关系。其中，等式左边第二项包括流体对流的影响。这样，如果浓度N的值(即N的导数已知)和流体中污染体的扩散系数D已知，则唯一的未知量就是速度v，即只有xi方向上的分量vi为未知量。

分子扩散系数D是物质的属性。在一定压力和温度下，D与另一种物质有关，污染体浓度N在不断变化，因此，应该即时测量。在这种情况下，应该将它的值从流体灰度图像中抽取出来。

引入跟踪物的流体灰度图像可形成具有不同灰度的像素的矩阵。一般来说，灰度的整数值从0(黑)到255(白)，它可以在0～1区间被归一化。如果引入流体的污染物是用外光源照明的，那么它的灰度值一般比流体的灰度值高。流体中污染体的浓度越高，在灰度图像中污染体的灰度值越高。因此，有:

式中:A为灰度图像中选定像素窗口的平均灰度值。文献［8］和文献［9］对此已有论述，式(4)中灰度值和污染体浓度的比值是线性的。选定的像素窗口的平均灰度可以用如下关系［2］得到:

式中:E(i，j)为l×m像素窗口中实测的(i，j)点单个像素的灰度强度。

浓度的时间导数∂N/∂t可以利用如图1所示的连续图像中固定窗口的平均灰度计算得到。

图1 流体连续图像窗口Fig.1 Continous image windows of the fluid

利用已知的两连续图像的时间间隔Δt(通过相机的图像采集频率得到)，则∂N/∂t项可以用式(6)估算:

式中:ΔA为固定窗口中相距Δt的两连续图像之间平均灰度的差值。污染体浓度在空间上的偏导数∂N/∂xi可以通过计算流体的单个图像获得。为了得到平均灰度A在空间上的差值，将现在的窗口向选定的xi方向移动，其示意图如图2所示。

图2 在窗口中移动单个图像示意图Fig.2 Diagram of moving single image in the window

空间导数可以用标准数字技术完成。式(7)为忽略截断误差的中心差分法［10］:

式(3)在两维系统中表示两个线性微分方程，它有两个未知数，即两个速度分量。求解这样一组方程需要明确的初始条件，非常繁琐和困难，因为流体的运动特性一般都未知。但是，在速度分量和其空间导数仍未知的情况下，这组微分方程可以转换为线性方程组。这组方程即使是超定的(即方程的个数多于未知数的个数)，也可以通过用数字的方法很容易地解出。与微分方程组相比，线性方程组的未知数的个数增加了一倍，这就意味着对两维系统至少需要四个一般方程。每一个方程需要两个流体连续的图像(因为要计算时间差)，即至少需要五个连续的流体图像才能得到流体的两个速度分量的信息。

上述研究方法是基于已知的物理关系，反映了借助机器视觉确定流体动力学特性的基本原理。它的有效性可以用传统的测量方法验证。为了达到这个目的，采用热丝风速计进行试验验证。

2 试验结构

测量来自垂直管空气射流速度场的测量装置如图3所示。

图3 测试部分的试验装置Fig.3 Experimental device of the test section

流体的雷诺数在管道的出口处进行测定，记为Re=1300。由于管道中的流体是层流，可以采用热线风速计法(hot-wire anemometer，HWA)和计算机辅助视觉法(computer-aided visualization，CAV)完成测量。在选定的测量点，用两分支热丝风速计测定速度的分布。两分支热丝风速计测量流向(y)和侧向(x)瞬时速度波的步骤与 Bruun［11］和 Jφrgenson［12］的研究步骤一致。在9×17网孔中选定153个测量点，每两个点之间的距离用定位表设置为5 mm。

为了用热丝风速计测量流体速度，本文采用DANTEC 55P62恒温热丝风速计传感器系统。该传感器用镀铂钨丝制成，长度和直径分别为1.25 mm、φ5 μm。风速计放大器的截止频率设置为10 kHz，风速计的工作温度为250℃。风道中的两分支热丝风速计传感器的位置要求能从两个方向测量瞬时速度。因此，风速计传感器采用PC控制的精确定位装置进行定位。热丝风速计的信号采用NI公司的16位数据采集板获得。数据采集时间为6 s，采集频率为50 kHz。

风速计的输出在传送到具有SCXI模块的A/D转换器之前，先用频率为10 kHz的低通四阶贝塞尔滤波器予以滤波。系统采用LabVIEW软件采集并存储数据。热丝风速计的校准是在专门校准风速计的测量台上完成的。为了与国王定律［11］保持一致，对第一和第二个热丝的常数设定如下:

此外，在校准和测量过程中还对温度进行了补偿。

温度的测量采用A级四线电阻温度计Pt100和Agilent 34970A数据采集器。根据国王定律方程式［11］和Jφrgenson方程［12］，从热丝风速计的输出计算得到实际的流向流速vy和侧向流速vx。测量也是用同样的设备完成的，且在测量之前进行了校准［13］。测量的总误差［14］是由工作点的选择和波动的不确定性、定位起始点的选择、风扇的转速、定位表对热丝温度传感器的定位、线性化、A/D的分辨率、数据采集、温度补偿、湿度以及热丝风速计的频率响应有限等因素引起的。风扇由一台分离变压器供电，便于与测量设备的供电电源分开。瞬时速度的整体测量误差估计为测量值的2.8%。

在CAV测量中，为了使流体可视化，在充分发展流中引入含有气化石蜡油的被动示踪烟气［15］。测量点与上述热丝风速计的测量点完全相同。通过流体入口将被动示踪体引入管道。相机放置于与风道的测量段垂直并且距离测量段表面1 m的位置。相机的图像采集频率为400 Hz，拍摄图像的快门速度为0.6 ms，所记录的图像具有8位灰度值、600×250像素的分辨率。本文采用具有直线导光的不间断无闪烁光源(Vega Velum DC 150 W)进行照明;利用LabVIEW软件包实现对相机的设置、图像的采集和保存。在每一个测量点采集1000个连续的图像，图像的个数受被动示踪体烟气的发生时间限制。被动示踪体瞬时浓度的总体测量误差估计为测量值的 3.8%［15］。

3 试验结果

试验中采集的图像序列用Matlab程序包进行处理。两个速度分量的测定与试验测量(热丝法)的情况一样，都是在同样的测量点、在x-y平面利用机器视觉测定的。根据文献［16］，空气(T=23℃、p=100 kPa)中的癸烷(C10H22)蒸汽 D=7.01 ×10－12m2/s。用 30 个连续图像序列来确定x-y平面上的流速场，序列中图像的数量保证了方程式(3)即使在被动示踪体稀少(测量平面x-y外沿)的情况下，对所有的测量点都有效。另一方面，当图像序列很大时，求取速度分量的平均值需要非常多的计算时间，有可能导致给定采集频率上的信息丢失。为了减小导数数值计算的截断误差，式(7)中步长Δxi应选得尽可能小。因此，步长接近于1个像素，它的大小取决于图像分辨率。

本文采用HWA和CAV两种方法测量流速场。测量结果显示，在大部分测量点，HWA法和CAV法的测量结果一致。

两种方法的测量结果表明，离开节流孔的射流流速减小，射流外边沿的流体沿径向朝外侧运动(射流沿轴向和径向散开)。尽管如此，两种方法所测得的速度矢量的方向和大小都有差别，尤其在射流的外边沿。为了正常工作，序列中的每个图像都要有足够的被动示踪体浓度值，而在这些位置，由于石蜡油蒸汽很少，影响了CAV法的效果。

对两种方法测得的速度分量分别从流向方向(y)和正常方向(x)进行定量分析。用x方向和y方向的归一化值来表示结果，这里 x*、y*∈［0，1］。x*=0 和 x*=1分别对应于图3最左边和最右边的测量点;类似地，y*=0和y*=1分别对应于距离管道节流孔最远和最近的测量点。

经归一化流向，当观测点y*的值分别为0.125、0.375、0.625 和 0.875 时，两种测量方法所测得的流速分量vy的值如图4所示。

图4 不同流向观测点流速分量vyFig.4 Velocity component vyin different directional observing points

通过比较两种测量结果可知，不管采用何种测量方法，y方向的流速具有相似的外形。但是，采用CAV法时vy的波动比HWA法的测量结果要大得多，在远离管道节流阀处表现得更为突出。这种现象很可能是流体的漩涡和特殊测量方法的数据采集频率共同作用的结果。当选定雷诺数，且观测点y*＞0.3时，管道中的流体形成了强烈的漩涡。

HWA法的数据采集频率大约为50 kHz，在每个点的测量历时约2 s。因此，HWA法获得的结果是这段时间的平均值，这个时间足够长，足以保证大漩涡的旋转对流体速度的平均值影响不大。对于同样的时间段，一方面CAV法大约需要700个图像序列，这将极大地影响处理/计算时间;另一方面，这种情况下CAV法的采集频率仅为HWA法采集频率的1/140。为了获得更精确的结果，应该增加图像序列的数量，但该方法只适用于流体平均速度(波动)偏差较大、持续时间足够长且被CAV法检测到的情况。即与HWA法相比，CAV法需要相机具有较高的采集频率，两者才有更好的可比性。

根据HWA法，图4中每一个流向观测点的vy(沿侧向距离x的计算值)平均值的相对差如表1所示。

表1 vy平均值的相对差Tab.1 Relative difference of vyaverage values

表1表明四个选定的流向观测点的平均值相对差大约为10%。因此，可以得出结论:与HWA法相比，即使在图像序列很短的情况下，CAV法也能很好地预测vy的值。

根据CAV法和HWA法，vx平均值的相对差比较如表2所示。

表2 vx平均值的相对差Tab.2 Relative difference of vxaverage values

表2中，vx平均值相对差比表1中vy平均值相对差稍大一些。其主要原因是:与vy的值相比，vx值更小(接近于0，最多到0.04 m/s)。由于相对差再大也不会超过15%，所以，对基本CAV法而言，这可以认为是较好的结果。

4 结束语

本文阐述了借助于机器视觉对流体动力学性能量化分析的基本原则。通过流体中污染体(被动示踪体)的连续图像序列，测量了流体的两维速度场。这是一种通过对流扩散方程把速度场和浓度场联系起来的方法。所需参数的计算是用图像处理和普通的数字技术来实现的。该方法本身是用传统的热丝风速计法与机器视觉法同时进行而予以验证的。

两种测试方法的结果对比表明:视觉法(即使最基本的形式)能够相当精确地计算出流体的大小和方向;两种方法在不同流向观测点的速度分量外形是相似的，两种结果的差别非常小。由此可以得出利用机器视觉测量流体速度场是可行的结论。

从现阶段来看，视觉法的精度劣于传统方法如HWA、PIV和LDA。因此，未来的研究不应该仅仅专注于数字计算领域(如步长大小、不同导数计算方法的应用对最终结果的影响等)，而是更应该专注于测量准备工作领域(如照明、三维流体结构、非定常性、序列中图像的个数等对结果误差的影响)。根据现有的算法和设备要求(合适的图像分辨率和采集频率)，视觉方法仅可应用于有很高雷诺数的湍流，具有一定的有效性和局限性。

［1］Sirok B，Potocar E，Novak M.Analysis of the flow kinematics behind a pulsating adaptive airfoil using computer-aided visualization［J］.Strojniski Vestnik-Journal of Mechanical Engineering，2000，46(6):330－341.

［2］Sirok B，Bajcar T，Dular M.Reverse flow phenomenon in a rotating diffuser［J］.Journal of Flow Visualization and Image Processing，2002，9(2＆3):193 －210.

［3］Sirok B，Blagojevic B，Bajcar T，et al.Simultaneous study of pressure pulsation and structural fluctuations of a cavitated vortex core in the draft tube of a Francis turbine［J］.Journal of Hydraulic Research，2003，41(4):541 －548.

［4］Bajcar T，Blagojevic B，Sirok B，et al.Influence of flow properties on a structure of a mineral wool primary layer［J］.Experimental Thermal and Fluid Science，2007，32(2):440 －449.

［5］Hu H，Saga T，Kobayashi T，et al.Analysis of a turbulent jet mixing flow by using a PIV-PLIF combined system［J］.Journal of Visualization，2004，7(1):33 －42.

［6］Reungoat D，Riviere N，Faure J P.3C PIV and PLIF measurement in turbulent mixing［J］.Journal of Visualization，2007，10(1):99 －110.

［7］McComb W D.The physics of fluid turbulence［M］.Oxford:Clarendon Press，1996.

［8］Aider J L，Westfreid J E.Visualizations and PDF of the fluctuations of a passive scalar in a turbulent Goertler flow［J］.Experimental and Numerical Flow Visualization，1995，218:123 －130.

［9］Simoens S，Ayrault M.Concentration flux measurements of a scalar quantity in turbulent flows［J］.Experiments in Fluids，1994，16(3 －4):273－281.

［10］Mathews J H，Kurtis D F.Numerical methods using Matlab［M］.4th Edition.London:Prentice-Hall，2004.

［11］Bruun H H.Hot-wire anemometry-principles and signal analysis［M］.New York:Oxford University Press，1995.

［12］Dantec Dynamics A/S.How to measure turbulence with hot-wire anemometers-a practical guide［EB/OL］.［2002 －02 －01］［2012 －03 － 02］.Http://www.dantecdynamics.com/Admin/Public/Download.aspx?file=files%2fsupport_and_download%2fresearch_and_education%2fpracticalguide.pdf.

［13］Coleman H W，Steele W G.Experimentation and Uncertainty Analysis for Engineers［M］.New York:John Wiley ＆ Sons Inc.，1999.

［14］Eberlinc M，Sirok B，Hocevar M，et al.Numerical and experimental investigation of axial fan with trailing edge self-induced blowing［J］.Journal of Engineering Research，2009，73(3):129 －138.

［15］Eberlinc M，Dular M，Sirok B，et al.Influence of blade deformation on integral characteristic of axial flow fan［J］.Journal of Mechnical Engineering，2008，54(3):159 －169.

［16］Won D，Shaw C Y.Determining coefficients for mass-transfer models for volatile organic compound emissions from architectural coatings［C］//Proceedings of the 7th International Conference of Healthy Buildings，2003:1－6.