CN102564508B - 基于视频图像的河流流量在线测验实现方法 - Google Patents

Abstract

本发明公开了一种基于视频图像的河流流量在线测验实现方法，属于非接触式明渠测流技术领域。该方法中涉及水流示踪系统、水面照明系统、视频流采集系统、摄像机标定系统、图像处理系统和流量估计系统六个部分，水流示踪系统在测流河段投放示踪粒子，摄像机拍摄水面视频流，然后按固定时间间隔取帧并数字化处理，图像处理系统根据摄像机标定系统提供的地面控制点坐标信息对数字图像进行处理，最后流量估计系统采用流速-面积法计算出河流流量。本发明在降低河流测流复杂度的同时能保证较高的测量精度，可实现极端环境下河流流量的动态监测。

Description

基于视频图像的河流流量在线测验实现方法

技术领域

本发明涉及一种河流流量测验方法，尤其涉及一种基于视频图像的河流流量在线测验实现方法，属于非接触式明渠测流技术领域。

背景技术

随着工农业生产和生活用水量的增加，水资源供需矛盾日益尖锐，导致河道断流现象及水污染灾害频发，已对沿河地区的经济、居民生活和生态环境等产生诸多不利影响。同时，污水流量的长期监测也是一项艰巨的任务。由于排污渠道环境恶劣，其湿度高、杂物多，并可能具有腐蚀性的气体和液体，会给设备和操作人员带来危险。因此，要求流量测验仪器能够以最小的维护投入进行长期的自动测量。对于平原灌区，在整个渠系的排水过程中，若精确了解到某时某处的流量，可为整个渠系的计划用水和水量调配等工作提供量化参考，长期观测也可为水利规划设计和科学研究提供数据积累，因此如何迅速、准确地测得河流流量显得至关重要。

目前用于河流流量在线测验的方法可分为以下两大类：

(1)建筑物测流法，采用堰槽量水设备(如巴歇尔槽)，让水流通过束缩的过水断面形成缓流到急流的过渡即临界流，保证稳定的水位-流量关系，然后采用非接触式水位计(如超声波水位计)测出水位，根据率定曲线计算出动态流量。该方法易于实现自动化，但需对已有河道进行改造，投资大；并且，率定曲线是一种全局的基于模型的曲线，不能够捕获局部的水力影响，对于非均匀流及紊流通常不能给出精确的结果。因此只适用于少数小型的人工渠道。

(2)流速-面积法，以测速垂线为界将河流断面划分为几个子区域，根据实测流速及水深计算各子区域的断面面积、平均流速及流量，然后将各个子区域的流量求和得到总的断面流量。根据传感器相对水流的位置关系，流速测量方法又可分为以下三类：

a)接触式。其中基于转子式流速仪的方法被认为是河流流量测验的标准方法。但由于流速仪测得的是点流速，历时长，需配合横跨河道的桥梁或缆道系统才能实现自动测验，因此采用铅鱼流速仪所测流速值的空间和时间分辨率都较差。此外，对于污水渠道，这种流速仪必须放置于污水内，因此不能可靠地用于长期监测。

b)微侵入式。典型的是目前较为先进的声学多普勒流速剖面仪法(ADCP)。座底式或水平式ADCP可在短时内直接测得断面子区域中各水层的平均流速，实现河流流量在线测验，但存在盲区且其对含沙量较为敏感，通常只适用于形状规则的人工渠道；此外，由于换能器必须浸没于水中，易受高速水流及漂浮物的冲击并受河床冲淤影响，导致仪器运行的不可靠甚至损毁。

c)非接触式。比较成熟的是电波流速仪法，它测量的是位于波束与水面交界处的一维表面流速，空间分辨率有限，并且当风或湍流在水面产生细小的波动时才能产生足够强的反射回波，因此不适用于水面平静的低流速河流。

可见，现有的流量测验手段在原理上存在先天不足，受现场条件下流态和环境复杂度的影响以及技术和经济因素的限制，使目前河流流量的实时监测中存在着许多不精确和不确定性问题。尤其是在洪水、决堤、大坝合龙等极端高流速环境下，以及针对流速极低而过水面积巨大的大尺度河流流量测验，目前还没有行之有效的测验手段。

发明内容

本发明针对现有河流测流技术存在的不足，而提出一种基于视频图像的河流流量在线测验实现方法。

在该方法中设置有如下硬件系统：

a.用于实现水面流动可视化的水流示踪系统；

b.用于拍摄和存储水面实时视频流的视频流采集系统；

c.用于给图像正射变换提供地面控制点及其坐标信息的摄像机标定系统；

d.用于处理数字图像的图像处理系统；

e.用于计算河流流量的流量估计系统；

该方法分为三个阶段：

(一)测站部署阶段，包括如下步骤：

(1)选择测流河段；

(2)获取测流断面水下地形；

(3)部署现场硬件系统；

(4)勘测地面控制点在世界坐标系下的平面坐标和高程坐标；

(5)计算用于图像畸变校正的DLT变换系数；

(二)流场测量阶段，包括如下步骤：

(1)启动水流示踪系统投放示踪粒子；

(2)启动视频流采集系统拍摄包含示踪粒子的水面视频流，然后按固定的时间间隔取帧并数字化处理；

(3)由图像处理系统处理数字图像：包括图像畸变校正、图像增强、流速矢量提取、流速场后处理，最后得到修正后的水面瞬时流速场；

(三)流量估计阶段，该阶段是先勘测测流断面水位值，然后流量估计系统采用流速-面积法计算出河流流量。

技术效果：

1.时空分辨率高。仪器测得的是水面瞬时二维流速矢量场，包含了水面流速分布的细节信息，而非传统仪器测量的点流速或线流速；本方法能在数秒内完成视频图像采集、流速场获取和流量估计，并能在常规和极端水文过程期间对难以到达的测点实施密集测量。

2.测量值可靠。当缺乏水面天然示踪物时，通过播撒人工示踪粒子的方式实现流动可视化，相比传统浮标法测流，粒子的跟随性更好，测量值更可靠。

