李 莹，周云龙，范振儒

(1.东北电力大学能源与动力工程学院，吉林吉林132012;2.国家开发银行内蒙古自治区分行，呼和浩特010010)

气固两相流动广泛存在于工业生产过程中，如煤粉燃烧、气力输送、烟尘排放等，自然界的沙尘暴、宇宙尘埃也属气固两相范畴。同时，工业上有许多场合需要研究两相流中固体颗粒的空间流动速度和方向，比如循环流化床锅炉内中的煤粉流速，炉内各局部空间的煤粉及灰颗粒的运动速度、方向。然而气固两相流动十分复杂，参数众多，现有的理论模型还不能完整地阐述其流动变化规律，那么首先要解决的是两相流体的参数的检测。在众多流动特征中流型及其流场、涡量场和空间速度分布等参数的准确测量会使人们更加深入地认识这种流动现象的物理本质，更重要的是，它有助于发展相关的理论数学模型，建立有效的理论计算方法，这类计算方法的建立，对工程设计人员发展更高效的工业设备非常重要。有关这方面的研究具有相当的难度。迄今为止，测量两相流中颗粒流速及流场的方法主要有激光多普勒测速法、分子示踪全场测速法、粒子示踪全场测速法和图像测速法等［1］。其中当今较流行的流体动态速度场检测手段有PIV和PTV技术。PIV法［2］计算某一诊断窗口内粒子的平均位移，一定程度上会产生平均效应、实验步骤复杂方法，计算繁琐。PTV法［3］适用于稀疏颗粒运动场合，如在输移现象、传送现象等尤其关注拉格朗日运动规律的情况下具有很好的发展。在两相流体流场及速度检测方面，周云龙等［4］应用高速摄像机对油气水三相流的水包油流型进行了拍摄，并跟踪气泡及油滴，实现对流速与流场的测定。朱佳琪等［5］提出的图像测速法是以直流灯作为正面光源，采用高速CCD拍摄流化床内的二维颗粒运动图像，基于图像的互相关算法计算得到循环流化床内颗粒二维全场速度。Antonio Busciglio等［6］利用图像分析技术，对流化床气固两相流的鼓泡床中气泡行为进行了研究，并利用拉格朗日测速法和欧拉测速法对鼓泡床中的流场和相关速度参数的进行了测定。周洁等［7］利用光信号互相关方法对气固两相流中固体颗粒一维流速的测量进行了研究。参考文献［8-12］是我国学者利用PIV技术检测流体流场和速度所取得的成果。

本文提出了一种两相流流场及速度测试的新方法。实验中，光源的位置、亮度恒定不变，摄像机镜头的拍摄位置及焦距不变，检测区域内图像灰度值变化是因两相流体运动而引起的，故可通过跟踪检测区域内光流场的变化来检测两相流体的流场。该方法对数字图像建立整体光流模型，计算出全场像素点的速度矢量，即得到两相流体的速度场。实验结果表明:这种方法可以准快速准确地反映两相流体的变化规律，具有深入研究的价值。

1 实 验

本实验装置是在东北电力大学的透明有机玻璃流化床装置系统上完成的，实验系统如图1所示。该实验装置主要包括两部分，即流体控制系统和图像采集系统。流体控制系统主要由空压机、风室、布风板、流化床主体和布袋分离器组成。实验工质采用空气和透明玻璃珠，空气经空压机升压和孔板流量计计量后进入风室，再通过布风板装置分散均匀的在流化床主体内垂直向上流动，气体-固体在流化床主体内有效混合并使固体颗粒流态化，最后进入布袋分离器，将空气分离出来并排入大气，分离后的固体颗粒保留在布袋内供循环使用。图像采集系统主要包括照明系统和高速摄影系统。由于高速摄影机对光线的亮度有较高的要求，照明系统的光源使用6400K色温的三基色光管，光线明亮无闪烁。由于两相流流型变化复杂，高速摄影系统采用瑞士WEINBERGER公司研发的 SpeedCam Visario系统，其最大分辨率为1536×1024，最大帧频达到 10000帧/s，能够清晰的抓拍各种流型的瞬变图像。在图像摄取过程中采用逆光照明，拍摄各种流态的阴影［13］。为了使光线分布均匀，可在观测区域后侧的有机玻璃板上，蒙上两层绘图用的硫酸纸，可获得满意的拍摄图像。

图1 流化床气固两相流实验系统

本实验装置中流化床主体尺寸为200 mm ×50 mm×2000 mm，布风板上布有3排共有29个风帽。每个风帽同一横截面上均匀的开有8个直径为2 mm的小孔。实验床料使用透明玻璃珠颗粒，颗粒的堆积密度为2.4～2.6 g/cm3，成圆率为80%，其筛分比见表1。

