CN102866260A - 非接触式河流表面流场成像量测方法 - Google Patents

Abstract

本发明公开了一种非接触式河流表面流场成像量测方法，属于明渠测流技术领域。所述方法首先通过在图像传感器前加装近红外滤光片实现近红外成像，增强水流示踪物与水面背景间的亮度对比；其次对图像进行空域高通滤波，抑制水面光学噪声；然后采用快速傅立叶互相关算法进行运动矢量估计，得到各分析区域的位移矢量；接下来采用一种基于时空联合滤波的方法重建时间平均流场；最后采用直接线性变换方法对流速矢量进行正射校正，实现流场定标及可视化。本发明能够稳定、可靠地用于复杂水面光学环境下的河流表面流场量测。

Description

非接触式河流表面流场成像量测方法

技术领域

本发明涉及一种明渠水流监测方法，尤其涉及一种非接触式河流表面流场成像量测方法，属于明渠测流技术领域。

背景技术

洪水期间精确的流速、流场、时序信息以及总径流变化速率的高分辨率估计，是河流水文学、河流地貌学以及河流生态学研究的重要依据。然而，高洪时期的水流往往具有含沙量高，混杂物多，流速快等特点，导致传统测流仪器无法正常施测。大尺度粒子图像测速（Large-Scale Particle Image Velocimetry，LSPIV）是一种安全、高效的全场流速测量技术。它采用模式识别的方法匹配、跟踪无畸变连续视频图像中的水面示踪物，实现流场的定量表示及可视化。作为实验室环境下的粒子图像测速（PIV）技术在大尺度现场环境下的扩展，不仅可用于常规条件下天然河道水流紊动特性的研究，其非接触特性更使之成为极端条件下河流流量测量少数可行的方法之一。

但是，极端条件下流体物化属性及测验环境复杂多变，特别是观测窗口内光强、水面反射及水下散射和天空散射等光线扰动的复杂光学环境下，水流示踪/漂浮物的尺度、旋转与遮挡等非刚体运动状态变化将导致运动矢量估计的困难，且目标时空分布不均或密度较低，极易导致水面流场重建错误。现有方法直接移植或借用光学成像检测和图像处理领域的现有技术，达不到水流示踪物可靠和稳定地连续检测的要求，使得LSPIV在野外难以推广应用。如何在极端测验条件下，解决水面弱小目标光学信息的连续采集和检测、运动矢量估计等关键技术，实现流场重建与定标是一项极具挑战性的研究课题。

已有的研究成果表明：水体对太阳电磁辐射的吸收在700-1600nm的近红外波段较强，1400nm和1900nm附近的吸收率甚至接近100%；对于河流、湖泊、湿地等天然水体，即使水很浅也能吸收绝大部分的红外辐射；而多数陆生植物以及滨海浅水的海草却总在近红外波段出现反射峰。因而在近红外图像中，水面背景与树叶、小木片等水面示踪物间必然存在较大的亮度差异。基于上述原理，若在近红外波段下成像可在一定程度上抑制水面光学噪声，较之目前现有的基于可见光成像的LSPIV方法，有利于提高运动矢量估计的峰值信噪比。

发明内容

本发明针对现有方法存在的不足，提供了一种非接触式河流表面流场成像量测方法。

该方法通过以下技术手段实现：

方法首先通过在图像传感器前加装近红外滤光片实现近红外成像，增强水流示踪物与水面背景间的亮度对比；其次对图像进行空域高通滤波，抑制水面光学噪声；然后采用快速傅立叶互相关算法进行运动矢量估计，得到各分析区域的位移矢量；接下来采用一种基于时空联合滤波的方法重建时间平均流场；最后采用直接线性变换方法对流速矢量进行正射校正，实现流场定标及可视化。

非接触式河流表面流场成像量测方法，分为图像采集、图像处理、流场定标及流场可视化四个步骤。其特征在于：图像采集和图像处理在以图像处理专用DSP为核心的智能相机上完成，仅上传时间平均流场信息及序列的首帧图像；流场定标及流场可视化在上位机上完成。其中，图像处理步骤可分为图像增强、运动矢量估计和时间平均流场重建三项任务。