3.抗毁性强。仪器不接触水体，不会因为高速水流的冲击和漂浮物的缠绕而损毁，维护成本低，可进行远程控制和长期监测。

4.操作安全。仪器可布设于岸边、桥上或空中自治运行，工作人员可通过计算机网络远程操控，无需涉水测量，极大保障了测流人员的生命安全。

5.原理直观。相比声学和雷达测流技术，基于视频图像的方法更加直观明了；数字图像易于理解、处理、存储和传输，因此还可用于系统自身的故障诊断以及现场的工情监测。

6.成本低廉，配置灵活。仪器采用了市面上成熟且通用的硬件产品，如CCTV摄像机、LED补光灯、图像采集卡和PC机等，系统易于获取和搭建，可根据具体应用灵活配置，并且随着新兴成像技术和图像处理技术的发展及应用，仪器的性能可得到进一步提升。

7.可测量极低流速。非接触式的测流方式不会对水流造成干扰，适用于浅水、宽断面、低流速的河流。

8.可测量极高流速。视频图像具有较高的分辨率和帧速率，覆盖了从几平方米到几百平方米的水面，可捕捉高速流动的水面目标。

附图说明

图1是本发明方法涉及的系统架构示意图，图中标号名称：1、水流示踪系统；1-1、缆道；1-2、物料斗；1-3、示踪粒子；2、水面照明系统；3、视频流采集系统中的摄像机；4、水面；5、河岸；6、摄像机标定系统中的地面控制点；7表示时间平均水面流速；8表示水面插值流速；9表示深度平均流速；10、流量估计(测流)断面；11、勘测的水下地形。

图2是本发明方法的流程图。

图3是基于直接线性变换的图像畸变校正算法示意图。

图4是基于时域中值滤波的背景差分图像增强算法示意图。

图5是基于相关跟踪的流速矢量提取算法示意图。

图6是基于空域中值滤波的流速场后处理算法示意图。

具体实施方式

下面对本发明方法作进一步说明。

本发明方法中涉及的硬件系统架构如图1所示，主要是由水流示踪系统、水面照明系统、视频流采集系统、摄像机标定系统、图像处理系统和流量估计系统六个部分组成，其中前四个部署在测验现场，后两个为监控中心计算机运行的处理和计算系统。

下面对上述六个系统进行分别介绍：

1.水流示踪系统，该系统由示踪粒子和粒子投放装置组成。理想条件下，标识流动的粒子应当具备以下两个条件：一是要能够精确跟随局部流体的运动，惯性及浮力是反应其是否符合流动特性的主要因素；二是要足够大使其能够被视频流采集系统的成像器件检测到。尽管在河流水面上可能存在一些天然示踪物，如轻质漂浮物、泡沫、冰块、悬浮质及细小的波纹，但它们的出现是随机且不稳定的，在低流速条件下的跟随性和可见性也参差不齐。本方法采用投放人工粒子的方式以在天然示踪物不足时进行补充，可获得比较理想的效果。具体选用一种廉价的用作包装填充材料的缓冲粒(又称泡泡粒)产品作为人工示踪粒子，它是由99％的玉米淀粉制成，是一种生物可降解的环保材料。粒子投放于水中数分钟后就会自行分解，不会对水体和生态环境造成污染；吸水后其密度与待测水流的密度相近，能满足跟随性的要求；其部分浸没于水中可避免风力干扰引起的复杂度；并且白色的粒子能够与深色的水体形成对比，便于检测。粒子投放装置包括悬索式缆道和物料斗，物料斗悬挂于悬索上，其经缆道控制在水面上方做往复运动，内盛的示踪粒子从底部漏出播撒于水面。

2.水面照明系统，该系统由光源及光强传感器组成。光源采用高亮卤素灯或LED补光灯，功率和发射角度应当和照射距离及镜头焦距相匹配，这样才能在视场内获得较为均匀的光照；光强传感器用于感知环境光的强弱，以实现控制照明系统自动开启或关闭。

3.视频流采集系统，该系统由测验现场的CCTV摄像机及监控中心的视频服务器组成。摄像机安装于云台上的防水盒内，可对其旋转角和俯仰角进行调节。摄像机拍摄PAL制式(720×576像素，25帧/秒)或NTSC制式(640×480像素，30帧/秒)的模拟视频流，通过同轴电缆上传至监控中心的视频服务器。视频服务器通过视频采集卡将模拟视频流进行数字化处理，一路经过压缩编码(MPEG-4或H.264格式)转换成数字视频流存储于硬盘阵列上，可由本地/远程的客户端访问，或通过视频解码盒解码转换为视频信号在测控屏上显示；另一路不经压缩直接存储为原始数字视频帧数据(RAW或BMP格式)，由图像服务器调用。

4.摄像机标定系统，该系统由至少六个地面控制点(GCP)组成，用于后期图像的几何校正。地面控制点是采用对测距激光具有强反射的铝合金材料制成的形状规则的十字形或正方形靶面。

5.图像处理系统，该系统运行于监控中心的图像服务器上。它首先从视频服务器上读取原始的数字视频帧数据，然后对图像进行畸变校正和图像增强的预处理，接下来采用模式识别的方法获得水面示踪物在无畸变连续视频帧中的位置，再将测得的位移除以连续帧的时间间隔得到每个分析区域的流速矢量，对视场中包含的各分析区域重复上述过程得到原始的瞬时流速场，最后利用邻域信息对流速矢量进行校正得到修正的瞬时流速场。

6.流量估计系统，该系统运行于监控中心的图像服务器上。它是利用水面流速场以及勘测的断面水下地形和水位值，基于流速-面积法估计出河流流量。图像服务器将得到的流速场及流量结果存储在本地硬盘上共客户端调用，并将信息叠加在水面的视频帧上通过VGA接口显示于测控屏。

本发明方法的流程如图2所示，总体分为测站部署阶段、流场测量阶段和流量估计阶段三个阶段，下面对每一阶段分别说明。

(一)测站部署阶段，包括如下步骤内容：

