CN103996171A - 基于时空图像的流体运动矢量估计方法 - Google Patents

Abstract

本发明公开了一种基于时空图像的流体运动矢量估计方法，属于流场测量技术领域。所述方法首先采用单台相机以固定时间间隔采集流体图像序列；然后采用像差修正的直接线性变换法对图像进行光学畸变校正，并沿测流方向设置一组测速线，每条测速线均以像素为横坐标、以时间为纵坐标合成一幅时空图像；接下来对边缘检测后的时空图像进行频谱分析，并采用离散逼近法求解频谱主方向得到与其垂直的纹理主方向；最后根据物方坐标和像方坐标间的映射关系求解测速线的物理尺度，并结合纹理主方向和时间间隔计算一维流速矢量的大小，实现时均流场定标。本发明能够用于河流水面等天然示踪水流的时均流场测量，特别适用于高洪期河流流量的在线监测。

Description

基于时空图像的流体运动矢量估计方法

技术领域

本发明涉及一种流体运动矢量估计方法，尤其涉及一种基于时空图像的流体运动矢量估计方法，属于流场测量技术领域。

背景技术

流体中存在着大量不具有固定质量并存在相对运动的流体微团，使得流体的运动形态比刚体和固体更为复杂。在过去的二十多年间，光学、电子及计算机技术的发展促进了粒子图像测速(PIV)这种非接触式瞬时全场流速测量技术的实现和应用。它通过粒子图像序列的分析和计算获得局部流体运动位移及速度的大小、方向、特征和分布情况，极大提高了实验室环境下各种复杂流动的测量能力。其中从粒子图像估计流体的运动矢量是PIV技术的核心和难点，方法的选择不仅依赖于硬件系统而且取决于待测流体的特点。

对于利用树叶、树枝等天然漂浮物以及漩滚、表面波等天然水面模式作为水流示踪物的大尺度粒子图像测速(LSPIV)应用而言，河流水面流场的运动矢量估计又面临着特有的挑战。首先，大多数天然水面模式产生于明渠紊流，它们没有明确的粒子特性，通常存在振荡、扩散及形变等非刚性运动，因而缺乏稳定的几何特征和统计特征，使得经典的粒子跟踪测速(PTV)方法难以直接适用。其次，现场条件下的待测流场面积往往覆盖数百至数千平方米，为满足流速测量的目标精度，图像分辨率需要达到百万至千万像素，远大于目前实验室中十万像素的数量级，使得基于直接互相关(DCC)等空域相关匹配的运动矢量估计方法计算量相当可观，难以满足流场、流量连续监测的实时性要求。再则，河流水面的流速分布在断面方向表现出从中泓区域向两岸递减的趋势，流速梯度加之图像透视畸变和水流示踪物分布不均的影响往往导致水面瞬时流场具有复杂的区域特征。采用固定大小的观测窗口会因为对局部流场的适应性较差而降低运动矢量估计的精度和可靠性。因此研究针对河流水面这类特殊流体的运动矢量估计方法具有重要的理论意义和应用价值。

发明内容

本发明针对现有方法存在的不足，提供了一种基于时空图像的流体运动矢量估计方法。

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

由于满足质量守恒定律的流体，其运动在短时内通常满足连续性的假设。对于跟随性良好的水流示踪物，其在三维时空域中的位置必然满足某种相关性。这种相关性在一维图像空间和一维序列时间组成的时空图像中表现为较为显著的纹理主方向，如果能够精确估计出时空图像的纹理主方向就能获得流体在指定空间方向的运动矢量。然而从空间域中可以看出，对于天然示踪的流体，水流示踪物的时空分布不均往往造成时空图像受到非规则纹理或噪声的干扰，影响纹理主方向的估计精度。图像的边缘极大程度上反应了图像纹理的方向，并且具有方向性规则纹理的图像在傅立叶频谱中的能量主要聚集在过原点的直线上，直线指向代表的频谱主方向与纹理主方向垂直。因此利用时空图像边缘信息的频谱整体估计纹理主方向可以有效抑制背景噪声，提高估计精度和运算效率。