所述的图像采集步骤通过在智能相机的单色CMOS图像传感器前加装近红外滤光片实现近红外成像，增强水流示踪物与水面背景间的亮度对比。

所述的图像增强任务采用空域高通滤波的方法抑制水面光学噪声。滤波采用如下7×7的卷积模板：

$H=\left[\begin{array}{ccccccc}-1& -1& -1& -3& -1& -1& -1\\ -1& -1& -3& 3& -3& -1& -1\\ -1& -3& 3& 6& 3& -3& -1\\ -3& 3& 6& 10& 6& 3& -3\\ -1& -3& 3& 6& 3& -3& -1\\ -1& -1& -3& 3& -3& -1& -1\\ -1& -1& -1& -3& -1& -1& -1\end{array}\right].$

所述的运动矢量估计任务采用快速傅立叶变换互相关算法（FFT-CC）计算各分析区域的位移矢量。

所述的时间平均流场重建任务采用一种基于时空联合滤波的方法重建时间平均流场。首先以矢量角度的全局统计为依据识别错误矢量，然后对瞬时位移场进行时间轴滤波更新时间平均位移场，最后采用空域中值滤波修正、平滑时间平均流场。

所述的流场定标步骤采用直接线性变换（DLT）的方法对时间平均流场中的位移矢量进行透视畸变校正并从图像坐标变换为世界坐标。其特征在于：1）不对所有图像中的每个像素进行坐标变换，仅完成首帧图像的透视畸变校正及位移矢量单位从像素到米的变换。2）计算DLT系数时采用一种人机交互的地面控制点（GCP）坐标提取方法：当智能相机上传首帧图像后用户手动选取目标点，然后系统在目标点的邻域内自动检测GCP的中心角点并拟合得到亚像素图像坐标。

所述的流场可视化步骤首先采用线性插值的方法对流场进行网格化处理，然后采用三线法绘制矢量箭头。

本发明具有以下有益效果：

1、抗干扰能力强。近红外波段成像，可有效抑制复杂的水面光学噪声、增强目标与背景间的亮度对比，从而提高水流示踪物运动矢量估计的准确性。仪器采用采用在单色CMOS图像传感器前加装近红外红外滤光片这种简单易行的光学滤波方法，使后续流程中原本复杂的图像增强算法得以省略或简化。

2、量测精度高。在时间平均流场重建方法中，时间轴滤波充分利用了示踪物在水面运动的时间冗余信息，而空域中值滤波对奇异矢量的修正和缺失矢量的插补十分有效。相比直接求平均等线性方法，这种基于非线性滤波方法更加适用于示踪物密度低、时空分布不均，且受环境干扰影响更大的天然河流。

3、实时性强。方法综合考虑了智能相机和上位机的处理能力及传输带宽，采用的算法占用更少的系统资源并满足实时性要求。

鉴于以上特点，本系统能够稳定、可靠地用于复杂水面光学环境下的河流表面流场量测。

附图说明

图1是本发明涉及的量测方法的工作流程图。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

本发明的量测方法的工作流程如图1所示。智能相机采用130万像素的单色CMOS芯片MT9M001（Aptina公司）作为图像传感器，并以图像处理专用DSP芯片TMS320DM642（TI公司）作为处理器；上位机为通用PC机。按照信息获取和处理的过程，工作流程可分为图像采集、图像处理、流场定标及流场可视化四个步骤。图像采集和图像处理在智能相机上完成，仅上传时间平均流场信息及序列的首帧图像，用于直观显示和流场定标中变换系数的求解。图像及流场数据的传输通过智能相机和上位机端分别基于TI NDK开放包和WinSocketAPI函数开发的网络通信模块实现。流场定标及流场可视化在上位机上完成。