表1 颗粒的筛分比

距布风板上方10～560 mm的空间区域，采集图像像素为1020×340，帧频为500帧/秒。大气温度为 28.6 ℃，风流量分别为 25 m3/h，40 m3/h，65 m3/h，90 m3/h，110 m3/h。拍摄在上述工况下，流化床气固两相流的运动情况。

2 MQD互相关算法

把两相流体运动流场的两幅图像看成是两个随机分布的数据集合，便可用相关函数来度量两幅图像的相关性，即相似程度。

图2(a)为t时刻参考目标像位置，图2(b)为t+△t时刻候补目标像位置，箭头指示参考目标在△t间隔内的位移。f(x，y)和g(x，y)分别为图2(a)中目标区和图2(b)中同一位置目标区的灰度分布函数。当图像用像素来构成时，取fij和gij为像素的灰度值，N和M为横纵方向像素数量。灰度分布函数f(x，y)和g(x，y)分别表示成fij和gij的集合。MQD互相关算法［14］是灰度分布相关法的改进算法，差别就在于它用灰度差的平方代替相关系数作为相似性分析的依据，即相关系数:

当Cfg值为最大时，参考目标与候补目标为同一粒子，这样就可确定粒子的位移及其速度。以下就是采用MQD互相关算法具体的实现过程。

(1)首先对两相流动实验中获得的连续两幅图像进行检测，保证两幅图像中的像素矩阵维数相同，并根据浓度将两幅图像划分均匀网格。确定搜索区域最大速度。

(2)选定图2(a)t时刻目标像中一个中心坐标为x，y大小为M×N的区域f，取出f的灰度矩阵F和时刻t+△t相同位置的相同大小的灰度矩阵G，进行一次式(1)的计算，并将此值作为相关矩阵C(m，n)的中心元素。

(3)改变 g(x，y)为 g(x*，y*)，再进行和上面一样的一次运算得到C(m，n)另一个元素，循环操作上述过程，直到计算出在搜索最大速度对应的位移区域内g(x，y)周围所有可能存在的相关元素，组成矩阵C(m，n)。

互相关MQD算法原理图

(4)对矩阵C(m，n)进行峰值检测，假定相关小区域内峰值在亚像元精度内分布符合高斯曲线，找到C(m，n)的最大值Cmax，将C矩阵的中心元素对应的位置和Cmax对应的位置相减，即可获得△t时间间隔内目标位移。

(5)利用目标位移与△t的比值就可得到目标运动速度。并可绘制出速度矢量图。

3 光流法分析原理

光流的概念是Gibson于1950年首先提出的。所谓光流是指图像中模式运动的速度。光流场是一种二维瞬时速度场，其中的二维速度矢量是检测区域中可见点的三维速度矢量在成像表面的投影。光流不仅包含了被观察物体的二维运动信息，而且携带着丰富的有关景物三维结构信息。在各种不同问题中，光流扮演着重要角色［15］。

3.1 运动场和光流场

物体点的运动除了导致其对应的图像点的运动外，同时也使在图像上的对应物体的亮度模式运动。或者说，人与被观察的物体发生相对运动的时候，被观察物体表面带光学特征部位的移动给人提供了运动以及结构的信息。当镜头与场景目标之间发生了相对运动的时候，所观察到的亮度模式运动就被成为光流，或者说，物体带光学特征的部位的移动投影到视网膜平面上就形成了光流。光流表达了图像的变化，包括了目标运动的信息，可以用来确定观察者对目标的运动情况。光流有三个要素:一是运动(速度场)，这是光流形成的必要条件;二是带光学特征的部位(例如有灰度的象素点)，他能携带信息;三是成像投影(从场景到图像平面)，因而能够被观察到。

采用光流分析的方法可以确定图像点上的运动方向和运动速率。基于光流的图像分析的直接目标是确定运动场。光流又不总是对应于实际的运动场，光流与运动场虽然有着密切的关系但是又不完全对应。场景中的目标运动导致图像中的亮度模式运动，而亮度模式的可见运动又产生了光流。在理想情况下光流与运动场相对应，但是在实际中也有不对应的时候。也就是说:运动产生光流，因而有光流就一定存在着运动，然而并不是说有了运动就一定能够产生光流。不过在绝大多数的情况下，光流与运动场都是相互对应的。所以在许多情况下我们可以根据光流与运动场的相互关系由图像的变化来估计相对运动［16］。

3.2 光流约束方程

光流场是指图像亮度模式的表观(或视在)运动。设I(x，y，t)为t时刻第k帧上图像点(x，y)的灰度，则在亮度恒常性的假设条件下结合k+1帧图像可得到光流基本约束方程

