帧差与快速密集光流结合的图像法测流研究

王剑平，朱 芮，张 果，何兴波，蔡如鹏

(1.昆明理工大学 信息工程与自动化学院，云南 昆明 650000；2.贵州省黔西南州水文水资源局，贵州 兴义 562400)

中国河流众多，水资源丰富，水资源的合理以及有效利用对我国的经济发展具有重要的意义。天然河流的流量监测是管理地表水资源的重要要素之一，高洪时期的流速和流量信息的获取也是山洪地质灾害防治非工程措施的重要组成部分。20世纪80年代以前，绝大部分国内外的水文监测站都使用基于转子式流速仪的流速–面积法来测量各个河道断面布置点的平均流速，并进行流量计算。但由于自然河流复杂多变的流动特性以及严峻的户外测量环境，点流速测量的代价高昂，导致流速仪等传统的接触式测量方法不能正常开展或者仪器无法正常施测。因此，发展非接触式测流方法是必要且迫切的。

近几年来，机器视觉及传感器技术的迅速发展推动非接触式测流技术取得巨大的突破，基于光学、图像以及雷达的各种测流技术都已经得到了快速地发展和应用，尤其是图像法测流技术，该技术因其简便、高效、安全等优点得到了普遍的关注。粒子图像测速（particle image velocimetry，PIV）技术发展于20世纪80年代，该技术在流场中布撒示踪粒子，以激光片光作为光源，用底片或CCD（charge coupled device，CCD）相机记录示踪粒子的图像，利用自相关法、互相关法或光学杨氏条纹法处理这些图像，从而获得流场的速度分布。Fujita等提出了大尺度粒子图像测速（large-scale particle image velocimetry，LSPIV）方法，改进PIV技术用于大尺度水面流场观测。该方法以波纹、植物碎片、泡沫等天然水面漂浮物和天然水面模式代替粒子作为示踪物，用普通的数字相机或视频录像机记录自然光光源下的河流图像，简化了硬件配置，还避免了示踪物对河流的污染。该方法自提出以来，被广泛地应用于水文站流量、流速的监测。但由于采用极易受到环境影响的天然水面漂浮物和诸如泡漩等不稳定的天然水面模式作为示踪物，该方法存在稳定性差的问题。且LSPIV方法在运动估计时采用空域相关匹配法，相关性计算的时间复杂度达O(

n

)，所需存储空间大且计算效率低。相比于LSPIV，Fujita等在2007年提出了基于时空图像的测流（spatio-temporal image velocimetry，STIV）方法。该方法在断面流量测量方面展现出了更好的性能，目前正在被广泛地研究，但该方法只能获得一维流速场，且对河流流态稳定性要求高。王万良等将深度学习用于流速测量，采用条件边界平衡生成对抗网络对流速图像进行生成，采用多特征融合的卷积分类网络对流速图像进行分类，此方法仅能获得流速区间，但为非接触式测流技术提供了一种全新的思路。Lin等结合无人机利用热水作为示踪剂，利用Lucas-Kanade光流法取得了优于常用跟踪方法的测流结果。但由于将热水作为大尺度天然河流的示踪剂并不具有普遍的可行性，该方法也存在局限性。Tauro等利用金字塔Lucas-Kanade光流法跟踪了FAST、Shi-Tomasi、SIFT、SURF、ORB5种特征点，并利用先验信息对特征点的跟踪轨迹进行滤波，只保留视场中与实际物体运动相关的轨迹，此方法在实验过程中取得了良好的结果，但特征点的选取是一个复杂过程，考虑天然河道实际工况的复杂性，在不同的现场环境下选取光流跟踪特征点的统一规则并没有制定，且先验信息的获取并不具有普遍性，故该方法在实际场景下稳定性较差。

由于天然环境下水流的复杂特性，主流的图像法测流技术依然存在着不稳定和耗时严重的问题。为此，提出一种结合帧间差分与快速密集光流（fast optical flow using dense inverse search grouping，DIS）的分组测流方法（frame difference-fast optical flow using dense inverse search-grouping，FD-DIS-G），以实现河流平均流速和断面流量的准确、实时、可持续监测。

