ABSTRACT

Space-time image velocimetry (STIV) is a time-averaged velocity measurement method which takes river surface images as the object of analysis, and detects the texture orientation of a generated space-time image to obtain the one-dimensional flow velocity. STIV method has special application potential in real-time monitoring of river surface velocity and cross-section discharge since it has the advantages of high spatial resolution, strong real-time performance, safety and convenience. Because this method is still in its infancy in China and the shortcomings of the current flow measurement methods, this study attempts to introduce STIV method to better serve the domestic hydrological measurement work. In order to verify the feasibility of the method and propose research needs, a comparative experiment was conducted in the Chongyang hydrological station with the impellor-type current meter. The result shows that the relative error of the flow velocity of STIV method was less than 16%. The feasibility of this method is preliminarily verified in this experiment.

Keywords:Flow Measurement, Image Processing, Spatial-Temporal Image, Comparative Experiment, STIV

基于河流表面时空图像识别的测流方法

赵浩源，陈华^*，刘维高，黄凯霖，刘炳义，刘德地，王俊

武汉大学水资源与水电工程科学国家重点实验室，湖北 武汉

收稿日期：2019年12月7日；录用日期：2019年12月20日；发布日期：2019年12月27日

摘 要

时空图像测速方法(STIV)是一种以河流表面图像为分析对象、通过检测合成时空图像的纹理主方向获取一维时均流速的测量方法。由于具有空间分辨率高、实时性强、安全简便等优点，STIV方法在河流水面流速、断面流量的实时监测中具有特别的应用潜力。鉴于该方法在国内的研究尚处于起步阶段以及现有测流方法的缺点，本研究尝试引进STIV方法以更好地服务于国内水文测验工作。为了验证该方法的可行性并提出研究需求，在崇阳水文站与旋桨式流速仪进行比测试验。结果表明，STIV方法所获取的表面流速的相对误差均控制在16%以内，初步说明了该方法的可行性。

关键词 :流速测量，图像处理，时空图像，比测试验，STIV

Copyright © 2020 by author(s) and Wuhan University.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

流量数据是水文研究最为基础的数据之一，也是一切水文研究的基础。由于天然河流复杂的流动特性以及恶劣的户外测验环境，天然河流的流量测量一直是水文测验中一项困难的工作。尤其是在当前洪水与干旱等极端情况频发的情况下，对于及时获取河流流速、径流量变化等关键信息的需求更为迫切，传统的流量测量方法由于自身种种缺点已经难以满足要求。

常见的流量测量方法包括流速仪法、浮标法以及多普勒流速剖面仪法(ADCP)等。流速仪法测流首先使用流速仪测量垂线的流速分布，再通过流速面积法计算流量。该方法是最基础也是目前应用最为广泛的测流方法，但是其费时费力。例如对于一个宽约300 m，设置有10条测速线的断面，即使采用最简单的一点法施测，一个测回也需要40~50分钟，而对于类似长江这样的大江大河，一个测回所需时间可达5个小时，甚至更久。此外，流速仪法属于接触式测流方法，难以保证水文测验人员自身的生命安全，尤其是在高洪时期；浮标法测流是另一种常用的测流方法，该方法通过测量出断面上水面浮标速度分布，再结合断面数据确定虚流量，最后乘以浮标系数得到断面流量 [1]。但是该方法在测量浮标流速、浮标系数的选取方面存在诸多不确定性，且该方法不适宜在发生特大洪水的情况下使用。正是由于上述缺点，该方法目前已很少采用；走航式ADCP于上世纪90年代初被长江水利委员会水文局率先引入到国内，并逐渐发展成为一种重要的测流方法。ADCP通过跟踪水体中散射体的运动来得到水流速度，与流速仪不同的是，ADCP可以一次测量垂线上的多点流速，再使用流速面积法得到断面流量 [2]。也正是由于ADCP的测量原理，导致了其在河流含沙量较高等情况下容易出现测量结果不稳定甚至严重失真的情况，此外高昂的价格也限制了ADCP在更大范围的推广使用。