基于上述原理，方法首先采用单台相机以固定时间间隔采集流体图像序列；然后采用像差修正的直接线性变换法对图像进行光学畸变校正，并沿测流方向设置一组测速线，每条测速线均以像素为横坐标、以时间为纵坐标合成一幅时空图像；接下来对边缘检测后的时空图像进行频谱分析，并采用离散逼近法求解频谱主方向得到与其垂直的纹理主方向；最后根据物方坐标和像方坐标间的映射关系求解测速线的物理尺度，并结合纹理主方向和时间间隔计算一维流速矢量的大小，实现时均流场定标。

所述的流体图像序列为顺序采集的M幅8位灰度图像，M为单次时均流场测量的时长T和每帧图像采集的时间间隔Δt之比；

所述的一组测速线由L条测速线l_u组成，i＝1，2，…，M表示图像在序列中的编号，j＝1，2，…，L表示测速线在图像中的编号，测速线的宽度为1个像素，长度N不小于水流示踪物在连续图像间运动位移的2倍；

所述的时空图像大小为M×N像素，共有L幅；

所述的边缘频谱分析，首先采用Canny算子对时空图像进行边缘检测得到二值化的边缘图像，然后依次进行二维离散傅立叶变换、中心平移、对数变换和灰度拉伸得到直角坐标系下的频谱图像；

所述的频谱主方向定义为经过频谱图像中心的直线中线积分峰值对应的直线角度；

所述的离散逼近法包括如下步骤：

1)以点(M/2，N/2)为原点、M和N中的短边为直径建立极坐标系；

2)在极坐标系中设置一条过原点且半径为R的搜索线，对于固定的角度γ计算线积分：

$F\left(\mathrm{\&gamma;}\right)=\underset{r=1}{\overset{R}{\mathrm{\&Sigma;}}}F(r,\mathrm{\&gamma;}),R=\left\{\begin{array}{cc}M/2,& M<N\\ N/2,& M\&GreaterEqual;N\end{array}---\left(1\right)\right.$

3)在0～180°的半圆空间内搜索F(γ)的最大值，首先以1°为步进进行粗定位，再以0.1°为步进在最大值附近±0.5°的范围内进行精搜索；

所述的物方坐标和像方坐标间的映射关系采用像差修正的直接线性变换公式描述如下：

