基于序贯图像矢量平均的时均流场重建方法

摘要

本发明公开了一种基于序贯图像矢量平均的时均流场重建方法，属于非接触式明渠测流技术领域。所述方法首先载入两帧连续图像作为一组图像对，并采用预先设定的掩膜提取图像对中的有效水面区域；其次采用快速傅立叶变换互相关算法对图像对中的所有分析区域IA进行运动矢量估计；然后采用全局角度直方图方法检测流动主方向并识别瞬时流场中的错误矢量；接下来进行流场时间滤波，舍弃当前错误的瞬时矢量，并用当前正确的瞬时矢量替代此前错误的时均矢量，用当前正确的瞬时矢量平均此前正确的时均矢量；最后进行错误矢量插值和流场平滑滤波得到重建的时均流场。本发明以时均流场的矢量正确率为依据控制求平均的进程，改善了流场重建的时间分辨率。

说明书

技术领域

本发明涉及一种二维时均流场重建方法，尤其涉及一种基于序贯图像矢量平均的时均流 场重建方法，属于非接触式明渠测流技术领域。

背景技术

天然河流中的水体在河槽中运动时受到断面形状、坡度、糙率、水深、弯道、风、气压、 潮汐等因素的影响而产生紊流。紊流内部水质点的瞬时流速在大小和方向上均随时间变化， 呈现出脉动现象，但在足够长的时段内其均值保持稳定。因此，对于面向断面流量监测的河 流水面成像测速系统而言，时均流场的重建比瞬时流场的获取更为重要。然而在天然河流中， 由于水面示踪物密度低、时空分布不均及水面光学噪声的影响，视场中的某些待测区域可能 会暂时性地缺乏明显的示踪物，导致无法估计出这些区域的瞬时矢量或出现错误矢量。

对于基于灰度相关匹配法的运动矢量估计，按照求平均操作在处理流程中所处的阶段及 作用对象的不同，可以分为图像平均、相关平均和矢量平均三种时均流场重建策略。图像平 均策略的优点是仅需计算一次相关测度，运算速度快；缺点是需要一次读入完整的图像序列， 占用大量存储空间；并且受到平均帧数的限制，对图像质量有较高的要求。相关平均策略的 优点是能够有效提高相关曲面的信噪比，降低出现错误矢量的概率；缺点是对相关曲面的直 流偏置噪声敏感，依赖于有效的图像背景抑制方法。矢量平均策略的优点是灵活、高效，可 以处理任意方法获得的网格化矢量场；由于是时均流场重建的最高层次，可以获得瞬时流场 以及特征级的中间参数，具有较高的容错性；缺点是当示踪物密度较低时，对相关曲面信噪 比低下引起的错误矢量敏感。以上策略都在一定程度上提高了位移估计的可靠性，但当错误 矢量个数超过一定比例（如30%）时，直接采用线性求平均的方法都将导致时均流场重建精 度过低。因此错误矢量的识别和修正对低信噪比条件下的时均流场重建具有重要意义。

发明内容

本发明针对现有方法存在的不足，提供了一种基于序贯图像矢量平均的时均流场重建方 法。

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

方法首先载入两帧连续图像作为一组图像对，并采用预先设定的掩膜提取图像对中的有 效水面区域；其次采用快速傅立叶变换互相关算法对图像对中的所有分析区域IA进行运动矢 量估计；然后采用全局角度直方图检测流动主方向并识别瞬时流场中的错误矢量；接下来进 行流场时间滤波，舍弃当前错误的瞬时矢量，并用当前正确的瞬时矢量替代此前错误的时均 矢量，用当前正确的瞬时矢量平均此前正确的时均矢量；最后进行错误矢量插值和流场平滑 滤波得到重建的时均流场。

所述的错误矢量识别方法，其实现过程如下：对于i时刻的瞬时流场，首先计算所有矢量 S_i(x,y)与x方向的夹角θ_i(x,y)，得到瞬时流场中矢量角度的极坐标分布：

θ_i(x,y)＝arctan(Δy_i/Δx_i)   (1)

其中，(x,y)表示矢量的图像坐标，Δx、Δy分别表示矢量在x和y方向的大小；然后以Δθ为 角度区间建立矢量方向在0-359°范围内的全局角度直方图；接下来依次搜索各角度区间，将 矢量数量最多的区间标记为流动主方向θ_M；最后以角度阈值θ_TH为判别窗口识别每个矢量的 类型，并建立矢量类型标志如下：

