面向频域时空图像测速的流速有效性识别及修正方法

1.本发明涉及图像法测流技术，特别是一种面向频域时空图像测速的流速有效性识别及修正方法。


背景技术：

2.时空图像测速法(stiv)是一种应用在平直河道条件下，把测速线设定为分析区域、根据合成时空图像的纹理主方向估计一维时均流速的测速方法。它利用水流示踪物在三维时空域中运动的连续性，采用平行于顺流方向的测速线作为分析区域，在图像空间和序列时间组成的时空图像中检测和示踪物运动相关的纹理方向，从而计算测速线上的流速值。由于具有空间分辨率高、实时性强的优点，在河流水面流速、流量的实时监测中具有特别的应用潜力。
3.频域时空图像测速法，通过将空域复杂的纹理主方向检测转换为频域中搜索图像频谱主方向的线性运算，能显著提高算法抗噪性，降低复杂度。然而在实际应用中，由于耀光、阴雨、遮挡等自然环境干扰，容易导致背景噪声过大，时空图像幅度谱的信噪比偏低，可能导致在幅度谱上检测的频谱主方向有误，进而计算出错误的流速值，也不利于后续的流场重建。因此有必要在幅度谱中对检测的频谱主方向加入重定位措施，并采用自适应置信度机制来获取置信度阈值，进而进行有效性判定以及对错误流速值进行修正，从而消除野外监测条件下环境扰动造成的随机噪声及粗大误差，提高有效数据率。


技术实现要素：