式中:Ix为点(x，y)在帧间隔△t内x方向和y方向的速度，即光流矢量。在上面方程中，Ix，Iy，It可以直接从图像中计算出来。

光流约束方程式(2)含两个未知量，求解光流(u，v)尚需加上其他的约束条件，假定光流在整个图像上的变化平滑，也就是加上平滑约束条件［17］。根据光流约束方程，光流误差为

其中x=(x，y)T。对于光滑变化的光流，其速度分量平方和积分为

将光滑性测度同加权微分约束测量组合起来，有

其中，α是控制平滑度的参数，α越大，则平滑度就越高，估计的精度也就越高。使用变分法将式(3)转化为一对偏微分方程

用有限差分方法将每个方程中的拉普拉斯算子折换成局部领域图像流矢量的加权和，并使用迭代方法求解这两个差分方程。

下面只考虑离散情况。在一点(i，j)及其4邻域上，根据光流约束方程，光流误差的离散量表示为

光流的平滑量也可由点(i，j)与它的4邻域点的光流值差分按如下公式计算:

则极小化函数为

E关于u和v的微分是

从上面两个方程便可以求出u和v。实际中，经常将求解u和v表示成迭代方程

其中，n为迭代次数。只要迭代次数k足够大，就可以得到比较稳定的光流矢量(un+1，vn+1)。

4 测量结果分析

按照上述实验的方法和步骤，可以获取流化床气固两相流4种典型流型的运动图像序列，如图3至图6所示，其中(a)、(b)是典型流型连续两帧的运动图像，(c)、(d)分别是用MQD法和光流法计算出的速度矢量图，(e)是用光流法得到的速度矢量来计算得出的等涡量场。实验中帧频设为125帧/秒，所以两幅图像的时间间隔为1/500=0.002 s，大气温度为28.6℃，四种流型对应的分量风流量分别为25 m3/h，40 m3/h，65 m3/h，80 m3/h。如图3所示，在鼓泡床中，气泡主要是向上运动，在运动的过程中，伴随着气泡的扩大，床层的高度也有变化。图4所示节涌床中，流化床内有多个气泡，气泡的尺寸也随高度的上升不断增大，此时床层相对较高，并伴有床层蠕动。从图5所示湍动床可以看出，流化床内已无明显的气泡和床层，气相和固相可以明显的区分，气相部分在断续中有连接，形状复杂多变。图6所示为贴壁流，气相和固相依旧可以明显的区分，只是流化床两侧壁面的固相流体是沿壁面向下流动，流化床中间底部的固相向上输送。

从图3至图6的(c)、(d)中可以看出，采用光流法能计算得出的观察区域内的速度矢量图，其计算结果与MQD法基本相似，同时也与实际的流体流动变化规律相吻合，由此可以说明，光流法的可行性。

图3 鼓泡床

图4 节涌床

图5 湍动床

通过对比分析两种方法，可知光流法比MQD法具有一定程度的优越性。实验是在CPU为Pentium4主频为2.93 GHz和内存为1 GB的计算机上进行的，MQD互相关法计算耗时平均在3 3s以上，光流法耗时平均在1.5 s以内。可见，用光流法计算流体运动的流场，耗时短，可为今后在线检测流场提供一个可行的方案。通过比较图3鼓泡床和图4节涌床的两种流场图可知，两种方法都能正确的描述气泡、床层的运动情况。鼓泡床中，气泡逸出、破灭后，床层的变化趋势是，中间下落，两边有突起的迹象。从图中可知，MQD法只表示出中间下落的部分，而光流法能还能体现出两侧微微突起变化趋势。节涌床中，最高处的大气泡内部并没有太多的变化，并且气泡准备逸出，床层的右侧并没有向下流动的趋势，这些都是MQD法描述的与实际不相符的信息。

运用同样的方法，可继续比较湍动床和贴壁流的两种流场矢量图的效果。最后可以得出，在检测图像中流体整体运动时，光流法所描述的流体的流场与实际情况更接近，较MQD互相关法有较高的敏感性和鲁棒性。

图3至图6中的(e)描述的是流体的等涡量场。涡线的任意一点的切线方向与在该点的涡量方向一致，所以等涡量线可以直观地反映其涡量场的方向。由于数据量有限且速度矢量未经插值，只能定性的观察，若要涡量做定量的计算，需要有足够的数据量和速度矢量的插值。在此，本文不做过多的叙述。

沿水平轴方向将流场速度矢量图等分成若干个区域，并统计出各区域内流体运动的平均速度，可得流化床内流体在垂直轴的运动速度的分布情况，如图7至图9所示。取流体垂直轴的上升速度为正方向。

