基于光流算法的粒子图像测速技术研究

王宏伟, 黄 湛

(中国航天空气动力技术研究院, 北京 100074)

基于光流算法的粒子图像测速技术研究

王宏伟, 黄 湛*

(中国航天空气动力技术研究院, 北京 100074)

光流测量技术作为一种新的空气动力学实验技术，以其像素级分辨率的矢量场测量优势获得广泛的应用。光流测量技术使用光流约束方程，配合平滑限定条件，可以进行速度场测量，获得高分辨率的全局矢量场。首先通过研究积分最小化光流测速理论和算法，采用C++编写光流速度测量程序，然后通过3种典型人工位移图像对光流计算程序进行验证，并将结果和标准位移分布进行比对分析，以指导如何在实际应用中获得高精度光流速度场，最后进行小型风洞后向台阶实验，利用高速相机拍摄示踪粒子图像，使用光流计算程序获得速度矢量场，同采用互相关算法的粒子图像测速计算结果进行比较，体现出光流计算方法像素级分辨率的矢量场测量优势。

光流；速度测量；粒子图像；运动估计；积分最小化

0 引 言

在空气动力学发展过程中，应用实验手段获得流体的空间运动结构，获取全场流动信息，是了解空气动力学诸多流动现象及机理的关键，同时这一需求也不断刺激空气动力学实验测试技术水平的提升。传统的流动测量方法大多采用单点测量或介入式测量手段，不能满足复杂流动情况以及全场速度测量的需求。

近年来，粒子图像测速技术利用拍摄跟随流场运动的粒子图像并通过相关计算获得速度场，并以其非介入式、全场流动速度测量、测量精度高及测速范围广等优点在全局速度测量方面具有应用优势。不过，粒子图像测速得到的速度矢量为其判读区域内粒群的平均速度，忽略了此判读区域内速度场的变化，这使得粒子图像测速方法的空间分辨率受到一定的限制，尤其对于速度梯度大、速度变化剧烈的区域，其速度场测量精度是受到制约的。

相机所记录的物体的图像，本质上是物体所反射的光的分布。当相机连续拍摄运动物体时，物体所反射的光就会在相机底片上形成一系列连续变化的图像，由这些图像所组成的连续变化的信息不断“流过”相机底片，如同光在流动，这是光流定义的物理基础。物体在光的照射下，其表面灰度会呈现出一定的分布，而所谓光流，就是指图像中灰度模式的运动速度。粒子图像测速技术中拍摄到的粒子图像显示的是粒子浓度的分布，具有一定浓度的运动粒子群通过其散射光在相机底片上成像，所以粒子图像测速技术中采用相关方法得到的速度场实际上是光流速度场(Optical Flow)，光流速度场能否代表真实速度场要取决于示踪物质能否充分地跟随流体运动。

Horn和Schunck早在20世纪80年代就开创性地进行了光流研究和计算[3]，Horn等认为同一运动物体引起的光流场应该是连续的、平滑的，即同一物体上相邻点的速度是相似的，那么其投影到图像上的光流变化也应该是平滑的，并提出了一种利用加在光流场上的附加速度平滑约束，即全局平滑约束将光流场的计算问题转化为一个变分问题，给出像素级分辨率的运动估计。随后，光流技术在发展过程中吸引了各个领域的科研人员进行研究[4-6]，也产生了很多其它解算方法，如图像相关测速法、基于能量的方法[7]、基于相位的方法[8]和神经动力学方法等，针对Horn等提出的微分法在计算较大像素位移时遇到的困难，Tianshu Liu等在其工作中提到[9]，Corpetti等人提出在时间步上采用积分形式的光流方程代替微分形式的变分方程[10]，Ruhnau等使用高斯滤镜生成二元图像金字塔的多分辨率策略[11]。这些算法有各自的优势，不过由于具有一定的适用范围和算法复杂度使其相对于微分法无法获得广泛应用。