4.发明目的：本发明的目的是提供一种面向频域时空图像测速的流速有效性识别及修正方法，从而消除野外监测条件下环境扰动造成的随机噪声及粗大误差，对后续的流场重建以及流量的计算提供可靠的保证。
5.技术方案：本发明所述的首先根据频域时空图像测速法在幅度谱中检测到的目标信号峰值及其对应方向，并由信号峰值的方向向幅度谱的最近边界进行搜索谷值以及对应方向，再取谷值点相对于信号峰值点的幅值对称点，通过对峰值对称性判断后重定位频谱主方向。其次将谷值作为背景噪声幅值计算置信度参数，并计算各条测速线的起点距以及相应的流速，接着根据自适应置信度机制获取置信度阈值来设置对应流速的有效性标志位。根据有效性标志位对错误流速识别并进行插值修正，最后统计流速场的有效数据率。一种面向频域时空图像测速的流速有效性识别及修正方法。包括：
6.(1)根据频域时空图像测速法得到幅度谱极坐标投影曲线
[0007][0008]
p
θ
表示在θ方向上的信号幅值，r表示频谱图上信号累加的半径，是固定值。 p(i,j,k,θ)表示频谱图上θ方向第k个且坐标为(i,j)的单个像素值。
[0009]
目标信号峰值ps及其对应主方向θ
max
由幅度谱极坐标投影曲线中对先验固定的有
效区间[θ
left
,θ
right
]搜索最大值获得。以当前遍历的幅值以及前后邻域内的10 个幅值取均值，在找到均值中最大的幅值作为检测到的目标信号峰值ps，对应的方向就是θ
max
。
[0010][0011]
谷值pn和方向θn按照目标信号峰值的检测方向向幅度谱最近的边界进行搜索获得。幅度谱最近的边界由目标信号峰值的方向在检测区间的位置确定，当检测的目标信号幅值靠近幅度谱极坐标投影曲线的右边界时；
[0012][0013]
当检测的目标信号幅值靠近左边界时，θ的取值在[θ
left
5,θ
max-5]。其次，根据对称关系，获得幅值对称点pn'。θn'是pn'对应的方向，通过幅度谱极坐标投影曲线遍历获得。
[0014]
pn'＝pnꢀꢀꢀ
(4)
[0015]
其次通过多组幅值点检测目标信号峰值ps邻域内的对称性，对称性的指标用 c表示，单组幅值点指标如下。
[0016][0017]
当c不超过对称性阈值c
thresh
时，判定为满足对称性，对主方向进行三点高斯拟合。(θn',pn')、(θ
max
,ps)以及(θn,pn)进行三点高斯拟合后的方向θm'。
[0018][0019]
得到θm'后，如果θm'在拟合区间内，则将拟合后的方向定为频谱主方向；否则，认为目标信号峰值ps对应的方向为频谱主方向。
[0020]
当c超过对称性阈值c
thresh
时，判定为不满足对称性，对主方向重新定位θ
m”。
[0021][0022]
θm如下。
[0023][0024]
(2)从将谷值pn作为背景均值噪声，计算信噪比。
[0025][0026]
测速线的起点距d
l
由物点到物平面光心垂足点的纵向距离和相机当前的起点距求得。流速值v
l
由频域时空图像测速法原理求得。
[0027]
(3)当所有测速线处理完毕后，考虑到自然环境波动造成的影响，通过
[0028]
自适应置信度机制来获取置信度阈值，thresh1和thresh2。具体表现为对所有测速线的置信度进行升序排列，根据排序后的置信度集合中按照测点的经验值划定一个有效流速数量比例阈值th，获取对应的置信度th1，并将其与一个经验性置信度阈值th2比较，获取自适应的置信度阈值。比较关系如下。
[0029][0030]
对有效性标志位f
l
进行赋值。f
l
表示第l条测速线的有效标志位。
[0031][0032]
当有效标志位为1的时候，代表流速值是正确的；当有效标志位为0的时候，代表流速值是不可靠的；当标志位为-1的时候，代表流速值是错误的。所有的有效性标志位都保存在集合s中，通过键值对key-value的方式进行存储，即测速线标号与对应的有效标志位进行匹配。
[0033]
s＝{l,f
l
}
ꢀꢀꢀ
(12)
[0034]
根据有效性标志位集合s定位到错误流速值对应的测速线标号le，以当前的错误流速矢量对应的测速线标号为基准，在集合s范围内向前以及向后寻找正确的流速矢量对应的测速线标号l
p
和ln，并根据已知的起点距并根据已知的起点距和正确流速值和正确流速值在d
l
和v
l
的坐标关系中采用线性插值修正。修正的结果如下。
[0035][0036]
(4)根据有效性标志位f
l
进行有效数据的统计，计算有效数据率pc。nc[0037]
表示根据有效性标志位统计的正确流速值的测速线数量，nc表示测速线的总数。
[0038][0039]
一种计算机存储介质，其上存储有计算机程序，该计算机程序被处理器执行时实现上述的一种面向频域时空图像测速的流速有效性识别及修正方法。
[0040]
一种计算机设备，包括储存器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述的一种面向频域时空图像测速的流速有效性识别及修正方法。
[0041]
有益效果：与现有技术相比，本发明具有如下优点：本发明增加了一种有效性识别和错误流速修正的方法，能够消除野外监测条件下环境扰动造成的随机噪声及粗大误差，对后续的流场重建以及流量的计算提供了可靠的保证。
附图说明
[0042]
图1是面向频域时空图像测速的流速有效性识别及修正方法的流程图。
[0043]
图2是频域时空图像测速法示意图。
[0044]
图3是倾斜视角下的中心透视投影测量模型示意图；其中图3a是x方向的示意图；图3b是剖面图。
[0045]
图4是有效性识别的结果图；其中图4a是正确情况的直方图；图4b是不可靠情况的直方图；图4c是错误情况的直方图。
[0046]
图5是进行插值修正前后的结果对比图；其中图5a是错误流速修正前的示意图；图5b是错误流速修正后的示意图。
具体实施方式
[0047]
下面结合附图对本发明的技术方案作进一步说明。
[0048]
如图1所示，本发明的一种面向频域时空图像测速的流速有效性识别及修正方法包括以下几个步骤：
[0049]
(1)首先在实验室对相机的内参进行标定，采用张正有法。通过不同姿态的棋盘格图像进行标定，计算内参矩阵k和畸变参数矩阵d：
[0050][0051]
d＝[k
1 k
2 p
1 p2]
ꢀꢀꢀ
(2)
[0052]
式中(c
x
,cy)为畸变图像的像主点坐标；f
x
和fy表示在像平面上的等效焦距； k1和k2表示径向畸变参数；p1和p2表示切向畸变参数；根据相机的像元尺寸s计算相机焦距f：
[0053][0054]
接着相机捕获30s视频，帧速率为25fps，获取750张原始图片，对每一幅图像进行畸变校正：
[0055][0056]
式中(x',y')和(x,y)分别表示有畸变和无畸变的坐标。
[0057]
然后通过测点提供的断面地形数据遍历获取水边线的范围，在畸变校正后的图像中按固定的步进设置测速线。
[0058]
(2)如图2所示，根据时空图像测速的流程，利用opencv提供的图片处理函数提取测速线上的像素，由空间序列和时间序列合成时空图像，进行图像补零，将时空图像扩充为正方形尺寸，进而对补零图像进行二维傅里叶变换，得到的频谱ftms如下：
[0059][0060]
式中(x,y)表示时空图像的坐标，(u,v)表示变换后的坐标，n表示补零后的尺寸。
[0061]
在频谱图中以图像中心为原点，按照0.1度为步进，统计每个方向θ上的半径为r的
范围内坐标为(i,j)的像素值的和p
θ
，由此获取幅度谱极坐标投影关系。
[0062][0063]
(3)在幅度谱极坐标投影曲线中检测目标信号峰值ps，目标信号峰值ps及其对应主方向θ
max
由幅度谱极坐标投影曲线中对先验固定的有效区间[θ
left
,θ
right
] 搜索最大值获得。以当前遍历的幅值以及前后邻域内的10个幅值取均值，在找到均值中最大的幅值作为检测到的目标信号峰值ps，对应的方向就是θ
max
。
[0064][0065]
其次计算谷值pn和方向θn，按照目标信号峰值的检测方向向幅度谱最近的边界进行搜索获得。幅度谱最近的边界由目标信号峰值的方向在检测区间的位置确定，当检测的目标信号幅值在幅度谱极坐标投影曲线的右侧时；
[0066][0067]
当检测的目标信号幅值靠近左边界时，θ的取值在[θ
left
5,θ
max-5]。其次，根据对称关系，获得幅值对称点pn'。θn'是pn'对应的方向，通过幅度谱极坐标投影曲线遍历获得。
[0068]
pn'＝pnꢀꢀꢀ
(9)
[0069]
其次通过多组幅值点检测目标信号峰值ps邻域内的对称性，对称性的指标用 c表示，单组幅值点指标如下。
[0070][0071]
当c不超过对称性阈值c
thresh
时，判定为满足对称性，对主方向进行三点高斯拟合。(θn',pn')、(θ
max
,ps)以及(θn,pn)进行三点高斯拟合后的方向θm'。
[0072][0073]
得到θm'后，如果θm'在拟合区间内，则将拟合后的方向定为频谱主方向；否则，认为目标信号峰值ps对应的方向为频谱主方向。
[0074]
当c超过对称性阈值c
thresh
时，判定为不满足对称性，对主方向重新定位θ
m”。
[0075][0076]
θm如下。
[0077][0078]
(4)根据目标信号幅值ps与谷值pn计算置信度。
[0079][0080]
其次计算测速线的起点距d
l
，由物点到物平面光心垂足点的纵向距离和相机当前的起点距关系获得。
[0081]
如图3a所示，根据摄影测量模型的原理建立倾斜视角下的中心透视投影测量模型，模型对应关系是物平面坐标(x,y)与像平面坐标(x,y)，o表示相机透镜平面的光心；o和o'分别表示光心在像平面与物平面的投影点；c为像平面延长线和通过光心的水平线的交点；h为光心到物平面间的垂直距离，c为对应的垂足点；相机的俯仰角α定义为相机主光轴和物平面间的夹角；s表示图像传感器的像元尺寸；f为相机焦距；m和n表示图像的尺寸；i和j表示像素的坐标。根据上述模型计算像平面到物平面对应的物象尺度因子δs。
[0082]
在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与投影线pjo的夹角分别为和φ。对于p
i,j
，满足以下关系：
[0083][0084]
由于φ＝∠p
i,j
opj，且：
[0085][0086]
进而得到：
[0087][0088]
同理，对p
i 1,j
有：
[0089][0090]
在y方向，假设像素p
i,j
位于像主点o的下方；p
i,j
及其相邻像素p
i,j 1
在物平面主纵线上的投影点分别为pj和p
j 1
，射线pjo和p
j 1
o与物平面p
j 1
c的夹角分别为β和γ。对于p
i,j
，其像平面主纵线上的投影点pj满足以下关系：
[0091][0092]
考虑到α＝∠coo、β＝∠copj得：
[0093][0094]
将p
i,j
在x和y方向上的二维物象尺度因子δs结合以上公式获取：
[0095][0096]
β表示相机当前俯仰角，可通过姿态传感器获取；h表示当前相机的水面高程，可由相机的高程hc与当前水位值a的差值得到。hc根据断面地形数据与相机到地面高度可得：
[0097]
h＝h
c-a
ꢀꢀꢀ
(22)
[0098]
再根据中心透视投影模型，如图3b所示，计算物平面上的物点到像平面光心点竖直方向上的纵向距离dy。
[0099][0100]
式中h表示相机的水面高程，同上；β和α分别表示相机的俯仰角和横滚角，可由姿态传感器获取；i和j表示图像坐标；m和n表示图像的尺寸；f表示相机焦距。
[0101]
当相机位于断面参考点的同岸时，测速线的起点距d
l
如下：
[0102]dl
＝camd dyꢀꢀꢀ
(24)
[0103]
当相机位于断面参考点的对岸时，测速线的起点距d
l
如下：
[0104]dl
＝camd-dyꢀꢀꢀ
(25)
[0105]
式中camd是相机的起点距，由相机架设位置确定。
[0106]
根据相机采集的帧速率与帧间隔的比值求得时间间隔δt；并结合步骤(3)中求得的θm和(5)中求得的物象尺度因子δs计算当前测速线上的流速值v
l
，表示第l条测速线对应的流速值。
[0107][0108]
将对应起点距d
l
以及上式计算出来的对应流速值v
l
并入集合u：
[0109]
u＝{l,d
l
,v
l
}
ꢀꢀꢀ
(27)
[0110]
判断测速线是否全部处理，如果未处理完则重复(2)～(4)的步骤，直到所有测速线都处理完毕。
[0111]
(5)当所有测速线处理完毕后，考虑到自然环境波动造成的影响，通过自适应置信度机制来获取置信度阈值，thresh1和thresh2。具体表现为对所有测速线的置信度进行升序排列，根据排序后的置信度集合中按照测点的经验值划定一个有效流速数量比例阈值th，取0.7，获取对应的置信度th1，置位1.1，并将其与一个经验性置信度阈值th2比较，这里取1.2获取自适应的置信度阈值。比较关系如下。
[0112][0113]
根据比较关系thresh2取值1.1。对有效性标志位f
l
进行赋值。f
l
表示第l条测速线的有效标志位。
[0114][0115]
当有效标志位为1的时候，代表流速值是正确的；当有效标志位为0的时候，代表流速值是不可靠的；当标志位为-1的时候，代表流速值是错误的。
[0116]
如图4a、4b、4c所示，分别代表根据以上识别方法识别出来的正确流速，不可靠流速以及错误流速的单例。
[0117]
所有的有效性标志位都保存在集合s中，通过键值对key-value的方式进行存储，即测速线标号与对应的有效标志位进行匹配。
[0118]
s＝{l,f
l
}
ꢀꢀꢀ
(30)
[0119]
(6)处理完所有的测速线后，根据集合u中的数据将所有的测速线以流速值为纵坐标，起点距为横坐标建立坐标系，如图4所示。关系如下：
[0120]
v(l)＝f(d
l
)
ꢀꢀꢀ
(31)
[0121]
其次，在有效性标志位集合s中遍历有效性标志位为-1的错误流速值对应的测速线标号le，以当前的错误流速矢量对应的测速线标号为基准，在集合s范围内向前以及向后寻找正确的流速矢量对应的测速线标号l
p
和ln，并根据已知的起点距并根据已知的起点距和正确流速值和正确流速值在d
l
和v
l
的坐标关系中采用线性插值修正。修正的结果如下。修正前后的效果如图5a、5b所示。
[0122][0123]
(7)根据有效性标志位f
l
进行有效数据的统计，计算有效数据率pc。nc[0124]
表示根据有效性标志位统计的正确流速值的测速线数量，本例中值为122。nc表示测速线的总数，本例中值为144，带入下式求得85％。
[0125]