$\left\{\begin{array}{c}x+{\mathrm{\&delta;}}_x=\frac{l_{1}X+l_{2}Y+l_{3}Z+l_4}{l_{9}X+l_{10}Y+l_{11}Z+1}\\ y+{\mathrm{\&delta;}}_y=\frac{l_{5}X+l_{6}Y+l_{7}Z+l_8}{l_{9}X+l_{10}Y+l_{11}Z+1}\end{array}\right.---\left(2\right)$

其中，(X，Y，Z)散表示物方坐标，(x，y)表示像方坐标，11个未知量l_i称为直接线性变换系数，光学系统的总像差(δ_x，δ_y)通过以下三种畸变分量的线性叠加来建模：

$\left\{\begin{array}{c}{\mathrm{\&delta;}}_x=\tilde{x}-x=k_1{x}_d(x_d^2+y_d^2)+k_2{x}_d{(x_d^2+y_d^2)}^2+p_1(3x_d^2+y_d^2)+2p_2{x}_d{y}_d+s_1(x_d^2+y_d^2)\\ {\mathrm{\&delta;}}_y=\tilde{y}-y=k_1{y}_d(x_d^2+y_d^2)+k_2{y}_d{(x_d^2+y_d^2)}^2+2p_1{x}_d{y}_d+p_2(x_d^2+3y_d^2)+s_2(x_d^2+y_d^2)\end{array}\right.---\left(3\right)$

其中，x_d、y_d分别为像点辐射距的水平和垂直分量，k₁、k₂为径向畸变系数，p₁、p₂为离心畸变系数，s₁、s₂为薄棱镜畸变系数。

所述的时均流场定标方法，首先将测速线上首尾像素的物理距离除以像素距离得到单位物理尺度s_j

$s_j=\frac{\sqrt{{(X_N-X_1)}^2+{(Y_N-Y_1)}^2}}{\sqrt{{(x_N-x_1)}^2+{(y_N-y_1)}^2}}---\left(4\right)$

接下来根据时空图像中位移、时间和流速间的三角关系，对时均流速矢量u_j进行定标：

$u_j=\frac{s_j}{\mathrm{\&Delta;t}}\tan {\mathrm{\&theta;}}_j---\left(5\right)$

反复执行边缘频谱分析及时均流场定标的步骤直到读入全部L幅时空图像，最后得到时间T内测得的一维时均流速场。

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

1、鲁棒性强。利用水流示踪物时空分布的连续性从时空图像中直接估计出一维时均流场，这种时空联合的整体估计方法可以有效克服示踪物非刚性运动引起的不确定性。

2、分辨率高。观测窗口为单像素的测速线，相比空域相关匹配法中的矩形窗口，有效改善了相机在倾斜视角下拍摄引起的图像远场空间分辨率降低，使测量结果更具有代表性。

3、实时性强。方法将流体一维运动矢量估计的问题转化为时空图像频谱主方向的检测问题来解决，由于无需进行穷尽搜索的空域相关运算，计算量仅为传统PIV方法的十分之一。

鉴于以上特点，本方法能够用于河流水面等天然示踪水流的时均流场测量，特别适用于高洪期河流流量的在线监测。

附图说明

图1是本发明涉及的运动矢量估计方法的工作流程图。

具体实施方式

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

本发明的运动矢量估计方法的工作流程如图1所示。共分为以下四个主要步骤：

(1)图像序列采集

首先设置单次时均流场测量的时长T和每帧图像采集的时间间隔Δt。T应当足够大以消除流速脉动效应引起的随机测量误差，同时应当足够小以满足测量的时间分辨率；Δt根据拍摄视场、图像分辨率和待测流进行选取，应当足够大以保证结果的离散化精度，同时应当足够小以保证水流示踪物在图像序列中的分布满足连续性。然后计算图像序列的长度M=T/Δt。最后采用单台相机以时间间隔Δt连续采集M幅8位灰度图像GRI_i，i=1，2，…，M表示图像在序列中的编号。

(2)时空图像合成

首先在图像中沿测速方向设置L条测速线l_ij，j=1，2，…，L表示测速线在图像中的编号；测速线的宽度为1个像素，长度N不小于水流示踪物在连续图像间运动位移的2倍。然后以图像序列的长度M为纵坐标、以测速线的长度N为横坐标，建立L幅大小为M×N像素的时空图像STI_i。接下来依次读入一幅采集的灰度图像GRI_i，并利用现场布设的控制点采用像差修正的直接线性变换法对图像进行光学畸变校正。最后将无失真图像GRI_i中每条测速线l_i上的像素载入对应的时空图像STI_i直到读入M幅灰度图像，得到合成的L幅时空图像。

(3)边缘频谱分析

首先依次读入一幅合成的时空图像STI_i，采用Canny算子对其进行边缘检测，得到二值化的边缘图像EDI_i。然后对EDI_i依次进行二维离散傅立叶变换、中心平移、对数变换和灰度拉伸得到直角坐标系下的频谱图像FSI_i接下来将过中心直线中FSI_i线积分峰值对应的直线角度定义为频谱主方向γ_i，采用离散逼近法搜索γ_i，包括如下步骤：

1)以点(M/2，N/2)为原点、M和N中的短边为直径建立极坐标系；

2)在极坐标系中设置一条过原点且半径为R的搜索线，对于固定的角度γ计算线积分：

$F\left(\mathrm{\&gamma;}\right)=\underset{r=1}{\overset{R}{\mathrm{\&Sigma;}}}F(r,\mathrm{\&gamma;}),R=\left\{\begin{array}{cc}M/2,& M<N\\ N/2,& M\&GreaterEqual;N\end{array}---\left(1\right)\right.$

3)在0～180°的半圆空间内搜索F(γ)的最大值，首先以1°为步进进行粗定位，再以0.1°为步进在最大值附近±0.5°的范围内进行精搜索。