本文引入光流计算方法，编写光流速度测量程序，通过人工位移图像对该程序进行验证和分析，以合理的测速范围和最优的加权值指导如何在实际应用中获得高精度光流速度场；进行小型风洞实验，以粒子图像为速度测量载体，使用光流计算方法获得像素级分辨率的速度矢量场。旨在通过这些研究，发展以光流技术为基础的全场流动测量技术，期望在提高空间分辨率和避免速度梯度影响、了解复杂流动的空间结构、获取更多的动力学信息方面获得进步，以此丰富实验流体力学测量技术。

1 基于粒子图像的光流速度场求解

Horn的光流场解算方法是微分法的基础，它以灰度守恒方程为基础，附加速度平滑约束条件，利用连续变化的图像上各点灰度的时间、空间梯度来计算每一点的速度矢量。

光流的灰度约束条件认为，图像上各点灰度随时间和空间变化满足如下关系：

式中：I(x,y,t)是t时刻(x,y)位置的图像灰度，I(x+Δx,y+Δy,t+Δt)为原(x,y)位置的物体运动到(x+Δx,y+Δy,t+Δt)位置的灰度，令u=Δx/Δt,v=Δy/Δt，将上式进行泰勒展开并忽略高阶项可以获得光流约束方程：

式中：Ix,Iy,It是灰度的二维空间和时间梯度。显然该方程是求解u,v2个未知数的不定方程，需要引入辅助约束方程进行求解。

速度梯度平方和判定量为

实际拍摄的图像灰度将会受到硬件条件和噪声的影响而偏离基本的光流约束条件，所以不能肯定地认为ξI恒为0，因此有必要找到合适的加权因子β，使得2个判定误差在二维空间上的积分为最小，即对

ξ2=∬

进行最小化。通过变分法建立每个因变量的欧拉特征方程：

2 光流速度场求解算法验证

基于微分法编写了光流速度场计算程序，拟用于风洞试验速度场测量。为了验证光流算法的精度，本文选取灰度图像，通过人工位移获得标准速度场，然后采用所编写的光流计算程序对该图像对进行处理，将光流计算结果与标准速度场进行比较。实际计算中，将采用1个单位的时间梯度，所得的速度值将直接反映图像像素位移值，同时，也可考察2张图像的像素位移对光流计算精度的影响。

首先考察第1组图像(见图1)，图像为600pixel×600pixel大小的风暴云图，通过处理将图像向右平移一个像素(Δx=1pixel)，通过光流速度场计算程序得到的速度矢量分布(Skip 16pixel×16pixel，下同)和流线如图2所示，可以看出，在1像素位移条件下，光流计算方法所得到的速度场分布十分均匀。

图1 灰度原图与平移一个像素后图像

图2 采用光流计算方法得到的速度矢量图和流线图

图3 横向速度分布和纵向速度分布分析

在图1中所示的纵向考察线和横向考察线(下同)上对光流速度场和标准速度进行分析如图3所示，可以看出, 在1像素位移条件下光流算法获得的速度值与标准速度几乎没有偏差。

第2组图像(见图4)采取了Δx=1 pixel和Δy=1 pixel的斜向移动，处理得到的速度矢量分布和流线分布如图5所示，速度矢量分布很均匀，流场平滑性很好，不过已开始出现小范围的速度偏差。

通过纵向考察线和横向考察线将所计算的速度场与标准速度进行比较发现(见图6), 在Δx=1 pixel和Δy=1 pixel的斜向运动条件下，所计算的速度值能够保证较高的准度分布，不过已开始出现一定量的抖动偏差。这种偏差很有可能与2个方向速度偏差的合成有关；同时也暗示随着连续图像间运动位移的增大，光流算法的计算误差也会随之增加。

图4 灰度原图与x,y各平移一个像素后图像

图5 采用光流计算方法得到的速度矢量图和流线图

图6 横向速度分布和纵向速度分布与标准速度比较

第3组图像(见图7)为中心旋转1°前后的2张图像，通过计算获得的速度矢量和流线分布图如图8所示，可以看出，图像中心的速度矢量分布较为均匀和平滑，边缘处的速度偏差较大，平滑性较差。