$F_i(x,y)=\left(\begin{array}{cc}0,& |{\theta }_i(x,y)-{\theta }_M|>{\theta }_{\mathrm{TH}}\\ 1,& |{\theta }_i(x,y)-{\theta }_M|\le {\theta }_{\mathrm{TH}}\end{array}\right)---\left(2\right)$

其中，1表示正确矢量，0表示错误矢量。

所述的流场时间滤波过程用下式描述为：

${\bar{S}}_i(x,y)=\left(\begin{array}{cc}\mathrm{NULL},& {\bar{F}}_{i-1}(x,y)=0,F_i(x,y)=0\\ S_i(x,y),& {\bar{F}}_{i-1}(x,y)=0,F{(x,y)}_i=1\\ {\bar{S}}_{i-1}(x,y),& {\bar{F}}_{i-1}(x,y)=1,F_i(x,y)=0\\ \frac{N_i(x,y)-1}{{\text{N}}_i(x,y)}{\bar{S}}_{i-1}(x,y)+\frac{1}{N_i(x,y)}{S}_i(x,y),& {\bar{F}}_{i-1}(x,y)=1,{F(x,y)}_i=1\end{array}\right)---\left(3\right)$

其中，分别为当前和此前时刻的时均矢量；S_i(x,y)为当前时刻的瞬时矢 量；错误矢量标记为“NULL”，以区别于零矢量；N_i(x,y)为参与当前矢量平均的正确矢量 个数；F_i(x,y)分别为此前时均矢量的类型标志和当前瞬时矢量的类型标志，对二 者进行或运算得到当前时均矢量的类型标志，如下：

${\bar{F}}_i(x,y)={\bar{F}}_{i-1}(x,y)|F_i(x,y)---\left(4\right)$

其中，“|”表示或运算符。遍历每个所有IA后统计得到正确和错误矢量的个数N₁和N₀，则 时均流场的矢量正确率为：

${\bar{P}}_i=N_1/(N_{\text{1}}+N_0)\times 100\%---\left(5\right)$

若小于一个阈值P_TH（如90%），则继续处理下一时刻的图像对，直到大于P_TH并且求平 均的帧数达到预设的阈值，得到一个包含错误矢量的时均流场和对应的矢量类型查找表；反 之则认为时均流场重建失效。

所述的错误矢量插值采用错误矢量邻域内正确矢量的均值对错误矢量进行插值：

${\bar{S}}_0(x,y)=\frac{1}{N_1(x,y)}\underset{n=-1}{\overset{1}{\Sigma }}\underset{m=-1}{\overset{1}{\Sigma }}\bar{F}(x+n,y+m)·\bar{S}(x+n,y+m)---\left(6\right)$

其中，N₁(x,y)表示错误矢量邻域内正确矢量的个数，m、n表示相邻矢量的偏移量， 对于3×3邻域有-1≤m≤1、-1≤n≤1。

所述的流场平滑滤波，采用时均流场和滤波器的空域卷积实现：

$\tilde{S}(x,y)=\bar{S}(x,y)*G(x,y)---\left(7\right)$

其中，“*”表示卷积运算，滤波器采用大小为3×3的归一化二维离散高斯窗，表示为：

$G(x,y)=\frac{1}{16}\left(\begin{array}{ccc}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right)---\left(8\right)$

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

方法利用了断面流速方向一致性分布的特点，采用基于非线性统计的全局角度直方图方 法检测流动主方向，较好地解决了示踪物密度低、时空分布不均及水面光学噪声影响下区域 性错误矢量的识别问题，改善了河流水面流场的估计精度。此外，方法充分利用了运动目标 时空分布的冗余信息，以时均流场的矢量正确率为依据控制求平均的进程，有效提高了流场 测量的时间分辨率。相比图像平均和相关平均策略，本方法在实时性和硬件开销上具有明显 优势，适用于河流水面流场、流量连续在线监测系统。

附图说明

图1是本发明涉及的时均流场重建方法的工作流程图。

具体实施方式

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

本发明的时均流场重建方法的工作流程如图1所示。共分为以下七个步骤：

（1）载入图像对，从相机采集的图像序列中载入两帧连续图像作为一组图像对。

（2）水面区域提取，采用预先设定的掩膜提取图像对中的有效水面区域。