伴随着信息技术的快速发展，基于图像、光学以及雷达的测流方法迅速发展，特别是基于图像的测流方法，因其简单、高效、安全而得到广泛关注。极坐标摄影浮标法是最早提出的基于图像的测流方法，该方法实质上是对传统浮标法的一种改进 [3]。其基本思想是采用摄影照片代替现场测量，通过人工标识或计算机识别的方式从照片中提取测流浮标的方位，进而求出运动轨迹上的水面流速。卫星遥感图像法是另一种基于图像的测量方法，随着高分辨率成像及合成孔径雷达(SAR)遥感技术的发展，卫星遥感图像被研究用于获取地面河流的形态及水位信息，并结合地面测量值和水文模型估计河流流量。该方法大致可分为利用河宽估计流量、利用水位估计流量以及利用多变量估计流量三种。目前卫星遥感图像法测流主要用于宽浅河流的流量估计，由于对地面信息和历史数据的依赖及过大的测量误差使之尚无法实用化，但可以为洪水、湿地、泥石流、堰塞湖等难以到达地区的灾害应急监测提供及时的先验信息 [4]。近十几年来，另一种通过分析河流表面图像变化的测流方法受到了广泛关注。该方法主要在河岸一侧架设相机，利用相机拍摄的水面视频图像来分析水流示踪物运动矢量的大小。作为一种非接触式的测流方法，它克服了传统接触式测流方法的缺点，可以简单、快速、高效、安全地测量流速和流量信息，尤其使得高洪时期的测量工作成为可能。此外，鉴于目前各水文站均配备有河流监控系统，这无疑为该方法的推广使用提供了良好的条件，同时也降低了使用成本。根据运动矢量估计方法的不同，该方法主要分为大尺度粒子图像测速法(LSPIV)和时空图像测速法(STIV)两类。LSPIV方法于20世纪90年代由Fujita等人提出 [5]，该方法是将实验室流体力学研究中的粒子图像测速(PIV)技术改进用于现场河流的水面流场观测及流量测量。它以天然漂浮物及涟漪等表面流动特征为水流示踪物，以自然光为主要光源，以数码相机或视频摄像机为图像采集设备，通过匹配分析区域内跟随表面水流运动的示踪物图像获得二维流速场 [6]。在国内，河海大学课题组以LSPIV为基本原理，建立了一套河流水面成像测速(RSIV)体系并开展了比测试验 [7]。然而，LSPIV也存在一些缺点：1) 测量依赖于水面示踪物的可见性，当采用天然示踪物时，由于示踪物密度低、时空分布不均，以及水面成像条件复杂等的影响，极易引起时均流场重建的失效或误差过大 [6] ；2) 分析区域大小的选择存在诸多不确定性。若分析区域尺寸过大，则会由于空间平均效应降低测量的空间分辨率 [8]。分析区域过小又容易因缺乏目标信息而导致出现误匹配 [9] ；3)所需存储空间大，计算效率较低 [10]。相比于LSPIV，Fujita等人于2007年提出的STIV方法在断面流量测量方面展现出更好的性能 [10]。尽管STIV方法仅能获得平行于流动方向的一维流速分布，但其空间分辨率能够达到单像素水平，并且算法效率是LSPIV的10倍以上，这已经足够满足实际流量测量工作的需求。Fujita等人先后在日本的多条河流使用该方法进行流量测量试验，均取得了较好效果 [10] [11] [12]。此后针对STIV方法拍摄河宽受限的问题，Fujita等成功使用不同拍摄俯角来加以解决 [13]，又使用高分辨率的红外线摄像机完成了基于STIV的夜间流量观测 [14]。目前对于STIV的研究与应用主要集中在日本，而在国内尚处于起步阶段。

针对STIV方法在测流方面的优良特性以及现有测流方法的缺点，本研究尝试引进STIV方法以更好地服务于国内水文测验工作。为了评价该方法的可行性并提出研究需求，于2019年6月14日在湖北省咸宁市崇阳水文站与旋桨式流速仪进行了流速比测试验。

2. STIV方法介绍

