河工模型试验中的DPIV技术及其应用0 引言
河工模型是探讨河流工程问题的有效研究方法之一，但由于河工模型尺寸大，观测的范围广，使用单点式速度测量仪器费时费力；对于动床模型或非定常流动模型，床面和边界形态均在不断变化之中，因此，在河工模型试验中采用全场实时测速技术十分必要，具有很大的科研和经济价值。粒子示踪的图像全场测速技术(Particle Image Velocimetry[1，2]和Particle Tracking Velocimetry[3，4])以及数字化粒子图像测速技术(Digital PIV[5]和Digital PTV)具有方便快捷的特点，在流体力学研究中得到了广泛应用。国内科研人员根据PIV和PTV的基本原理，已成功开发出应用于大型河工模型表面流场测量的粒子示踪测速系统：清华大学研制的DPTV系统[6]、中国科学院力学研究所研制的DPIV系统等。
DPIV系统的核心是对流动图像进行处理，得到示踪粒子代表的流体的速度，这一过程称为速度提取。速度提取一般采用互相关算法和二次傅立叶变换[7，8，9]进行，需优良的硬件设备(运算速度高的计算机和浮点运算加速器件)才能满足实时测量的要求。本文考虑到河工模型中示踪粒子的分布特点以及自然光照明的条件等，采用了三种改进的速度提取算法，缩短了速度提取时间，在普通Pentium系列微机由软件实现河工模型试验流场一定分辨率下的实时测量。
1 系统的结构原理及工作参数
如图1所示，用于河工模型试验表面流场测量的DPIV系统包括一台微型计算机、图像卡、CCD摄像机、软件以及其它的监视设备。利用自然光或灯光照明，不需特定光源。位于模型上方的单台摄像机可以覆盖10m2～30m2左右的区域，系统可配置三十台工业CCD摄像机，对河工模型几百平分米至上千平方米的流动区域进行测量。
在河工模型试验中，定常流动和非定常流动对多台CCD摄像机信号的同步性要求不同。对于定常流动，则通过计算机同时选中所有图像卡上的第一个摄像机进行图像采集，采集完成之后，再切换到各个图像卡上的下一个摄像机，直到所有摄像机均完成采集，这样不同图像卡上相对应的摄像机得到的信号是基本同步的(仅相差计算机的几条指令所用的时间)，同一图像卡上各个摄像机得到的信号两两之间相差一个固定的时间(两帧图像所用时间0.08s、CCD切换时间、等待帧同步时间与软件选择的延时时间之和)，这个时间在0.08s～1.0s左右，而对于同步性要求高的非定常流动，则利用图像卡的一些特殊功能(如RGB分色输入，信号转化同步进行等)使同一图像卡上的各个摄像机信号保持完全同步。 | 
| 图1 用于河工模型的DPIV系统结构简图 | Scheme of DPIV used in physical modal | |
若镜头为理想镜头，且CCD摄像机轴线垂直于被测量平面，边缘与中心拥有相同的缩放倍率，不会产生失真；若CCD镜头轴线与被测量平面不垂直，会导致梯形失真，严重失真时需调整CCD摄像机的角度，轻度的失真则利用软件调整不同区域的缩放倍率；一般的镜头并非理想镜头，可能导致枕形失真、梯形失真等，严重时图形扭曲，增大测量系统的误差，为了克服CCD像机导致的失真，系统在使用前对每个CCD与镜头进行了搭配检验，保证水平线和竖直线在几个易产生失真的特定位置(中心和边缘)的误差均不超过一个像素；并且软件中配备了图像修正模块，利用栅格进行检验，修正由镜头造成的失真，若图像出现严重的扭曲变形，则通过更换镜头的办法予以解决。
系统采用了工业CCD作为图像摄取器件，因此性能受到CCD空间分辨率(768×576像素)和时间分辨率(两*能参数，若需要测量更快的速度，则可能通过0.02s的场间隔，将可测量的速度范围提高到20.0m/s，若镜头覆盖的面积发生变化，则“测量的速度范围”和“示踪粒子尺寸”的变化将和覆盖区域的边长的变化成正比。
表1 测量河工模型流场的DPIV系统的性能参数 |
Capability of DPIV in flow velocity measurement of physical modal |
|
摄像头数 | 6台～30台 |
|
判读区域尺寸 | (32×32)、(48×48)、(64×64)、(96×96) |
测量的速度范围 | 0.03m/s～10.0m/s |
两帧间延时时间 | 0s最高速～4.0s最低速 |
示踪粒子尺寸 | 16mm |
单头速度矢量数 | 96×96的判读区域：48～6400 |
| 64×64的判读区域：108～6400 |
| 48×48的判读区域：192～6400 |
| 32×32的判读区域：432～6400 |
单矢量判读时间 | 软件实现＜0.01s(96×96像素) |
|
2 DPIV速度提取算法
DPIV技术是利用CCD像机或扫描仪等直接得到记录粒子运动的图像序列，通过用相继两帧数字图像中相对应的小区域进行速度提取，得到各个小区域上粒子平均位移的大小和方向，也就是小区域中心处的运动速度。用img1(i,j)和img2(i,j)分别表示t0和t0+Δt时刻相继两帧数字图像中相对应的小区域，若判读区域的中心为(x,y)，判读区域的宽度为d，则i和j的取值范围为：i∈(x-d/2,x+d/2］,j∈(y-d/2,y+d/2］。假设在判读区域内各粒子及粒子团的速度大小相等，方向相同，若以(x,y)点为中心的判读区域的粒子的位移为(dx,dy)，则有
img2(i,j)=img1(i-dx,j-dy) | (1) |
速度提取就是根据图像函数img1(i,j)和img2(i,j)求解dx和dy的过程。
对PIV底片的速度提取主要采用二次傅立叶交换法[7]和图像自相关法[8]，对DPIV图像进行速度提取趋向于采用规范的二次傅立叶变换法和图像互相关算法[9]。二次傅立叶变换法和图像相关算法同样适合于河工模型流场测量中的速度提取，但根据河工模型的特点，本文提出了三种算法，提高了速度提取的效率，在一定的分辨率前提下达到了对河工模型表面流场的实时测量和记录的要求。由于河工模型中图像质量较差，因此在速度提取前进行了二值化处理。
2.1 区域内图像块相关及其误差
由于河工模型试验中示踪粒子数较少，在每一个速度提取区域内(32×32像素或其它尺寸)粒子数一般只有几粒至十几粒，根据此特点，对图像互相关算法进行的简化，得到区域内图像块相关。首先在t0时刻的图像的每一个判读区域内，多次选择包含超过70%粒子的一块区域(图2.1)，并在t0+Δt时刻的图像同一判读区域内通过相关计算寻找和所选择的小区域相关性最好的大小完全相同的图像块(图2.2)，利用这两个块之间的位移确定该求解区域内平均速度的大小和方向(图2.3)。
| 1-2.粒子图像 3.块比较的结果 4.速度矢量 | 图2 区域内图像块相关求解过程 | Calculating course of image piece correlation | | 
| 图3 块相关速度提取的误差 | Calculating error of piece correlation | |
由于相继两帧数字图像相隔时间Δt，导致第一幅图像上的部分粒子离开了判读区域(如图2.1右下侧粒子)，而在第二幅图像上又有部分新的粒子进入判读区域(如图2.2左上侧粒子)，这些粒子对图像互相关算法有一定的影响，使得判读精度降低，但对区域内图像块相关判读结果影响较小，整个系统的误差接近于模拟图像向数字图像转化时粒子位置的误差。依据河工模型流场的特点，随机产生了若干对粒子图像(32×32)，其中粒子数为1～50，粒子位移1～20像素，求解了块内粒子百分数变化时的速度，给出了图3所示的错误速度占总数的百分数。可以看出，选取粒子数低于70%时，误差较大，超过90%时，随着选取粒子百分数数的增加，求解误差增大，这是因为包含的不相关粒子数开始增加。粒子百分数为100%的区域内图像块相关即为图像互相关算法。
2.2 一维平均(压缩)法及其误差
将N×N二维图像img(i,j)平均为两个一维图像Ix(i)和Iy(j)
| (2) |
由 (i)和 (i)利用相关或傅立叶分析得到x方向位移dx，由 (j)和 (j)求得y方向位移dy。 以32×32的判读区域，产生了若干对粒子图像，粒子的位移为1～20像素，对具有不同粒子数的判读区域进行了求解，结果的误差见图4，随着判读区域内粒子数目的增加，错误结果占总数的百分比逐步提高。由于未考虑判读区域内粒子进出边界带来的影响，对结果进行了一次修正(松弛算法)，判读结果的错误率有一定程度的下降，但随粒子数的增加错误率提高的趋势和未修正时一致。 | 
| 图4 一维平均法误差 | Velocity calculating error of one dimension compress method | |
2.3 傅立叶变换相位谱分析
对图像函数img１(i,j)和img2(i,j)分别进行傅立叶变换
| (3) |
式中
，用R和θ表示F的幅度谱和连续相位谱，根据R可以判断图像的相关性。若R1=R2，则图像完全相关，而粒子的位移信息用连续相位谱之差θ2(m,n)-θ1(m,n)来表示。傅立叶变换得到的相位谱通常是缠绕相位(在-π～π之间)，需要经过相位展开计算得到连续相位谱，才能用于粒子位移的计算。
相关算法和傅立叶变换(幅度谱分析和相位谱分析)均为较精确的分析手段，而区域内图像块相关和一维平均(压缩)法是近似算法，尤其是一维平均(压缩)法目前只成功应用于河工模型的二值化图像，本文利用研制的DPIV系统在水槽内拍摄了一个漩涡运动的图像，得到了如图5所示的速度矢量分布(96×96像素的判读区域，无重叠)，图6给出了利用各种求解算法对该图像进行判读所用的全部时间。
| 图5 漩涡流动的速度场 | Velocity field of vortex | | 
| 图6 五种求解方式的比较 | Comparison of five methods | |
速度提取的精度和速率与判读区域有关系：判读精度与区域大小成反比；提取速率与2N2～N4(N为区域尺度)成反比。
河工模型试验中会出现在判读区域内没有粒子的情况，影响判读的精度。因此在图像采集时，连续记录了四幅图像，在未得到正确速度矢量的情况下，可以利用其它相关结果进行补充，而且软件可以自动对得到的速度矢量进行判断，并利用连续性方程对错误矢量进行修正。
3 DPIV技术在河工模型中的应用实例 目前，该技术已在长江科学院和武汉水利电力大学得到应用。在长江科学院，主要应用于南水北调穿黄工程和三峡船闸冲放水的模型试验研究中；在武汉水利电力大学工程泥沙实验室，已完成的项目有：长荆铁路桥模型试验、三峡下游芦家河浅滩段模型关键技术研究、阳逻电厂取水口河段模型试验研究和大别山火电厂取水保证率分析等具体研究工作。图7为武汉水利电力大学提供的三峡下游芦家河浅滩段模型中的一幅流场图。 | 
| 图7 芦家河浅滩段模型试验的流场图 | Flow field of physical model for the ford of Lujiahe | |
4 结论
利用普通CCD摄像机镜头配合图像卡以及微型计算机构成的DPIV系统，可以实现河工模型无接触快速测量速度的需要。
根据河工模型中粒子分布的特点，提出了三种DPIV速度提取的改进算法，提高了效率，达到了在一定分辨率前提下对河工模型流场的实时测量和记录的要求。
参考文献
[1]Adrian R J., Multi point optical measurements of simultaneous vectors in unsteady flow a review,Int.J.Heat & Fluid Flow,127(1986):127.[2]Adrian R J., Particle imaging techniques for experimental fluid mechanics,Annu.Rev.Fluid Mech.,23(1991):261.
[3]Oimotakis P.E., Debussy F.D., et al.Particle streak velocity field measurement in a two dimensions mixing layer.Phys.Fluids,24(1983):283.
[4]Khalighi B., Study of the intake tumble motion by flow visualization and particle tracking velocimetry,Exp.Fluids,7(1989):230.
[5]Willert C.E.,Gharib M.,Digital particle velocimetry,Exp.Fluids,10(1991):181.
[6]王兴奎等.图像处理技术在河工模型试验流场量测中的应用.泥沙研究，1996，(4)：21.
[7]许宏庆等.应用PIV技术测量二维瞬时流场.空气动力学报.11，1993，(11)：409。
[8]Meinhart C.D., Adrian R.J., et al., A parallel digital processor system for PIV,Meas.Sci.Tech., 4(1993):619.
[9]Dabiri D., Gharib M.,Digital partical image velocimetry:the method & implementation.Exp.Fluids,11(1991):77.