最后根据纹理主方向与频谱主方向垂直的关系计算时空图像STI_i的纹理主方向θ_i，如下：

(4)时均流场定标

首先根据以下物方坐标(X，Y，Z)和像方坐标(x，y)间的映射关系计算测速线l_i上的首尾像素(x₁，y₁)和(x_N，y_N)的物方坐标(X₁，Y₁，Z₁)和(X_N，Y_N，Z_N)：

$\left\{\begin{array}{c}x+{\mathrm{\&delta;}}_x=\frac{l_{1}X+l_{2}Y+l_{3}Z+l_4}{l_{9}X+l_{10}Y+l_{11}Z+1}\\ y+{\mathrm{\&delta;}}_y=\frac{l_{5}X+l_{6}Y+l_{7}Z+l_8}{l_{9}X+l_{10}Y+l_{11}Z+1}\end{array}\right.---\left(3\right)$

其中，11个未知量l_i称为直接线性变换系数，光学系统的总像差(δ_x，δ_y)通过以下三种畸变分量的线性叠加来建模：

$\left\{\begin{array}{c}{\mathrm{\&delta;}}_x=\tilde{x}-x=k_1{x}_d(x_d^2+y_d^2)+k_2{x}_d{(x_d^2+y_d^2)}^2+p_1(3x_d^2+y_d^2)+2p_2{x}_d{y}_d+s_1(x_d^2+y_d^2)\\ {\mathrm{\&delta;}}_y=\tilde{y}-y=k_1{y}_d(x_d^2+y_d^2)+k_2{y}_d{(x_d^2+y_d^2)}^2+2p_1{x}_d{y}_d+p_2(x_d^2+3y_d^2)+s_2(x_d^2+y_d^2)\end{array}\right.---\left(4\right)$

其中，x_d、y_d分别为像点辐射距的水平和垂直分量，k₁、k₂为径向畸变系数，p₁、p₂为离心畸变系数，s₁、s₂为薄棱镜畸变系数。

然后将测速线上首尾像素的物理距离除以像素距离得到单位物理尺度s_j：

$s_j=\frac{\sqrt{{(X_N-X_1)}^2+{(Y_N-Y_1)}^2}}{\sqrt{{(x_N-x_1)}^2+{(y_N-y_1)}^2}}---\left(5\right)$

接下来根据时空图像中位移、时间和流速间的三角关系，对时均流速矢量u_j进行定标：

$u_j=\frac{s_j}{\mathrm{\&Delta;t}}\tan {\mathrm{\&theta;}}_j---\left(6\right)$

反复执行步骤(3)和(4)直到读入全部L幅时空图像，最后得到时间T内测得的一维时均流速场。

Claims (5)

1.一种基于时空图像的流体运动矢量估计方法，分为图像序列采集、时空图像合成、边缘频谱分析及时均流场定标四个步骤,其特征在于：方法首先采用单台相机以固定时间间隔采集流体图像序列；然后采用像差修正的直接线性变换法对图像进行光学畸变校正，并沿测流方向设置一组测速线，每条测速线均以像素为横坐标、以时间为纵坐标合成一幅时空图像；接下来对边缘检测后的时空图像进行频谱分析，并采用离散逼近法求解频谱主方向得到与其垂直的纹理主方向；最后根据物方坐标和像方坐标间的映射关系求解测速线的物理尺度，并结合纹理主方向和时间间隔计算一维流速矢量的大小，实现时均流场定标；

所述的流体图像序列为顺序采集的M幅8位灰度图像，M为单次时均流场测量的时长T和每帧图像采集的时间间隔Δt之比；

所述的一组测速线由L条测速线l_u组成，i＝1，2，…，M表示图像在序列中的编号，j=1，2，…，L表示测速线在图像中的编号，测速线的宽度为1个像素，长度N不小于水流示踪物在连续图像间运动位移的2倍；

所述的时空图像大小为M×N像素，共有L幅；

所述的频谱主方向定义为经过频谱图像中心的直线中线积分峰值对应的直线角度。

2.根据权利要求1所述的基于时空图像的流体运动矢量估计方法，其特征在于：所述的边缘频谱分析首先采用Canny算子对时空图像进行边缘检测得到二值化的边缘图像，然后依次进行二维离散傅立叶变换、中心平移、对数变换和灰度拉伸得到直角坐标系下的频谱图像。

3.根据权利要求1所述的基于时空图像的流体运动矢量估计方法，其特征在于：所述的离散逼近法包括如下步骤：

1)以点(M/2，N/2)为原点、M和N中的短边为直径建立极坐标系；

2)在极坐标系中设置一条过原点且半径为R的搜索线，对于固定的角度γ计算线积分：一的搜索线，对于固定的角度计算线积分：