2.1. 原理简介

在忽略风的影响下，诸如涟漪、波纹等河流表面流动特征是随水流一起运动的，因此可以认为其运动速度近似等于河流表面流速，这些表面流动特征的运动又会导致河流表面灰度发生变化，综上所述，河流表面灰度的变化大小可以反映河流表面流速的大小。我们可以在所拍摄的河流视频中沿水流方向设置一系列测速线，逐帧提取每条测速线的灰度信息以合成该测速线的时空图像。由于灰度的变化，在每幅时空图像中会呈现出带状纹理，带状纹理与竖直方向所夹的角度(纹理角)即反映了表面流速信息。

2.2. 时空图像生成

时空图像包含着河流表面灰度的变化信息，通过对时空图像中纹理角的识别可以得到表面流速信息。下面以一个实例来说明时空图像的生成：图1是现场拍摄的视频图像，图中红色箭头代表水流方向，红线代表沿水流方向设置的测速线，测速线长度为200像素，宽度为1像素。根据测速线的像素坐标位置，逐帧提取测速线的灰度，并按照从上往下的顺序进行排列，即可得到如图2所示的时空图像。

如图2所示，时空图像的横坐标代表测速线的长度，纵坐标代表视频拍摄时间。例如，在上面的例子中，测速线长度为200像素，则时空图像的横向宽度为200像素；视频的持续时间为15 s，每秒30帧，则时空图像的纵向宽度为450像素。

图1. 测速线设置

图2. 测速线的时空图像

2.3. 纹理角识别

在已经得到时空图像的情况下，需要采用一种合适的图像识别方法来检测时空图像中纹理角的大小。在目前的STIV方法中，选取灰度梯度张量法来对纹理角α进行角度识别，计算公式如下 [15] ：

$\tan 2\alpha =\frac{2G_{xt}}{G_{tt}-G_{xx}}$ (1)

在(1)式中：

$G_{xx}={\displaystyle {\iint }_S\frac{\partial I\left(x,t\right)}{\partial x}\frac{\partial I\left(x,t\right)}{\partial x}\text{d}x\text{d}x}$ (2)

$G_{xt}={\displaystyle {\iint }_S\frac{\partial I\left(x,t\right)}{\partial x}\frac{\partial I\left(x,t\right)}{\partial t}\text{d}x\text{d}t}$ (3)

$G_{tt}={\displaystyle {\iint }_S\frac{\partial I\left(x,t\right)}{\partial t}}\frac{\partial I\left(x,t\right)}{\partial t}\text{d}t\text{d}t$ (4)

式中：I(x, t)表示时空图像中灰度的大小，S表示积分区域。 $\frac{\partial I\left(x,t\right)}{\partial x}$ 与 $\frac{\partial I\left(x,t\right)}{\partial t}$ 表示灰度的梯度，反映灰度沿x方向与t方向的变化情况，在本研究中采用Sobel算子H_x与H_t求解灰度梯度：

$H_x=\left[\begin{array}{ccc}-3& 0& 3\\ -10& 0& 10\\ -3& 0& 3\end{array}\right]，H_t=\left[\begin{array}{ccc}-3& -10& -3\\ 0& 0& 0\\ 3& 10& 3\end{array}\right]$ (5)

在这里引入一个相关性参数C，以检查时空图像纹理的清晰程度：

$C=\frac{\sqrt{{\left(G_{xx}-G_{tt}\right)}^2+4G_{xt}^2}}{G_{xx}+G_{tt}}$ (6)

参数C的取值范围为0至1，参数C的值越大，表明时空图像的纹理越清晰。

为了尽可能排除时空图像中噪声对计算结果的干扰，我们将时空图像分成若干个小部分，先由(1)~(5)式求出各部分角度，再对各部分的角度值求平均，得到最终的纹理角大小，过程如图3所示。

用以下公式求解最终的纹理角大小：

$\bar{\alpha }=\frac{{\displaystyle \underset{i}{\sum }{\alpha }_i\cdot C\left({\alpha }_i\right)}}{{\displaystyle \underset{i}{\sum }C\left({\alpha }_i\right)}}$ (7)