图像采集。测量任务开始后，智能相机首先根据上位机传递的参数对CMOS图像传感器进行初始化的参数配置。配置为连续捕获模式。由于TMS320DM642的视频口不支持FRAME_VALID信号，因此通过记录每行的水平消隐脉冲数和每帧的垂直消隐脉冲数来实现帧同步。当LINE_VALID有效时，TMS320DM642在PIXCLK的同步下并行读入像素数据D[9:0]，并暂存于视频端口内部集成的FIFO中；当FIFO达到溢出状态时触发产生一个EDMA事件，通过EDMA通道将数据送入外扩的SDRAM中，供DSP调用进行图像处理。

图像处理。图像处理可分为图像增强、运动矢量估计和时间平均流场重建三项任务：

1.图像增强。近红外成像可以在很大程度上增强示踪物和水面的对比度，但却无法彻底去除某些强噪声，特别是倒影。在示踪物密度低的区域这些静态模式占主导，将无法估计出运动矢量，因此需要进行图像增强处理实现背景抑制。在利用近红外成像系统获取的水面图像中，含噪水面主要是大面积缓慢变化的低频部分，而示踪物则表现为与背景不相关的孤立亮斑，其边缘和角点包含了大量高频信息。对此，采用空域高通滤波的方法抑制背景。

图像的空域高通滤波可描述为输入图像I中的一系列滑动窗口W_ij与高通滤波模板H的卷积运算：

$G_{\mathrm{ij}}=W_{\mathrm{ij}}*H=\underset{m=1}{\overset{S}{\mathrm{\Σ}}}\underset{n=1}{\overset{S}{\mathrm{\Σ}}}I(i-S/2+m-1,j-S/2+n-1)\×(m,n),---\left(1\right)$

i=1,2,…,M；j=1,2,…,N

其中，(i，j)表示滑动窗口的中心坐标，M、N分别为图像的长度和宽度，W_ij与H具有相同的尺寸S。图像滤波过程分4步完成：

（1）将窗口在图中滑动，其中心与图像中待求像素的位置重合；

（2）将模板中的系数与滑动窗口中的对应元素相乘；

（3）对所有乘积项求和并赋值给图像中对应模板中心的像素；

（4）当滑动窗口遍历图像中的每个像素后就得到高通滤波图像。

高通滤波模板的设计是影响滤波性能的关键。考虑到水面示踪物在图像中的大小通常在1×1至5×5个像素之间的特点，本发明设计了一种7×7的高通滤波模板，如下所示：

$H=\left[\begin{array}{ccccccc}-1& -1& -1& -3& -1& -1& -1\\ -1& -1& -3& 3& -3& -1& -1\\ -1& -3& 3& 6& 3& -3& -1\\ -3& 3& 6& 10& 6& 3& -3\\ -1& -3& 3& 6& 3& -3& -1\\ -1& -1& -3& 3& -3& -1& -1\\ -1& -1& -1& -3& -1& -1& -1\end{array}\right].$

模板H的中心权值10最大，用以提高弱目标的灰度值；高权值6和3按菱形阶梯状分布，用于扩展小目标的面积；低权值-3用于削弱孤立噪声点，并和高权值3形成强边缘；低权值-1使滤波后的背景灰度变得更均匀；高低权值之和相抵消可有效地去除低频分量。

2.运动矢量估计。采用基于区域匹配的快速傅立叶变换互相关算法（FFT-CC）估计运动矢量。其算法流程如下：

（1）选取分析区域。读入t=t₀时刻的图像A和t=t₀+dt时刻的图像B作为运动矢量估计的一组图像对。在两图中分别选取以点I₀(m,n)和I₁(m,n)为中心、大小均为M×N个像素的子区域，称为分析区域（IA）和匹配区域（MA），每个IA对应一个位移矢量

M、N的取值应为2的整数倍，并需兼顾流场的空间分辨率、示踪物的位移及区域内的信息量等因素。为保证计算精确度，通常取示踪物位移大小的2倍以上。