通过在前述2个方向上将光流计算的速度矢量和标准速度进行比较(见图9)，可以明显发现：图像位移量不超过2个像素部分的计算结果与标准速度结果十分吻合，偏差很小；图像位移量超过2个像素部分的计算结果则出现了明显的偏差。这很好地说明了所得到的速度场分布中为何中心区域速度分布均匀，平滑性好，而边缘区域速度偏差较大，平滑性较差，同时也在很大程度上对如何在实际测量应用中保证光流计算获得高精度的结果进行了指导。

图7 灰度原图与旋转1度后的图像

图8 采用光流计算方法得到的速度矢量图和流线图

图9 横向速度分布和纵向速度分布与标准速度比较

3 加权因子最优值估计

β2作为加权因子表征了实验和数值离散噪音的相对大小。在实际计算中β2通过调整速度梯度平方和比重来影响计算结果的平滑性，即通过选取不同的加权值获得不同平滑程度的流场速度分布。下面将以旋转灰度图像为计算基础，在不同β2值下对光流计算结果中纵向考察线速度分布与标准速度分布进行比较(见图10～14)。

由图中数据可以看出，随β2增大，纵向考察线上速度场的分布能够更好地吻合标准速度分布。在实际应用粒子图像进行光流速度场计算时，β2相对较大，说明噪声较大，这与粒子图像不同于激光诱导荧光这类图像有关，粒子图像的灰度分布存在一定的间断性，需要大的加权值来限定误差。

图10 β2=1时纵向考察线速度分布与标准速度分布比较

图11 β2=4时纵向考察线速度分布与标准速度分布比较

图12 β2=7时纵向考察线速度分布与标准速度分布比较

图13 β2=10时纵向考察线速度分布与标准速度分布比较

图14 β2=13时纵向考察线速度分布与标准速度分布比较

为了进一步分析β2的选取对速度场特性影响，根据不同β2值解算出的光流速度场与标准速度分布之间的相关系数来表现β2对解吻合的影响(见图15)。根据图中数据结果可以看出，解的准确性在β2增长初期会迅速增加，在10～18之间达到平稳阶段并取得极值，超过20以后，反而会下降，因此选取适当的β2值将对光流算法获得的速度场的准确性有很大帮助。不过这里β2是与灰度的时间和空间梯度间隔相关的常数，实际应用中应结合具体参数进行选取计算。

图15 相关系数随β2变化曲线

4 光流速度场测量试验与结果

本文针对光流测速算法对流场的要求，进行了后向台阶内流场涡结构观测实验。

4.1 光流速度场测量试验

实验的难点在于如何保证实验条件满足光流流场计算的约束方程。根据前述分析，应用如上微分法对光流流场进行解算时，要求连续拍摄的每张图像具有相同的光照条件，相邻图像对应像素点的位移不至太大，最好保证在2个像素以内，从而进行有效的光流速度场计算。

实验在实验室自行搭建的小型低速风洞系统内进行，整个系统由风洞主体、流动控制系统、示踪粒子发生和播发系统、激光照明及附属光路系统、图像采集系统和计算机处理系统组成。实验布局如图16所示。

图16 光流测速实验布局

风洞主体通过其前段可控转速的涡轮风扇电机控制流动速度，通过调节变频器可以将流经中部实验段的流动速度控制在1～20m/s的范围内，风洞主体通过其尾部的蜂窝整流器吸入混有示踪粒子的气流并进行整流，保证来流的稳定条件，并与3个示踪粒子播撒管共同工作使示踪粒子与吸入气流混合均匀。为保证恒定的光照条件，实验采用氩离子激光器输出的连续激光通过扩束镜和片光源为流场提供连续的片光照明；为保证相邻图像间位移不至太大，实验采用高速相机进行图像采集，观测流场前先对标尺进行拍摄用以标定空间尺度，采集结束后传输到计算机上存储并处理,现场设备及布置照片如图17～20所示。

图17 实验光路及实验段

图18 粒子播撒架

图19 粒子发生器

4.2 粒子图像光流场分析与图像均一化