1 FD-DIS-G测流算法

河流的流动反映在图像上为图像灰度的表面运动，本方法利用密集光流法跟踪像素点的运动，不依赖于天然水面模式作为示踪物，有利于增强方法的稳定性；且只选取测速点邻域的块区域进行快速密集光流计算，保证了算法的实时性；在光流计算之前进行了运动显著性图计算，通过捕捉水面不显著的运动，进一步保障了算法的准确性。

基于FD-DIS-G的分组测流方法主要包括以下步骤：计算运动显著性图、计算帧间光流位移、分组处理奇异值以及计算表面流速。算法流程图如图1所示。

图1 算法流程图Fig. 1 Algorithm flowchart

1.1 运动显著性图

对车辆或行人等运动目标进行运动检测时，目标和背景之间通常很容易就被区分开来，但在天然河流的视频图像中水面是作为整体一起运动的。图2为河流表面图和运动显著性图。

图2 河流表面图与运动显著性图Fig. 2 River surface map and motion saliency map

1.2 基于密集逆向搜索的快速光流估计

由第1.1节得到运动显著性图，根据显著性图中断面线的起点加上起点距

b

可得测速点

P

的坐标，并以

P

为中心点确定当前帧（第

t

帧）的块区域

G

，下一帧的块区域

G

，块区域大小为(δ×δ)， δ为块区域的边长，由相邻两个测速点之间的像素距离决定。于是对两帧块区域之间位移 Δ

d

的求解可转化为对两块区域之间的密集光流位移

A

进行计算，计算过程如下：Step1：以 ρ为尺度金字塔下采样系数，分别对

G

和

G

构建尺度金字塔，如图3所示。定义模板

T

为金字塔在第

s

层块区域

G

里的子区域，模板大小为(θ×θ)。用均匀的网格块初始化

G

且每个网格块的大小和模板

T

相同，如图4所示，并用式（2）初始化第

k

个网格块

I

中心点的光流位移(

u

,

v

)。

图3 块区域金字塔Fig. 3 Block area pyramid

图4 块区域网格化Fig. 4 Block area meshing

式中：(

x

,

y

)为网格块

I

中心点的坐标；为

A

的元素，包含了

x

、

y

方向的位移；

A

为第

s

+1层的光流位移。

Step4：定义能量函数（8）对当前尺度产生的密集光流进行优化：

1.3 分组处理奇异数据

视频拍摄过程中常常会出现丢帧、无效帧或光照突变的情况，这会导致算法在进行帧间计算时产生奇异数据，从而影响整体计算结果，降低算法的准确率。因此，奇异数据的处理是必要的。本研究方法设计了一种分组处理奇异数据的方法对运动显著性图之间得到的光流位移 Δ

d

进行处理，使得算法更具有鲁棒性。视频分组图如图5所示。

图5 视频分组处理图Fig. 5 Video grouping processing diagram

式中，

s

为

M

个数据的标准差。

Z

分数以标准差为单位衡量某一原始数据偏离平均值的距离，

Z

分数的值就表示了原始数据与平均值之间相差的标准差个数。本算法定义，与最佳估计值的偏差超过2倍标准差的测量值称为奇异值，因此当|

Z

|＞2时，该数据被认为是奇异值，应该舍掉该数据。

1.4 地面标定与透视投影变换

采用透视投影变换的方法将流场中的运动矢量从图像坐标系转换到世界坐标系，从而获得两帧之间的真实位移，并除以时间间隔得到以“m·s”为单位的河流表面流速，再由流速–面积法求得平均流速和断面流量。

综上，可以通过上述变换将图像中的帧间位移转化为实际的帧间位移，从而获得河流表面流速。

1.5 表面流速估测与流量计算

由第1.4节所述的变换可以得到实际的帧间位移，并除以时间间隔得到以“m·s”为单位的表面流速。设得到测速点

i

的表面流速为

v

（

i

=1,2,···,

n

；

n

为断面线上测速点的个数），测速点

i

-1到

i

对应的过水断面的面积为

s

，测速点

i

-1到

i

之间的平均表面流速为：

最后通过河流的流量和面积可以得到整条河流的平均流速

v

¯：

