1、基于序贯矢量平均的河流水面时均流场重建方法张振 1，陈红 2，刘海韵 1，石爱业 1（1. 河海大学计算机与信息学院，江苏 南京 211100；2. 河海大学水利水电学院，江苏 南京 210098）摘要：大尺度粒子图像测速（LSPIV）是一种新兴的瞬时全场流速测量技术。针对其在天然河流中由于示踪物密度低、时空分布不均及水面光学噪声干扰造成的流场估计精度过低问题，通过分析三种时均流场重建策略的敏感性，提出了一种基于序贯矢量平均的时均流场重建方法。方法利用断面流速方向一致性分布的特点，采用基于非线性统计的全局角度直方图方法检测流动主方向，较好地解决了区域性错误矢量的识别问题，改善了河流水面流场的估
2、计精度。此外，利用示踪物时空分布的冗余信息，以时均流场的矢量正确率为依据控制求平均的进程，有效提高了流场测量的时间分辨率。关键词：大尺度粒子图像测速；河流水面；时均流场；平均策略中图分类号：TH815 文献标识码：A 国家标准学科分类代码：510.40A Sequential Vector Averaging Method for reconstructing time-averaged flow field of River SurfaceZhang Zhen1, Chen Hong2, Liu Haiyun1, Shi Aiye1（ 1. College of Computer and I
3、nformation Engineering, Hohai University, Nanjing 211100, Nanjing;2. College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, Nanjing）Abstract： Large Scale Image Velocimetry (LSPIV) is a promising technique for measuring whole-field velocities instantaneously. Howev
4、er it suffers low flow field estimation precision in the application of natural rivers due to low seeding density, non-uniform tracer distribution and free-surface optical noises. By analyzing sensitivity of three time-averaged flow field reconstruction strategies, a method based on sequential vecto
5、r averaging is proposed. It uses a nonlinear statistic-based global angle histogram method to detect the main flow direction according to the consistent distribution of cross-section velocity directions. Experiment results show that it is well suited for regional error vector identification in insta
6、ntaneous flow field with low correct rate. Make use of the temporal-spatial redundancy information of flow tracers, the averaging process is controlled according to the correct rate of time-averaged flow field, which efficiently improves the temporal resolution of flow field.Key words：Large Scale Im
7、age Velocimetry; river surface; time-averaged flow field; averaging strategy1 引言天然河流的水文测验，特别是高洪期流速、流场及流量等信息的获取是洪水预报与防治的重要措施，以及河流水文学、河流动力学等研究的科学基础 1。然而，天然河流中的水体在河槽中运动时受到断面形状、坡度、糙率、水深、弯道、风、气压、潮汐等因素的影响而产生紊流。紊流内部水质点的瞬时流速在大小和方向上均随时间变化，呈现出脉动现象，但在足够长的时段内其均值保持稳定 2。因此，对于河流断面流量监测等应用而言，时均流场的重建比瞬时流场的获取更为重要。作者简介
8、：张振（1985-），男，江苏武进人，讲师，博士，主要从事光电成像与多传感器系统、大尺度粒子图像测速研究。大尺度粒子图像测速（LSPIV）作为一种新兴的瞬时全场流速测量技术，不仅可用于常规条件下明渠紊动特性和时均特性的研究，其非接触特性更具有极端条件下河道水流监测的应用潜力 3。目前，LSPIV中主要采用基于灰度相关匹配的运动矢量估计方法 4-6。它通过图像匹配估计流场中局部粒子群的运动矢量，直接得到网格化的二维流速场。具有物理意义明确、实现简单、抗干扰能力强的特点，在穷尽搜索策略下可以达到很高的匹配精度。然而，对于以天然漂浮物和水面模式为水流示踪物的河流水面成像测速而言，受水流示踪物密度低、
9、时空分布不均及水面光学噪声的影响，视场中的某些待测区域可能会暂时性地缺乏明显的示踪物，导致无法估计出这些区域的瞬时矢量或出现错误矢量 7-9。瞬时流场的矢量正确率通常仅为30%左右，远低于实验室中采用人工示踪粒子的PIV系统，使得传统基于局部流体运动连续性假设的局部校验方法难以用于错误矢量的识别，若直接采用线性求平均的方法将导致流场估计误差过大 10-12。针对河流水面成像测速的特点，本文将首先分析三种经典的时均流场重建策略对图像质量和图像数量的敏感性；然后利用河流流场中流速分布的特点和目标运动的时空冗余信息，提出一种基于序贯矢量平均的时均流场重建方法；最后通过一组实测图像序列评价方法的性能。
10、(a) 图像平均策略 (b) 相关平均策略 (c) 矢量平均策略图 1 三种时均流场重建策略Fig.1 Three kinds of time average flow field reconstruction strategy2 时均流场重建策略及敏感性分析对于基于灰度相关匹配的运动矢量估计方法，按照求平均操作在处理流程中所处的阶段及作用对象的不同，可以将后续的时均流场重建策略分为图像平均、相关平均和矢量平均三种 13，如图1所示。2.1 测试图像序列为分析三种策略的敏感性，我们采用研制的近红外智能相机采集了六安市横排头淠河总干渠的河流水面图像序列 14。如图2所示，横排头闸下游测流河段的顺
11、直长度约600m，为人工河槽，左、右岸均为混凝土护坡；河床系粘土，较稳定；河段上游建有一座五孔进水闸，闸孔出流在水面产生顺流而下的漩滚模式和白色泡沫。实验断面位于进水闸下游550m 处，河宽约74 m；相机架设于左岸距水面高程约为5m的水泥高台上，视场覆盖完整的断面，以1s为时间间隔拍摄 pixel的8bit灰度图像。12804从图3中可以看出，由于视场尺度较大，水面细小的白色泡沫几乎已不可见，有效的示踪物以漩滚模式为主，还夹杂着少数漂浮物。图像中部和远场分别受水面耀光和植被倒影的影响，形成片状亮区和条带区的干扰模式。这里选取序列中的30帧，并将其中 pixel的单个分析区域（IA）图像对作为
12、测试子256图。图像增强采用基于视觉感受野的自适应背景抑制方法 9，运动矢量估计采用基于快速哈特利变换互相关（FHT-CC）的方法 6。29组图像对的相关曲面信噪比（SNR）和瞬时位移表明：序列中的漩滚模式由于存在自身形变运动，在较长的时间间隔下容易引起图像对中目标模式的杂乱，导致相关曲面的SNR较低，仅分布在46之间，并且瞬时位移表现出强烈的脉动现象，在方向上存在明显偏差的错误矢量为9个，约占矢量总数的30%。图 2 测流河段的遥感视图及现场视图 图 3 河流水面图像及测试子图Fig.2 Remote sensing and field view of the measured sectio
13、n Fig3 River-surface image and the test subgraph2.2 图像平均策略图像平均操作处于运动矢量估计中的图像预处理阶段，是一种像素级的时均流场重建策略。它首先将幅图像 分为奇数组 和偶数组 ，分别对每组图像中对应像素的灰度值求平均得到1对时均图像，然Nfofef后计算该图像对的相关测度得到1个时均相关曲面，最后通过峰值检测得到时均矢量。(a) 10 帧 (b) 20 帧 (c) 30 帧(a) 10 frames (b) 20 frames (c) 30 frames图 4 图像平均策略下不同帧数得到的时均图像Fig.4 Time-averaged
14、images for different frame numbers with image-averaging strategy(a) 10 帧 (b) 20 帧 (c) 30 帧(a) 10 frames (b) 20 frames (c) 30 frames图 5 图像平均策略下不同帧数得到的时均相关曲面Fig.5 Time-averaged correlation surfaces for different frame numbers with image-averaging strategy图像平均策略的相关测度可以表示为：(1),xyoexycdfxyfd其中， 和 。将上式展开，
15、得到：1312/oNfff 24/e Nfff(2)1416132335551214161, NxyNNNfffcd fffff 可以看出，式中仅 个对角项表示有效图像对间的连续相关，是对位移估计有用的目标项，以 为级数/2N N增加；而其他 个非对角项表示非连续图像间的随机相关，构成了相关曲面的噪声项，以 为14 2级数增加。因此，噪声项的增加速度高于目标项将使得图像平均策略中相关曲面SNR随帧数增加而降低。从图4给出的时均图像中可以看出：由于水流运动，序列中稠密分布的水面模式在时均图像中被拓展为片状纹理，在一定程度上提高了IA中示踪物的密度。但由于目标相对于背景的面积和灰度值均较小，图像的
16、SNR 较低。随着平均帧数的增加，目标强度的衰减加剧。在图5的时均相关曲面中，当 时，噪声10N峰已经超过目标峰在相关曲面中占主导，导致匹配失效。2.3 相关平均策略相关平均操作处于运动矢量估计中的图像特征提取阶段，是一种特征级的时均流场重建策略。它首先将 幅图像依次划分成 组图像对，然后分别计算每组图像对的相关测度，并对 个瞬时相关曲面N1N 1求平均得到1个时均相关曲面，最后通过峰值检测得到时均矢量。相关平均策略的相关测度可以表示为：(3)1,xyiixycdfxyfd其中， 。根据线性运算的交换律，可将求平均算子放在积分运算内，即：1,2iN(4) 1 12341,xyiixy Ncdf
17、fdfffN 可以看出，式中仅包括式(2)的对角项，即对位移估计有用的目标项，以 为级数增加。因此，相关平均策略可以得到更高信噪比的时均相关曲面。(a) 10 帧 (b) 20 帧 (c) 30 帧图 6 相关平均策略下不同帧数得到的时均相关曲面Fig.6 Time-averaged correlation surfaces for different frame numbers with correlation-averaging strategy从图6给出的时均相关曲面中可以看出：相关噪声的方差随着帧数的增加明显降低，而目标峰更加尖锐，说明当瞬时相关曲面的SNR较高时，相关平均可以有效滤除
18、随机噪声。但零位移附近的噪声峰在相关平均过程中也被增强了，说明这种由图像背景泄漏引起的直流偏置噪声并不能通过求平均来消除，当其累加速度超过目标峰时（ ）匹配就失效。30N2.4 矢量平均策略矢量平均操作处于运动矢量估计中的位移矢量输出阶段，是一种决策级的时均流场重建策略。它首先将 幅图像依次划分成 组图像对，然后分别计算每组图像对的相关测度得到 个瞬时相关曲面，N1 1N最后对峰值检测后输出的 个瞬时矢量求平均得到时均矢量。矢量平均策略的原理最为直观，可用瞬时矢量的算术平均表示为：(5)11,NNi iiiSxySxyxy从表1给出的结果可以看出：矢量平均得到的时均矢量未出现明显的方向错误，说
19、明矢量平均对方向的容错性较好；但相比图像平均和相位平均策略，在大小上存在不同程度的估值过低现象。这是由于错误矢量被带入平均操作造成的，并且错误比例越高、平均时间越长，退化作用越严重。因此，矢量平均策略要求所有的瞬时矢量必须都可靠才能获得准确的时均矢量。表2给出了剔除错误矢量后的相关平均和矢量平均结果。其中位移估计的可靠性和稳定性均得到了明显改善。综上分析，图像平均策略的优点是仅需计算一次相关测度，运算速度快；缺点是需要一次读入完整的图像序列，占用大量存储空间；并且受到平均帧数的限制，对图像质量有较高的要求。相关平均策略的优点是能够有效提高相关曲面的信噪比，降低出现错误矢量的概率；缺点是对相关曲
20、面的直流偏置噪声敏感，依赖于有效的图像背景抑制方法。矢量平均策略的优点是灵活、高效，可以处理任意方法获得的网格化矢量场，由于是时均流场重建的最高层次，可以获得瞬时流场以及特征级的中间参数，具有较高的容错性；缺点是当示踪物密度较低时，对相关曲面信噪比低下引起的错误矢量敏感。表 1 三种平均策略下不同帧数得到的时均相关曲面信噪比和位移Table 1 SNRs and displacements of time-averaged correlation surfaces for different frame numbers with three averaging-strategies策略 图像平
21、均 相关平均 矢量平均帧数 SNR x y SNR x y SNR x y10 3.5379 -0.1458 0.4778 4.2142 -54.3038 4.1978 4.7153 -38.5417 4.520020 3.7369 -0.0476 0.1104 4.3347 -39.6791 4.2494 4.8924 -39.2058 3.601230 5.0387 -0.0155 0.0966 4.2819 -0.1042 0.3135 4.8302 -34.7993 3.2022表 2 剔除错误矢量后相关平均和矢量平均策略的对比Table 2 Comparison of correla
22、tion-averaging and vector-averaging strategies after eliminating the error vectors策略 相关平均 矢量平均帧数 SNR x y SNR x y10 4.5442 -46.0585 4.3114 4.8716 -49.5478 4.139920 4.6432 -46.4824 4.3808 4.8016 -50.6635 4.179930 4.4217 -40.3267 4.3458 4.7414 -48.2347 4.29983 序贯矢量平均方法这里以矢量平均策略为基础，提出一种基于序贯矢量平均的时均流场重建方法
23、。具体实现流程如图7所示，主要包含以下三个关键步骤。3.1 错误矢量识别在采用灰度相关匹配法得到的瞬时位移矢量场中不可避免地会产生流动特征异于相邻矢量的错误矢量。目前在PIV技术中应用最为广泛的错误矢量识别方法主要有两类 15, 16：一类是基于局部流体运动连续性假设的局部校验方法。方法假设流体在流动过程中其运动状态满足连续性，只要被检测矢量与邻域内其他矢量在大小或方向上的均值或中值之差在容许的误差范围内，则认为当前矢量为正确矢量。另一类是基于全局流体运动一致性假设的全局校验方法。方法利用瞬时流场中矢量大小或方向的全局特征作为错误矢量的识别依据，一般适用于流动一致性较好的均匀流场。图 7 序贯
24、矢量平均方法的流程图Fig.7 Flow chart of the sequential vector-averaging method在天然河流中流速脉动往往导致位移大小随时间的变幅可达30%以上，而受边壁效应的影响，断面流速的空间分布呈现出中泓向两岸递减的趋势。这种时空分布的变化性使得基于矢量大小的全局校验方法失效。相比之下，明渠水流多为方向一致性较好的单向流，在实验中发现，80%以上大小错误的矢量，其方向也具有明显差异。鉴于此，本文采用矢量方向作为错误矢量识别的依据，设计了一种基于非线性统计的全局角度直方图算子求解瞬时流场的流动主方向。具体实现方式是：对于 时刻的瞬时位移场，首先计算i各
25、矢量 与 方向的夹角：,iSxy(6),arctn/i iixyyx得到图8所示的矢量极坐标分布；然后以 为角度区间统计矢量方向在0-359范围内的全局直方图，如图9所示；接下来依次搜索各角度区间，将对应矢量数量最多的区间标记为流动主方向 （红色区域）；M最后以角度阈值 为判别窗口（绿色区域）识别每个矢量的类型，并建立矢量类型标志如下：TH(7)0,1iMTHi xyFxy其中，1表示正确矢量，0表示错误矢量。这种方法的好处是无需流动主方向的先验知识。3.2 流场时间滤波对时均流场重建而言，平均时间是一个重要的敏感因素：时间过长会使结果失去代表性，降低时间分辨率；时间过短会降低测量的精确度，不
26、足以消除随机误差。因此，理想的重建策略应当在保证时均流场精确度和可靠性的前提下，尽可能地提高时间分辨率。MM图 8 矢量极坐标分布 图 9 矢量角度全局直方图Fig.8 Vector distribution under polar coordinates Fig.9 Global histogram of vector angles为达到该目的，本文采用了一种序贯式的矢量平均策略，或称为时间轴滤波。基本思想是：舍弃当前错误的瞬时矢量，并用当前正确的瞬时矢量替代此前的错误矢量，用当前正确的瞬时矢量平均此前正确的时均矢量。滤波过程可用下式描述为：(8)11 1NUL, ,0,1, , , 1,i
27、 iii i iii ii ii FxyxySxySxy Sxyxyxy 其中， 、 分别为当前和此前时刻的时均矢量； 为当前时刻的瞬时矢量；错误矢量,iSxy1,ixy ,i标记为“NULL”； 为参与当前矢量平均的正确矢量个数； 、 分别为此前时均矢N 1iFxy,i量的类型标志和当前瞬时矢量的类型标志，对二者进行“或”运算得到当前时均矢量的类型标志：(9)1,|,iiiFxyxy遍历每个IA后统计得到正确和错误矢量的个数 和 。则时均流场的矢量正确率为：N0(10)1/iP若 小于一个阈值 （如90% ），则继续处理下一时刻的图像对，直到 大于 并且达到预设的平均帧iPTHP iPTH数
28、。如果遍历所有帧后仍未达到指定的正确率，则认为时均流场重建失效。由于这种序贯式的矢量平均策略以时均流场的矢量正确率为依据控制求平均的进程，可以最大限度地保证时间分辨率。3.3 流场空间滤波流场时间滤波后得到一个包含错误矢量的时均流场和对应的矢量类型查找表。为保证流场的完整性，采用错误矢量邻域内正确矢量的均值对空缺矢量进行插值：(11)11, ,FnmSxyFxnySxnymN其中， 为33邻域内正确矢量的个数，邻域大小可根据流场的空间分辨率选取。对于时均流场正,TNxy确率较高的情况，可以增加一个空间平滑滤波过程，使有限采样时间内的离散化流场变得更平滑，更接近真实流场的连续性。时均流场和滤波器
29、的空域卷积表示为：(12),xyxySG其中滤波函数采用归一化的二维离散高斯窗，其33形式可表示为：(13)1246G由于矢量插值仅利用邻域内的正确矢量，可以避免错误矢量带来的误差；而平滑滤波是对流场的全局处理，平滑的程度由窗口大小决定，可能对正确矢量产生退化作用，因此对于平均帧数较少或时均流场正确率小于预设值的情况，应谨慎使用。4 方法的性能评价这里采用图3序列的完整断面图像测试时均流场重建方法的有效性。其中，有效水面区域的大小为pixel；运动矢量估计采用固定窗口的FHT-CC 方法，窗口尺度设为 pixel， 和 方向的128075 128xy平移步进分别设为64和32pixel ，共有
30、380个IA ；用于建立全局角度直方图的角度区间设为10，角度阈值设为15，时均流场正确率的阈值设为100%，以便观察序列完整的重建过程。图10给出了瞬时流场及时均流场的矢量正确率随时均流场重建过程的变化曲线。可以看出，29幅瞬时流场的矢量正确率分布在1434%之间，远低于实验室中采用人工粒子的示踪条件。在第1幅瞬时流场图（图11a）中，错误矢量（黑色箭头）主要分布在图像近场的低流速区域和远场的水面倒影区域。随着重建过程的进行，水面图像中原本缺乏足够目标信息的区域逐渐获得了正确的位移矢量估计结果。在时间轴滤波的作用下，时均流场中的错误矢量逐渐被正确矢量（白色箭头）替代，时均流场的正确率从27%
31、逐渐提升至第29幅的80%（图11b）。在重建过程的初始阶段，矢量正确率随平均帧数的增多迅速提升，提升速率和示踪物的运动速度有关。但受强噪声的影响，某些区域依然无法获得正确矢量，导致正确率曲线在75%附近趋于平缓。采用局部邻域均值对错误矢量进行插值（图11c），得到的流场中部分和大小错误矢量相关的粗大误差并没有很快被平均掉。对插值后的流场进行一次空间平滑滤波后（图11d），流场的连续性得到了明显改善，符合断面流速分布规律。图 10 瞬时流场及时均流场的矢量正确率Fig.10 Correct ratios of vectors in instantaneous and time-averaged
32、 flow fields(a) 第 1 幅 (b) 第 29 幅(a) 1st frame (b) 29th frame(c) 矢量插值 (d) 空间平滑(c) vector interpolated (d) spatial smoothed图 11 时均流场重建过程Fig. 11 Reconstruction process of time-averaged flow field5 总结针对天然河流中水流示踪物密度低、时空分布不均及水面光学噪声影响导致时均流场估计精度过低的问题，提出了一种序贯矢量平均的时均流场重建方法。方法利用断面流速方向一致性分布的特点，采用基于非线性统计的全局角度直方图
33、方法检测流动主方向，较好地解决了瞬时流场矢量正确率较低情况下区域性错误矢量的识别问题；利用运动目标时空分布的冗余信息，以时均流场的矢量正确率为依据控制求平均的进程，仅需顺序采集20帧左右的图像就可以达到80%以上的矢量正确率，并且单帧执行时间在0.5s以内，改善了流场测量的时间分辨率。相比图像平均和相关平均策略，本方法在实时性和硬件开销上具有明显优势。需要注意的是，全局角度校验法能够适应瞬时流场正确率在30%左右的情况，但当示踪物密度极低、流场正确率小于10%时，估计流动主方向的可靠性将大幅降低。因此有必要在错误矢量识别中进一步考虑局部流速矢量的大小和主流方向，以适应流态较复杂的情况。参考文献
34、1 徐立中, 张振, 严锡君, 王慧斌, 王鑫. 非接触式明渠水流监测技术的发展现状J. 水利信息化 ,2013, (3): 37-44. (XULizhong, ZHANG Zhen, YAN Xijun, WANG Huibin, WANG Xin. Advances of Non-contact Instruments and Techniques for Open-channel Flow Measurements J. Water Resources Informatization, 2013, (3): 37-44.(in Chinese)2 赵志贡. 水文测验学M. 郑州: 黄河
35、水利出版社, 2005. (ZHAO Zhigong. Hydrometry M. Zhengzhou: the Yellow River water conservancy press, 2005. (in Chinese)3 张振，徐枫，王鑫，徐立中 . 河流水面成像测速研究进展J. 仪器仪表学报, 2015, 36(7): 1441-1450. (ZHANG Zhen, XU Feng，WANG Xin，XU Lizhong. ResearchProgress on River Surface Imaging Velocimetry J. Chinese Journal of Scien
36、tific Instrument, 2015, 36(7): 1441-1450.(in Chinese)4 I. Fujita, M. Muste, A. Kruger. Large-scale particle image velocimetry for flow analysis in hydraulic engineering applicationsJ. Journal of hydraulic Research, 1998, 36(3): 397-414.5 J. Le Coz, A. Hauet, G. Pierrefeu, G. Dramais, B. Camenen. Per
37、formance of image-based velocimetry (LSPIV) applied to flash-flood discharge measurements in Mediterranean riversJ. Journal of Hydrology, 2010, 394(1): 42-52.6 严锡君, 张振, 陈哲, 吕莉, 王慧斌. 基于FHT-CC 的流场图像自适应运动矢量估计方法J. 仪器仪表学报, 2014, 35(01): 50-58. (YAN Xijun, ZHANG Zhen, CHEN Zhe, LV Li, WANG Huibin. FHT-CC-
38、based adaptive motion vector estimation method for flow field image J. Chinese Journal of Scientific Instrument, 2014, 35(01): 50-58. (in Chinese)7 张振 , 徐立中, 王慧斌. 河流水面成像测速中的水流示踪物综述J. 水利水电科技进展. 2014, 34(03): 81-88. (ZHANG Zhen, XU Lizhong, WANG Huibin. Review of natural flow tracers for river surface
39、 imaging velocimetry J. Advances in Science and Technology of Water Resources, 2014, 34(03): 81-88. (in Chinese)8 Zhang Zhen, Wang Xin, Fan Tanghuai, et al. River surface target enhancement and background suppression for unseeded LSPIV J. Flow Measurement and Instrumentation, 2013, 30: 99-111.9 张振,
40、陈哲, 吕莉, 王鑫, 徐立中. 基于视觉感受野的自适应背景抑制方法J. 仪器仪表学报, 2014, 35(01): 191-199. (ZHANG Zhen, CHEN Zhe, LV Li, WANG Xin, XU Lizhong. Adaptive background suppression method based on visual receptive field J. Chinese Journal of Scientific Instrument, 2014, 35(01): 191-199. (in Chinese)10 I. Fujita, S. Aya. Refinem
41、ent of LSPIV technique for monitoring river surface flowsC, ASCE 2000 Joint Conference. on Water Resources Engineering and Water Resources Planning & Management, Minneapolis: 2000.11 M. Jodeau, A. Hauet, A. Paquier, J. Le Coz, G. Dramais. Application and evaluation of LS-PIV technique for the monito
42、ring of river surface velocities in high flow conditionsJ. Flow Measurement and Instrumentation, 2008, 19(2): 117-127.12 A. J. Bechle, C. H. Wu, W. C. Liu, N. Kimura. Development and application of an automated river-estuary discharge imaging systemJ. Journal of Hydraulic Engineering, 2012, 138(4): 327-339.