1、选择测流河段。

具体实施时按照以下条件选取测流河段：(a)河槽的底坡、断面形状、糙率等因素比较稳定，槽段内的断面规则、水流平顺、流速分布对称均匀；(b)河段内无巨大块石、丛生水草阻水，无巨大漩涡、乱流等现象，当天然河道达不到以上要求时，必须进行人工整治使其符合测流条件；(c)宜选择河面较窄、通视条件好、横断面较单一、受风浪影响较小的河段，这样有利于示踪粒子的投放，并能减小图像远场分辨率不足对测流精度的影响。有条件的测点还可利用桥梁、堰闸布置测验。

2、勘测断面水下地形。

人工渠道的断面结构通常是已知的，而天然河道的水下地形可通过专用的声学仪器(如多波束声纳或声学多普勒流速剖面仪)直接测得。河流的水下地形勘测也可在测流时或测流后进行，前提是在该时间段内河床和边坡不会发生显著的变化。

3、部署现场仪器系统。

水流示踪系统部署于紧邻测流断面的上游；水面照明系统部署于河流测流断面的两岸，光源可与视频流采集系统中的摄像机布置于同一位置，也可根据实际情况分散部署；视频流采集系统中的CCTV摄像机架设于岸边制高点(可以是高层建筑、电线塔或临时搭建的高杆)上或横跨断面的桥梁上，调整云台的角度使其光轴尽量垂直于水流方向，且成像平面尽量平行于水面；由于在水面上放置控制点较为困难，在布设摄像机标定系统时，可将地面控制点平行放置于待测河段两岸河面最高水位以上的某一水平面处，并且尽量均匀而分散地覆盖完整的摄像机视场区域。

4、勘测地面控制点坐标。

采用全站仪、差分全球定位系统(DGPS)或激光测距仪测量各地面控制点(GCP)在世界坐标系O_W-X_WY_WZ_W下的平面坐标

以及高程坐标

i＝1，L，N，这里N为GCP的数量。

5、计算变换系数。

拍摄一幅视频帧，提取出GCP的图像坐标，并基于直接线性变换的方法求解出用于图像畸变校正的DLT变换系数，供镜头畸变校正和图像正射校正使用。

至此完成测站的前期准备工作。

(二)流场测量阶段，包括如下步骤内容：

1、控制照明并投放粒子。

启动水面照明系统，光强传感器根据环境光强判断是否开启补光光源。启动粒子投放装置，往复运动的物料斗将内盛的示踪粒子均匀地播撒于水面，粒子群顺流而下进入待测区域。

2、拍摄视频流并取帧。

摄像机拍摄包含示踪粒子的水面视频流，由视频服务器取帧并数字化为8位灰度图像。

3、处理数字图像。

数字视频图像的处理流程如下：

(1)图像畸变校正。CCTV摄像机的光学系统一般为由多片透镜组成的光学镜头，按理想化的小孔成像原理建模时存在着模型误差，引起包括径向畸变、切向畸变和薄棱镜畸变在内的二维图像的非线性畸变。此外，为了使图像包含较大的水流区域(通常是一个完整的测验断面)，摄像机通常处于一个倾斜于流动平面的视角下进行拍摄，由此引入了图像的透视畸变。本方法采用了一种基于直接线性变换的图像畸变校正算法，首先通过建立理想针孔摄像机模型、平面间变高单应关系以及镜头非线性畸变模型，建立镜头畸变模型下的图像坐标和2D世界坐标间的映射关系，接着采用最小二乘法求解变换系数，包括DLT正变换系数、非线性畸变系数和DLT逆变换系数，最后进行镜头畸变校正和图像正射校正以获得真实世界中的度量信息，并用双线性内插法对整数坐标位置上的像素灰度进行插值，以便在图像上进行显示。

(2)图像增强。通过观察可发现，水面背景像素的分布通常是均衡的，即单模态的，并且示踪粒子相对于整个水面场景是一种运动连续的小目标。但相比传统实验室中使用的粒子图像测速技术，现场环境下的测量条件更为复杂：由于河流的开放性，不可避免地会在一定程度上受到外界环境的影响，如强风降雨在水面产生的扰动(动态)、阳光在水面产生的镜面反射(准静态)以及树木或建筑的倒影(静态)等。如果在图像预处理阶段不能合理地消除这些不利的水面模式，将会对后续的流速矢量提取环节造成影响，因此需要采用图像增强算法对水面灰度图像进行预处理。本方法结合水面图像的特点，设计了一种基于时域中值滤波的背景差分图像增强算法，背景差分可充分利用图像中变化部分的信息，使相关曲面的峰值锐化，而时域中值滤波背景估计的速度较快，保证所估计的背景像素可取到各帧输入图像中稳定时间最长的像素值。因此该算法能较好地抑制多种不利的水面模式，突出前景的运动粒子，从而达到图像增强的效果。

(3)流速矢量提取。流速矢量提取算法是本方法的关键。在粒子密度中等、模式细节丰富且分布较为均匀的情况下，河流水面示踪粒子及其它水面模式的相对位置在相邻帧间隔内不会发生显著变化。由此，本方法依据模式识别的原理进行目标跟踪，具体采用了一种基于相关跟踪的流速提取算法，它将归一化的互相关系数作为模板图像与匹配图像灰度分布的相似性度量，通过对相关曲面中的峰值位置进行搜索及亚像素定位，确定第一帧图像中的分析区域相对于第二帧图像中对应窗口的最大似然位移。由于该算法利用图像的统计特征而不是粒子的几何特征进行匹配，单个示踪粒子的尺寸并不重要，因此对图像质量要求不高，可在低信噪比的条件下正常工作，能适应变化较复杂的流体运动目标。

(4)流速场后处理。受示踪粒子浓度低和分布不均的影响以及外界环境的干扰，视场中的某些待测区域可能会暂时性地缺乏明显的示踪物，导致无法估计出这些区域(标记为NULL的分析区域)的瞬时表面流速场(矢量缺失)，或导致出现伪矢量(幅值或方向的偏差较大)。因此，需要采用一些准则以及全局或局部的滤波过程来减少错误矢量的个数，本方法采用空域中值滤波的方法修正伪矢量，并对缺失矢量进行插补。