从图7可以看出，流化床内不同流型上升速度的空间分布情况。在曲线中可以看出在鼓泡床和节涌床中，气泡上升缓慢，整个流化床只有在有气泡的地方有上升速度，其他地方的上升速度近似为零。在湍动床中，流体在床中间区域运动较剧烈，两侧运动较缓慢，贴壁流则相反。

如图8所示，流化床内不同流型下落速度的空间分布情况。在曲线中可以看出在鼓泡床中，在气泡逸出后，床层回落，同时也就具有了下落速度。在节涌床中，下落速度主要源自于气泡上移后，尾部的脱落漩涡。湍动床和贴壁流下落速度与上升速度分布曲线近似。

在统计流体垂直方向的上升速度和下落速度的代数和基础上，得到了不同流型垂直轴的平均速度的空间分布曲线，如图9所示。此图可以描述流化床中气固两相流不同流型的在垂直轴上的总体运动情况，该方法为今后，进一步定量研究不同流型的运动情况提供了参考。

图7 流化床内流体上升速度分布图

图8 流化床内流体下落速度分布图

图9 流化床内流体平均速度分布图

5 结 论

(1)利用高速摄影系统直接检测流化床气固两相流的流场，发挥了高速摄像机的分辨率与帧频高的优点，同时解决了PIV技术的平均效应问题和实验步骤繁琐的弊端，是一种方便、有效的检测方法。

(2)通过跟踪检测区域内光流场的变化来检测两相流体的速度场，较MQD相关法检测出的结果，光流法测得的速度矢量场场更能真实的反映流体的变化情况，并且耗时短。

(3)利用光流法测得的速度矢量绘制4种典型流型的等涡量场，通过统计获得了流型在垂直轴上的上升和下落的平均速度以及总体速度的空间分布曲线，分别为定性和定量分析流化床内气固两相流的运动机理提供参考。

(4)实验证明，高速摄影系统结合光流法检测流化床气固两相流的流场及其相关参数的方法是可靠的，能够准快速准确地反映流化床内气固两相流动的变化情况。这种非接触的流场检测方法在实际的工业生产中有较强的适用性，为在线检测两相流体的参数提供支持。

［1］康琦，申功忻.全场测速技术进展［J］.力学进展，1997，27(l):106-120.

［2］J Westereel.Fundamentals of digital particle image velocimetry［J］.Measurement Science and Technology，1997，8(12):1379-1392.

［3］A.Clarke.The application of particle tracking velocimetry and flow visualization to curtain coating［J］.Chemical Engineering Science，1995，50(15):2397-2407.

［4］周云龙，李洪伟，范振儒.基于PTV法对油气水三相流流场的测定［J］.化工学报，2008，59(10):2505-2510.

［5］朱佳琪，马增益，严建华，等.图像法用于循环流化床颗粒二维速度场可视化的实验研究［J］.电站系统工程，2005，21(5):32-34.

［6］Antonio Busciglio，Giuseppa Vella，Giorgio Micale，et al.Analysis of the bubbling behavior of 2D gas fluidized beds.Part I.Digital Image A-nalysis Technique［J］.Chemical Engineering Journal，2008，140(1-3):398 – 413.

［7］周洁，袁镇福，岑可法，等.光信号互相关测量两相流中颗粒流动速度的研究［J］.中国电机工程学报，2003，23(1):185-188.

［8］虞建，李荣先，周力行.射流稀疏颗粒图像PIV算法应用［J］.化工学报，2005，56(7):1206-1208.

［9］翁文国，廖光煊，王喜世.基于互相关的DPIV图像诊断方法研究［J］.实验力学，1999，14(3):323-329.

［10］李水清，严建华，张志霄，等.基于DPIV技术回转圆筒内颗粒流场可视化研究［J］.中国电机工程学报，2002，22(4):56-60.

［11］石惠娴.循环流化床环核结构速度分布PIV测试［J］.水动力学研究与进展(A辑)，2006，21(1):8-12.

［12］孙立志.基于多尺度的光流算法在PIV中的应用［D］.大连:大连理工大学，2006.

［13］Ha Y J，Ltu Z C，Hanratty T J.A backlighted imaging technique for particle size measurements in two phase flows［J］.Experiments in Fluids，1998，25(3):226-232.

［14］Gui L，Merzkirch W.A comparative study of t he MQD method and several correlation-based PIV evaluation algorithms［J］.Experiments in Fluids，2000，28:36-44.

［15］王睿，张广军，阎鹏.基于光流分层方法的平面3D运动估测［J］.光学技术，2007，33(1):102-105.

［16］张广军.机器视觉［M］.北京:科学出版社，2005.

［17］Berthold K.P.Horn，Brian G.Schunck.Determining optical flow［J］.Artificial Intelligence，1981，17(1-3):185-203.