$F\left(\mathrm{\&gamma;}\right)=\underset{r=1}{\overset{R}{\mathrm{\&Sigma;}}}F(r,\mathrm{\&gamma;}),R=\left\{\begin{array}{cc}M/2,& M<N\\ N/2,& M\&GreaterEqual;N\end{array}---\left(1\right)\right.$

3)在0～180°的半圆空间内搜索F(γ)的最大值，首先以1°为步进进行粗定位，再以0.1°为步进在最大值附近±0.5°的范围内进行精搜索。

4.根据权利要求1所述的基于时空图像的流体运动矢量估计方法，其特征在于：所述的物方坐标和像方坐标间的映射关系采用像差修正的直接线性变换公式描述如下：

$\left\{\begin{array}{c}x+{\mathrm{\&delta;}}_x=\frac{l_{1}X+l_{2}Y+l_{3}Z+l_4}{l_{9}X+l_{10}Y+l_{11}Z+1}\\ y+{\mathrm{\&delta;}}_y=\frac{l_{5}X+l_{6}Y+l_{7}Z+l_8}{l_{9}X+l_{10}Y+l_{11}Z+1}\end{array}\right.---\left(2\right)$

其中，(X，Y，Z)表示物方坐标，(x，y)表示像方坐标，11个未知量l_i称为直接线性变换系数，光学系统的总像差(δ_x，δ_y)通过以下三种畸变分量的线性叠加来建模：

$\left\{\begin{array}{c}{\mathrm{\&delta;}}_x=\tilde{x}-x=k_1{x}_d(x_d^2+y_d^2)+k_2{x}_d{(x_d^2+y_d^2)}^2+p_1(3x_d^2+y_d^2)+2p_2{x}_d{y}_d+s_1(x_d^2+y_d^2)\\ {\mathrm{\&delta;}}_y=\tilde{y}-y=k_1{y}_d(x_d^2+y_d^2)+k_2{y}_d{(x_d^2+y_d^2)}^2+2p_1{x}_d{y}_d+p_2(x_d^2+3y_d^2)+s_2(x_d^2+y_d^2)\end{array}\right.---\left(3\right)$

其中，x_d、y_d分别为像点辐射距的水平和垂直分量，k₁、k₂为径向畸变系数，p₁、p₂为离心畸变系数，s₁、s₂为薄棱镜畸变系数。

5.根据权利要求1所述的基于时空图像的流体运动矢量估计方法，其特征在于：所述的时均流场定标方法为首先将测速线上首尾像素的物理距离除以像素距离得到单位物理尺度s_j：

$s_j=\frac{\sqrt{{(X_N-X_1)}^2+{(Y_N-Y_1)}^2}}{\sqrt{{(x_N-x_1)}^2+{(y_N-y_1)}^2}}---\left(4\right)$

接下来根据时空图像中位移、时间和流速间的三角关系，对时均流速矢量u_j进行定标：

$u_j=\frac{s_j}{\mathrm{\&Delta;t}}\tan {\mathrm{\&theta;}}_j---\left(5\right)$

反复执行边缘频谱分析及时均流场定标的步骤直到读入全部L幅时空图像，最后得到时间T内测得的一维时均流速场。