高速相机拍摄的后向台阶内流场粒子图像如图21所示，大量连续拍摄的粒子图像记录了连续变化的光流流场信息，虚线框内为计算区域。为了对实际拍摄示踪粒子图像的光流场特性进行验证，将实际拍摄的光流场图像进行旋转获得标准位移场，如图22所示；采用与前述光流解算方法验证相同的比较方式对计算获得的速度场进行考察，如图23和24所示；并分析了速度分布与标准速度分布间的相对误差和绝对误差，如图25所示。

图20 实验现场

图21 后向台阶流场粒子图像及流场计算区域

图22 使用示踪粒子图像旋转1°前后图像及验证计算区域

图23 用光流算法对旋转前后粒子图像进行计算的速度矢量与流线图

分析可以看出，在示踪粒子图像上，当像素位移小于一定值时，计算获得的速度分布能够很好的与标准值吻合；横向考察线和纵向考察线上速度分布绝对误差基本处于5%线以下，不过由于靠近0位移位置时，速度绝对值较小，故而相对误差较大；在±1像素位置绝对误差和相对误差最小，都达到了1%左右，该位置±0.5像素邻域是速度相对精度和绝对精度都比较高的区域。

图24 横向速度分布和纵向速度分布与标准速度比较

图25 考察线上速度分布的绝对误差和相对误差

由于激光经过柱透镜制造的片光存在一定的扩张角度，在激光片光传播过程中，片光的宽度始终在变大，所以激光片光在展开区域的能量分布随时空变化也会有一定的差异。此外由于实验段的有机玻璃制造工艺以及存在一定的磨损，其中存在许多细微的气泡或表面有轻微划痕，当激光入射进有机玻璃遇到这些气泡或划痕时，激光就要被削减一部分能量，导致射入流场里的激光片光含有大量的条纹，即激光片光的能量分布是不连续的，自然所得的灰度场里也含有大量的条纹。

如果还要考虑噪声的话，可以在实验前或实验后获得一个噪声灰度场。获得噪声灰度场的一种方法是，在不加示踪物质的情况下，将摄像机对准拍摄区域，连续拍摄一段时间，将得到的图像进行平均，就可以得到噪声灰度场INoise。这时相对灰度场的表达如式(9)所示。

4.3 试验结果与分析

经过光流算法计算得到的真实流场速度分布如图26所示，结果显示最大位移为2.2像素，表明区域内速度计算结果基本位于可信区间内。

由于所观测流场的绝对速度较低(<1m/s),后向台阶内流场结构不具有典型性，所以将计算区域内的光流处理结果同采用互相关算法的PIV结果(见图27)进行比较。图28为PIV显示的流场结构，同图29的光流算法流场结构比较可以发现，在同样像素分辨率下，光流算法的结果在速度大小、方向和等速度线分布上基本与PIV计算结果基本一致，大尺度涡结构和位置也与其基本无差异，但光流计算结果中可以观测到更小尺度的涡结构，等速度线更加平滑。

图26 光流计算结果(Skip 8pixel×8pixel)

图27 PIV计算结果(8pixel×8pixel判读区)

比较两者流线图(图28和29)可以看出，在同样像素分辨率下，PIV算法的计算结果中较大尺度的涡结构可以显示出来，但是对于比较小尺度的涡结构显示却不明显，而光流算法的计算结果中小尺度涡结构被明显地显示出来。

为进一步分析，将两者小尺度涡结构区域图像放大(图30和31)，可以看出，光流计算结果在每个像素上都获得一个速度矢量，可以显示更加精细的流场结构，而PIV算法由于采用了8pixel×8pixel的判读小区，其流场结构的精细显示受到限制。

图28 PIV计算获得的流场结构

图29 光流计算获得的流场结构

图30 PIV结果中小尺度结构区域放大

图31 光流结果中小尺度结构区域放大

5 结论与展望

通过对光流算法研究与编程实现，成功搭建了可用于光流测速的实验平台，进行了小型风洞内的光流速度场测量实验，比较光流算法和PIV对拍摄的粒子图像处理得到的速度场，可以得出结论：在本实验流动条件下，光流算法可以在像素级分辨率情况下获得优于PIV获得的速度矢量场，适合用于复杂流动的速度场测量。不过基于微分法的光流解算方法在计算较大位移区间时存在缺陷，未来将进一步通过算法研究和实验手段进行改进。