（2）求解相关系数。首先将IA和MA对应的两个二维实矩阵f(m,n)和g(m,n)扩充成M×N的复数型矩阵，即每个像素由整型的实部和整型的虚部组成，其中虚部为0。接着进行二维快速傅立叶变换（FFT2），采用TI公司针对C64x系列DSP优化设计的快速傅立叶变换库函数DSP_fft 32x32()实现。先逐行进行k点一维行变换，然后对所得矩阵进行转置，再逐列进行k点一维列变换，得到频域形式：

$\left\{\begin{array}{c}F(u,v)=\mathrm{FFT}2\left(f(m,n)\right)\\ G(u.v)=\mathrm{FFT}2\left(g(m,n)\right)\end{array}\right.,---\left(2\right)$

则功率谱即为频域互相关函数：

Φ(u，v)=F(u，v)G^*(u，v)。                   （3）

其中，G^*(u,v)为G(u,v)的共轭，矩阵的点乘采用DSP_mul32()函数实现。最后进行二维傅立叶逆变换（IFFT2），采用DSP_ifft 32x32()函数依次对行和列进行一维逆变换，得到空域互相关函数：

φ(m，n)=IFFT2(Φ(u,v))。                    （4）

所形成的相关曲面反映了互相关系数在X和Y方向上的分布情况。

（3）检测相关峰值。对相关曲面中各点及其邻域进行遍历搜索得到峰值φ_max(m,n)=φ(m_p,n_p)，其对应坐标(m_p,n_p)即为分析区域内模式在X和Y方向上相对位移的估计值。

（4）确定绝对位移。在计算相关函数时，实际是把有限长度的序列看作周期性序列中的一个周期来处理，也就是把有限长度为M×N的序列延拓到整个变量空间，形成周期性序列。这样，所得到的相关函数φ(m，n)也是周期性的，即：

φ(m，n)=φ(m±M，n±N)。             （5）

因此，φ(m，n)取得最大值对应于模式正向移动和反向移动这两种情况。由于步骤（1）中假设示踪物的位移小于IA大小1/2，则绝对位移可以按以下方式确定：

$\mathrm{\Δx}\left\{\begin{array}{cc}{m}_p,& m_p\≤[M/2]\\ m_p-M,& M_p>[M/2]\end{array}\right.$

$\mathrm{\Δy}=\left\{\begin{array}{cc}{n}_p,& n_p\≤[N/2]\\ n_p-N,& n_p>[N/2]\end{array}\right..$

(6)

将分析区域按照预设的步进或重叠度在图像中漫游，可得到图像对A和B的瞬时位移场

其中单个位移矢量

可用其分析区域的中心点坐标(m，n)和矢量长度D_m，n(Δx，Δy)唯一表示。

3.时间平均流场重建。

本发明结合河流表面流场的先验信息及流速的时空冗余信息，采用矢量角度的全局统计作为错误矢量识别的依据，设计了一种基于时空联合滤波的时间平均流场重建方法，具体步骤如下：

（1）矢量类型识别。首先计算瞬时位移场中各矢量

与X方向的夹角

并建立方向矩阵，然后以10°为角度区间统计矢量方向在0-359°范围内的全局直方图，将对应矢量数量最多的区间标记为流场的主方向θ_M，最后根据角度阈值θ_TH（如15°）形成的判别窗口依次识别每个矢量的类型，并建立矢量类型查找表T如下：

$T=\left\{\begin{array}{cc}0,& |\mathrm{\θ}-{\mathrm{\θ}}_M|>{\mathrm{\θ}}_{\mathrm{TH}}\\ 1,& |\mathrm{\θ}-{\mathrm{\θ}}_M|\≤{\mathrm{\θ}}_{\mathrm{TH}}\end{array}\right.,---\left(7\right)$

其中，1表示正确矢量，0表示错误矢量，其个数分别用N1、N0表示。则流场的矢量正确率为P=N1/(N1+N0)，若P大于一个可接受的低阈值P_THl（如50%），则认为是有效的瞬时流场，参与后续处理，否则丢弃。

（2）时间轴滤波。按步骤（1）计算得到t₀和t₁时刻两个有效流场的查找表T₀和T₁后，对瞬时位移场

和

进行时间轴滤波得到更新的时间平均位移场

其原则如下：

${\stackrel{\→}{S}}_{\mathrm{TA}}=\left\{\begin{array}{cc}{\stackrel{\→}{S}}_0,& T_0=1{,T}_1=0\\ {\stackrel{\→}{S}}_1,& T_0=0,T_1=1,\\ ({\stackrel{\→}{S}}_0+{\stackrel{\→}{S}}_1)/2,& T_0=1,T_1=1\end{array}\right.---\left(8\right)$

其中，

为更新后的时间平均位移矢量。然后对T₀和T₁进行“或”运算得到新的时间平均查找表T_TA并统计

的正确率P_TA。若P_TA小于一个高阈值P_THh（如90%），则读入下一时刻的图像对(I_t,I_t+1)、估计瞬时位移矢量场并重复步骤（1）和（2），直到P_TA大于P_THh或处理完一次测量的所有N个瞬时流场。

（3）空域中值滤波。若时间轴滤波后的流场依然存在错误矢量，则采用3×3的中值滤波模板分别选取IA邻域内正确的位移分量，按大小排序后取其中值，就此得到时空联合滤波后重建的时间平均流场，并以网络帧的形式发送给上位机。

流场定标。为获得正射视角下的矢量场，本发明采用直接线性变换（DLT）的方法对时间平均流场中的位移矢量进行透视畸变校正并从图像坐标变换为世界坐标。和现有方法在运动矢量估计前首先对所有图像中的每个像素进行坐标变换不同，本发明仅完成首帧图像的透视畸变校正及位移矢量单位从像素到米的变换，大大减少了高分辨率图像坐标变换所需的计算量。其中，为计算DLT系数，采用一种人机交互的地面控制点（GCP）坐标提取方法：当智能相机上传首帧图像后用户手动选取目标点，然后系统在目标点的邻域内自动检测GCP的中心角点并拟合得到亚像素图像坐标。

流场可视化。首先采用线性插值的方法对流场进行网格化处理，然后采用三线法绘制矢量箭头。

至此完成一次河流表面流场的量测。

Claims (1)

1.一种非接触式河流表面流场成像量测方法，分为图像采集、图像处理、流场定标及流场可视化四个步骤, 其特征在于：图像采集和图像处理在以图像处理专用DSP为核心的智能相机上完成，仅上传时间平均流场信息及序列的首帧图像；流场定标及流场可视化在上位机上完成，其中，图像处理步骤可分为图像增强、运动矢量估计和时间平均流场重建三项任务;

所述的图像采集步骤通过在智能相机的单色CMOS图像传感器前加装近红外滤光片实现近红外成像，增强水流示踪物与水面背景间的亮度对比;

所述的图像增强任务采用空域高通滤波的方法抑制水面光学噪声;滤波采用如下7×7的卷积模板：

，

所述的运动矢量估计任务采用快速傅立叶变换互相关算法（FFT-CC）计算各分析区域的位移矢量;

所述的时间平均流场重建任务采用一种基于时空联合滤波的方法重建时间平均流场;首先以矢量角度的全局统计为依据识别错误矢量，然后对瞬时位移场进行时间轴滤波更新时间平均位移场，最后采用空域中值滤波修正、平滑时间平均流场；

所述的流场定标步骤采用直接线性变换（DLT）的方法对时间平均流场中的位移矢量进行透视畸变校正并从图像坐标变换为世界坐标;1）不对所有图像中的每个像素进行坐标变换，仅完成首帧图像的透视畸变校正及位移矢量单位从像素到米的变换;2）计算DLT系数时采用一种人机交互的地面控制点（GCP）坐标提取方法：当智能相机上传首帧图像后用户手动选取目标点，然后系统在目标点的邻域内自动检测GCP的中心角点并拟合得到亚像素图像坐标;

所述的流场可视化步骤首先采用线性插值的方法对流场进行网格化处理，然后采用三线法绘制矢量箭头。

Images