式(7)的思想是：将各部分的相关性 $C\left({\alpha }_i\right)$ 作为各部分角度值的权重，对各部分加权平均求最终的纹理角 $\bar{\alpha }$。这样，条纹更加清晰的部分将获得更大权重，而因噪声影响使得条纹不清晰的部分权重较小，使得求解结果更加准确。

图3. 纹理角求解过程

2.4. 相机标定与测速线长度求解

得到时空图像的纹理角后，根据测速线的实际长度即可求得测速线的流速。测速线的长度可以根据相平面坐标(x, y)与实际空间直角坐标(X, Y, Z)的关系求得 [16] ：

$\{\begin{array}{l}x=x_p+f\frac{r_{11}\left(X-X_p\right)+r_{12}\left(Y-Y_p\right)+r_{13}\left(Z-Z_p\right)}{r_{31}\left(X-X_p\right)+r_{32}\left(Y-Y_p\right)+r_{33}\left(Z-Z_p\right)}+\Delta x\hfill \\ y=y_p+f\frac{r_{21}\left(X-X_p\right)+r_{22}\left(Y-Y_p\right)+r_{23}\left(Z-Z_p\right)}{r_{31}\left(X-X_p\right)+r_{32}\left(Y-Y_p\right)+r_{33}\left(Z-Z_p\right)}+\Delta y\hfill \end{array}$ (8)

式中，(x_p , y_p)表示像主点的相平面坐标，(X_p, Y_p, Z_p)表示相机的实际空间直角坐标，f代表焦距，(Δx, Δy)代表镜头畸变校正数，r_ij (i, j = 1~3)是两个坐标的转换系数。

式(8)中含有12个未知参数，因此我们至少需要知道6个点的相平面坐标与实际空间坐标才能求解上述未知数(每个点可列出2个方程)，为此我们需要在现场设定6个或者6个以上地面标定牌，用全站仪测出它们相对于测站点的实际空间直角坐标，其相平面坐标则可直接从视频图像中获得。在求得上述12个未知数后，我们就可以根据(8)式并结合水位值算出时空图像上任意一点的实际空间直角坐标，从而求得所设置测速线的实际长度。

2.5. 表面流速与流量计算

我们假设在实际空间直角坐标系中，涟漪等表面流动特征在时间T内沿测速线运动的距离为L，相应地，在相平面坐标系下则对应于在k帧内运动了i个像素，则测速线上的流速大小为：

$v=\frac{L}{T}=\frac{i\cdot S_x}{k\cdot S_t}=\tan \bar{\alpha }\cdot \frac{S_x}{S_t}=\tan \bar{\alpha }\cdot S_x\cdot \text{fps}$ (9)

式中： $\bar{\alpha }$ 表示纹理角大小，S_x表示每个像素所代表的实际距离(单位为：m/像素)，fps表示相机的帧速率(单位为：帧/s)。

最后可根据流速面积法，按照下面的公式计算流量：

$Q={\displaystyle \underset{i}{\sum }\eta \cdot V_i\cdot A_i}$ (10)

式中： $\eta$ 表示表面流速系数，V_i表示区间表面流速，A_i表示区间面积。

3. 比测试验

3.1. 试验介绍

目前在水文测验中，将流速仪法测量得到的流量结果视为流量真值，因此流速仪法是比较其他测流方法的标准 [2]。为了验证上述基于图像识别的测流方法的可行性与准确性，于2019年6月14日在崇阳水文站与旋桨式流速仪进行了流速比测试验。流速仪测量与视频拍摄在无风时段同步进行，以满足STIV使用条件且保证二者测量的是同一时段水流。流速仪在起点距分别为70 m，80 m，90 m，100 m，110 m，130 m和140 m的垂线利用两点法测量垂线流速，同时测量位于水面下50 cm处的流速，将其近似视为表面流速。根据视频中流速仪的位置，在相同起点距处设置长度相等的测速线，按照第2部分的方法计算表面流速与流量，并比较二者测量结果。

3.2. 试验布置与设备

崇阳水文站位于湖北省咸宁市崇阳县陆水河干流畔，隶属于长江水利委员会水文局，是国家级基本水文站。测流实验地面标定布置如图4中红框所示。

图4. 地面标定点布设

试验器材如表1所示。

表1. 试验器材清单

3.3. 试验数据

试验采集的7个地面标定点的实际空间直角坐标与像平面坐标如表2所示。

表2. 标定点坐标

由流速仪测量得到的表面流速与垂线平均流速如表3所示，其中垂线平均流速是由两点法算出。

表3. 流速仪测量结果

3.4. 试验结果

3.4.1. 表面流速

图像法测流的计算结果如表4所示。

表4. 图像法测流结果

3.4.2. 断面流量

如图5所示，利用流速面积法求解断面流量(岸边流速系数取0.7)，流速仪法流量计算表如表5所示。

表5. 流速仪法流量计算表

由表5可知，通过流速仪法计算得到的断面流量为59.78 m³/s。

图5. 断面示意图

图6. 流量计算结果

根据崇阳站多年观测经验，分别取表面流速系数为0.70~0.75，按照流速面积法计算断面流量。以取表面流速系数为0.7为例，断面示意图如图5所示，图像法流量计算表如表6所示。

表6. 图像法流量计算表(表面流速系数为0.7)

由表6可知，当表面流速系数取0.7时，通过图像法计算得到的断面流量为57.87 m³/s。按照同样的方法计算取其他流速系数时得到的断面流量，结果如图6所示。

3.5. 结果讨论

整理上述计算结果，可得表7。

表7. 计算结果

由表7可以看出，由图像法计算得到的表面流速均大于流速仪法实测值，相对误差基本在15%以内。事实上，为了消除水流脉动等因素的影响，流速仪在测量水面一点流速时，转子的实际入水深度至少要在水下50 cm以上才能保证测量值的稳定，因此流速仪测得的并不是真实的表面流速，而是水下50 cm左右处的流速，只是将其近似认为是表面流速，所以由图像法得到的表面流速大于流速仪实测值。从流量的计算结果来看，如图6所示，通过图像法得到的流量均在合理范围内，相对误差均不超过±5%，说明只要根据水文站特性适当选取表面流速系数即可得到令人满意的流量计算结果。

4. 结语

针对当前常见流量测量方法存在的问题，本文尝试引进STIV方法以更好地服务于流量测量工作。该方法认为在忽略风的影响下，诸如涟漪等河流表面流动特征是随水流一起运动的，因此可以认为其运动速度近似等于河流表面流速，这些表面流动特征的运动又会导致河流表面灰度发生变化，在逐帧连续提取的时空图像中形成条状纹理。纹理角度的大小就反映了表面流速的大小，所以可以通过识别纹理角的大小来计算表面流速。为了验证该方法的可行性，于崇阳水文站与旋桨式流速仪进行了比测试验，结果表明：图像法获得的表面流速的相对误差均控制在16%以内，由表面流速计算得到的流量也在合理范围内，初步说明了该方法的可行性。

未来的工作将结合现有结果对STIV方法进行完善，包括：1) 考虑风速对STIV方法的影响；2) 完善图像处理方法，提高示踪物特征，以增强STIV方法在诸如下雨、光照变化、倒影等复杂水面成像条件下的适应性；3) 对测速线长度、视频帧数等因素进行敏感性分析及不确定度评估，以便为STIV方法的应用提供更好指导。

致谢

感谢长江水利委员会水文局对该研究的支持，为该研究提供试验场地与试验数据。

基金项目

感谢国家自然科学基金重点项目51539009对本研究的支持。

文章引用

赵浩源,陈 华,刘维高,黄凯霖,刘炳义,刘德地,王 俊. 基于河流表面时空图像识别的测流方法
Application of Flow Measurement Method Based on Space-Time Image Recognition of River Surface[J]. 水资源研究, 2020, 09(01): 1-11. https://doi.org/10.12677/JWRR.2020.91001

参考文献