至此获得修正后的水面瞬时流速场。

(三)流量估计阶段。

河流流量是单位时间内通过测量断面的水流体积，即流量＝断面平均流速×过水断面面积。因此，本方法利用河流水下地形、水面流速场、水下纵向流速分布及水位值，基于流速-面积法估计河流流量。流速-面积法的原理是先以测速垂线作为划分界线，将河流断面划分为若干个子断面，再根据实测流速及水深计算各子断面的过水面积、平均流速及流量，然后将各个子断面的流量求和得到总的断面流量。从本质上看它依赖于河流的断面流速分布，在本方法中，水面流速场反映了水面横向流速分布的规律，而建立水面流速与深度平均流速间关系的流速系数则反映了纵向流速分布的规律。流量估计方法流程如下：

1、勘测河流测流断面水位。在确认河流水下地形参数有效(即在测流期间断面的水下地形未发生显著变化)的基础上，通过超声波水位计或水尺勘测实时水位值l。

2、计算时间平均流速场。对于一段长度为N帧的视频序列共有N-1个图像对及其对应的瞬时流速场。对分析区域i的瞬时流速矢量V_i在时域进行信息融合，获得时间平均流速矢量U_i，i＝1，2，L，M，M为时间平均流速矢量总数，简单的融合方法可采用算术平均：

$U_i=\frac{\underset{N-1}{\mathrm{\Σ}}{V}_i}{N-1}$

然后，将所选勘测断面中分析区域i的时间平均流速矢量U_i投影到流动主方向X，得到标量的时间平均水面流速U_i。

3、计算水面插值流速。将测流断面划分为以n+1条测速垂线为边界的n个子断面，则子断面S_v的边界为垂线v和v+1。由于测速垂线v不一定对应于分析区域i，因此对垂线相邻分析区域的时间平均水面流速U_i进行线性插值，求得测速垂线的水面插值流速U′_v。

4、计算深度平均流速。用水面插值流速U′_v计算垂线的深度平均流速

假设每条测速垂线处纵向流速剖面的形状都相同，则深度平均流速其中α为流速系数，它反映了纵向流速分布的规律，其取值依赖于剖面形状和河床的特性，在对数律的纵向流速分布和典型河床糙率的条件下，通常取α＝0.85。

5、计算子断面流量。根据经典的流速-面积法计算子断面S_v的流量：

$D_v=A_v\left(l\right)\·{\stackrel{\‾}{U}}_{\mathrm{sv}}$

其中：是子断面的平均流速，假定其大小为子断面两条边界垂线深度平均流速的均值，即A_v(l)是水位值为l时子断面S_v的过水面积，

h_v和h_v+1是子断面S_v两条边界垂线处的水深，由水下地形及水位值l共同给出；w_v是子断面S_v的宽度，即子断面S_v两条边界垂线之间的水平距离。

6、计算断面流量。将各个子断面的流量求和得到河流的总断面流量：