式中，

S

为断面总面积。

2 试验结果与分析

2.1 试验设计

在水文测验中，流速仪法得到的测量结果普遍被认为是真值。因此本算法将流速仪法得到的结果作为比测标准，以垂线平均流速、平均流速和流量作为衡量指标，对某水文站经过断面修整的河道进行算法比测实验，并选取了贵州省黔西南州保基水文站的天然河道进行算法普遍适用性验证实验。

2.2 比测实验

比测实验在某水文站经过断面修整的河道进行了比测试验，选取LSPIV方法、STIV方法、DIS方法（未进行运动显著性处理）以及未进行分组奇异值处理的方法（frame difference-fast optical flow using dense inverse search，FD-DIS）测得的流速和流量做为对比对象，并分析了5种算法的运行时间。实验设备处理器的型号为i5-7300CPU，操作系统为Windows10(x64)，程序运行环境为Python3.7和OpenCV4.2。根据该水文站多年实测经验，得岸边系数为0.80，水面流速系数为0.82。实验地面标定点布置如图6所示。左岸布置

A

、

B

两个点，右岸布置

D

、

C

两个点，起点与

D

点重合，终点与

A

点重合，以

DA

为断面线，断面数据如表1所示。

图6 比测实验的地面标定点布设示意图Fig. 6 Layout diagram of ground mark in comparison test

表1 比测实验的断面数据
Tab. 1 Cross section data of the comparison test

注：断面总面积为9.70 m。

起点距b/m深度/m 0.10.59 0.50.70 1.00.81 1.50.82 2.00.76 2.50.86 3.00.93 3.50.90 4.00.96 4.51.00 5.00.99 5.50.94 6.01.04 6.50.95 7.00.86 7.50.78 8.00.77 8.50.70 9.00.92 9.50.85 10.00.81 10.50.74 11.00.55 11.50.42 11.90.39

为了控制变量单一，于同一时段进行流速仪测量与视频拍摄以保证两者是测的同一段时段的水流。测量期间水流稳定、天气晴朗、测量流速仪型号为LJ-20。利用流速仪法在起点距分别为2、3、4、5、6、7、8、9、10、11 m的垂线上测得垂线平均流速，结果如表2所示，并计算得到了河流平均流速为0.47 m/s，断面流量为4.56 m/s。根据流速仪的位置，在相同起点距位置确定测速点，如图6中的小圆点所示，利用本文所介绍的方法测量测速点的表面流速，并按照河流流量测验规范计算平均流速和流量，将本文方法(FD-DIS-G)的测量结果与流速仪法、LSPIV方法、STIV方法、DIS方法以及FD-DIS方法的测量结果进行比较。

表2 比测实验的流速仪测量结果
Tab. 2 Current meter measurement results of the experiment were compared

起点距b/m 垂线平均流速/(m·s-1)部分平均流速/(m·s-1)部分面积/m2部分流量/(m3·s-1)0（岸边）0——0～2—0.291.450.42 2 0.36——2～3—0.400.850.34 3 0.43——3～4—0.390.920.36 4 0.35——4～5—0.400.990.40 5 0.45——5～6—0.600.980.59 6 0.76——6～7—0.720.950.68 7 0.68——7～8—0.640.800.51 8 0.61——8～9—0.540.770.42 9 0.48——9～10—0.430.860.37 100.38——10～11—0.440.710.31 110.49——11～11.9—0.390.420.16 11.9（岸边）0——

2.2.1 运动显著性图

实验视频拍摄时长为10 s，共切分成325帧，捕获视频的原始图像分辨率为1 920×1 080 像素/英寸，为捕捉细微的水面运动，先对325帧图片进行运动显著性计算。在利用视频图像进行表面流速测量时，大部分方法对图像的预处理是对图像进行灰度化，但在低流速下的河流中如果只对图像进行灰度化处理会难以突出很多微小的水面运动。本算法在进行运动矢量估计前先进行运动显著性计算，分别对帧差图中大于和小于设定阈值的像素值进行了缩放，认为大于设定阈值的像素为细微的运动区域对其进行放大，从而达到了运动显著性计算的效果，突出了微小的水面运动。图7为获取运动显著性图过程。从图7（d）中看出进行运动显著性计算后水面运动更为清晰，且捕捉到了很多细微的水面运动，如7（d）中方框框出的部分。

图7 获取运动显著性图过程示意Fig. 7 Obtaining motion saliency map

2.2.2 密集光流场

由于LSPIV方法使用空域或频域相关匹配来获得运动矢量时相关性的计算耗时长，本算法提出利用DIS方法进行运动矢量的估计，在相邻两帧运动显著性图中选取块区域计算密集光流位移。由于没有在全局区域进行密集光流的计算，且在计算密集光流时运用了逆向搜索的法则，使得本方法不仅可以得到亚像素精度的运动矢量，并在运算速度上也得到了保证。将视频图像中得到的密集光流可视化，如图8所示。

图8 密集光流可视化Fig. 8 Dense optical flow visualization

2.2.3 奇异数据处理

设计了一种分组处理奇异值的方法，对相邻4帧中较大、较小的光流估计进行剔除，并保留中间值，并对所有组保留的光流估计值按照其是否满足独立同分布分别进行奇异值处理，使得由于不明显的光照变化或视频缺帧丢帧导致的奇异光流估计被舍弃。如图9所示，奇异值（图9中红色的点）剔除后每个测量点的帧间像素位移分布更集中，从而降低了由于视频拍摄过程中操作不当而产生的误差，故本方法测得的值更稳定，更接近于流速仪真值，因此提高了本方法的精度。

图9 奇异数据处理Fig. 9 Abnormal data processing

2.2.4 结果对比

FD-DIS-G方法、LSPIV方法、STIV方法、DIS方法以及FD-DIS方法的测流结果如表3所示。将FDDIS-G方法得到的垂线平均流速、平均流速和流量分别与流速仪法测量结果进行对比得到结果如图10（a）、（b）、（c）所示。将表3中5种方法的测量结果与流速仪真值相比可计算出5种方法的绝对误差（表4）和相对误差（表5）。5种方法的耗时对比结果如图10（d）所示。

表3 5种方法测量结果对比
Tab. 3 Comparison of measurement results of five methods

垂线平均流速/(m·s–1)平均流速/(m·s–1) 流量/(m3·s–1)b=2 m b=3 m b=4 m b=5 m b=6 m b=7 m b=8 m b=9 m b=10 m b=11 m流速仪0.360.430.350.450.760.680.610.480.380.490.474.56 LSPIV0.230.010.450.580.380.410.570.620.370.620.393.77 STIV0.380.460.030.450.130.790.720.720.610.400.434.21 DIS0.280.270.470.610.480.490.590.690.550.340.454.32 FD-DIS0.380.440.480.540.480.530.580.780.720.770.525.08 FD-DIS-G0.330.410.450.520.460.510.540.750.680.640.494.75方法

表5 5种方法的相对误差对比
Tab. 5 Relative error comparison of five methods

方法垂线平均流速相对误差/%平均流速相对误差/%流量相对误差/%b=2 mb=3 mb=4 mb=5 mb=6 mb=7 mb=8 mb=9 mb=10 mb=11 m LSPIV34.9397.7028.6729.2150.2140.157.0128.283.6126.8617.0017.30 STIV4.606.7091.000.8082.8016.8018.049.6061.5017.7012.807.70 DIS22.2037.2034.3035.6036.8027.903.3043.8044.7030.604.305.30 FD-DIS5.602.3037.1020.0036.8022.104.9062.5089.5057.1010.6011.40 FD-DIS-G8.304.7028.6015.6039.5025.0011.5056.3079.0030.604.304.20

图10 5种方法测量结果Fig. 10 Five methods of measurement results

分析垂线平均流速可得（图10（a））：在靠近起点的7个点中，本文提出的FD-DIS-G算法和未经过分组奇异值处理的FD-DIS方法所得结果是最稳定且接近于流速仪真值的，相比与LSPIV方法和STIV方法具有良好的一致性；但在靠近终点的3个点中，奇异值处理就体现出了其作用，尤其是最靠近终点的测点提高了0.13 m/s的精度。从表4、5可以看出，该条河流流速较小，在起点距5 m处，即使是小于0.1 m/s的绝对误差也会表现为15%左右的相对误差。而在靠近终点的3个点中，5种方法的相对误差都偏大，这是由于靠近终点的地方距离镜头较远，镜头畸变和靠近岸边表面流态较混乱造成。总体来看，由于捕捉到了细微的水面运动并进行了奇异值处理，本文算法稳定性得到了提高，使得测量整体精度也得到了提高。

表4 5种方法的绝对误差对比
Tab. 4 Absolute error comparison of five methods

方法垂线平均流速绝对误差/(m·s–1)平均流速绝对误差/(m·s–1)流量绝对误差/(m3·s–1)b=2 m b=3 m b=4 m b=5 m b=6 m b=7 m b=8 m b=9 m b=10 m b=11 m LSPIV0.130.420.100.130.380.270.040.140.010.130.080.79 STIV0.020.030.320.000.630.110.110.240.230.090.060.35 DIS0.080.160.120.160.280.190.020.210.170.150.020.24 FD-DIS0.020.010.130.090.280.150.030.300.340.280.050.52 FD-DIS-G0.030.020.100.070.300.170.070.270.300.150.020.19

分析平均流速的相对误差可得（表5）：本算法得到的平均流速是5种方法里面最精确的，与流速仪真值相比相对误差仅为4.3%，本算法既稳定的接近于流速仪真值且与流速仪真值的绝对误差小，相比LSPIV方法和STIV方法，本文算法分别提高了12.7%和8.5%。

分析流量的相对误差可知（表5）：与流速仪测量的流量真值相比，本方法的相对误差仅为4.2%，符合水文站的监测要求，其精度比LSPIV方法和STIV方法分别提高了13.1%和3.5%，与没有计算运动显著性图的DIS方法和未经过奇异值处理的FD-DIS方法相比，精度分别提高了1.1%和7.2%，说明通过运动显著性计算捕捉到的微小运动和奇异值的处理对本算法精度的提高都是至关重要的。

由图10（d）可知，本研究方法耗时74.337 6 s，相比于广泛使用的LSPIV方法提高了13倍，虽比没有计算运动性显著性图的DIS方法和未经过奇异值处理的FD-DIS方法耗时要长，但相比算法提高的精度，60 s左右的耗时增加是在水文站实时监测的可接受范围内的。

2.3 普适性验证实验

在贵州省黔西南州保基水文站进行了算法普遍适用性验证试验，根据该水文站多年实测，岸边系数为0.70，水面流速系数为0.88。试验地面标定点布置如图11所示。

图11 普适性验证实验的地面标定点布设示意图Fig. 11 The layout diagram of ground mark for universal validation experiment

左岸布置

A

、

B

2个点，右岸布置

D

、

C

2个点，断面线起点为

E

，断面线终点为

F

，断面数据如表6所示。于同一时段进行流速仪测量与视频拍摄。测量期间水流稳定、天气阴、测量流速仪型号为LS25-3A。利用流速仪法在起点距分别为7、12、17、22、27、32、37和42 m的垂线上测得垂线平均流速，结果如表7所示，并计算得到了河流平均流速为1.57 m/s，断面流量149 m/s。根据流速仪的位置，在相同起点距的位置确定测速点，如图11中的小圆点所示。利用FD-DIS-G方法和DIS方法测量测速点的表面流速并计算垂线平均流速、流量及平均流速可以得到表8，将表8的结果与表7的流速仪测量结果进行对比分析和计算可以得到两种方法与流速仪真值相比的绝对误差（表9）和相对误差（表10）。

表6 普适性验证实验的断面数据
Tab. 6 Sectional data of the universal validation experiment

注：断面总面积为94.9 m。

起点距b/m深度/m 1.81.37 7.01.90 12.01.75 17.01.72 22.01.84 27.01.88 32.02.12 37.02.33 42.02.47 48.41.35

表7 普适性验证实验的流速仪测量结果
Tab. 7 Current meter measurement results of the universal validation experiment

部分面起点距b/m 垂线平均流速 /(m·s–1)部分平均流速/(m·s–1)积/m2部分流量/(m3·s–1)0（岸边）0——0～7—0.9310.109.39 7 1.33——7～12—1.608.9714.40 121.88——12～17—1.948.5916.70 172.01——17～22—2.148.8619.00 222.27——22～27—2.039.1818.60 271.79——27～32—1.7210.0017.20 321.64——32～37—1.6811.1018.60 371.73——37～42—1.5812.3019.40 421.43——42～51.4—1.0015.8015.80 51.40——

表9 流速仪方法、FD-DIS-G和DIS方法的绝对误差分析
Tab. 9 Analysis between Flow meter, FD-DIS-G and DIS

方法垂线平均流速绝对误差/(m·s–1)平均流速绝对误差/(m·s–1)流量绝对误差/(m3·s–1)b=7 mb=12 mb=17 mb=22 mb=27 mb=32 mb=37 mb=42 m FD-DIS-G0.400.400.190.510.060.350.190.470.044 DIS0.960.930.690.620.120.010.340.040.3836

分析表8、9的垂线平均流速可以看出，当河流水面宽大于50 m，河流平均流速接近2 m/s时，FD-DISG方法除了靠近左岸和右岸的两个测量点相对误差超过30%，其余测量点的相对误差都在23%以内，相比于未进行运动显著性计算的DIS方法，其稳定性和精确性得到了提升。FD-DIS-G方法在靠近左岸和右岸的测量点测得的值比流速仪真值分别偏小0.40 m/s和偏大0.47 m/s，这是因为靠近岸边的测量点受到了岸边水草等杂物和镜头畸变的双重影响。分析表10的平均流速和流量可知，FD-DIS-G方法所计算出的流量和平均流速与流速仪测量所得到的结果相比，相对误差分别仅为2.7%和2.5%，相比于未进行运动显著性计算的DIS方法分别提高了21.5%和21.7%。结合比测实验和普适性验证实验的结果可知，本文提出的方法适用于不同环境和不同规模的河流流量、平均流速的测量，具有较好稳定性和广泛适用性。

表8 流速仪方法、FD-DIS-G和DIS方法的对比
Tab. 8 Comparison between Flow meter, FD-DIS-G and DIS

垂线平均流速/(m·s–1)平均流速/(m·s–1) 流量/(m3·s–1)b=2 m b=3 m b=4 m b=5 m b=6 m b=7 m b=8 m b=9 m b=10 m b=11 m流速仪1.331.882.012.271.791.641.731.431.571491.331.88 DIS0.370.951.321.651.671.651.391.391.191130.370.95 FD-DIS-G0.931.481.821.761.851.991.541.901.531450.931.48方法

表10 流速仪方法、FD-DIS-G和DIS方法的相对误差分析
Tab. 10 Analysis between Flow meter, FD-DIS-G and DIS

方法垂线平均流速相对误差/%平均流速相对误差/%流量相对误差/%b=7 mb=12 mb=17 mb=22 mb=27 mb=32 mb=37 mb=42 m FD-DIS-G30.121.49.522.33.121.511.132.92.52.7 DIS71.849.434.327.47.00.819.52.624.224.2

3 结 论

本文研究了一种基于帧间差分与DIS快速光流结合的分组图像法(FD-DIS-G)测流方法。相比于LSPIV和STIV，利用帧差法计算运动显著性图捕捉不显著的运动区域，并利用DIS快速密集光流法计算运动显著性图上的密集流场，还分组对数据进行奇异值处理，提高了算法精度。实验结果表明，本文提出的方法计算出的河流平均流速和断面流量与流速仪测流结果相比相对误差小于5%，并且算法运行时间相比于相关性计算耗时严重的LSPIV方法提高了约13倍。本文提出的方法能满足实际工程测量要求，是一种设备和人力成本低且可连续的河流流量监测方法。

本文未考虑风场对流速的影响，在低流速情况下风场影响较大，未来工作将考虑风场进行建模及分析计算。同时目前视频测流工作中对算法的有效评价指标研究还不充分，未来将研究相关评价指标并建立基准评价数据集，结合深度学习的目标跟踪方法，搭建端到端的图像法测流神经网络框架，以建立更加科学合理的计算模型。