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

1.本发明属于图像法测流技术领域，涉及一种时空图像测速法，尤其涉及一种基于最大流速先验信息的频域时空图像测速法。


背景技术：

2.时空图像测速方法(stiv)是一种以测速线为分析区域、通过检测合成时空图像的纹理主方向估计一维时均流速测量方法。它利用水流示踪物在三维时空域中运动的连续性，采用平行于顺流方向的测速线作为分析区域，在图像空间和序列时间组成的时空图像中检测和示踪物运动相关的纹理方向特征，直接估计指定空间方向上的时均运动矢量。由于具有空间分辨率高、实时性强的优点，在河流水面流速、流量的实时监测中具有特别的应用潜力。
3.频域时空图像测速法，通过将空域复杂的纹理主方向检测转换为频域中搜索图像频谱主方向的线性运算，能显著提高算法抗噪性，降低复杂度。然而在实际应用中，阴影、耀光、倒影、障碍物遮挡、风浪、雨雾等环境干扰因素容易导致背景噪声过大，时空图像幅度谱的信噪比偏低，使得时空图像纹理主方向的检测出现粗大误差，从而引起流速测量失效。
4.所以，需要一个新的技术方案来解决这个问题。


技术实现要素：

5.发明目的：为了克服现有技术中存在的不足，提供一种基于最大流速先验信息的频域时空图像测速法，能够根据每条测速线的实际流速范围确定最佳搜索区间，进而抑制不相关噪声对有用信号检测的干扰。
6.技术方案：为实现上述目的，本发明提供一种基于最大流速先验信息的频域时空图像测速方法，包括如下步骤：
7.s1：假设水深与表面流速成正比例关系，根据已知的断面地形和水位信息计算断面最大水深和测速线水深；
8.s2：基于断面最大水深和测速线水深，利用当前水位下断面最大流速先验值计算测速线的最大流速先验值；
9.s3：基于变高水面摄影测量模型计算测速线的尺度变换因子，并根据频域时空图像测速法的逆变换确定时空图像频谱主方向的搜索区间；
10.s4：在搜索区间内搜索幅度谱极坐标投影的极值，将对应角度作为测速线的频域主方向估计值，并根据频域时空图像测速法的正变换计算测速线的流速值。
11.进一步地，所述步骤s1中测速线水深的计算方式为：
12.已知断面任意两点地形(x0，h0)、(x1，h1)，根据线性地形插函数得到两点间单条测速线上的断面地形h
l
：
[0013][0014]
以地面为基准点的水位值b和单点地形插值h
l
计算得到测速线水深d
l
：
[0015]
d
l
＝b
‑
h
l
ꢀꢀꢀꢀ
(2) 进一步地，所述步骤s2中测速线的最大流速先验值的计算方式为：
[0016][0017]
其中，v
max
为测速线的最大流速先验值；d
smax
为断面最大水深，断面最大水深为实测的已知值，通过实际输入的水位值和断面最低点相减得到；v
fmax
为当前水位下断面最大流速先验值，不同水位下的断面最大流速先验值根据经验或实测值确定。
[0018]
进一步地，所述步骤s3中测速线的尺度变换因子的计算公式为：
[0019][0020]
其中，j、n、f分别为图片尺寸、像素大小、相机焦距，可以通过相机的内参标定得到，图像传感器的像元尺寸s可由内参标定结果计算得到；相机俯仰角α通过相机的外参标定校准；摄像机水面高程h等于水位值b和摄像机距地面高度hc之和：
[0021]
h＝h
c
b
ꢀꢀꢀꢀ
(5)
[0022]
进一步地，所述步骤s3中基于变高水面摄影测量模型计算测速线的尺度变换因子的具体过程为：
[0023]
a1：在实验室对摄像机进行内参标定：
[0024]
在实验室中标定相机内参矩阵k和畸变参数矩阵d如下：
[0025][0026]
d＝[k
1 k
2 p
1 p2]
ꢀꢀꢀꢀ
(7)
[0027]
其中(c
x
，c
y
)表示畸变图像的像主点坐标，c
y
表示像主点的纵坐标，f
x
、f
y
分别表示摄像机在像平面x轴和y轴方向上的等效焦距，k1表示第一径向畸变系数，p1表示第一切向畸变系数，k2表示第二径向畸变系数，p2表示第二切向畸变系数，根据图像传感器的像元尺寸s得到相机焦距f：
[0028]
f＝(f
x
f
y
)
·
s/2
ꢀꢀꢀ
(8)
[0029]
畸变参数用于图像的非线性畸变校正：
[0030][0031]
其中，(x
′
，y
′
)和(x，y)分别为畸变和无畸变的相机坐标，它们与对应的图像坐标 (u
′
，v
′
)和(u，v)间满足：
[0032][0033][0034]
以上三式建立了无畸变图像坐标到畸变图像坐标的变换关系；
[0035]
a2：进行相机的外参标定：
[0036]
将摄像机架设到岸边，读取集成在摄像机上的三轴传感器数值得到俯仰角α＝，横
滚角ω，由于相机和激光测距仪的中心并不完全重合，畸变校正后需要继续进行偏心角的检校，偏心角定义为相机和测距仪的三个旋转角，方位角置零，检校俯仰角α以及横滚角ω：
[0037][0038]
其中，e
x
、e
y
表示所有网格在x和y方向尺寸的均方根误差：
[0039][0040]
m、l表示棋盘格内角点的数量，(x
i，j
，y
i，j
)表示其物方坐标，d
x
、d
y
表示网格尺寸的真值；
[0041]
a3：建立倾斜视角下的中心透视投影成像模型：
[0042]
像平面和物平面坐标系分别用(x，y)，(x，y)表示；0为透镜平面的光心，o、o
′
分别为其在像平面和物平面上的投影点；c为像平面延长线和通过光心的水平线的交点；h为光心到物平面间的垂直距离，c为对应的垂足点；相机的俯仰角α定义为相机主光轴和物平面间的夹角；s表示图像传感器的像元尺寸；f为相机焦距；(m，n)表示图像的大小；下标(i，j)表示像素的坐标；
[0043]
像平面坐标(x，y)与图像坐标(i，j)的关系可由以下公式表示：
[0044][0045]
对成像模型进行合理地简化：实际布设岸基式系统时，可以分别通过对齐测量断面及调节云台水平将方位角和横滚角的值置零，因此可以仅考虑相机主光轴平行于断面方向仅存在俯仰角的情况；
[0046]
物像尺度变换因子s可由x方向或y方向对应换算关系计算得出，这里以x方向为例计算；在x方向，假设像素p
i，j
位于像主点o的左侧，p
i，j
及其相邻像素p
i 1，j
在物平面上的投影点分别为p
i，j
和p
i 1，j
，射线p
i，j
o和p
i 1，j
o与投影线p
j
o的夹角分别为和φ，对于p
i，j
和p
i 1，j
，分别有：
[0047][0048][0049]
在y方向，假设像素p
i，j
位于像主点o的下方，p
i，j
及其相邻像素p
i，j 1
在物平面主纵线上的投影点分别为p
j
和p
j 1
，射线p
j
o和p
j 1
o与物平面p
j 1
c的夹角分别为β和γ，对于p
i，j
，其像平面主纵线上的投影点力满足以下三角关系：
[0050][0051]
由于α＝∠coo、β＝∠cop
j
，代入上式得：
[0052]
[0053]
将像素p
i，j
在x方向上的二维物像尺度因子s用其物点p
i，j
和相邻像素对应物点间的距离来描述，即：
[0054][0055]
结合以上各式，x方向的物像尺度变换因子s可由以下公式表示：
[0056][0057]
其中，j、n、s、f分别为图片尺寸，像素大小，像元尺寸，相机焦距，是上述相机内参标定求出的数值；相机俯仰角α由上述相机的外参标定校准；摄像机水面高程h等于水位值b和摄像机距地面高度h
c
之和：
[0058]
h＝h
c
b
ꢀꢀꢀꢀ
(21)
[0059]
进一步地，所述步骤s3中时空图像频谱主方向的搜索区间的确定方法为：
[0060]
所述搜索区间的取值和光流速度的方向有关，当光流速度为正时，频域搜索区间为 [90
°
θ
min
，θ
max
]；当光流速度为负时，频域搜索区间为[180
°‑
θ
max
，90
°‑
θ
min
]；
[0061]
其中，θ
max
为上边界，用于滤除时空图像中风浪、紊流等水平方向的背景噪声，采用以下公式求出：
[0062][0063]
θ
min
为下边界，用于滤除时空图像中阴影、耀光、遮挡等垂直方向的背景噪声，取值为定值。
[0064]
进一步地，所述步骤s4具体为：
[0065]
所述幅度谱极坐标投影的极值是指将傅里叶变换后的频谱信息映射到极坐标下，沿角度及半径方向投影遍历搜索区间，计算出投影的极值；
[0066]
以图像中心为原点o、n/2为半径设置搜索线，以δθ＝0.1
°
为步进、0～180
°
为区间建立极坐标系，统计搜索测速线上的幅度分布f(t，θ)，按下式对幅度谱的平方项进行累加得到能量
‑
角度分布：
[0067][0068]
在搜索区间[θ
min
，θ
max
]搜索p(θ)极值，将对应的频域主方向θ代入以下频域时空图像测速法的正变换公式，可得该条测速线上的流速值v：
[0069][0070]
完成单条测速线的表面流速测量；δt为时间间隔，计算公式为：
[0071][0072]
有益效果：本发明与现有技术相比，能够根据每条测速线的实际流速范围确定最佳搜索区间，进而抑制不相关噪声对有用信号检测的干扰，可以有效滤除野外复杂光照、水流及气象条件下倒影、耀光、紊流、雨雾、遮挡等环境噪声对时空图像测速法的干扰，从而抑制测量粗大误差，显著改善低流速及高流速测量的稳定性和可靠性。
附图说明
[0073]
图1是本发明方法的流程图；
[0074]
图2是测量系统布设示意图；
[0075]
图3是断面插值地形示意图；
[0076]
图4是变高平面摄影测量模型示意图，其中图4.1是x方向的立体视图，图4.2是剖面视图；
[0077]
图5是噪声的能量区间图；
[0078]
图6是设置搜索区间前后的搜索结果对比图，其中图6.1是受到干扰时的错误结果，图6.2是消除了粗大误差的正确结果；
[0079]
图7是设置搜索区间前后的测速结果对比图，其中图7.1是使用本发明方法之前的测速结果，图7.2是使用了本发明方法之后的测速结果。
具体实施方式
[0080]
下面结合附图和具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本技术所附权利要求所限定的范围。
[0081]
本发明提供一种基于最大流速先验信息的频域时空图像测速方法，如图1所示，其包括如下步骤：
[0082]
s1：假设水深与表面流速成正比例关系，根据已知的断面地形和水位信息计算断面最大水深和测速线水深：
[0083]
s2：基于断面最大水深和测速线水深，利用当前水位下断面最大流速先验值计算测速线的最大流速先验值：
[0084]
s3：基于变高水面摄影测量模型计算测速线的尺度变换因子，并根据频域时空图像测速法的逆变换确定时空图像频谱主方向的搜索区间：
[0085]
s4：在搜索区间内搜索幅度谱极坐标投影的极值，将对应角度作为测速线的频域主方向估计值，并根据频域时空图像测速法的正变换计算测速线的流速值。
[0086]
本实施例中以某水文站为例，将上述方法进行实例应用，具体如下：
[0087]
(1)参照图3，以海平面为基准，已知断面任意两点地形(x0，h0)、(x1，h1)，根据线性地形插函数得到两点间单条测速线上的断面地形h
l
：
[0088][0089]
以海平面为基准点的水位值b和单点地形插值h
l
计算得到测速线水深d
l
：
[0090]
d
l
＝b
‑
h
l
ꢀꢀꢀꢀ
(27)
[0091]
根据断面最大水深d
smax
和测速线水深d
l
，建立比例关系，将断面最大流速值v
fmax
转化为测速线的实际流速最大值v
max
：
[0092][0093]
本实施例中按以上思路，已知断面地形由水文站提供，以海平面为基准，每隔五米测量一次地形。取(85.0，983.10)，(90.0，983.00)，插值地形为983.02m。输入水位为 989.04m，断面最低点为982.60m，最大水深为6.44m，测速线水深为6.02m。此水位值对应水
位级下断面最大流速为3.5m/s，计算可得测速线最大流速为3.27m/s。
[0094]
(2)通过视频帧速率和采样帧间隔计算时间间隔δt：
[0095][0096]
本实施例中实际测量时采样帧间隔为1，帧速率为25fps，δt取0.04s。
[0097]
(3)在实验室对摄像机进行内参标定：
[0098]
相机的成像原理建立在多个坐标系的变换之上，即摄像机坐标系、图像物理坐坐标系、图像像素坐标系和世界坐标系。根据这四个坐标系之间的关系，现实场景下的任意一点都可变换至其中的任一坐标系中。根据这四个坐标系之间的关系，现实场景下的任意一点都可变换至其中的任一坐标系中，具体如下：
[0099]
step 1：世界坐标系中的一点通过旋转和平移变换转换至摄像机坐标系；
[0100]
step 2：通过投影变换可将step 1中得到的相机坐标点变换为图像物理坐标系中的像点p’；
[0101]
step 3：由于图像物理坐标系和图像像素坐标系同处一个平面，像点p’可通过缩放和平移变换最终转换为像素坐标系中的点p(μ，v)。
[0102]
以上坐标系的转换过程可通过一个矩阵来表示：
[0103][0104]
其中：s为像元尺寸，为旋转平移变换矩阵，为投影变换矩阵，f为焦距，为缩放平移变换矩阵。分别表示单位距离上横、纵坐标的像素个数。
[0105]
标定后得到像元尺寸s＝0.0013，内参矩阵k和畸变参数矩阵d如下：
[0106][0107]
d＝[k
1 k
2 p
1 p2]＝[
‑
0.4074 0.2169
ꢀ‑
0.0007
ꢀ‑
0.0004]
ꢀꢀꢀꢀꢀ
(32)
[0108]
其中，c
x
表示像主点的横坐标，c
y
表示像主点的纵坐标，f
x
表示摄像机在像平面x 轴上的等效焦距，f
y
表示摄像机在像平面y轴上的等效焦距，相机焦距f根据图像传感器的像元尺寸s计算得到：
[0109]
f＝(f
x
f
y
)
·
s/2＝3.7447mm
ꢀꢀꢀꢀ
(33)
[0110]
k1表示第一径向畸变系数，p1表示第一切向畸变系数，k2表示第二径向畸变系数， p2表示第二切向畸变系数，畸变参数用于图像的非线性畸变校正：
[0111]
[0112]
其中，(x
′
，y
′
)和(x，y)分别为畸变和无畸变的相机坐标，它们与对应的图像坐标(u
′
，v
′
) 和(u，v)间满足：
[0113][0114][0115]
以上三式建立了无畸变图像坐标到畸变图像坐标的变换关系。
[0116]
(4)将摄像机架设到岸边，如图2所示，读取集成在摄像机上的三轴传感器数值得到俯仰角α＝15.718140，横滚角ω＝1.026672。由于相机和激光测距仪的中心并不完全重合，畸变校正后需要继续进行偏心角的检校。偏心角定义为相机和测距仪的三个旋转角，方位角置零，检校俯仰角α以及横滚角ω。
[0117][0118]
其中，e
x
、e
y
表示所有网格在x和y方向尺寸的均方根误差。
[0119][0120]
其中，m、l表示棋盘格内角点的数量，(x
i，j
，y
i，j
)表示其物方坐标，d
x
、d
y
表示网格尺寸的真值。
[0121]
标定前可通过角点检测法检测角点，进而可对棋盘格每行每列相邻角点进行距离，得到未检校时角点距离均方根误差的最小值。
[0122]
粗检校时可先控制俯仰角和横滚角在正负1
°
内变化，重复上述检测步骤，求出均方根误差最小时对应的俯仰角α和横滚角ω。
[0123]
粗检校完成后进行进一步精确计算，以0.1
°
为基准，取粗检校所得角度相邻两点，通过三点高斯曲线拟合得到曲线最低点即为最终标定所得结果，最终计算得到俯仰角α＝19.818140和横滚角ω＝
‑
0.009788。
[0124]
(5)建立倾斜视角下的中心透视投影成像模型：像平面和物平面坐标系分别用(x， y)，(x，y)表示；o为透镜平面的光心，o、o
′
分别为其在像平面和物平面上的投影点；c 为像平面延长线和通过光心的水平线的交点；h为光心到物平面间的垂直距离，c为对应的垂足点；相机的俯仰角α定义为相机主光轴和物平面间的夹角；s表示图像传感器的像元尺寸；f为相机焦距；(m，n)表示图像的大小；下标(i，j)表示像素的坐标。
[0125]
像平面坐标(x，y)与图像坐标(i，j)的关系可由以下公式表示：
[0126][0127]
对成像模型进行合理地简化：实际布设岸基式系统时，可以分别通过对齐测量断面及调节云台水平将方位角和横滚角的值置零，因此可以仅考虑相机主光轴平行于断面方向仅存在俯仰角的情况。
[0128]
本实施例获取到如图4所示的变高平面摄影测量模型示意图，其中，图4.1是x方向
的立体视图，图4.2是剖面视图。
[0129]
物像尺度变换因子s可由x方向或y方向对应换算关系计算得出，本实施例以x方向为例计算。参照图4.1，在x方向，假设像素p
i，j
位于像主点o的左侧。p
i，j
及其相邻像素p
i 1，j
在物平面上的投影点分别为p
i，j
和p
i 1，j
，射线p
i，j
o和p
i 1，j
o与投影线 p
j
o的夹角分别为和φ。对于p
i，j
和p
i 1，j
，分别有：
[0130][0131][0132]
在y方向，假设像素p
i，j
位于像主点o的下方。p
i，j
及其相邻像素p
i，j 1
在物平面主纵线上的投影点分别为p
j
和p
j 1
，射线p
j
o和p
j 1
o与物平面p
j 1
c的夹角分别为β和γ。对于p
i，j
，其像平面主纵线上的投影点p
j
满足以下三角关系：
[0133][0134]
由于α＝∠coo、β＝∠∠cop
j
，代入上式得：
[0135][0136]
将像素p
i，j
在x方向上的二维物像尺度因子s用其物点p
i，j
和相邻像素对应物点间的距离来描述，即：
[0137][0138]
结合以上各式，x方向的物像尺度变换因子s可由以下公式表示：
[0139][0140]
其中，j、n、s、f分别为图片尺寸，像素大小，像元尺寸，相机焦距，是上述相机内参标定求出的数值。相机俯仰角α由上述相机的外参标定校准。摄像机水面高程h等于水位值b和摄像机距地面高度h
c
之和：
[0141]
h＝h
c
b
ꢀꢀꢀ
(46)
[0142]
实际测量中摄像机水面高程计算得到18.760022m，代入计算s＝0.13024。
[0143]
(6)根据频域时空图像测速法的逆变换确定时空图像频谱主方向的搜索区间，上边界可由测速线流速先验值和尺度变换因子确定：
[0144][0145]
计算得到θ
max
约为176
°
。
[0146]
对时空图像进行二维傅里叶变换，用对数幅值显示变换结果，将变换后的频谱中心从四角移到中心，一、三象限对调，二、四象限对调。此时直流分量集中到频谱中心。在幅度谱中对幅度谱的平方项进行累加，对幅度谱进行能量
‑
角度分布计算，图5的结果显示低频噪声主要分布在90
°±5°
的范围内，所以θ
min
的取值为5
°
。
[0147]
当光流速度为正，频谱主方向在90
°
～180
°
范围，此时搜索区间的下边界取值为 90
°
θ
min
，搜索区间为[90
°
θ
min
，θ
max
]；当光流速度为负，频谱主方向在0
°
～90
°
范围，此时搜
索区间的下边界θ
min
的值为90
°‑
θ
min
，搜索区间为[180
°‑
θ
max
，90
°‑
θ
min
]。
[0148]
此时搜索区间为[95
°
，176
°
]，为了直观的体现出效果，本实施例将设置搜索区间前后的结果进行对比，具体如图6所示，其中图6.1是受到干扰时的错误结果，图6.2是消除了粗大误差的正确结果，从图6的对比在中可明显看出上边界附近的误差被隔离在区间以外。
[0149]
(7)以图像中心为原点o、n/2为半径设置搜索线，以δθ＝0.1
°
为步进、0
°
～180
°
为区间建立极坐标系，统计搜索测速线上的幅度分布f(r，θ)，按下式对幅度谱的平方项进行累加得到能量
‑
角度分布，按下式对幅度谱的平方项进行累加得到能量
‑
角度分布：
[0150][0151]
在搜索区间[95
°
，176
°
]搜索p(θ)极值，将对应的频域主方向θ代入以下频域时空图像测速法的正变换公式，可得该条测速线上的流速值v：
[0152][0153]
为了验证本发明方法的效果，本实施例将现有频域时空图像测速法应用于相同水文站，测量到单条测速线的表面流速为2.396m/s，图7是设置搜索区间前后的测速结果对比图，其中，图7.1是使用本发明方法之前的测速结果，图7.2是使用了本发明方法之后的测速结果，经过对比可以明显看出本方法能够有效的消除远岸的粗大误差，并且准确识别出位于近岸的漩涡，使得测量结果更精确。