$Q=\underset{v=1}{\overset{n}{\mathrm{\Σ}}}{D}_v=\underset{v=1}{\overset{n}{\mathrm{\Σ}}}{A}_v\left(l\right)\·{\stackrel{\‾}{U}}_{\mathrm{sv}}=\underset{v=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{2}(h_v+h_{v+1})\×\frac{1}{2}({\stackrel{\‾}{U}}_v+{\stackrel{\‾}{U}}_{v+1})w_v$

至此完成河流流量估计。

下面结合附图对本方法中涉及的部分算法展开说明。

基于直接线性变换的图像畸变校正算法如图3所示。图像畸变校正实质上是建立3D世界坐标系中的点在2D世界坐标系中的坐标与它在镜头畸变模型下的图像坐标之间的关系。校正原理建立在含畸变的针孔摄像机模型基础上，校正过程包括镜头畸变校正和图像正射校正，具体内容如下：

I、理想针孔摄像机模型的建立

CCTV摄像机的成像系统建立在理想针孔摄像机模型的基础上，模型和图像畸变校正算法共包含5种坐标系：反映真实物理空间的3D世界坐标系O_W-X_WY_WZ_W、反映摄像机方位的摄像机坐标系O_C-xyz、反映CCD/CMOS感光器件所在平面的像平面坐标系O_I-X_IY_I、反映采集的视频帧在计算机上显示的图像坐标系O_S-X_SY_S，以及正射校正后水面对应的2D世界坐标系O_w-X_wY_w。

其中，摄像机坐标系O_C-xyz的坐标原点为透镜的光学中心O_C，z轴与光轴重合；像平面坐标系O_I-X_IY_I的中心位于光轴与像平面的交点O_I，X_I、Y_I平行于x、y轴；图像坐标系O_S-X_SY_S的坐标原点O_S与2D世界坐标系O_w-X_wY_w的坐标原点O_w重合。

对于3D世界坐标系中的点P_W(X_W，Y_W，Z_W)，其在摄像机坐标系中的坐标为P_C(x，y，z)，在理想针孔模型下的像平面坐标为P_I(X_I，Y_I)，在镜头畸变模型下的像平面坐标为P_I(X_I，Y_I)，在理想针孔模型下的图像坐标为P_S(X_S，Y_S)，在镜头畸变模型下的图像坐标为P_S(x_S，y_S)，在理想针孔模型下的2D世界坐标为P_w(X_w，Y_w)，在镜头畸变模型下的2D世界坐标为P_w(x_w，y_w)。

II、平面间变高单应关系的建立

在理想针孔摄像机模型中，摄像机中心、图像点和目标点间存在共线的约束，图像和世界坐标系间的坐标空间变换可由共线方程推导得到如下线性关系：

$X_S=\frac{l_1{X}_W+l_2{Y}_W+l_3{Z}_W+l_4}{l_9{X}_W+l_{10}{Y}_W+l_{11}{Z}_W+1}$ (1)

$Y_S=\frac{l_5{X}_W+l_6{Y}_W+l_7{Z}_W+l_8}{l_9{X}_W+l_{10}{Y}_W+l_{11}{Z}_W+1}$

该变换即为直接线性变换(DLT)的正变换表达式，l_i称作DLT正变换系数，i＝1，2，K，11。通常情况下，可以把水面近似看作一个平面，Z_W＝0的直接线性变换描述了它和图像平面间的单应关系。

由于在水面上放置控制点较为困难，在布设标定系统时可将GCP平行放置于河面最高水位以上高程坐标为

的某一平面处，

可表示为：

$Z_W^G=d_1{X}_W^G+d_2{Y}_G^W+d_3---\left(2\right)$

其中，d₁、d₂分别为高程相对

方向变化的系数，d₃(t)是t时刻GCP所在平面与水平面间的高度差，水位值由自记式水位计测得。

假设河流表面的倾斜很小，则可将上式简化为

将其代入直接线性变换公式中可得：

$X_S^G=\frac{l_1{X}_W^G+l_2{Y}_W^G+l_3{d}_3\left(t\right)+l_4}{l_9{X}_W^G+l_{10}{Y}_W^G+l_{11}{d}_3\left(t\right)+1}$ (3)

$Y_S^G=\frac{l_5{X}_W^G+l_6{Y}_W^G+l_7{d}_3\left(t\right)+l_8}{l_9{X}_W^G+l_{10}{Y}_W^G+l_{11}{d}_3\left(t\right)+1}$

可见，结合水位变化进行修正的直接线性变换的方程组共有11个待求的未知量。

III、镜头非线性畸变模型的建立

为消除或减轻非线性畸变，在像平面坐标系中建立以下非线性畸变模型：

$X_I=x_I+k_1{x}_I(x_I^2+y_I^2)+k_2{x}_I{(x_I^2+y_I^2)}^2+p_1(3x_I^2+y_I^2)+2p_2{x}_I{y}_I+s_1(x_I^2+y_I^2)$ (4)

$Y_I=y_I+k_1{y}_I(x_I^2+y_I^2)+k_2{y}_I{(x_I^2+y_I^2)}^2+2p_1{x}_I{y}_I+p_2(3x_I^2+y_I^2)+s_2(x_I^2+y_I^2)$

式中k₁、k₂为径向畸变系数，p₁、p₂为切向畸变系数，s₁、s₂为薄棱镜畸变系数。根据实际情况可对该畸变模型进行简化，例如可只需要考虑径向畸变。因为引入过多的非线性参数有可能会引起解的不稳定性。

此外，畸变像平面坐标(x_I，y_I)与畸变图像坐标(x_S，y_S)的关系可表示为：

$x_I=(x_S-C_x)S_x^{-1}{d}_x$ (5)

y_I＝(y_S-C_y)d_y

式中：C_x、C_y为镜头光轴与像平面的交点在图像坐标系X_S、Y_S轴上的坐标值，S_x为垂直比例因子，d_x为CCD在x方向相邻像敏单元中心距，d_y为CCD在y方向相邻像敏单元中心距。

同理，理想像平面坐标(X_I，Y_I)与理想图像坐标(X_S，Y_S)的关系可表示为：

$X_I=(X_S-C_x)S_x^{-1}{d}_x$ (6)

Y_I＝(Y_S-C_y)d_y

对应的逆关系式为：

$X_S=X_I{S}_x{d}_x^{-1}+C_x$ (7)

$Y_S=Y_I{d}_y^{-1}+C_y$

设

则最终待求的系数为(C_x，C_y，a₀，k₁，k₂，p₁，p₂，s₁，s₂)。

IV、变换系数求解

1.求解DLT正变换系数

在测站部署阶段勘测N≥6个共面的GCP，将N个平面坐标值

及其高程坐标值代入式(1)直接线性变换方程组的矩阵形式：

WL＝I                    (8)

式中：

$W=[2N\×11]=\left[\begin{array}{ccccccccccc}{X}_W^{G1}& Y_W^{G1}& d_3& 1& 0& 0& 0& 0& -X_S^{G1}{X}_W^{G1}& -X_S^{G1}{Y}_W^{G1}& -X_S^{G1}{d}_3\\ X_W^{G2}& Y_W^{G2}& d_3& 1& 0& 0& 0& 0& -X_S^{G2}{X}_W^{G2}& -X_S^{G2}{Y}_W^{G2}& -X_S^{G2}{d}_3\\ M& M& M& M& M& M& M& M& M& M& M\\ X_W^{\mathrm{GN}}& Y_W^{\mathrm{GN}}& d_3& 1& 0& 0& 0& 0& -X_S^{\mathrm{GN}}{X}_W^{\mathrm{GN}}& -X_S^{\mathrm{GN}}{Y}_W^{\mathrm{GN}}& -X_S^{\mathrm{GN}}{d}_3\\ 0& 0& 0& 0& X_W^{G1}& Y_W^{G1}& d_3& 1& -Y_S^{G1}{X}_W^{G1}& -Y_S^{G1}{Y}_W^{G1}& -Y_S^{G1}{d}_3\\ 0& 0& 0& 0& X_W^{G2}& Y_W^{G2}& d_3& 1& -Y_S^{G2}{X}_W^{G2}& -Y_S^{G2}{Y}_W^{G2}& -Y_S^{G2}{d}_3\\ M& M& M& M& M& M& M& M& M& M& M\\ 0& 0& 0& 0& X_W^{\mathrm{GN}}& Y_W^{\mathrm{GN}}& d_3& 1& -Y_S^{\mathrm{GN}}{X}_W^{\mathrm{GN}}& -Y_S^{\mathrm{GN}}{Y}_W^{\mathrm{GN}}& -Y_S^{\mathrm{GN}}{d}_3\end{array}\right]$

L＝[11×1]＝[l₁，l₂，l₃，l₄，l₅，l₆，l₇，l₈，l₉，l₁₀，l₁₁]^T

$I=[2N\×1]={[X_S^{G1},X_S^{G2},L,X_S^{\mathrm{GN}},Y_S^{G1},Y_S^{G2},L,Y_S^{\mathrm{GN}}]}^T$

用最小二乘法求解如下方程可得DLT正变换系数L：

L＝(W^TW)^-1W^TI。                (9)

2.求解非线性畸变系数

由于摄像机镜头存在畸变，以上求得的矩阵L不是准确值。因此构建目标函数：

$F=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}[{(X_S^{\mathrm{Gi}}-\stackrel{\‾}{X_S^{\mathrm{Gi}}})}^2+{(Y_S^{\mathrm{Gi}}-\stackrel{\‾}{Y_S^{\mathrm{Gi}}})}^2]---\left(10\right)$

