CN113804916B - 一种基于最大流速先验信息的频域时空图像测速法 - Google Patents

一种基于最大流速先验信息的频域时空图像测速法 Download PDF

Info

Publication number
CN113804916B
CN113804916B CN202111096392.2A CN202111096392A CN113804916B CN 113804916 B CN113804916 B CN 113804916B CN 202111096392 A CN202111096392 A CN 202111096392A CN 113804916 B CN113804916 B CN 113804916B
Authority
CN
China
Prior art keywords
image
camera
value
line
velocity
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202111096392.2A
Other languages
English (en)
Other versions
CN113804916A (zh
Inventor
张振
董睿
袁章
周冰倩
杨家伟
张婷
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN202111096392.2A priority Critical patent/CN113804916B/zh
Publication of CN113804916A publication Critical patent/CN113804916A/zh
Application granted granted Critical
Publication of CN113804916B publication Critical patent/CN113804916B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01PMEASURING LINEAR OR ANGULAR SPEED, ACCELERATION, DECELERATION, OR SHOCK; INDICATING PRESENCE, ABSENCE, OR DIRECTION, OF MOVEMENT
    • G01P5/00Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft
    • G01P5/18Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft by measuring the time taken to traverse a fixed distance
    • G01P5/20Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft by measuring the time taken to traverse a fixed distance using particles entrained by a fluid stream

Landscapes

  • Engineering & Computer Science (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种基于最大流速先验信息的频域时空图像测速方法,包括:假设水深与表面流速成正比例关系,计算断面最大水深和测速线水深;利用当前水位下断面最大流速先验值计算测速线的最大流速先验值;计算测速线的尺度变换因子,并根据频域时空图像测速法的逆变换确定时空图像频谱主方向的搜索区间;在搜索区间内搜索幅度谱极坐标投影的极值,将对应角度作为测速线的频域主方向估计值,并根据频域时空图像测速法的正变换计算测速线的流速值。本发明可以有效滤除野外复杂光照、水流及气象条件下倒影、耀光、紊流、雨雾、遮挡等环境噪声对时空图像测速法的干扰,从而抑制测量粗大误差,显著改善低流速及高流速测量的稳定性和可靠性。

Description

一种基于最大流速先验信息的频域时空图像测速法
技术领域
本发明属于图像法测流技术领域,涉及一种时空图像测速法,尤其涉及一种基于最大流速先验信息的频域时空图像测速法。
背景技术
时空图像测速方法(STIV)是一种以测速线为分析区域、通过检测合成时空图像的纹理主方向估计一维时均流速测量方法。它利用水流示踪物在三维时空域中运动的连续性,采用平行于顺流方向的测速线作为分析区域,在图像空间和序列时间组成的时空图像中检测和示踪物运动相关的纹理方向特征,直接估计指定空间方向上的时均运动矢量。由于具有空间分辨率高、实时性强的优点,在河流水面流速、流量的实时监测中具有特别的应用潜力。
频域时空图像测速法,通过将空域复杂的纹理主方向检测转换为频域中搜索图像频谱主方向的线性运算,能显著提高算法抗噪性,降低复杂度。然而在实际应用中,阴影、耀光、倒影、障碍物遮挡、风浪、雨雾等环境干扰因素容易导致背景噪声过大,时空图像幅度谱的信噪比偏低,使得时空图像纹理主方向的检测出现粗大误差,从而引起流速测量失效。
所以,需要一个新的技术方案来解决这个问题。
发明内容
发明目的:为了克服现有技术中存在的不足,提供一种基于最大流速先验信息的频域时空图像测速法,能够根据每条测速线的实际流速范围确定最佳搜索区间,进而抑制不相关噪声对有用信号检测的干扰。
技术方案:为实现上述目的,本发明提供一种基于最大流速先验信息的频域时空图像测速方法,包括如下步骤:
S1:假设水深与表面流速成正比例关系,根据已知的断面地形和水位信息计算断面最大水深和测速线水深;
S2:基于断面最大水深和测速线水深,利用当前水位下断面最大流速先验值计算测速线的最大流速先验值;
S3:基于变高水面摄影测量模型计算测速线的尺度变换因子,并根据频域时空图像测速法的逆变换确定时空图像频谱主方向的搜索区间;
S4:在搜索区间内搜索幅度谱极坐标投影的极值,将对应角度作为测速线的频域主方向估计值,并根据频域时空图像测速法的正变换计算测速线的流速值。
进一步地,所述步骤S1中测速线水深的计算方式为:
已知断面任意两点地形(x0,h0)、(x1,h1),根据线性地形插函数得到两点间单条测速线上的断面地形hl
Figure BDA0003267386900000021
以地面为基准点的水位值B和单点地形插值hl计算得到测速线水深Dl
Dl=B-hl (2)
进一步地,所述步骤S2中测速线的最大流速先验值的计算方式为:
Figure BDA0003267386900000022
其中,Vmax为测速线的最大流速先验值;Dsmax为断面最大水深,断面最大水深为实测的已知值,通过实际输入的水位值和断面最低点相减得到;Vfmax为当前水位下断面最大流速先验值,不同水位下的断面最大流速先验值根据经验或实测值确定。
进一步地,所述步骤S3中测速线的尺度变换因子的计算公式为:
Figure BDA0003267386900000023
其中,j、n、f分别为图片尺寸、像素大小、相机焦距,可以通过相机的内参标定得到,图像传感器的像元尺寸s可由内参标定结果计算得到;相机俯仰角α通过相机的外参标定校准;摄像机水面高程H等于水位值B和摄像机距地面高度Hc之和:
H=Hc+B (5)
进一步地,所述步骤S3中基于变高水面摄影测量模型计算测速线的尺度变换因子的具体过程为:
A1:在实验室对摄像机进行内参标定:
在实验室中标定相机内参矩阵K和畸变参数矩阵D如下:
Figure BDA0003267386900000024
D=[k1 k2 p1 p2] (7)
其中(Cx,Cy)表示畸变图像的像主点坐标,Cy表示像主点的纵坐标,fx、fy分别表示摄像机在像平面x轴和y轴方向上的等效焦距,k1表示第一径向畸变系数,p1表示第一切向畸变系数,k2表示第二径向畸变系数,p2表示第二切向畸变系数,根据图像传感器的像元尺寸s得到相机焦距f:
f=(fx+fy)·s/2 (8)
畸变参数用于图像的非线性畸变校正:
Figure BDA0003267386900000031
其中,(x′,y′)和(x,y)分别为畸变和无畸变的相机坐标,它们与对应的图像坐标(u′,v′)和(u,v)间满足:
Figure BDA0003267386900000032
Figure BDA0003267386900000033
以上三式建立了无畸变图像坐标到畸变图像坐标的变换关系;
A2:进行相机的外参标定:
将摄像机架设到岸边,读取集成在摄像机上的三轴传感器数值得到俯仰角α=,横滚角ω,由于相机和激光测距仪的中心并不完全重合,畸变校正后需要继续进行偏心角的检校,偏心角定义为相机和测距仪的三个旋转角,方位角置零,检校俯仰角α以及横滚角ω:
Figure BDA0003267386900000034
其中,EX、EY表示所有网格在X和Y方向尺寸的均方根误差:
Figure BDA0003267386900000035
M、L表示棋盘格内角点的数量,(Xi,j,Yi,j)表示其物方坐标,DX、DY表示网格尺寸的真值;
A3:建立倾斜视角下的中心透视投影成像模型:
像平面和物平面坐标系分别用(x,y),(X,Y)表示;0为透镜平面的光心,o、O′分别为其在像平面和物平面上的投影点;c为像平面延长线和通过光心的水平线的交点;H为光心到物平面间的垂直距离,C为对应的垂足点;相机的俯仰角α定义为相机主光轴和物平面间的夹角;s表示图像传感器的像元尺寸;f为相机焦距;(m,n)表示图像的大小;下标(i,j)表示像素的坐标;
像平面坐标(x,y)与图像坐标(i,j)的关系可由以下公式表示:
Figure BDA0003267386900000041
对成像模型进行合理地简化:实际布设岸基式系统时,可以分别通过对齐测量断面及调节云台水平将方位角和横滚角的值置零,因此可以仅考虑相机主光轴平行于断面方向仅存在俯仰角的情况;
物像尺度变换因子S可由x方向或y方向对应换算关系计算得出,这里以x方向为例计算;在x方向,假设像素pi,j位于像主点o的左侧,pi,j及其相邻像素pi+1,j在物平面上的投影点分别为Pi,j和Pi+1,j,射线Pi,jO和Pi+1,jO与投影线PjO的夹角分别为
Figure 1
和φ,对于pi,j和pi+1,j,分别有:
Figure BDA0003267386900000042
Figure BDA0003267386900000043
在y方向,假设像素pi,j位于像主点o的下方,pi,j及其相邻像素pi,j+1在物平面主纵线上的投影点分别为Pj和Pj+1,射线PjO和Pj+1O与物平面Pj+1C的夹角分别为β和γ,对于pi,j,其像平面主纵线上的投影点力满足以下三角关系:
Figure BDA0003267386900000044
由于α=∠cOo、β=∠cOpj,代入上式得:
Figure BDA0003267386900000045
将像素pi,j在x方向上的二维物像尺度因子S用其物点Pi,j和相邻像素对应物点间的距离来描述,即:
Figure BDA0003267386900000046
结合以上各式,x方向的物像尺度变换因子S可由以下公式表示:
Figure BDA0003267386900000047
其中,j、n、s、f分别为图片尺寸,像素大小,像元尺寸,相机焦距,是上述相机内参标定求出的数值;相机俯仰角α由上述相机的外参标定校准;摄像机水面高程H等于水位值B和摄像机距地面高度HC之和:
H=Hc+B (21)
进一步地,所述步骤S3中时空图像频谱主方向的搜索区间的确定方法为:
所述搜索区间的取值和光流速度的方向有关,当光流速度为正时,频域搜索区间为 [90°+θmin,θmax];当光流速度为负时,频域搜索区间为[180°-θmax,90°-θmin];
其中,θmax为上边界,用于滤除时空图像中风浪、紊流等水平方向的背景噪声,采用以下公式求出:
Figure BDA0003267386900000051
θmin为下边界,用于滤除时空图像中阴影、耀光、遮挡等垂直方向的背景噪声,取值为定值。
进一步地,所述步骤S4具体为:
所述幅度谱极坐标投影的极值是指将傅里叶变换后的频谱信息映射到极坐标下,沿角度及半径方向投影遍历搜索区间,计算出投影的极值;
以图像中心为原点O、N/2为半径设置搜索线,以Δθ=0.1°为步进、0~180°为区间建立极坐标系,统计搜索测速线上的幅度分布F(t,θ),按下式对幅度谱的平方项进行累加得到能量-角度分布:
Figure BDA0003267386900000052
在搜索区间[θmin,θmax]搜索P(θ)极值,将对应的频域主方向θ代入以下频域时空图像测速法的正变换公式,可得该条测速线上的流速值V:
Figure BDA0003267386900000053
完成单条测速线的表面流速测量;Δt为时间间隔,计算公式为:
Figure BDA0003267386900000054
有益效果:本发明与现有技术相比,能够根据每条测速线的实际流速范围确定最佳搜索区间,进而抑制不相关噪声对有用信号检测的干扰,可以有效滤除野外复杂光照、水流及气象条件下倒影、耀光、紊流、雨雾、遮挡等环境噪声对时空图像测速法的干扰,从而抑制测量粗大误差,显著改善低流速及高流速测量的稳定性和可靠性。
附图说明
图1是本发明方法的流程图;
图2是测量系统布设示意图;
图3是断面插值地形示意图;
图4是变高平面摄影测量模型示意图,其中图4.1是x方向的立体视图,图4.2是剖面视图;
图5是噪声的能量区间图;
图6是设置搜索区间前后的搜索结果对比图,其中图6.1是受到干扰时的错误结果,图6.2是消除了粗大误差的正确结果;
图7是设置搜索区间前后的测速结果对比图,其中图7.1是使用本发明方法之前的测速结果,图7.2是使用了本发明方法之后的测速结果。
具体实施方式
下面结合附图和具体实施例,进一步阐明本发明,应理解这些实施例仅用于说明本发明而不用于限制本发明的范围,在阅读了本发明之后,本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
本发明提供一种基于最大流速先验信息的频域时空图像测速方法,如图1所示,其包括如下步骤:
S1:假设水深与表面流速成正比例关系,根据已知的断面地形和水位信息计算断面最大水深和测速线水深:
S2:基于断面最大水深和测速线水深,利用当前水位下断面最大流速先验值计算测速线的最大流速先验值:
S3:基于变高水面摄影测量模型计算测速线的尺度变换因子,并根据频域时空图像测速法的逆变换确定时空图像频谱主方向的搜索区间:
S4:在搜索区间内搜索幅度谱极坐标投影的极值,将对应角度作为测速线的频域主方向估计值,并根据频域时空图像测速法的正变换计算测速线的流速值。
本实施例中以某水文站为例,将上述方法进行实例应用,具体如下:
(1)参照图3,以海平面为基准,已知断面任意两点地形(x0,h0)、(x1,h1),根据线性地形插函数得到两点间单条测速线上的断面地形hl
Figure BDA0003267386900000061
以海平面为基准点的水位值B和单点地形插值hl计算得到测速线水深Dl
Dl=B-hl (27)
根据断面最大水深Dsmax和测速线水深Dl,建立比例关系,将断面最大流速值Vfmax转化为测速线的实际流速最大值Vmax
Figure BDA0003267386900000062
本实施例中按以上思路,已知断面地形由水文站提供,以海平面为基准,每隔五米测量一次地形。取(85.0,983.10),(90.0,983.00),插值地形为983.02m。输入水位为989.04m,断面最低点为982.60m,最大水深为6.44m,测速线水深为6.02m。此水位值对应水位级下断面最大流速为3.5m/s,计算可得测速线最大流速为3.27m/s。
(2)通过视频帧速率和采样帧间隔计算时间间隔Δt:
Figure BDA0003267386900000071
本实施例中实际测量时采样帧间隔为1,帧速率为25fps,Δt取0.04s。
(3)在实验室对摄像机进行内参标定:
相机的成像原理建立在多个坐标系的变换之上,即摄像机坐标系、图像物理坐坐标系、图像像素坐标系和世界坐标系。根据这四个坐标系之间的关系,现实场景下的任意一点都可变换至其中的任一坐标系中。根据这四个坐标系之间的关系,现实场景下的任意一点都可变换至其中的任一坐标系中,具体如下:
Step 1:世界坐标系中的一点通过旋转和平移变换转换至摄像机坐标系;
Step 2:通过投影变换可将Step 1中得到的相机坐标点变换为图像物理坐标系中的像点p’;
Step 3:由于图像物理坐标系和图像像素坐标系同处一个平面,像点p’可通过缩放和平移变换最终转换为像素坐标系中的点p(μ,v)。
以上坐标系的转换过程可通过一个矩阵来表示:
Figure BDA0003267386900000072
其中:s为像元尺寸,
Figure BDA0003267386900000073
为旋转平移变换矩阵,
Figure BDA0003267386900000074
为投影变换矩阵,f为焦距,
Figure BDA0003267386900000075
为缩放平移变换矩阵。
Figure BDA0003267386900000076
分别表示单位距离上横、纵坐标的像素个数。
标定后得到像元尺寸s=0.0013,内参矩阵K和畸变参数矩阵D如下:
Figure BDA0003267386900000077
D=[k1 k2 p1 p2]=[-0.4074 0.2169 -0.0007 -0.0004] (32)
其中,Cx表示像主点的横坐标,Cy表示像主点的纵坐标,fx表示摄像机在像平面x轴上的等效焦距,fy表示摄像机在像平面y轴上的等效焦距,相机焦距f根据图像传感器的像元尺寸s计算得到:
f=(fx+fy)·s/2=3.7447mm (33)
k1表示第一径向畸变系数,p1表示第一切向畸变系数,k2表示第二径向畸变系数,p2表示第二切向畸变系数,畸变参数用于图像的非线性畸变校正:
Figure BDA0003267386900000081
其中,(x′,y′)和(x,y)分别为畸变和无畸变的相机坐标,它们与对应的图像坐标(u′,v′) 和(u,v)间满足:
Figure BDA0003267386900000082
Figure BDA0003267386900000083
以上三式建立了无畸变图像坐标到畸变图像坐标的变换关系。
(4)将摄像机架设到岸边,如图2所示,读取集成在摄像机上的三轴传感器数值得到俯仰角α=15.718140,横滚角ω=1.026672。由于相机和激光测距仪的中心并不完全重合,畸变校正后需要继续进行偏心角的检校。偏心角定义为相机和测距仪的三个旋转角,方位角置零,检校俯仰角α以及横滚角ω。
Figure BDA0003267386900000084
其中,EX、EY表示所有网格在X和Y方向尺寸的均方根误差。
Figure BDA0003267386900000085
其中,M、L表示棋盘格内角点的数量,(Xi,j,Yi,j)表示其物方坐标,DX、DY表示网格尺寸的真值。
标定前可通过角点检测法检测角点,进而可对棋盘格每行每列相邻角点进行距离,得到未检校时角点距离均方根误差的最小值。
粗检校时可先控制俯仰角和横滚角在正负1°内变化,重复上述检测步骤,求出均方根误差最小时对应的俯仰角α和横滚角ω。
粗检校完成后进行进一步精确计算,以0.1°为基准,取粗检校所得角度相邻两点,通过三点高斯曲线拟合得到曲线最低点即为最终标定所得结果,最终计算得到俯仰角α=19.818140和横滚角ω=-0.009788。
(5)建立倾斜视角下的中心透视投影成像模型:像平面和物平面坐标系分别用(x,y),(X,Y)表示;O为透镜平面的光心,o、O′分别为其在像平面和物平面上的投影点;c 为像平面延长线和通过光心的水平线的交点;H为光心到物平面间的垂直距离,C为对应的垂足点;相机的俯仰角α定义为相机主光轴和物平面间的夹角;s表示图像传感器的像元尺寸;f为相机焦距;(m,n)表示图像的大小;下标(i,j)表示像素的坐标。
像平面坐标(x,y)与图像坐标(i,j)的关系可由以下公式表示:
Figure BDA0003267386900000091
对成像模型进行合理地简化:实际布设岸基式系统时,可以分别通过对齐测量断面及调节云台水平将方位角和横滚角的值置零,因此可以仅考虑相机主光轴平行于断面方向仅存在俯仰角的情况。
本实施例获取到如图4所示的变高平面摄影测量模型示意图,其中,图4.1是x方向的立体视图,图4.2是剖面视图。
物像尺度变换因子S可由x方向或y方向对应换算关系计算得出,本实施例以x方向为例计算。参照图4.1,在x方向,假设像素pi,j位于像主点o的左侧。pi,j及其相邻像素pi+1,j在物平面上的投影点分别为Pi,j和Pi+1,j,射线Pi,jO和Pi+1,jO与投影线 PjO的夹角分别为
Figure BDA0003267386900000092
和φ。对于pi,j和pi+1,j,分别有:
Figure BDA0003267386900000093
Figure BDA0003267386900000094
在y方向,假设像素pi,j位于像主点o的下方。pi,j及其相邻像素pi,j+1在物平面主纵线上的投影点分别为Pj和Pj+1,射线PjO和Pj+1O与物平面Pj+1C的夹角分别为β和γ。对于pi,j,其像平面主纵线上的投影点pj满足以下三角关系:
Figure BDA0003267386900000095
由于α=∠cOo、β=∠∠cOpj,代入上式得:
Figure BDA0003267386900000101
将像素pi,j在x方向上的二维物像尺度因子S用其物点Pi,j和相邻像素对应物点间的距离来描述,即:
Figure BDA0003267386900000102
结合以上各式,x方向的物像尺度变换因子S可由以下公式表示:
Figure BDA0003267386900000103
其中,j、n、s、f分别为图片尺寸,像素大小,像元尺寸,相机焦距,是上述相机内参标定求出的数值。相机俯仰角α由上述相机的外参标定校准。摄像机水面高程H等于水位值B和摄像机距地面高度HC之和:
H=Hc+B (46)
实际测量中摄像机水面高程计算得到18.760022m,代入计算S=0.13024。
(6)根据频域时空图像测速法的逆变换确定时空图像频谱主方向的搜索区间,上边界可由测速线流速先验值和尺度变换因子确定:
Figure BDA0003267386900000104
计算得到θmax约为176°。
对时空图像进行二维傅里叶变换,用对数幅值显示变换结果,将变换后的频谱中心从四角移到中心,一、三象限对调,二、四象限对调。此时直流分量集中到频谱中心。在幅度谱中对幅度谱的平方项进行累加,对幅度谱进行能量-角度分布计算,图5的结果显示低频噪声主要分布在90°±5°的范围内,所以θmin的取值为5°。
当光流速度为正,频谱主方向在90°~180°范围,此时搜索区间的下边界取值为90°+θmin,搜索区间为[90°+θmin,θmax];当光流速度为负,频谱主方向在0°~90°范围,此时搜索区间的下边界θmin的值为90°-θmin,搜索区间为[180°-θmax,90°-θmin]。
此时搜索区间为[95°,176°],为了直观的体现出效果,本实施例将设置搜索区间前后的结果进行对比,具体如图6所示,其中图6.1是受到干扰时的错误结果,图6.2是消除了粗大误差的正确结果,从图6的对比在中可明显看出上边界附近的误差被隔离在区间以外。
(7)以图像中心为原点O、N/2为半径设置搜索线,以Δθ=0.1°为步进、0°~180°为区间建立极坐标系,统计搜索测速线上的幅度分布F(r,θ),按下式对幅度谱的平方项进行累加得到能量-角度分布,按下式对幅度谱的平方项进行累加得到能量-角度分布:
Figure BDA0003267386900000111
在搜索区间[95°,176°]搜索P(θ)极值,将对应的频域主方向θ代入以下频域时空图像测速法的正变换公式,可得该条测速线上的流速值V:
Figure BDA0003267386900000112
为了验证本发明方法的效果,本实施例将现有频域时空图像测速法应用于相同水文站,测量到单条测速线的表面流速为2.396m/s,图7是设置搜索区间前后的测速结果对比图,其中,图7.1是使用本发明方法之前的测速结果,图7.2是使用了本发明方法之后的测速结果,经过对比可以明显看出本方法能够有效的消除远岸的粗大误差,并且准确识别出位于近岸的漩涡,使得测量结果更精确。

Claims (3)

1.一种基于最大流速先验信息的频域时空图像测速方法,其特征在于,包括如下步骤:
S1:假设水深与表面流速成正比例关系,根据已知的断面地形和水位信息计算断面最大水深和测速线水深;
S2:基于断面最大水深和测速线水深,利用当前水位下断面最大流速先验值计算测速线的最大流速先验值;
S3:基于变高水面摄影测量模型计算测速线的尺度变换因子,并根据频域时空图像测速法的逆变换确定时空图像频谱主方向的搜索区间;
S4:在搜索区间内搜索幅度谱极坐标投影的极值,将对应角度作为测速线的频域主方向估计值,并根据频域时空图像测速法的正变换计算测速线的流速值;
所述步骤S1中测速线水深的计算方式为:
已知断面任意两点地形(x0,h0)、(x1,h1),根据线性地形插值函数得到两点间单条测速线上的断面地形hl
Figure FDA0003783382080000011
以地面为基准点的水位值B和单点地形插值hl计算得到测速线水深Dl
Dl=B-hl (2);
所述步骤S2中测速线的最大流速先验值的计算方式为:
Figure FDA0003783382080000012
其中,Vmax为测速线的最大流速先验值,Dsmax为断面最大水深,Vfmax为当前水位下断面最大流速先验值;
所述步骤S3中测速线的尺度变换因子的计算公式为:
Figure FDA0003783382080000013
其中,j、n、f分别为图片尺寸、像素大小、相机焦距,通过相机的内参标定得到,图像传感器的像元尺寸s可由内参标定结果计算得到;相机俯仰角α通过相机的外参标定校准;摄像机水面高程H等于水位值B和摄像机距地面高度HC之和:
H=Hc+B (5);
所述步骤S3中基于变高水面摄影测量模型计算测速线的尺度变换因子的具体过程为:
A1:在实验室对摄像机进行内参标定:
在实验室中标定相机内参矩阵K和畸变参数矩阵D如下:
Figure FDA0003783382080000021
D=[k1 k2 p1 p2] (7)
其中(Cx,Cy)表示畸变图像的像主点坐标,Cy表示像主点的纵坐标,fx、fy分别表示摄像机在像平面x轴和y轴方向上的等效焦距,k1表示第一径向畸变系数,p1表示第一切向畸变系数,k2表示第二径向畸变系数,p2表示第二切向畸变系数,根据图像传感器的像元尺寸s得到相机焦距f:
f=(fx+fy)·s/2 (8)
畸变参数用于图像的非线性畸变校正:
Figure FDA0003783382080000022
其中,(x',y')和(x,y)分别为畸变和无畸变的相机坐标,它们与对应的图像坐标(u',v')和(u,v)间满足:
Figure FDA0003783382080000023
Figure FDA0003783382080000024
以上三式建立了无畸变图像坐标到畸变图像坐标的变换关系;
A2:进行相机的外参标定:
将摄像机架设到岸边,读取集成在摄像机上的三轴传感器数值得到俯仰角α,横滚角ω,由于相机和激光测距仪的中心并不完全重合,畸变校正后需要继续进行偏心角的检校,偏心角定义为相机和测距仪的三个旋转角,方位角置零,检校俯仰角α以及横滚角ω:
Figure FDA0003783382080000025
其中,EX、EY表示所有网格在X和Y方向尺寸的均方根误差:
Figure FDA0003783382080000031
M、L表示棋盘格内角点的数量,(Xi,j,Yi,j)表示其物方坐标,DX、DY表示网格尺寸的真值;
A3:建立倾斜视角下的中心透视投影成像模型:
像平面和物平面坐标系分别用(x,y),(X,Y)表示;O为透镜平面的光心,o、O'分别为其在像平面和物平面上的投影点;c为像平面延长线和通过光心的水平线的交点;H为光心到物平面间的垂直距离,C为对应的垂足点;相机的俯仰角α定义为相机主光轴和物平面间的夹角;s表示图像传感器的像元尺寸;f为相机焦距;(m,n)表示图像的大小;下标(i,j)表示像素的坐标;
像平面坐标(x,y)与图像坐标(i,j)的关系可由以下公式表示:
Figure FDA0003783382080000032
对成像模型进行简化;
物像尺度变换因子S可由x方向或y方向对应换算关系计算得出,这里以x方向为例计算;在x方向,假设像素pi,j位于像主点o的左侧,pi,j及其相邻像素pi+1,j在物平面上的投影点分别为Pi,j和Pi+1,j,射线Pi,jO和Pi+1,jO与投影线PjO的夹角分别为
Figure FDA0003783382080000037
和φ,对于pi,j和pi+1,j,分别有:
Figure FDA0003783382080000033
Figure FDA0003783382080000034
在y方向,假设像素pi,j位于像主点o的下方,pi,j及其相邻像素pi,j+1在物平面主纵线上的投影点分别为Pj和Pj+1,射线PjO和Pj+1O与物平面Pj+1C的夹角分别为β和γ,对于pi,j,其像平面主纵线上的投影点pj满足以下三角关系:
Figure FDA0003783382080000035
由于α=∠cOo、β=∠cOpj,代入上式得:
Figure FDA0003783382080000036
将像素pi,j在x方向上的二维物像尺度因子S用其物点Pi,j和相邻像素对应物点间的距离来描述,即:
Figure FDA0003783382080000041
结合以上各式,x方向的物像尺度变换因子S可由以下公式表示:
Figure FDA0003783382080000042
其中,j、n、s、f分别为图片尺寸,像素大小,像元尺寸,相机焦距,是上述相机内参标定求出的数值;相机俯仰角α由上述相机的外参标定校准;摄像机水面高程H等于水位值B和摄像机距地面高度HC之和:
H=Hc+B (21);
所述步骤S3中时空图像频谱主方向的搜索区间的确定方法为:
所述搜索区间的取值和光流速度的方向有关,当光流速度为正时,频域搜索区间为[90°+θminmax];当光流速度为负时,频域搜索区间为[180°-θmax,90°-θmin];
其中,θmax为上边界,用于滤除时空图像水平方向的背景噪声,采用以下公式求出:
Figure FDA0003783382080000043
θmin为下边界,用于滤除时空图像垂直方向的背景噪声,取值为定值。
2.根据权利要求1所述的一种基于最大流速先验信息的频域时空图像测速方法,其特征在于,所述步骤S4具体为:
以图像中心为原点O、N/2为半径设置搜索线,以Δθ=0.1°为步进、0~180°为区间建立极坐标系,统计搜索测速线上的幅度分布F(r,θ),按下式对幅度谱的平方项进行累加得到能量-角度分布:
Figure FDA0003783382080000044
在搜索区间[θminmax]搜索P(θ)极值,将对应的频域主方向θ代入以下频域时空图像测速法的正变换公式,可得该条测速线上的流速值V:
Figure FDA0003783382080000045
完成单条测速线的表面流速测量。
3.根据权利要求2所述的一种基于最大流速先验信息的频域时空图像测速方法,其特征在于,所述步骤S4中Δt为时间间隔,计算公式为:
Figure FDA0003783382080000051
CN202111096392.2A 2021-09-17 2021-09-17 一种基于最大流速先验信息的频域时空图像测速法 Active CN113804916B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111096392.2A CN113804916B (zh) 2021-09-17 2021-09-17 一种基于最大流速先验信息的频域时空图像测速法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111096392.2A CN113804916B (zh) 2021-09-17 2021-09-17 一种基于最大流速先验信息的频域时空图像测速法

Publications (2)

Publication Number Publication Date
CN113804916A CN113804916A (zh) 2021-12-17
CN113804916B true CN113804916B (zh) 2022-09-30

Family

ID=78895923

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111096392.2A Active CN113804916B (zh) 2021-09-17 2021-09-17 一种基于最大流速先验信息的频域时空图像测速法

Country Status (1)

Country Link
CN (1) CN113804916B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114898577B (zh) * 2022-07-13 2022-09-20 环球数科集团有限公司 一种用于高峰期道路管理的道路智能管理系统与方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008058264A (ja) * 2006-09-04 2008-03-13 Kobe Univ 実河川を対象とした流速観測装置、流速観測方法,および流速観測プログラム
CN104777327B (zh) * 2015-03-17 2018-03-20 河海大学 基于激光辅助标定的时空图像测速系统和方法
CN106092061B (zh) * 2016-05-31 2018-07-24 河海大学 基于倾斜视角下透镜成像模型的河流水面流场定标方法
CN107632168B (zh) * 2017-09-27 2020-03-31 杭州开闳环境科技有限公司 一种河道二维流速及流量测量系统及方法
CN110057295B (zh) * 2019-04-08 2020-12-25 河海大学 一种免像控的单目视觉平面距离测量方法
CN110135056B (zh) * 2019-05-14 2020-06-30 珠江水利委员会珠江水利科学研究院 一种河道内桥墩有效阻水比的快速自动分析方法

Also Published As

Publication number Publication date
CN113804916A (zh) 2021-12-17

Similar Documents

Publication Publication Date Title
CN105823416B (zh) 多相机测量物体的方法和装置
Fujita et al. Development of a non‐intrusive and efficient flow monitoring technique: The space‐time image velocimetry (STIV)
CN103512579B (zh) 一种基于热红外摄像机和激光测距仪的地图构建方法
CN113359097B (zh) 一种毫米波雷达和相机联合标定的方法
CN113819974B (zh) 一种免水尺的河流水位视觉测量方法
CN105526906B (zh) 大角度动态高精度激光测角方法
CN102168954A (zh) 基于单目摄像机的深度、深度场及物体大小的测量方法
CN113223075A (zh) 一种基于双目摄像头式的船舶高度测量系统及方法
CN108775872A (zh) 基于自动变焦实时图像处理的桥梁挠度检测方法
CN109489566A (zh) 锂电池隔膜材料分切宽度检测方法、检测系统及装置
CN113804916B (zh) 一种基于最大流速先验信息的频域时空图像测速法
CN108710148A (zh) 三维倾角域稳相叠前深度偏移方法和装置
CN105491315A (zh) 一种投影仪伽马校正方法
Huang et al. Measurement method and recent progress of vision-based deflection measurement of bridges: A technical review
Zhen et al. Design and evaluation of an FFT-based space-time image velocimetry (STIV) for time-averaged velocity measurement
CN108507564B (zh) 一种基于点扩散函数拟合的星敏感器质心定位方法
CN102798380A (zh) 线阵图像中目标运动参数的测量方法
Zhang et al. Freight train gauge-exceeding detection based on three-dimensional stereo vision measurement
CN117092621A (zh) 基于光线追踪校正的高光谱图像-点云立体配准方法
CN113592934B (zh) 一种基于单目相机的目标深度与高度测量方法及装置
CN111398625B (zh) 一种物理模型试验中的测速方法
CN115100446A (zh) 一种sar与可见光遥感图像匹配的相似性度量方法
Cao et al. River surface velocity estimation using optical flow velocimetry improved with attention mechanism and position encoding
CN114565642A (zh) 基于无人机热红外影像序列的光流流速场测量方法及系统
CN114397476A (zh) 面向频域时空图像测速的流速有效性识别及修正方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant