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

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

王剑平, 朱芮, 张果, 等. 帧差与快速密集光流结合的图像法测流研究 [J]. 工程科学与技术, 2022, 54(4): 195-207. doi: 10.15961/j.jsuese.202101046
引用本文: 王剑平, 朱芮, 张果, 等. 帧差与快速密集光流结合的图像法测流研究 [J]. 工程科学与技术, 2022, 54(4): 195-207. doi: 10.15961/j.jsuese.202101046
WANG Jianping, ZHU Rui, ZHANG Guo, et al. Image Flow Measurement Based on the Combination of Frame Difference and Fast and Dense Optical Flow [J]. Advanced Engineering Sciences, 2022, 54(4): 195-207. doi: 10.15961/j.jsuese.202101046
Citation: WANG Jianping, ZHU Rui, ZHANG Guo, et al. Image Flow Measurement Based on the Combination of Frame Difference and Fast and Dense Optical Flow [J]. Advanced Engineering Sciences, 2022, 54(4): 195-207. doi: 10.15961/j.jsuese.202101046

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

基金项目: 国家重点研发计划资助项目(2017YFB0306405);国家自然科学基金资助项目(61364008);云南省基础研究计划重点资助项目(202101AS070016)
详细信息
    • 收稿日期:  2021-10-18
    • 网络出版时间:  2022-07-06 10:57:50
  • 作者简介:

    王剑平(1975—),男,副教授,博士. 研究方向:机器视觉与智能控制. E-mail:kmustwjp@126.com

  • 中图分类号: TV123

Image Flow Measurement Based on the Combination of Frame Difference and Fast and Dense Optical Flow

  • 摘要: 图像法测流技术因其简便、高效、安全等优点得到了普遍的关注,开始应用于国内外水文站。目前主流的图像法测流技术主要有大尺度粒子图像测速(large-scale particle image velocimetry,LSPIV)和时空图像测速(spatio-temporal image velocimetry,STIV)等,LSPIV方法采用天然粒子作为示踪物,利用空域互相关法处理图像获得流场矢量,STIV方法根据河流的主流向设定测速线,对视频图像中测速线的灰度进行采样形成时空图像,利用时空图像中的纹理角求得流速。LSPIV方法依赖示踪物的可见性,存在稳定性差等缺点,STIV方法存在对断面流态稳定性要求高和仅能测量1维流速等缺点。本文提出一种结合帧间差分与快速密集光流的分组测流方法(frame difference-fast optical flow using dense inverse search-grouping,FD-DIS-G),利用帧差法计算运动显著性图,通过捕捉细微的水面运动处理河流运动在视频中表现不明显的问题,减少对天然示踪物的依赖,使用快速密集光流法(fast optical flow using dense inverse search,DIS)计算运动显著性图中小块区域之间的密集光流位移,提高流量测量的精度和时效性,有效克服流态不稳定的情况。同时,设计一种分组处理奇异值的方法,提高了算法的整体准确性,增强了算法的稳定性。将流速仪测量得到的垂线平均流速、平均流速以及断面流量作为比测标准,利用水文站所拍得的天然河道水流视频进行比测实验,实验结果表明,相比于目前广泛使用的图像测流方法,本文方法在平均流速和断面流量上的精度有明显的提升,垂线平均流速测量的稳定性有显著的增强且实时性好。

     

    Abstract: Due to its simplicity, efficiency and safety, image velocimetry has gained widespread attention and is beginning to be used in domestic and international hydrographic stations. At present, the mainstream image-based flow measurement techniques mainly include large-scale particle image velocimetry (LSPIV) and spatio-temporal image velocimetry (STIV), etc. While the LSPIV method uses natural particles as tracers and uses the spatial correlation method to process the image to obtain the flow field vector, the STIV method sets the velocity line according to the main flow direction of the river, samples the grey scale of the velocity line in the video image to form the spatio-temporal image, and uses the texture angle in the spatio-temporal image to obtain the flow velocity. Generally, the LSPIV method has disadvantages such as poor stability due to its reliance on the visibility of the tracer, and the STIV method has disadvantages such as high requirements for the stability of the cross-sectional flow regime and the ability to measure only one-dimensional flow velocities. In this paper, a grouping flow measurement method based on the combination of inter-frame difference and fast optical flow using dense inverse search (FD-DIS-G) was proposed. In the method, the frame difference method was used to calculate a kinematic saliency map that captures subtle water surface motions to deal with the problem that river motions are not evident in the video, which reduces the reliance on natural tracers. Then, the fast optical flow is used to calculate the dense optical flow displacements between small areas of the kinematic prominence map, improving the accuracy and timeliness of flow measurements and effectively overcoming the flow instability. Meanwhile, a method of grouping processing abnormal values was designed, which improves the overall accuracy of the algorithm and enhances the stability of the algorithm. The results of the experiments showed that the accuracy of the mean flow velocity and cross-sectional flow rate is significantly improved compared with the widely-used image flow measurement methods, and the proposed mean flow velocity measurement methods achieves real-time performance and its stability of is significantly enhanced.

     

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

    近几年来,机器视觉及传感器技术的迅速发展推动非接触式测流技术取得巨大的突破[2],基于光学、图像以及雷达的各种测流技术都已经得到了快速地发展和应用,尤其是图像法测流技术,该技术因其简便、高效、安全等优点得到了普遍的关注[5]。粒子图像测速(particle image velocimetry,PIV)技术[6-7]发展于20世纪80年代,该技术在流场中布撒示踪粒子,以激光片光作为光源,用底片或CCD(charge coupled device,CCD)相机记录示踪粒子的图像,利用自相关法、互相关法或光学杨氏条纹法处理这些图像,从而获得流场的速度分布[8-9]。Fujita等[10]提出了大尺度粒子图像测速(large-scale particle image velocimetry,LSPIV)方法,改进PIV技术用于大尺度水面流场观测。该方法以波纹、植物碎片、泡沫等天然水面漂浮物和天然水面模式代替粒子作为示踪物,用普通的数字相机或视频录像机记录自然光光源下的河流图像[11],简化了硬件配置,还避免了示踪物对河流的污染。该方法自提出以来,被广泛地应用于水文站流量、流速的监测[12-13]。但由于采用极易受到环境影响的天然水面漂浮物和诸如泡漩等不稳定的天然水面模式作为示踪物,该方法存在稳定性差的问题。且LSPIV方法在运动估计时采用空域相关匹配法,相关性计算的时间复杂度达O(n4),所需存储空间大且计算效率低。相比于LSPIV,Fujita等[14-15]在2007年提出了基于时空图像的测流(spatio-temporal image velocimetry,STIV)方法。该方法在断面流量测量方面展现出了更好的性能,目前正在被广泛地研究,但该方法只能获得一维流速场,且对河流流态稳定性要求高[16]。王万良等[17]将深度学习用于流速测量,采用条件边界平衡生成对抗网络对流速图像进行生成,采用多特征融合的卷积分类网络对流速图像进行分类,此方法仅能获得流速区间,但为非接触式测流技术提供了一种全新的思路。Lin等[18]结合无人机利用热水作为示踪剂,利用Lucas-Kanade光流法[19]取得了优于常用跟踪方法的测流结果。但由于将热水作为大尺度天然河流的示踪剂并不具有普遍的可行性,该方法也存在局限性。Tauro等[20]利用金字塔Lucas-Kanade光流法[19]跟踪了FAST、Shi-Tomasi、SIFT、SURF、ORB5种特征点,并利用先验信息对特征点的跟踪轨迹进行滤波,只保留视场中与实际物体运动相关的轨迹,此方法在实验过程中取得了良好的结果,但特征点的选取是一个复杂过程,考虑天然河道实际工况的复杂性,在不同的现场环境下选取光流跟踪特征点的统一规则并没有制定,且先验信息的获取并不具有普遍性,故该方法在实际场景下稳定性较差。

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

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

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

    图  1  算法流程图
    Fig.  1  Algorithm flowchart
    下载: 全尺寸图片

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

    图  2  河流表面图与运动显著性图
    Fig.  2  River surface map and motion saliency map
    下载: 全尺寸图片

    在没有任何示踪剂的水面且时间间隔很短时,细微的水面运动是很难察觉的,如图2(a)(b)所示。这种情况下,水面像素点的跟踪是困难且不准确的,从而导致光流估计并不能很好地反应河流表面的速度。因此,本方法提出在运动显著性图上计算光流的方法,利用连续2帧图像中波浪运动造成的差异放大运动区域,生成运动显著性图,在运动显著性图上进行光流估计。定义视频中当前帧图像的灰度矩阵为 ${{\boldsymbol{F}}_t}$ ${f}_{i,j}^t$ 为矩阵 ${{\boldsymbol{F}}_t}$ 的元素,前一帧图像的灰度矩阵为 ${{\boldsymbol{F}}_{{{{t} - 1}}}}$ ${f}_{i,j}^{t - 1}$ 为矩阵 ${{\boldsymbol{F}}_{t - 1}}$ 的元素,t时刻的运动显著性图为 ${{\boldsymbol{D}}_t}$ ${d}_{i,j}^t$ ${{\boldsymbol{D}}_t}$ 元素,可以表示为:

    $$ {d}_{i,j}^t = \left\{ {\begin{array}{*{20}{l}} {{k_1}| {f_{i,j}^t - f_{i,j}^{t - 1}} |,}&{| {f_{i,j}^t - f_{i,j}^{t - 1}} | \ge h}; \\ {{k_2}| {f_{i,j}^t - f_{i,j}^{t - 1}} |,}&{| {f_{i,j}^t - f_{i,j}^{t - 1}} | \lt h} \end{array}} \right. $$ (1)

    式中: $h$ 为预设定的阈值; $ {k_1} $ $ {k_2} $ 为预设定的系数,通过给 $ {k_1} $ 设定比1大的值,给 $ {k_2} $ 设定比1小的值,可突出运动区域,如图2(d)所示。为后续更准确地光流估计提供了基础。

    由第1.1节得到运动显著性图,根据显著性图中断面线的起点加上起点距b可得测速点P的坐标,并以P为中心点确定当前帧(第t帧)的块区域 ${{\boldsymbol{G}}_t}$ ,下一帧的块区域 ${{\boldsymbol{G}}_{t + 1}}$ ,块区域大小为 $ (\delta \times \delta ) $ $ \delta $ 为块区域的边长,由相邻两个测速点之间的像素距离决定。于是对两帧块区域之间位移 $ \Delta d $ 的求解可转化为对两块区域之间的密集光流位移 ${\boldsymbol{A}}$ 进行计算,计算过程如下:

    Step1:以 $\;{\rho}$ 为尺度金字塔下采样系数,分别对 ${{\boldsymbol{G}}^t}$ ${{\boldsymbol{G}}^{t + 1}}$ 构建尺度金字塔,如图3所示。定义模板 ${{\boldsymbol{T}}_s}$ 为金字塔在第 $ s $ 层块区域 ${\boldsymbol{G}}^{t,s}$ 里的子区域,模板大小为 $ ({\theta _s} \times {\theta _s}) $ 。用均匀的网格块初始化 ${\boldsymbol{G}}^{t,s}$ 且每个网格块的大小和模板 $ {T_s} $ 相同,如图4所示,并用式(2)初始化第k个网格块 ${\boldsymbol{I}}^k$ 中心点的光流位移 $( {u_{k,init}}, {v_{k,init}} )$

    图  3  块区域金字塔
    Fig.  3  Block area pyramid
    下载: 全尺寸图片
    图  4  块区域网格化
    Fig.  4  Block area meshing
    下载: 全尺寸图片
    $$ ( {u_{k,init}}, {v_{k,init}} ) ={\boldsymbol{a}}_{(x_k/\rho, y_k/\rho)}^{s+1}\cdot\rho $$ (2)

    式中: $(x_k,y_k)$ 为网格块 ${\boldsymbol{I}}^k$ 中心点的坐标; ${\boldsymbol{a}}_{(x_k/\rho, y_k/\rho)}^{s+1} $ ${{\boldsymbol{A}}^{s + 1}} $ 的元素,包含了xy方向的位移; ${{\boldsymbol{A}}^{s + 1}}$ 为第 $ s + 1 $ 层的光流位移。

    Step2:求解每个网格块 ${\boldsymbol{I}}^k$ 中心点的光流位移 $( {u_{k}}, {v_{k}} )$ ,在区域 ${\boldsymbol{G}}^{{t + 1},s}$ 中搜索与 ${\boldsymbol{G}}^{t,s}$ 中网格块 ${\boldsymbol{ I}}^k$ 最匹配的子窗口,即找到 $(u_k,v_k)$ 使得网格块 ${\boldsymbol{I}}^k$ ${\boldsymbol{G}}^{{t + 1},s}$ 在子窗口上的灰度值的差异平方和最小化,如式(3)所示:

    $$ (u_k,v_k)= {{\rm{argmin}}} {\sum\limits_x {\left[ {{\boldsymbol{g}}_{(x'+u_k,y'+v_k)}^{t+1,s} - {\boldsymbol{i}}_{(x', y')}^k} \right]} ^2} $$ (3)

    式中, $(x', y') $ 表示块区域 ${\boldsymbol{G}}^{{t + 1},s} $ 里每个像素点的坐标, ${{\boldsymbol{g}}_{(x'+u_k,y'+v_k)}^{t+1,s}}$ ${\boldsymbol{G}}^{{t + 1},s} $ 的元素, $ {\boldsymbol{i}}_{(x', y')}^k $ ${\boldsymbol{I}}^k $ 的元素。式(3)要成立需要引入逆Lucas-Kanade方法[23],对 $ {u_i} $ 进行迭代优化,其中包括:

    ①通过对块区域 ${\boldsymbol{G}}^{{t + 1},s}$ 进行像素值提取和双线性插值找到一个更新向量 $(\Delta u_k,\Delta v_k)$ ,如式(4)所示:

    $$ (\Delta u_k,\Delta v_k) = {{\rm{argmin}}} \sum\limits_{(x',y')} {{{\left[ { {\boldsymbol{i}}_{(x'+\Delta u_k, y'+\Delta v_k)}^k - {\boldsymbol{g}}_{(x'+u_k,y'+v_k)}^{t+1,s}} \right]}^2}} $$ (4)

    ②更新 $(u_k,v_k) = (u_k-\Delta u_k,v_k-\Delta v_k)$ ;为了增强鲁棒性,对于更新的 $(u_k,v_k)$ ,若有 $(u_k,v_k)$ $(u_{k,init},v_{k,init})$ 差值的二范数 ${\| {({u_{i,init}},{v_{i,init}}) - (u_k,v_k)} \|_2}$ 超过了模板 $ {T_s} $ 的大小 $ \theta $ ,则将 $(u_k,v_k)$ 都重新设定为初始化的光流 $(u_{k,init},v_{k,init})$

    Step3:块区域 ${\boldsymbol{G}}^{t,s}$ 中每个网格块中心点的光流 $(u_k,v_k)$ 求解完毕后可得第 $ s $ 层块区域 ${\boldsymbol{G}}^{t,s}$ 的密集光流位移 ${{\boldsymbol{A}}^s}$ ,如式(5)所示:

    $$ {\boldsymbol{a}}_{(x'', y'')}^s = \frac{1}{Z}\sum\limits_{k=1}^{{N_s}} {\frac{{{\lambda _{(x'', y'')}^k}}}{{\max ( {1,{{\| {{{\boldsymbol{d}}_{(x'',y'')}^k}} \|}_2}} )}}} (u_k,v_k) $$ (5)

    式中: $(x'',y'') $ 为块区域 ${\boldsymbol{G}}^{t,s} $ 中像素点的坐标, ${\boldsymbol{a}}_{(x'', y'')}^s $ ${\boldsymbol{A}}^s $ 的元素, ${\lambda _{(x'', y'')}^k}$ 为网格块 ${\boldsymbol{I}}^k$ 与以 $(x'',y'')$ 为中心、和网格同大小区域的重叠度,当且仅当二者完全重叠时 ${\lambda _{(x'', y'')}^k} = 1$ $ {N_s} $ 为网格数目; ${{\boldsymbol{d}}}_{(x'', y'')}^k$ ${\boldsymbol{D}}^k $ 的元素; ${\boldsymbol{D}}^k$ 为网格块 ${\boldsymbol{I}}^k$ 和块区域 ${\boldsymbol{G}}^{{t + 1},s}$ 的像素值差矩阵; $ Z $ 为对块区域 ${\boldsymbol{G}}^{t,s} $ 内所有重叠度 ${\lambda _{(x'', y'')}^k}$ 的标准化,且有:

    $$ {{\boldsymbol{d}}}_{(x'', y'')}^k ={\boldsymbol{g}}_{(x''+u_k,y''+v_k)}^{t+1,s} - {\boldsymbol{i}}_{(x'', y'')}^k $$ (6)
    $$ Z = \frac{{\displaystyle\sum\limits_{k=1}^{{N_s}} {{\lambda _{i,x}}} }}{{\max \left( {1,{{\left\| {{{{\boldsymbol{d}}}_{(x'',y'')}^k}} \right\|}_2}} \right)}} $$ (7)

    式中, $ \max () $ 为取最大值的函数, $ {\left\| \;\; \right\|_2} $ 为二范数。

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

    $$ E(U) = \int_\varOmega \alpha \omega \left( {{E_{\rm{I}}^2}} \right) + \beta \omega \left( {{E_{\rm{G}}^2}} \right) + \gamma \omega \left( {{E_{\rm{S}}^2}} \right){\rm{d}}x $$ (8)

    式中: $\varOmega$ 为图像区域; ${E_{\rm{I}}}$ , ${E_{\rm{G}}}$ , ${E_{\rm{S}}}$ 分别为图像域上的强度数据项、梯度数据项和平滑度数据项[21],对每一项都进行正则化; $\omega \left( {{E^2}} \right) = \sqrt {{E^2} + {\varepsilon ^2}}$ (取 $ \varepsilon = 0.001 $ ); $ \alpha $ $ \;\beta $ $ \gamma $ 为正则化系数。

    Step5:经过每一个尺度的迭代计算后,最终求得密集光流位移 ${\boldsymbol{A}}$ ,即为块区域 ${{\boldsymbol{G}}^t}$ 中每个像素点到块区域 ${{\boldsymbol{G}}^{t + 1}}$ 中每个像素点的光流,则块区域之间的位移 $ \Delta d $ 可表示为:

    $$ \Delta d = \frac{{\displaystyle\sum \left\|{\boldsymbol{a}}_{(x,y)}\right\|_2}}{{\delta \times \delta }} $$ (9)

    式中, $(x,y) $ 为块区域 ${{\boldsymbol{G}}^t} $ 中像素点的坐标, ${\boldsymbol{a}}_{(x,y)} $ ${{\boldsymbol{A}}}$ 的元素, $ \delta \times \delta $ 为块区域的大小。

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

    图  5  视频分组处理图
    Fig.  5  Video grouping processing diagram
    下载: 全尺寸图片

    首先,将由运动显著性图计算得到的光流位移分为M组,每组有由4帧相邻运动显著图得到的3个光流位移数据 $\{ \Delta {d_{i,1}},\Delta {d_{i,2}}, \Delta {d_{i,3}} \}$ ,其中 $ i $ 为分组序号,然后将3个数据中的最大值 $ \Delta {d_{i,\max }} $ 最小值 $ \Delta {d_{i,\min }} $ 去掉,取中间值 $\Delta {d_{i,{\rm{mid}}}}$ 作为这组数据的平均帧间位移 $\Delta {d_{i,{\rm{avg}}}}$ ,最后对M个平均帧间位移数据进行去掉奇异值的处理,处理过程为:

    Step1:若 ${\Delta {d_{i,{\text{ mid }}}}}$ 是独立同分布的,根据中心极限定理[24]有,对于任何 ${\Delta {d_{i,{\text{ mid }}}}}$ ,服从正态分布,将其转化为标准正态分布,则分布密度函数 $f\left( {\Delta {d_{i,{\text{ mid }}}}}- \mu \right)$ 为:

    $$ f\left( {\Delta {d_{i,{\text{ mid }}}}} - \mu\right) = \frac{1}{{\sigma \sqrt {2\text{π} } }}\exp \left( { - \frac{{{{\left( {\Delta {d_{i,{\rm{mid}}}} - \mu } \right)}^2}}}{{2{\sigma ^2}}}} \right) $$ (10)

    式中,μM ${\Delta {d_{i,{{{\rm{ mid}} }}}}}$ 的期望, $\sigma$ M ${\Delta {d_{i,{{ {\rm{mid}} }}}}}$ 标准差。由分布密度函数 $f\left( {\Delta {d_{i,{{ {\rm{mid}} }}}}}- \mu \right)$ 的定义可知, ${\Delta {d_{i,{\rm{mid}}}} }$ 取值在对称区间 $[-\eta,\eta]$ 内的概率为该区间上密度函数的积分,即:

    $$ \begin{aligned}[b]& P\left( { - \eta \le \Delta {d_{i,{\rm{mid}}}}- \mu \le \eta} \right) = \\&\qquad \frac{1}{{\sigma \sqrt {2\text{π} } }}\int_{ - \eta}^{\eta} {{{\rm{e}}^{ - \frac{{{{\left( {\Delta {d_{i,{\rm{mid}}}} - \mu } \right)}^2}}}{{2{\sigma ^2}}}}}} {\rm{d}}\left( {\Delta {d_{i,{\rm{mid}}}} - \mu } \right) =\\&\qquad \frac{2}{{\sigma \sqrt {2\text{π} } }}\int_0^{\eta} {{{\rm{e}}^{ - \frac{{{{\left( {\Delta {d_{i,{\rm{mid}}}} - \mu } \right)}^2}}}{{2{\sigma ^2}}}}}} {\rm{d}}\left( {\Delta {d_{i,{\rm{mid}}}} - \mu } \right)\\[-15pt] \end{aligned} $$ (11)

    $\delta= \dfrac{\Delta {d_{i,{\rm{mid}}}} - \mu } {\sigma }$ ,有

    $$ \begin{aligned}[b] & P\left( { - \eta \le \Delta {d_{i,{\rm{mid}}}} - \mu\le \eta} \right) = \\&\qquad \frac{2}{{\sqrt {2\text{π} } }}\int_0^{\frac{\eta}{\sigma }} {{{\rm{e}}^{ - \frac{{{\delta ^2}}}{2}}}} {\rm{d}}\delta = 2{\varPhi} \left( {\frac{\eta}{\sigma }} \right) \end{aligned} $$ (12)

    式中: ${\varPhi (\;)}$ 为正态分布函数,对于给定误差区间 $[- \eta, \eta ]$ ,由概率积分查表得当 $\dfrac{\eta}{\sigma } = 2$ 时,P=0.9545[24],即 $\Delta {d_{i,{\rm{mid}}}}$ 取值在对称区间 $[ \mu - 2\sigma, \mu + 2\sigma ]$ 的概率为0.9545。又因为在1次试验中概率小于5%的事件被认为几乎是不可能发生的,称为小概率事件。所以,可认为落在区间 $[ \mu - 2\sigma, \mu + 2\sigma ]$ 外的数据为奇异值并舍掉。

    Step2:若这M个数据不满足独立同分布,则将M组数据的平均值 ${\Delta {\bar{d}_{i,{\rm{mid}}}}}$ 作为最佳估计值,和约定真值,引入Z分数:

    $$ Z = \frac{{\Delta {d_{i,{\rm{mid}}}} - {\Delta {\bar{d}_{i,{\rm{mid}}}}} }}{s} $$ (13)

    式中, $s$ M个数据的标准差。Z分数以标准差为单位衡量某一原始数据偏离平均值的距离,Z分数的值就表示了原始数据与平均值之间相差的标准差个数。本算法定义,与最佳估计值的偏差超过2倍标准差的测量值称为奇异值,因此当 $|{{Z}}| > 2$ 时,该数据被认为是奇异值,应该舍掉该数据。

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

    在河流两岸选取4个标记点ABCD,并近似的认为4个标记点处于同一平面且与水平面共面。令ABCD4个点两两之间的距离为 ${L_{AB}}$ ${L_{BC}}$ ${L_{CD}}$ ${L_{DA}}$ ${L_{AC}}$ ${L_{BD}}$ ,以A为原点,AB所在的方向为x轴的正方向,由勾股定理求得4个标记点在实际平面中的坐标:

    $$ {\left\{ {\begin{array}{l} {x_A} = 0,{y_A} = 0; \\ {x_B} = L_{AB},{y_A} = 0; \\ {x_C} = \dfrac{{{L^2_{AB}} + {L^2_{AC}} - {L^2_{BC}}}}{{2L_{AB}}},{y_C} = \sqrt {L^2_{AC} - \dfrac{{L^2_{AB} + L^2_{AC} - {L^2_{BC}}}}{{2L_{AB}}}} ,\\ \qquad \cos (\angle {{BAC}}) \lt 0 ;\\ {x_C} = - \dfrac{{{L^2_{AB}} + {L^2_{AC}} - {L^2_{BC}}}}{{2L_{AB}}},{y_C} = \sqrt {{L^2_{AC}} - \dfrac{{{L^2_{AB}} + {L^2_{AC}} - {L^2_{BC}}}}{{2L_{AB}}}} ,\\ \qquad\cos (\angle {{BAC}}) \ge 0 ;\\ {x_D} = \dfrac{{L^2_{BD} - L^2_{DA} - L^2_{AB}}}{{2L_{AB}}},{y_D} = \sqrt {L^2_{DA} - \dfrac{{L^2_{BD} - L^2_{DA} - L^2_{AB}}}{{2L_{AB}}}} ,\\ \qquad\cos (\angle {{BAD}}) \lt 0 ;\\ {x_D} = - \dfrac{{L^2_{BD} - L^2_{DA} - L^2_{AB}}}{{2L_{AB}}},{y_D} = \sqrt {L^2_{DA} - \dfrac{{L^2_{BD} - L^2_{DA} - L^2_{AB}}}{{2L_{AB}}}} ,\\ \qquad \cos (\angle {{BAD}}) \ge 0 \end{array}} \right. }$$ (14)

    将4个标记点的图上坐标和实际坐标一一对应,由透视投影变换[25]定理求出单应矩阵 ${\boldsymbol M} = \left[ {\begin{array}{*{20}{l}} {{m_{11}}}&{{m_{12}}}&{{m_{13}}} \\ {{m_{21}}}&{{m_{22}}}&{{m_{23}}} \\ {{m_{31}}}&{{m_{32}}}&{{m_{33}}} \end{array}} \right]$ ,令图像平面内点R的坐标为 $ (u,v) $ ,其在世界坐标系中对应的坐标为 $ (x,y,z) $ ,则有对应关系:

    $$ \left[ {\begin{array}{*{20}{l}} x \\ y \\ {\textit{z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {{m_{11}}}&{{m_{12}}}&{{m_{13}}} \\ {{m_{21}}}&{{m_{22}}}&{{m_{23}}} \\ {{m_{31}}}&{{m_{32}}}&{{m_{33}}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} u \\ v \\ 1 \end{array}} \right] $$ (15)
    $$ \left\{ {\begin{array}{*{20}{l}} {x = {m_{11}}u + {m_{12}}v + {m_{13}}}, \\ {y = {m_{21}}u + {m_{22}}v + {m_{23}}}, \\ {{\textit{z}} = {m_{31}}u + {m_{32}}v + {m_{33}}} \end{array}} \right. $$ (16)

    由式(15),(16)可得经过变换后世界坐标系中物理平面中的坐标为: ${x}^{\prime }=\dfrac{x}{{\textit{z}}},{y}^{\prime }=\dfrac{y}{{\textit{z}}}$ ,即:

    $$ \left\{ {\begin{array}{{l}} {{x^\prime } = \dfrac{{{m_{11}}u + {m_{12}}v + {m_{13}}}}{{{m_{31}}u + {m_{32}}v + {m_{33}}}}} ,\\ {{y^\prime } = \dfrac{{{m_{21}}u + {m_{22}}v + {m_{23}}}}{{{m_{31}}u + {m_{32}}v + {m_{33}}}}} \end{array}} \right. $$ (17)

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

    由第1.4节所述的变换可以得到实际的帧间位移,并除以时间间隔得到以“m·s–1”为单位的表面流速。设得到测速点 $ i $ 的表面流速为 $ {v_i} $ $ i = 1,2,\cdots,n $ n为断面线上测速点的个数),测速点 $ i - 1 $ $ i $ 对应的过水断面的面积为 $ {s_i} $ ,测速点 $ i - 1 $ $ i $ 之间的平均表面流速为 $ {\bar v_i} $

    $$ {\bar v_i} = \frac{{{v_{i - 1}} + {v_i}}}{2} $$ (18)
    $$ {s_i} = \frac{{{d_{i - 1}} + {d_i}}}{2}{l_i} $$ (19)

    式中, $ {d_i} $ 为测速点 $ i $ 对应的水深, $ {l_i} $ $ i $ 点和 $ i - 1 $ 点之间的间隔距离。根据水流条件选取水面流速系数k=0.78~0.88,利用平均表面流速乘以系数即可得到该点的垂线平均流速 $ v_i^\prime = k{\bar v_i} $ 。根据河流流量测验规范[26]可得过水断面 $ {s_i} $ 的流量为 $ {q_i} = v_i^\prime {s_i} $ ,则河流的总流量 $ Q $ 为:

    $$ Q = \sum\limits_{i = 1}^n {{q_i}} = \sum\limits_{i = 1}^n {v_i'} {s_i} $$ (20)

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

    $$ \bar v = \frac{Q}{S} $$ (21)

    式中, $ S $ 为断面总面积。

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

    比测实验在某水文站经过断面修整的河道进行了比测试验,选取LSPIV方法[10]、STIV方法[16]、DIS方法[21](未进行运动显著性处理)以及未进行分组奇异值处理的方法(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所示。左岸布置AB两个点,右岸布置DC两个点,起点与D点重合,终点与A点重合,以DA为断面线,断面数据如表1所示。

    表  1  比测实验的断面数据
    Table  1  Cross section data of the comparison test
    起点距b/m 深度/m
    0.1 0.59
    0.5 0.70
    1.0 0.81
    1.5 0.82
    2.0 0.76
    2.5 0.86
    3.0 0.93
    3.5 0.90
    4.0 0.96
    4.5 1.00
    5.0 0.99
    5.5 0.94
    6.0 1.04
    6.5 0.95
    7.0 0.86
    7.5 0.78
    8.0 0.77
    8.5 0.70
    9.0 0.92
    9.5 0.85
    10.0 0.81
    10.5 0.74
    11.0 0.55
    11.5 0.42
    11.9 0.39
      注:断面总面积为9.70 m2
    图  6  比测实验的地面标定点布设示意图
    Fig.  6  Layout diagram of ground mark in comparison test
    下载: 全尺寸图片

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

    表  2  比测实验的流速仪测量结果
    Table  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.29 1.45 0.42
    2 0.36
    2~3 0.40 0.85 0.34
    3 0.43
    3~4 0.39 0.92 0.36
    4 0.35
    4~5 0.40 0.99 0.40
    5 0.45
    5~6 0.60 0.98 0.59
    6 0.76
    6~7 0.72 0.95 0.68
    7 0.68
    7~8 0.64 0.80 0.51
    8 0.61
    8~9 0.54 0.77 0.42
    9 0.48
    9~10 0.43 0.86 0.37
    10 0.38
    10~11 0.44 0.71 0.31
    11 0.49
    11~11.9 0.39 0.42 0.16
    11.9(岸边) 0
    2.2.1   运动显著性图

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

    图  7  获取运动显著性图过程示意
    Fig.  7  Obtaining motion saliency map
    下载: 全尺寸图片
    2.2.2   密集光流场

    由于LSPIV方法[10]使用空域或频域相关匹配来获得运动矢量时相关性的计算耗时长,本算法提出利用DIS方法[21]进行运动矢量的估计,在相邻两帧运动显著性图中选取块区域计算密集光流位移。由于没有在全局区域进行密集光流的计算,且在计算密集光流时运用了逆向搜索的法则,使得本方法不仅可以得到亚像素精度的运动矢量,并在运算速度上也得到了保证。将视频图像中得到的密集光流可视化,如图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方法[10]、STIV方法[16]、DIS方法以及FD-DIS方法的测流结果如表3所示。将FD-DIS-G方法得到的垂线平均流速、平均流速和流量分别与流速仪法测量结果进行对比得到结果如图10(a)(b)(c)所示。将表3中5种方法的测量结果与流速仪真值相比可计算出5种方法的绝对误差(表4)和相对误差(表5)。5种方法的耗时对比结果如图10(d)所示。

    表  3  5种方法测量结果对比
    Table  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.36 0.43 0.35 0.45 0.76 0.68 0.61 0.48 0.38 0.49 0.47 4.56
    LSPIV 0.23 0.01 0.45 0.58 0.38 0.41 0.57 0.62 0.37 0.62 0.39 3.77
    STIV 0.38 0.46 0.03 0.45 0.13 0.79 0.72 0.72 0.61 0.40 0.43 4.21
    DIS 0.28 0.27 0.47 0.61 0.48 0.49 0.59 0.69 0.55 0.34 0.45 4.32
    FD-DIS 0.38 0.44 0.48 0.54 0.48 0.53 0.58 0.78 0.72 0.77 0.52 5.08
    FD-DIS-G 0.33 0.41 0.45 0.52 0.46 0.51 0.54 0.75 0.68 0.64 0.49 4.75
    表  4  5种方法的绝对误差对比
    Table  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
    LSPIV 0.13 0.42 0.10 0.13 0.38 0.27 0.04 0.14 0.01 0.13 0.08 0.79
    STIV 0.02 0.03 0.32 0.00 0.63 0.11 0.11 0.24 0.23 0.09 0.06 0.35
    DIS 0.08 0.16 0.12 0.16 0.28 0.19 0.02 0.21 0.17 0.15 0.02 0.24
    FD-DIS 0.02 0.01 0.13 0.09 0.28 0.15 0.03 0.30 0.34 0.28 0.05 0.52
    FD-DIS-G 0.03 0.02 0.10 0.07 0.30 0.17 0.07 0.27 0.30 0.15 0.02 0.19
    表  5  5种方法的相对误差对比
    Table  5  Relative error comparison of five methods
    方法 垂线平均流速相对误差/% 平均流速相
    对误差/%
    流量相对
    误差/%
    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
    LSPIV 34.93 97.70 28.67 29.21 50.21 40.15 7.01 28.28 3.61 26.86 17.00 17.30
    STIV 4.60 6.70 91.00 0.80 82.80 16.80 18.0 49.60 61.50 17.70 12.80 7.70
    DIS 22.20 37.20 34.30 35.60 36.80 27.90 3.30 43.80 44.70 30.60 4.30 5.30
    FD-DIS 5.60 2.30 37.10 20.00 36.80 22.10 4.90 62.50 89.50 57.10 10.60 11.40
    FD-DIS-G 8.30 4.70 28.60 15.60 39.50 25.00 11.50 56.30 79.00 30.60 4.30 4.20
    图  10  5种方法测量结果
    Fig.  10  Five methods of measurement results
    下载: 全尺寸图片

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

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

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

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

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

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

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

    表  6  普适性验证实验的断面数据
    Table  6  Sectional data of the universal validation experiment
    起点距b/m 深度/m
    1.8 1.37
    7.0 1.90
    12.0 1.75
    17.0 1.72
    22.0 1.84
    27.0 1.88
    32.0 2.12
    37.0 2.33
    42.0 2.47
    48.4 1.35
      注:断面总面积为94.9 m2
    表  7  普适性验证实验的流速仪测量结果
    Table  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.93 10.10 9.39
    7 1.33
    7~12 1.60 8.97 14.40
    12 1.88
    12~17 1.94 8.59 16.70
    17 2.01
    17~22 2.14 8.86 19.00
    22 2.27
    22~27 2.03 9.18 18.60
    27 1.79
    27~32 1.72 10.00 17.20
    32 1.64
    32~37 1.68 11.10 18.60
    37 1.73
    37~42 1.58 12.30 19.40
    42 1.43
    42~51.4 1.00 15.80 15.80
    51.4 0
    表  8  流速仪方法、FD-DIS-G和DIS方法的对比
    Table  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.33 1.88 2.01 2.27 1.79 1.64 1.73 1.43 1.57 149 1.33 1.88
    DIS 0.37 0.95 1.32 1.65 1.67 1.65 1.39 1.39 1.19 113 0.37 0.95
    FD-DIS-G 0.93 1.48 1.82 1.76 1.85 1.99 1.54 1.90 1.53 145 0.93 1.48
    表  9  流速仪方法、FD-DIS-G和DIS方法的绝对误差分析
    Table  9  Analysis between Flow meter, FD-DIS-G and DIS
    方法 垂线平均流速绝对误差/(m·s–1) 平均流速绝对
    误差/(m·s–1)
    流量绝对
    误差/(m3·s–1)
    b=7 m b=12 m b=17 m b=22 m b=27 m b=32 m b=37 m b=42 m
    FD-DIS-G 0.40 0.40 0.19 0.51 0.06 0.35 0.19 0.47 0.04 4
    DIS 0.96 0.93 0.69 0.62 0.12 0.01 0.34 0.04 0.38 36
    表  10  流速仪方法、FD-DIS-G和DIS方法的相对误差分析
    Table  10  Analysis between Flow meter, FD-DIS-G and DIS
    方法 垂线平均流速相对误差/% 平均流速相
    对误差/%
    流量相对
    误差/%
    b=7 m b=12 m b=17 m b=22 m b=27 m b=32 m b=37 m b=42 m
    FD-DIS-G 30.1 21.4 9.5 22.3 3.1 21.5 11.1 32.9 2.5 2.7
    DIS 71.8 49.4 34.3 27.4 7.0 0.8 19.5 2.6 24.2 24.2

    分析表89的垂线平均流速可以看出,当河流水面宽大于50 m,河流平均流速接近2 m/s时,FD-DIS-G方法除了靠近左岸和右岸的两个测量点相对误差超过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%。结合比测实验和普适性验证实验的结果可知,本文提出的方法适用于不同环境和不同规模的河流流量、平均流速的测量,具有较好稳定性和广泛适用性。

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

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

  • 图  1   算法流程图

    Fig.  1   Algorithm flowchart

    下载: 全尺寸图片

    图  2   河流表面图与运动显著性图

    Fig.  2   River surface map and motion saliency map

    下载: 全尺寸图片

    图  3   块区域金字塔

    Fig.  3   Block area pyramid

    下载: 全尺寸图片

    图  4   块区域网格化

    Fig.  4   Block area meshing

    下载: 全尺寸图片

    图  5   视频分组处理图

    Fig.  5   Video grouping processing diagram

    下载: 全尺寸图片

    图  6   比测实验的地面标定点布设示意图

    Fig.  6   Layout diagram of ground mark in comparison test

    下载: 全尺寸图片

    图  7   获取运动显著性图过程示意

    Fig.  7   Obtaining motion saliency map

    下载: 全尺寸图片

    图  8   密集光流可视化

    Fig.  8   Dense optical flow visualization

    下载: 全尺寸图片

    图  9   奇异数据处理

    Fig.  9   Abnormal data processing

    下载: 全尺寸图片

    图  10   5种方法测量结果

    Fig.  10   Five methods of measurement results

    下载: 全尺寸图片

    图  11   普适性验证实验的地面标定点布设示意图

    Fig.  11   The layout diagram of ground mark for universal validation experiment

    下载: 全尺寸图片

    表  1   比测实验的断面数据

    Table  1   Cross section data of the comparison test

    起点距b/m 深度/m
    0.1 0.59
    0.5 0.70
    1.0 0.81
    1.5 0.82
    2.0 0.76
    2.5 0.86
    3.0 0.93
    3.5 0.90
    4.0 0.96
    4.5 1.00
    5.0 0.99
    5.5 0.94
    6.0 1.04
    6.5 0.95
    7.0 0.86
    7.5 0.78
    8.0 0.77
    8.5 0.70
    9.0 0.92
    9.5 0.85
    10.0 0.81
    10.5 0.74
    11.0 0.55
    11.5 0.42
    11.9 0.39
      注:断面总面积为9.70 m2

    表  2   比测实验的流速仪测量结果

    Table  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.29 1.45 0.42
    2 0.36
    2~3 0.40 0.85 0.34
    3 0.43
    3~4 0.39 0.92 0.36
    4 0.35
    4~5 0.40 0.99 0.40
    5 0.45
    5~6 0.60 0.98 0.59
    6 0.76
    6~7 0.72 0.95 0.68
    7 0.68
    7~8 0.64 0.80 0.51
    8 0.61
    8~9 0.54 0.77 0.42
    9 0.48
    9~10 0.43 0.86 0.37
    10 0.38
    10~11 0.44 0.71 0.31
    11 0.49
    11~11.9 0.39 0.42 0.16
    11.9(岸边) 0

    表  3   5种方法测量结果对比

    Table  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.36 0.43 0.35 0.45 0.76 0.68 0.61 0.48 0.38 0.49 0.47 4.56
    LSPIV 0.23 0.01 0.45 0.58 0.38 0.41 0.57 0.62 0.37 0.62 0.39 3.77
    STIV 0.38 0.46 0.03 0.45 0.13 0.79 0.72 0.72 0.61 0.40 0.43 4.21
    DIS 0.28 0.27 0.47 0.61 0.48 0.49 0.59 0.69 0.55 0.34 0.45 4.32
    FD-DIS 0.38 0.44 0.48 0.54 0.48 0.53 0.58 0.78 0.72 0.77 0.52 5.08
    FD-DIS-G 0.33 0.41 0.45 0.52 0.46 0.51 0.54 0.75 0.68 0.64 0.49 4.75

    表  4   5种方法的绝对误差对比

    Table  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
    LSPIV 0.13 0.42 0.10 0.13 0.38 0.27 0.04 0.14 0.01 0.13 0.08 0.79
    STIV 0.02 0.03 0.32 0.00 0.63 0.11 0.11 0.24 0.23 0.09 0.06 0.35
    DIS 0.08 0.16 0.12 0.16 0.28 0.19 0.02 0.21 0.17 0.15 0.02 0.24
    FD-DIS 0.02 0.01 0.13 0.09 0.28 0.15 0.03 0.30 0.34 0.28 0.05 0.52
    FD-DIS-G 0.03 0.02 0.10 0.07 0.30 0.17 0.07 0.27 0.30 0.15 0.02 0.19

    表  5   5种方法的相对误差对比

    Table  5   Relative error comparison of five methods

    方法 垂线平均流速相对误差/% 平均流速相
    对误差/%
    流量相对
    误差/%
    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
    LSPIV 34.93 97.70 28.67 29.21 50.21 40.15 7.01 28.28 3.61 26.86 17.00 17.30
    STIV 4.60 6.70 91.00 0.80 82.80 16.80 18.0 49.60 61.50 17.70 12.80 7.70
    DIS 22.20 37.20 34.30 35.60 36.80 27.90 3.30 43.80 44.70 30.60 4.30 5.30
    FD-DIS 5.60 2.30 37.10 20.00 36.80 22.10 4.90 62.50 89.50 57.10 10.60 11.40
    FD-DIS-G 8.30 4.70 28.60 15.60 39.50 25.00 11.50 56.30 79.00 30.60 4.30 4.20

    表  6   普适性验证实验的断面数据

    Table  6   Sectional data of the universal validation experiment

    起点距b/m 深度/m
    1.8 1.37
    7.0 1.90
    12.0 1.75
    17.0 1.72
    22.0 1.84
    27.0 1.88
    32.0 2.12
    37.0 2.33
    42.0 2.47
    48.4 1.35
      注:断面总面积为94.9 m2

    表  7   普适性验证实验的流速仪测量结果

    Table  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.93 10.10 9.39
    7 1.33
    7~12 1.60 8.97 14.40
    12 1.88
    12~17 1.94 8.59 16.70
    17 2.01
    17~22 2.14 8.86 19.00
    22 2.27
    22~27 2.03 9.18 18.60
    27 1.79
    27~32 1.72 10.00 17.20
    32 1.64
    32~37 1.68 11.10 18.60
    37 1.73
    37~42 1.58 12.30 19.40
    42 1.43
    42~51.4 1.00 15.80 15.80
    51.4 0

    表  8   流速仪方法、FD-DIS-G和DIS方法的对比

    Table  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.33 1.88 2.01 2.27 1.79 1.64 1.73 1.43 1.57 149 1.33 1.88
    DIS 0.37 0.95 1.32 1.65 1.67 1.65 1.39 1.39 1.19 113 0.37 0.95
    FD-DIS-G 0.93 1.48 1.82 1.76 1.85 1.99 1.54 1.90 1.53 145 0.93 1.48

    表  9   流速仪方法、FD-DIS-G和DIS方法的绝对误差分析

    Table  9   Analysis between Flow meter, FD-DIS-G and DIS

    方法 垂线平均流速绝对误差/(m·s–1) 平均流速绝对
    误差/(m·s–1)
    流量绝对
    误差/(m3·s–1)
    b=7 m b=12 m b=17 m b=22 m b=27 m b=32 m b=37 m b=42 m
    FD-DIS-G 0.40 0.40 0.19 0.51 0.06 0.35 0.19 0.47 0.04 4
    DIS 0.96 0.93 0.69 0.62 0.12 0.01 0.34 0.04 0.38 36

    表  10   流速仪方法、FD-DIS-G和DIS方法的相对误差分析

    Table  10   Analysis between Flow meter, FD-DIS-G and DIS

    方法 垂线平均流速相对误差/% 平均流速相
    对误差/%
    流量相对
    误差/%
    b=7 m b=12 m b=17 m b=22 m b=27 m b=32 m b=37 m b=42 m
    FD-DIS-G 30.1 21.4 9.5 22.3 3.1 21.5 11.1 32.9 2.5 2.7
    DIS 71.8 49.4 34.3 27.4 7.0 0.8 19.5 2.6 24.2 24.2
  • [1] 邱江潮,刘丙军,杨子博,等.基于水位流量关系的流量估算不确定性分析[J].水科学进展,2020,31(2):214–223. doi: 10.14042/j.cnki.32.1309.2020.02.007

    Qiu Jiangchao,Liu Bingjun,Yang Zibo,et al.Uncertainty analysis of estimated discharge based on stage-discharge rating curves[J].Advances in Water Science,2020,31(2):214–223 doi: 10.14042/j.cnki.32.1309.2020.02.007
    [2] 徐立中,张振,严锡君,等.非接触式明渠水流监测技术的发展现状[J].水利信息化,2013(3):37–44. doi: 10.19364/j.1674-9405.2013.03.010

    Xu Lizhong,Zhang Zhen,Yan Xijun,et al.Advances of non-contact instruments and techniques for open-channel flow measurements[J].Water Resources Informatization,2013(3):37–44 doi: 10.19364/j.1674-9405.2013.03.010
    [3] 王文娥,廖伟,漆力健.宽窄相间河道水流紊动特性试验研究[J].水科学进展,2020,31(3):394–403. doi: 10.14042/j.cnki.32.1309.2020.03.009

    Wang Wene,Liao Wei,Qi Lijian.Experiment of turbulent characteristics of flow in wide-and-narrow channels[J].Advances in Water Science,2020,31(3):394–403 doi: 10.14042/j.cnki.32.1309.2020.03.009
    [4] 张振,周扬,李旭睿,等.图像法测流系统开发与应用[J].水利信息化,2018(3):7–13. doi: 10.19364/j.1674-9405.2018.03.002

    Zhang Zhen,Zhou Yang,Li Xurui,et al.Development and application of an image-based flow measurement system[J].Water Resources Informatization,2018(3):7–13 doi: 10.19364/j.1674-9405.2018.03.002
    [5] 赵浩源,陈华,刘维高,等.基于河流表面时空图像识别的测流方法[J].水资源研究,2020(1):1–11.

    Zhao Haoyuan,Chen Hua,Liu Weigao,et al.Application of flow measurement method based on space-time image recognition of river surface[J].Journal of Water Resources Research,2020(1):1–11
    [6] 于长东,毕晓君,韩阳,等.基于轻量化深度学习模型的粒子图像测速研究[J].光学学报,2020,40(7):0720001. doi: 10.3788/AOS202040.0720001

    Yu Changdong,Bi Xiaojun,Han Yang,et al.Particle image velocimetry based on a lightweight deep learning model[J].Acta Optica Sinica,2020,40(7):0720001 doi: 10.3788/AOS202040.0720001
    [7] 胡文斌,马志敏,田猛,等.多光谱成像的粒子图像测速[J].光谱学与光谱分析,2018,38(7):2038–2043. doi: 10.3964/j.issn.1000-0593(2018)07-2038-06

    Hu Wenbin,Ma Zhimin,Tian Meng,et al.Multispectral-based particle image velocimetry[J].Spectroscopy and Spectral Analysis,2018,38(7):2038–2043 doi: 10.3964/j.issn.1000-0593(2018)07-2038-06
    [8] Adrian R J.Multi-point optical measurements of simultaneous vectors in unsteady flow—A review[J].International Journal of Heat and Fluid Flow,1986,7(2):127–145. doi: 10.1016/0142-727X(86)90062-7[LinkOut
    [9] Adrian R J.Particle-imaging techniques for experimental fluid mechanics[J].Annual Review of Fluid Mechanics,1991,23:261–304. doi: 10.1146/annurev.fl.23.010191.001401[LinkOut
    [10] Fujita I,Muste M,Kruger A.Large-scale particle image velocimetry for flow analysis in hydraulic engineering applications[J].Journal of Hydraulic Research,1998,36(3):397–414. doi: 10.1080/00221689809498626[LinkOut
    [11] 张振,徐枫,王鑫,等.河流水面成像测速研究进展[J].仪器仪表学报,2015,36(7):1441–1450. doi: 10.19650/j.cnki.cjsi.2015.07.001

    Zhang Zhen,Xu Feng,Wang Xin,et al.Research progress on river surface imaging velocimetry[J].Chinese Journal of Scientific Instrument,2015,36(7):1441–1450 doi: 10.19650/j.cnki.cjsi.2015.07.001
    [12] Lewis Q W,Lindroth E M,Rhoads B L.Integrating unmanned aerial systems and LSPIV for rapid,cost-effective stream gauging[J].Journal of Hydrology,2018,560:230–246. doi: 10.1016/j.jhydrol.2018.03.008[LinkOut
    [13] Akbarpour F,Fathi-Moghadam M,Schneider J.Application of LSPIV to measure supercritical flow in steep channels with low relative submergence[J].Flow Measurement and Instrumentation,2020,72:101718. doi: 10.1016/j.flowmeasinst.2020.101718[LinkOut
    [14] Watanabe K,Fujita I,Iguchi M,et al.Improving accuracy and robustness of space-time image velocimetry (STIV) with deep learning[J].Water,2021,13(15):2079. doi: 10.3390/w13152079[LinkOut
    [15] Fujita I,Notoya Y,Tani K,et al.Efficient and accurate estimation of water surface velocity in STIV[J].Environmental Fluid Mechanics,2019,19(5):1363–1378. doi: 10.1007/s10652-018-9651-3[LinkOut
    [16] Fujita I,Watanabe H,Tsubaki R.Development of a non-intrusive and efficient flow monitoring technique:The space-time image velocimetry (STIV)[J].International Journal of River Basin Management,2007,5(2):105–114. doi: 10.1080/15715124.2007.9635310[LinkOut
    [17] 王万良,杨胜兰,赵燕伟,等.基于条件边界平衡生成对抗网络的河流表面流速估测[J].浙江大学学报(工学版),2019,53(11):2118–2128. doi: 10.3785/j.issn.1008-973X.2019.11.009

    Wang Wanliang,Yang Shenglan,Zhao Yanwei,et al.Estimation of river surface flow velocity based on conditional boundary equilibrium generative adversarial network[J].Journal of Zhejiang University (Engineering Science),2019,53(11):2118–2128 doi: 10.3785/j.issn.1008-973X.2019.11.009
    [18] Lin D,Grundmann J,Eltner A.Evaluating image tracking approaches for surface velocimetry with thermal tracers[J].Water Resources Research,2019,55(4):3122–3136. doi: 10.1029/2018wr024507[LinkOut
    [19] Bouguet J.Pyramidal implementation of the Lucas kanade feature tracker.Intel Corporation,Microprocessor Research Labs.[EB/OL].[2021–10–01].http://robots.stanford. edu/cs223b04/algo_affine_tracking.pdf
    [20] Tauro F,Tosi F,Mattoccia S,et al.Optical tracking velocimetry (OTV):Leveraging optical flow and trajectory-based filtering for surface streamflow observations[J].Remote Sensing,2018,10(12):2010. doi: 10.3390/rs10122010[LinkOut
    [21] Kroeger T,Timofte R,Dai Dengxin,et al.Fast optical flow using dense inverse search[C]//European Conference on Computer Vision.Cham:Springer,2016:471–488.
    [22] 冯全,张彦洪,赵晓刚.基于机器视觉的河水表面流速估计[J].农业工程学报,2018,34(19):140–146. doi: 10.11975/j.issn.1002-6819.2018.19.018

    Feng Quan,Zhang Yanhong,Zhao Xiaogang.Estimation of surface speed of river flow based on machine vision[J].Transactions of the Chinese Society of Agricultural Engineering,2018,34(19):140–146 doi: 10.11975/j.issn.1002-6819.2018.19.018
    [23] Baker S,Matthews I.Lucas-kanade 20 years on:A unifying framework[J].International Journal of Computer Vision,2004,56(3):221–255. doi: 10.1023/b:visi.0000011205.11775.fd[LinkOut
    [24] 吴石林,张玘.误差分析与数据处理[M].北京:清华大学出版社,2010.
    [25] 拉斐尔·C.冈萨雷斯,理查德·E.伍兹.数字图像处理[M].阮秋琦,阮宇智.2版.北京:电子工业出版社,2011.
    [26] 中华人民共和国住房和城乡建设部.河流流量测验规范:GB 50179—2015[S].北京:中国计划出版社,2016.
    [27] 王俊,刘东生,陈松生.河流流量测验误差的理论与实践[M].武汉:长江出版社,2017.
图(11)  /  表(10)

本文结构

    /

    返回文章
    返回