式中：

表示第i个地面控制点的图像坐标值，其是由地面控制点的畸变图像坐标及系数(C_x，C_y，a₀，k₁，k₂，p₁，p₂，s₁，s₂)计算得出；

表示假设L为某一初值时，根据式(1)和

计算得到的第i个地面控制点的图像坐标理论值。

如果(C_x，C_y，a₀，k₁，k₂，p₁，p₂，s₁，s₂)和L选取准确，

将逼近于

根据摄像机参数d_x＝4.8/811＝0.00592mm/pixel、d_y＝3.6/508＝0.00707mm/pixel选取初值：

C_x＝320.000、C_y＝240.000、S_x＝1.000，畸变系数k₁，k₂，p₁，p₂，s₁，s₂的初值均为0。理论上，矩阵L的初值可为任意，但考虑迭代的收敛和迭代时间，L值应尽量接近真值。

将参数初值和共面点对坐标

代入式(10)，利用无约束最小二乘法求解，当式(10)取值最小，即：

$\min \left\{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\right[{(X_S^{\mathrm{Gi}}-\stackrel{\‾}{X_S^{\mathrm{Gi}}})}^2+{(Y_S^{\mathrm{Gi}}-\stackrel{\‾}{Y_S^{\mathrm{Gi}}})}^2\left]\right\}$

此时(C_x，C_y，a₀，k₁，k₂，p₁，p₂，s₁，s₂)和L的值就是所要求的解。

3.求解DLT逆变换系数

畸变校正后的水面对应于2D世界坐标系O_w-X_wY_w的平面。将图像坐标(X_S，Y_S)映射到2D世界坐标(X_w，Y_w)的逆关系式为：

$X_w=\frac{p_1{X}_S+p_2{Y}_S+p_3}{p_7{X}_S+p_8{Y}_S+p_9}$ (11)

$Y_w=\frac{p_4{X}_S+p_5{Y}_S+p_6}{p_7{X}_S+p_8{Y}_S+p_9}$

其中，逆变换系数P可由正变换系数L求得，其关系如下：

p₁＝d₃(l₆l₁₁-l₇l₁₀)+(l₆-l₈l₁₀)

p₂＝-[d₃(l₂l₁₁-l₃l₁₀)+(l₂-l₄l₁₀)]

p₃＝d₃(l₂l₇-l₃l₆)+(l₂l₈-l₄l₆)

p₄＝-[d₃(l₅l₁₁-l₇l₉)+(l₅-l₈l₉)]

p₅＝d₃(l₁l₁₁-l₃l₉)+(l₁-l₄l₉)

p₆＝-[d₃(l₁l₇-l₃l₅)+(l₁-l₄l₅)]

p₇＝l₅l₁₀-l₆l₉

p₈＝-(l₁l₆-l₂l₉)

p₉＝l₁l₆-l₂l₅

上述变换系数的求解步骤在测站部署阶段完成。

V、畸变图像校正

校正过程在流场测量阶段进行，包括镜头畸变校正和图像正射校正两部分，分别实现(x_S，y_S)到(X_S，Y_S)的映射，以及(X_S，Y_S)到(X_w，Y_w)的映射。算法实现步骤如下：

1.镜头畸变校正。首先进行畸变像平面坐标的变换，通过式(5)将畸变图像坐标(x_S，y_S)映射到畸变像平面坐标(x_I，y_I)；然后进行像平面坐标的畸变校正，引入预先计算好的畸变系数k₁，k₂，p₁，p₂，s₁，s₂，通过式(4)的非线性畸变模型将(x_I，y_I)表示为(X_I，Y_I)；接下来进行理想像平面坐标的变换，通过式(7)将理想像平面坐标(X_I，Y_I)映射到理想图像坐标(X_S，Y_S)。

2.图像正射校正。首先进行逆DLT变换，引入预先计算好的DLT逆变换系数P，通过式(11)将图像坐标(X_S，Y_S)映射到2D世界坐标(X_w，Y_w)；然后进行图像灰度插值，由于以上求得的2D世界坐标(X_w，Y_w)可能是非整数值，若要在图像上进行显示，需根据定义在整数坐标上的像素灰度值对这些位置进行插值，本方法采用的是双线性内插法，它利用周围4个邻点的灰度值在两个方向上作线性内插以得到待采样点的灰度值，即根据待采样点与相邻点的距离确定相应的权值计算出待采样点的灰度值。相比最邻近插值法，双线性内插法考虑了待采样点周围4个直接邻点对待采样点的影响，因此基本克服了前者灰度不连续的缺点。最后采用像素填充的方法将输出像素逐次映射到新的图像坐标系下，如果一个输出像素被映射到四个输入像素之间，则其灰度值由灰度级插值决定。

至此完成图像的畸变校正。

基于时域中值滤波背景估计的背景差分算法如图4所示。算法流程如下：

1.读取时间间隔为dt的N帧图像I_k作为输入，k＝1，2，L，N。dt在取值时尽量使粒子在相邻帧内运动的位移大于粒径的一半，N的取值要保证所估计的背景像素取到各帧输入图像中稳定时间最长的像素值，通常取N＝10～100即可获得较好的效果；

2.若输入图像受到视频流采集系统的电流暗噪声或编码噪声的干扰而降质时，可采用高斯滤波对图像进行平滑处理；