（3）运动矢量估计，采用快速傅立叶变换互相关算法对图像对中的所有分析区域IA进 行运动矢量估计，得到当前图像对的瞬时流场；

（4）错误矢量识别，对于i时刻的瞬时流场，首先计算所有矢量S_i(x,y)与x方向的夹角 θ_i(x,y)，得到瞬时流场中矢量角度的极坐标分布：

θ_i(x,y)＝arctan(Δy_i/Δx_i)   (1)

其中，(x,y)表示矢量的图像坐标，Δx、Δy分别表示矢量在x和y方向的大小；然后以Δθ为 角度区间建立矢量方向在0-359°范围内的全局角度直方图；接下来依次搜索各角度区间，将 矢量数量最多的区间标记为流动主方向θ_M；最后以角度阈值θ_TH为判别窗口识别每个矢量的 类型，并建立矢量类型标志如下：

$F_i(x,y)=\left(\begin{array}{cc}0,& |{\theta }_i(x,y)-{\theta }_M|>{\theta }_{\mathrm{TH}}\\ 1,& |{\theta }_i(x,y)-{\theta }_M|\le {\theta }_{\mathrm{TH}}\end{array}\right)---\left(2\right)$

其中，1表示正确矢量，0表示错误矢量。

（5）流场时间滤波，基本思想是：舍弃当前错误的瞬时矢量，并用当前正确的瞬时矢量 替代此前错误的时均矢量，用当前正确的瞬时矢量平均此前正确的时均矢量，描述为：

${\bar{S}}_i(x,y)=\left(\begin{array}{cc}\mathrm{NULL},& {\bar{F}}_{i-1}(x,y)=0,F_i(x,y)=0\\ S_i(x,y),& {\bar{F}}_{i-1}(x,y)=0,F{(x,y)}_i=1\\ {\bar{S}}_{i-1}(x,y),& {\bar{F}}_{i-1}(x,y)=1,F_i(x,y)=0\\ \frac{N_i(x,y)-1}{N_i(x,y)}{\bar{S}}_{i-1}(x,y)+\frac{1}{N_i(x,y)}{S}_i(x,y),& {\bar{F}}_{i-1}(x,y)=1,F{(x,y)}_i=1\end{array}\right)---\left(3\right)$

其中，分别为当前和此前时刻的时均矢量；S_i(x,y)为当前时刻的瞬时矢 量；错误矢量标记为“NULL”，以区别于零矢量；N_i(x,y)为参与当前矢量平均的正确矢量 个数；F_i(x，y)分别为此前时均矢量的类型标志和当前瞬时矢量的类型标志，对二 者进行或运算得到当前时均矢量的类型标志，如下：

${\bar{F}}_i(x,y)={\bar{F}}_{i-1}(x,y)|F_i(x,y)---\left(4\right)$

其中，“|”表示或运算符。遍历每个所有IA后统计得到正确和错误矢量的个数N₁和N₀，则 时均流场的矢量正确率为：

${\bar{P}}_i=N_1/(N_1+N_0)\times 100\%---\left(5\right)$

若小于一个阈值P_TH（如90%），则继续处理下一时刻的图像对，直到大于P_TH并且求平 均的帧数达到预设的阈值，得到一个包含错误矢量的时均流场和对应的矢量类型查找表；反 之则认为时均流场重建失效。

（6）错误矢量插值，采用错误矢量邻域内正确矢量的均值对错误矢量进行插值：

${\bar{S}}_0(x,y)=\frac{1}{N_1(x,y)}\underset{n=-1}{\overset{1}{\Sigma }}\underset{m=-1}{\overset{1}{\Sigma }}\bar{F}(x+n,y+m)·\bar{S}(x+n,y+m)---\left(6\right)$

其中，N₁(x,y)表示错误矢量邻域内正确矢量的个数，m、n表示相邻矢量的偏移量， 对于3×3邻域有-1≤m≤1、-1≤n≤1。由于矢量插值仅利用邻域内的正确矢量，避免了错 误矢量带来的误差。

（7）流场平滑滤波，采用时均流场和滤波器的空域卷积实现：

$\tilde{S}(x,y)=\bar{S}(x,y)*G(x,y)---\left(7\right)$

其中，“*”表示卷积运算，滤波器采用大小为3×3的归一化二维离散高斯窗，表示为：

$G(x,y)=\frac{1}{16}\left(\begin{array}{ccc}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right)---\left(8\right)$

至此完成一次河流水面二维时均流场重建过程。