[1] Lang D B. Laser doppler velocity and vorticity measure- ments in turbulent shear layers[D]. California:Caltech, Pasadena, 1985.

[2] Adrian R J. Particle-imaging techniques for experimental mechanics[J]. Annual Review of Fluid Mechanics, 1991, 23(1): 261-304.

[3] Horn B K P, Schunck B G. Determining optical flow[J]. Artificial Intelligence, 1981, 17: 185-203.

[4] Liu T, Shen L. Determination of velocity and skin friction fields from images by solving projected motion equations[C]. 22nd International Congress on Instrumenta- tion in Aerospace Simulation Facilities (ICIASF), International Congress on Instrumentation in Aerospace Simulation Facilities, Pacific Grove, CA, June 2007.

[5] Liu T, Sullivan J P. Pressure and temperature sensitive paints experimental fluid mechanics[M]. Berlin: Springer-Verlag, 2004.

[6] Tikhonov A N, Arsenin V Y. Solutions of ill-posed problems[M]. Winston, Washington D C, 1977.

[7] Heeger D J. Model for the extraction of image flow[J]. Optical Society of America, 1987, A4: 1455-1471.

[8] Fleet D J, Jespon A D. Computation of component image velocity from local phase information[J]. International Journal of Computer Vision, 1990, 5: 77- 104.

[9] Liu Tianshu, Shen Lixin. Fluid flow and optical flow[J]. J. Fluid Mech, 2008, 614: 253-291.

[10] Corpetti T, Heitz D, Arroyo G, et al. Fluid experimental flow estimation based on an optical flow scheme[J]. Exps Fluids, 2006, 40(1): 80-97.

[11] Ruhnau P, Kohlberger T, Schnorr C, et al. Variational optical flow estimationfor particle image velocimetry[J]. Exps Fluids, 2005, 38(1): 21-32.

[12] Zhan Huang, Hong Wei Wang, et al. The research of particle image velocimetry based on optical flow[C]. 16th International Symposium on Flow Visualization, Okinawa, Japan, 2004.

(编辑：杨 娟)

Research on particle image velocimetry based on optical flow

Wang Hongwei, Huang Zhan*
(China Academy of Aerospace Aerodynamics, Beijing 100074, China)

As a new aerodynamics experimental technique, optical flow test technique gradually attracts more and more attention due to its advantages in vector field measurement with pixel scale resolution and strong smoothness ability. By means of scalar constraint equation combined with smoothness constraint condition, optical flow test technique can measure global velocity vector field with high space resolution. In this paper, the integration minimization optical flow velocimetry theory and algorithm were firstly studied and programmed by C++; then, in order to verify the accuracy of the optical flow algorithm, three groups of grayscale images shifted by given displacements were used for processing by the optical flow algorithm program and the result can be used to guide how to obtain high-precision optical flow calculation result in the actual measurement applications; at last, backward step verification experiment, in which tracer particle images were acquired by high speed camera and velocity vector field was calculated by optical flow algorithm, was completed in a small wind tunnel. The result calculated by optical flow algorithm was compared with the result calculated by PIV′s cross-correlation algorithm, which shows that optical flow test technique possesses the advantages of pixel scale resolution.

optical flow;velcocity measurement;particle image;motion estimation;integration minimization

1672-9897(2015)03-0068-08

10.11729/syltlx20140115

2014-10-11;

2014-12-02

国家重点基础研究发展计划(2014CB744801)

WangHW,HuangZ.Researchonparticleimagevelocimetrybasedonopticalflow.JournalofExperimentsinFluidMechanics, 2015, 29(3): 68-75. 王宏伟, 黄 湛. 基于光流算法的粒子图像测速技术研究. 实验流体力学, 2015, 29(3): 68-75.

O355

A

王宏伟(1987-)，男，黑龙江讷河人，助理工程师。研究方向：实验流体力学。通信地址：北京市丰台区云岗西路17号(100074)。E-mail: www_whw_cn@163.com

*通信作者 E-mail: xunfang05@sohu.com