3.对图像I_k中的像素I_k(i，j)按灰度大小进行排序，并取中值作为背景图像B对应像素B(i，j)的灰度，即：B(i，j)＝Median[I₁(i，j)，I₂(i，j)，L，I_N(i，j)]，对每个像素执行此操作可得到估计的背景图像B；

4.将输入图像I_k与背景图像B做差分得到对应的前景图像F_k，即F_k＝I_k-B。对于差分得到的负值将其置零。

5.若前景图像中运动粒子和背景的对比不强，会使图像的可视效果较差，并可能引起后续流速提取算法中主次相关峰的差异较小导致错误匹配概率增加，此时可进一步采用直方图均衡化的方法进行灰度拉伸。

由于该算法是通过对若干帧图像在时间轴上的滤波来估计当前背景图像，因此场景背景在空间上的色彩分布复杂与否对估计结果影响不大，可以消除静态背景，并对准静态和动态扰动起到抑制作用。

基于相关跟踪的流速矢量提取算法如图5所示。算法流程如下：

1.分析区域划分。选取t＝t₀时刻的图像A和t＝t₀+dt时刻的图像B作为流速矢量提取的一组图像对。将图像A划分成大小为I_x×I_y的分析区域(Interrogation Area，IA)，每个IA对应一个流速矢量V。IA大小的选取需考虑两方面因素：第一，它应当足够小以保留所需的流动细节，尺度小于分析区域尺寸的流场会在分析过程中丢失；第二，它也应当足够大以包含可识别的示踪粒子或模式，变化信息过少会导致匹配的失效。

2.搜索区域选取。对于图像A中以点a(m，n)为中心的分析区域IA(m，n)，在图像B中以对应点b(m，n)为中心选取一个大小为S_x×S_y的搜索区域(Search Area，SA)记为SA(m，n)。SA大小的选取需考虑两方面因素：第一，它应当尽量小以减少计算量，SA过大会产生冗余的匹配；第二，它应当足够大以覆盖两帧间隔内示踪粒子的运动范围，若IA中的大部分粒子运动出SA则会导致匹配的失效。对于边界处的IA，将用图像的边缘对其SA的大小进行裁剪。

3.相关系数求解。在图像B的SA中选取和IA大小相同的匹配区域(Matching Area，MA)，MA是一种滑动窗口，从SA的左上方开始以K(K∈N，K≥1)个像素的水平和垂直步进向右下方滑动，以点b(u，v)为中心的匹配区域记为MA(u，v)。按如下定义的公式求解IA(m，n)和MA(u，v)的归一化互相关系数：

$R[\mathrm{IA}(m,n),\mathrm{MA}(u,v)]=\frac{\underset{i=1}{\overset{\mathrm{Ix}}{\mathrm{\Σ}}}\underset{j=1}{\overset{\mathrm{Iy}}{\mathrm{\Σ}}}{a}_{m,n}(i,j)\×b_{u,v}(i,j)}{\{{\underset{i=1}{\overset{\mathrm{Ix}}{\mathrm{\Σ}}}\underset{j=1}{\overset{\mathrm{Iy}}{\mathrm{\Σ}}}{\left[a_{m,n}(i,j)\right]}^2\underset{i=1}{\overset{\mathrm{Ix}}{\mathrm{\Σ}}}\underset{j=1}{\overset{\mathrm{Iy}}{\mathrm{\Σ}}}{\left[b_{u,v}(i,j)\right]}^2\}}^{1/2}}$

其中：a_m，n(i，j)和b_u，v(i，j)分别为IA(m，n)和MA(u，v)中的点(i，j)处的灰度值；互相关系数满足|R|≤1，越接近于1表示两块区域的相似性越高。

4.相关峰值搜索。对SA(m，n)中的所有MA计算归一化互相关系数后，将得到一个以x、y方向位移为坐标的相关曲面。与IA(m，n)相似性最高的MA(u，v)对应相关曲面的最高峰，而相似性较低的MA，如MA(u，v)，则对应于曲面的次高峰或其它区域。对相关曲面中各点遍历搜索后得到最高相关峰对应的整数坐标点b(u，v)。若最高峰对应的互相关系数小于预先设定的相似性阈值R_TH(如R_TH＝02)，则认为在该搜索区域中没有合适的匹配区域，将此IA标记为NULL(取值为“空”)，不参与后续步骤5～7的处理。

5.亚像素拟合。由于数字图像记录的是离散的灰度信息，因此MA的滑动是以整像素为单位的，MA(u，v)的中心像素点b(u，v)并非最佳匹配点，位移的最大误差可达0.5个像素。在大尺度场景下，该误差将直接影响流速测量结果。本方法将相关峰分布区域的曲面看作高斯曲面，通过3点高斯拟合的方法利用相关曲面的峰值邻域计算出具有亚像素精度的峰值坐标b(u₀，v₀)。计算公式如下：

$\left\{\begin{array}{c}{u}_0=u+\frac{\ln R(u-1,v)-\ln R(u+1,v)}{2\ln R(u-1,v)-4\ln R(u,v)+2\ln R(u+1,v)}\\ v_0=v+\frac{\ln R(u,v-1)-\ln R(u,v+1)}{2\ln R(u,v-1)-4\ln R(u,v)+2\ln R(u,v+1)}\end{array}\right.$

其中：R(u，v)是相关曲面峰值，R(u-1，v)、R(u+1，v)、R(u，v-1)、R(u，v+1)表示峰值附近4邻域的互相关系数。

6.流速矢量计算。最高相关峰坐标点b(u₀，v₀)在图像A中的对应点为a(u₀，v₀)，从点a(m，n)到点a(u₀，v₀)的位移即为分析区域IA(m，n)的最大似然位移：

D(m，n)＝(dx，dy)_m，n＝(u₀-m，v₀-n)

其中：d_x、d_y分别为x和y方向的位移分量。将位移除以两帧间的时间间隔dt就得到了流速矢量：

$V(m,n)=\frac{D(m,n)}{\mathrm{dt}}=\frac{{(d_x,d_y)}_{m,n}}{\mathrm{dt}}={(v_x,v_y)}_{m,n}=(\frac{u_0-m}{\mathrm{dt}},\frac{v_0-n}{\mathrm{dt}})$

其中：v_x、v_y分别为x和y方向的流速分量。

7.流速矢量绘制。为便于直观的显示流场，以IA(m，n)的中心a(m，n)为起点，L_x＝λ_xv_x、L_y＝λ_yv_y为x和y方向的长度(其中λ_x、λ_y为缩放比例)，采用三线法将矢量箭头绘制于图像A中。对于成像平面中流动区域内(可用掩膜标识出)的所有IA重复步骤2～7可得到t＝t₀时刻一个图像对的原始流速矢量场，它包含了在成像区域中每个分析区域内求得的瞬时流速矢量。

基于空域中值滤波的流速场后处理算法如图6所示。算法流程如下：

1.从原始流速矢量场中选取要处理的分析区域；

2.以分析区域为中心选取3×3(或5×5)的邻域建立x方向和y方向的滤波窗口：

$x=\left[\begin{array}{ccc}{x}_{11}& x_{12}& x_{13}\\ x_{21}& x_{22}& x_{23}\\ x_{31}& x_{32}& x_{33}\end{array}\right]$ $y=\left[\begin{array}{ccc}{y}_{11}& y_{12}& y_{13}\\ y_{21}& y_{22}& y_{23}\\ y_{31}& y_{32}& y_{33}\end{array}\right]$

其中位于顶点的窗口大小为2×2，位于边缘的窗口大小为2×3；

3.将滤波窗口中各区域的流速矢量分解到x方向和y方向；

4.采用冒泡法对x方向和y方向的流速分量按大小进行排序，其中标记为NULL的矢量缺失区域(如图中的x₁₃、x₂₃、x₃₃及y₁₃、y₂₃、y₃₃)不参与排序；

5.计算序列的中值作为分析区域的新流速分量：

x′₂₂＝Median(x)＝Median(x₁₁，x₁₂，x₂₁，x₂₂，x₃₁，x₃₂)＝(x₂₁+x₃₁)/2

y′₂₂＝Median(y)＝Median(y₁₁，y₁₂，y₂₁，y₂₂，y₃₁，y₃₂)＝(y₂₁+y₃₂)/2

6.按照新流速分量更新分析区域内的流速矢量。对图像中的所有分析区域重复步骤1～6得到修正后的瞬时流速场，并绘制于增强前的图像A。

Claims (6)

1.一种基于视频图像的河流流量在线测验实现方法，其特征在于：

该方法中设置有如下硬件系统：

a.用于实现水面流动可视化的水流示踪系统；

b.用于拍摄和存储水面实时视频流的视频流采集系统；

c.用于给图像正射变换提供地面控制点及其坐标信息的摄像机标定系统；

d.用于处理数字图像的图像处理系统；

e.用于计算河流流量的流量估计系统；

该方法分为三个阶段：

（一）测站部署阶段，包括如下步骤：

（1）选择测流河段；

（2）获取测流断面水下地形；

（3）部署现场硬件系统；

（4）勘测地面控制点在世界坐标系下的平面坐标和高程坐标；

（5）计算用于图像畸变校正的DLT变换系数；

（二）流场测量阶段，包括如下步骤：

（1）启动水流示踪系统投放示踪粒子；

（2）启动视频流采集系统拍摄包含示踪粒子的水面视频流，然后按固定的时间间隔取帧并数字化处理；

（3）由图像处理系统处理数字图像：包括图像畸变校正、图像增强、流速矢量提取、流速场后处理，最后得到修正后的水面瞬时流速场；

（三）流量估计阶段，该阶段是先勘测测流断面水位值，然后流量估计系统采用流速-面积法计算出河流流量，计算步骤如下：

（1）计算时间平均流速场：

对第i个分析区域的瞬时流速矢量V_i在时域进行信息融合，得到时间平均流速矢量

将

投影到流动主方向X，得到标量的时间平均水面流速U_i，i=1,2,...,M，M为时间平均流速矢量总数；

（2）计算水面插值流速：

根据流速-面积法，将测流断面划分为以n+1条测速垂线为边界的n个子断面，对测速垂线相邻分析区域的时间平均水面流速U_i进行线性插值，得到测速垂线的水面插值流速U′_v；

（3）计算深度平均流速

其中α为流速系数；

（4）计算子断面流量：

计算子断面S_v的流量

其中

是子断面的平均流速，

的大小为子断面两条边界垂线深度平均流速的均值，即

A_v(l)是水位值为l时子断面S_v的过水面积；

（5）将各个子断面的流量求和，得到所测河流的总断面流量，即求得河流流量。

2.根据权利要求1所述的基于视频图像的河流流量在线测验实现方法，其特征在于：该方法中还设置有由光源和光强传感器组成的水面照明系统，所述水面照明系统部署于测流断面的两岸。

3.根据权利要求1所述的基于视频图像的河流流量在线测验实现方法，其特征在于：所述阶段（一）的步骤（3）的内容如下：

水流示踪系统部署于紧邻测流断面的上游，架设视频流采集系统中的CCTV摄像机，使其光轴垂直于水流方向，将摄像机标定系统中的地面控制点平行设置于测流河段两岸最高水位以上一水平面处，地面控制点覆盖完整的摄像机视场区域。

4.根据权利要求1或3所述的基于视频图像的河流流量在线测验实现方法，其特征在于：所述摄像机标定系统是由至少六个地面控制点组成。

5.根据权利要求1所述的基于视频图像的河流流量在线测验实现方法，其特征在于：所述阶段（一）的步骤（5）的方法为：

利用视频流采集系统拍摄一幅视频帧，提取出地面控制点的图像坐标，基于直接线性变换的方法求解出用于图像畸变校正的DLT变换系数。

6.根据权利要求1、3或5所述的基于视频图像的河流流量在线测验实现方法，其特征在于：所述地面控制点是采用铝合金材料制成的十字形或正方形靶面。

Images