CN111798386B - River channel flow velocity measurement method based on edge identification and maximum sequence density estimation - Google Patents

River channel flow velocity measurement method based on edge identification and maximum sequence density estimation Download PDF

Info

Publication number
CN111798386B
CN111798386B CN202010588422.0A CN202010588422A CN111798386B CN 111798386 B CN111798386 B CN 111798386B CN 202010588422 A CN202010588422 A CN 202010588422A CN 111798386 B CN111798386 B CN 111798386B
Authority
CN
China
Prior art keywords
image
river
calibration
camera
pixel
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
CN202010588422.0A
Other languages
Chinese (zh)
Other versions
CN111798386A (en
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.)
Wuhan Dashuiyun Technology Co ltd
Original Assignee
Wuhan Dashuiyun Technology Co ltd
Wuhan University WHU
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 Wuhan Dashuiyun Technology Co ltd, Wuhan University WHU filed Critical Wuhan Dashuiyun Technology Co ltd
Priority to CN202010588422.0A priority Critical patent/CN111798386B/en
Publication of CN111798386A publication Critical patent/CN111798386A/en
Application granted granted Critical
Publication of CN111798386B publication Critical patent/CN111798386B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/40—Analysis of texture
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/40—Image enhancement or restoration using histogram techniques
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/13—Edge detection
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/80—Analysis of captured images to determine intrinsic or extrinsic camera parameters, i.e. camera calibration
    • G06T7/85—Stereo camera calibration
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/90—Determination of colour characteristics

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Operations Research (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Computing Systems (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种基于边缘识别与最大序列密度估计的河道流速测量方法，首先利用摄像头排水河道视频图像，根据设定标定点进行坐标转换标定，之后在视频中绘制测速线，通过测速线生成对应的时空图像，通过直方图均衡化来增强其纹理特征，之后通过Canny算子对均衡化之后的时空图像进行边缘识别，识别出闭合的多边形轮廓；然后求出时空图像的测速线的斜率；据斜率结合河面波纹纹理像素点在固定时间内的运动距离以及摄像机的帧率，可以计算出测速线的平均速度，对多根测速线进行加权平均即可得到河流的表面流速。本发明提供了一种无接触式河道流量测量方式，大大减少河道流量测量安全隐患，提高了测量精度。

Figure 202010588422

The invention discloses a river flow velocity measurement method based on edge identification and maximum sequence density estimation. First, a video image of a drainage river channel is used to perform coordinate conversion and calibration according to a set calibration point, and then a velocity measurement line is drawn in the video, which is generated by the velocity measurement line. Corresponding spatio-temporal image, its texture features are enhanced by histogram equalization, and then edge recognition is performed on the equalized spatio-temporal image by the Canny operator to identify the closed polygonal contour; then the slope of the velocity measurement line of the spatio-temporal image is obtained; According to the slope, combined with the moving distance of the ripple texture pixels on the river surface in a fixed time and the frame rate of the camera, the average speed of the speed line can be calculated, and the surface flow velocity of the river can be obtained by weighted average of multiple speed lines. The invention provides a non-contact river flow measurement method, which greatly reduces the hidden danger of river flow measurement and improves the measurement accuracy.

Figure 202010588422

Description

River channel flow velocity measurement method based on edge identification and maximum sequence density estimation
Technical Field
The invention belongs to the field of flow velocity measurement, relates to a river channel flow velocity measurement method, and particularly relates to a river channel flow velocity measurement method based on edge identification and maximum sequence density estimation.
Background
Hydrology tests are an important national foundation. The flow velocity measurement is taken as an important component of a hydrological test, and plays an important role in the aspects of hydraulic engineering planning, design, planning and management of water resources, flood prevention and drought control and the like. Along with the improvement of the informatization degree, the requirements on flow velocity measurement gradually get closer to automation and intellectualization. This requires that flow rate data be collected without the presence of a pilot for an extended period of time, video streaming providing such a means. Meanwhile, the flood always brings danger to the flow measurement work of hydrologic staff when the flood comes, and how to rapidly measure the flow velocity in a non-contact mode is very important.
For non-contact flow measurement, an ADCP flow velocity meter based on an acoustic principle is represented, which obtains a flow velocity by transmitting sound waves into water and performing doppler effect analysis on received scattered back signals. And an Image recognition-based method, wherein Adrian proposed a PIV algorithm (Particle Image Velocimetry) in 1991, and later Fujita made a further improvement on the basis of the PIV in the 90's of 20 th century, proposed an LSPIV (Large-Scale PIV) and successfully applied to the flow measurement of Yodo river.
At present, the method for measuring the flow velocity of the river channel is mainly contact flow measurement, and the non-contact method for measuring the flow velocity of the river channel is mainly a flow velocity meter. The principle of flow measurement is
Figure BDA0002554624130000011
Where N is the rotational speed of the rotor, T is the duration of the flow measurement, C is the instrument constant, and v is the flow rate. Often in the basic hydrological stationThe rotor current meter such as a rotary cup type and a rotary propeller type is adopted to test the flow velocity, a cableway system is needed to be matched, the influence of weather is easily caused, and the risk to the personal safety of testers is brought to a certain degree when flood comes. For a river with the width of about 300m and 10 speed measuring lines, the time of about 50 minutes is required for measuring the flow by only one point method.
In the current flow measurement method, the contact measurement has the problems of higher requirement on the matching environment of the flow measurement and long time consumption, and in the non-contact flow measurement, such as ADCP, the requirement on the environment of the flow measurement is relatively higher, and the safety of flow measurement personnel can be threatened when extreme events occur, so that the image flow measurement has important significance as a flow measurement means which has low environmental requirement and is free from the real-time monitoring of field personnel.
Disclosure of Invention
The invention mainly aims to provide a flow velocity measuring method based on edge identification and maximum sequence density estimation, and aims to solve the problems that the conventional river flow velocity measuring method is long in time consumption and cannot be normally used under the condition of mountain torrents disasters.
In order to achieve the purpose, the technical scheme adopted by the invention is as follows:
a river channel flow velocity measurement method based on edge identification and maximum sequence density estimation is characterized by comprising the following steps:
step 1, erecting a camera and setting a site calibration point: arranging cameras and calibration points on two sides of a river bank to be detected;
step 2, measuring coordinates of the calibration points on the spot and measuring the size of the cross section of the river channel: calibrating the calibration point erected in the step 1 by using a total station to obtain the world coordinate of the calibration point, measuring the attitude parameter of the camera by using an auxiliary measuring instrument to obtain the pixel coordinate of the calibration point, and measuring the size of the cross section of the river channel by using the total station;
step 3, calibrating the camera: obtaining the pixel coordinates and the space coordinates of the calibration points after the step 2, establishing the relationship between two-dimensional pixel points and three-dimensional coordinate points through DLT linear transformation, analyzing a DLT linear equation according to the world coordinates of the site calibration points and the corresponding world coordinates, and finding the optimal solution to obtain the internal and external parameters calibrated by the camera;
(xi p,yi p)＝f(xi 3D,yi 3D,zi 3D)
step 4, shooting a river video and generating a space-time image: shooting a river through an erected camera, drawing a speed measuring line in a video after the shooting is finished, and generating a corresponding space-time image through the speed measuring line;
step 5, space-time image preprocessing: the texture effect of the generated space-time image is not obvious, the texture characteristics of the space-time image are enhanced through histogram equalization, then the edge of the space-time image after equalization is identified through a Canny operator, a closed polygon outline is identified, and the identified sample is p ═ { p ═ p { (p } { (p) }iIn which p isiFor each polygon end point coordinate set
Figure BDA0002554624130000021
Step 6, slope of the local closed area: in order to obtain the slope of a speed measurement line of a spatio-temporal image, ridge regression analysis is performed on edge point rows constituting a closed polygon, and a linear slope sequence W of regression is obtained as { W }1,w2...wn}；
Step 7, the inclination angle of the space-time pattern texture: on the basis of the step 6, carrying out estimation of an optimal subinterval on the obtained sequence W, and taking a sequence average value in the optimal subinterval as an estimation value of the flow speed;
step 8, river surface flow rate: according to the slope estimation value calculated in the previous step, the average speed of the speed measurement lines can be calculated by combining the movement distance of the river surface ripple texture pixel points in a fixed time and the frame rate of the camera, and the surface flow velocity of the river can be obtained by carrying out weighted average on a plurality of speed measurement lines.
Further, in step 2, a camera is installed at a suitable position of the river channel, at least 6 ground calibration boards are prepared, and the calibration hanging boards are arranged on two sides of the river channel, so that the space coordinates of the calibration boards are measured by using a total station under the condition that the positions of the calibration boards are not changed during calibration.
Further, in step 3, in order to determine a corresponding relationship between coordinates of a pixel point and three-dimensional coordinates, and determine parameters inside and outside a camera model, a DLT linear transformation method is introduced, the coordinates of the pixel point in an image captured by the camera are recorded as (u, v), the coordinates of a corresponding point in a three-dimensional space are recorded as (x, y, z), and undetermined parameters calibrated by the camera are recorded as liThe transformation formula is as follows:
Figure BDA0002554624130000022
and solving the equation through a preset calibration point to complete the transformation from the three-dimensional coordinate point of the image sequence to the coordinate of the image pixel point.
Further, in the step 4, a specific method for generating the spatio-temporal image is as follows:
the method comprises the steps of drawing the change of brightness of pixel points on a flow velocity line along the time on an image through the flow velocity line drawn along the flowing direction in an obliquely shot video, wherein the horizontal axis represents the length of the flow velocity line, and the vertical axis represents the time axis, so that the generated image is called a space-time image.
Further, in the step 5, the specific steps are as follows:
step 5.1, the image collected by the camera is three-channel RGB, and for convenient processing, the image is converted into gray processing, and the calculation mode is as follows:
Figure BDA0002554624130000031
step 5.2, histogram equalization, wherein the pixel point with the gray level of i has niA probability of
Figure BDA0002554624130000032
Setting s and r as the gray levels of the image before and after transformation, performing normalization processing, and introducing a transformation functionThe number T is as follows
Figure BDA0002554624130000033
The equalized pixel values can be calculated from the histogram before transformation;
step 5.3, Canny edge identification, firstly, filtering the image, and processing by adopting Gaussian filtering; secondly, solving the gradient direction and the gradient amplitude by using a Sobel operator, wherein the specific method is to convolute the pixel matrix and the operator matrix of the picture:
Figure BDA0002554624130000034
where F is the pixel matrix of the image and the identified closed polygon is used as the basis for subsequent maximum sequence density estimation.
Further, the specific steps in the step 6 are as follows:
step 6.1, for each p identified in step 5iPerforming ridge regression analysis, i.e.
Figure BDA0002554624130000035
Wherein Y is { Y ═ Y1,y2..yn},X＝[x1，x2，...，xn]T w＝[w1,w2....,wn]T；
Step 6.2, the formula in the step 6.1 is iteratively solved by adopting gradient descent, namely
Figure BDA0002554624130000036
In the formula, eta is a step length parameter, and w obtained by solvingiThe tangent value of the angle representing the texture of the block can be considered, where
Figure BDA0002554624130000041
Further, the specific steps in step 7 are as follows:
w by point column regression analysis of each contour of step 6iSlope sequence of composition w ═ (w)1,w2,..,wn) Let ε be the allowable error, let the center value of the window be
Figure BDA0002554624130000042
Counting intervals by sliding windows
Figure BDA0002554624130000043
The number of internal samples, considering the interval covering the maximum number of samples N
Figure BDA0002554624130000044
Using the mean value in this interval
Figure BDA0002554624130000045
As an estimate of the tangent of the texture angle, wherein
Figure BDA0002554624130000046
Further, the specific steps in step 8 are as follows:
in the spatio-temporal image, the spatial-temporal image,
Figure BDA0002554624130000047
the maximum sequence density estimated value calculated in the step 7
Figure BDA0002554624130000048
The calculated flow velocity is that the distance of the moving characteristic quantity moving along the velocity measurement line in the time T is S under the assumption of three-dimensional space, and meanwhile, the pixel point moves by i pixels in the k frame, so that the flow velocity of the velocity measurement line is S
Figure BDA0002554624130000049
Wherein fps is a photographic cameraThe frame rate of (a) is higher than the frame rate of (b),
Figure BDA00025546241300000410
is the tangent of the slope of the spatio-temporal image, PxIs the spatial distance represented by the pixel; for the velocity v of different velocity linesi(i ═ 1,2.. n), use
Figure BDA00025546241300000411
To calculate the average flow velocity of the river surface, whereiAnd selecting the weight coefficient according to the corresponding specification.
The invention has the beneficial effects that:
the invention can carry out non-contact flow measurement on the river, creates a new measurement mode, almost achieves real-time measurement after establishing the database, effectively avoids potential safety hazards caused by contact measurement to operators, and is particularly suitable for occasions with torrential rivers.
Drawings
Fig. 1 is a schematic view of the arrangement of the camera and the index point on two sides of the river channel.
Fig. 2 is a schematic diagram of drawing a speed measurement line in a shooting video according to the present invention.
Fig. 3 is an STIV image generated in an embodiment of the present invention.
Fig. 4 is an STIV image after histogram equalization in an embodiment of the invention.
FIG. 5 shows the edge recognition result of the Canny operator according to the embodiment of the present invention.
FIG. 6 is a diagram illustrating maximum sequence density estimation according to an embodiment of the present invention.
FIG. 7 is an overall block diagram of an algorithm in an embodiment of the present invention.
Detailed Description
The invention is illustrated below with reference to the accompanying drawings, and as shown in fig. 1 to 7, a river channel flow velocity measurement method based on edge identification and maximum sequence density estimation includes the following steps:
step 1, erecting a camera and setting a site calibration point. Selecting proper positions on both sides of the river bank, erecting and fixing cameras, setting calibration points in order to obtain the relationship between the pixel coordinates and the three-dimensional coordinates during the primary measurement, calibrating by using calibration plates as calibration points on both sides of a river bank, uniformly arranging site calibration points on both sides of the river, measuring world coordinates and camera attitude parameters of the calibration points by using a total station and other auxiliary measuring instruments, particularly installing a camera at a proper position of the river channel for better acquiring calibration data, preparing at least 6 ground calibration plates, calibration boards are arranged on two sides of a river channel, as shown in fig. 1, the left bank of the river channel is provided with calibration boards 1,2, 3 and 4, the right bank is provided with calibration boards 5, 6 and 7, seven calibration boards are provided in total, and a camera 8 is fixed on the right bank, so that the total station is used for measuring the space coordinates of the calibration boards under the condition that the positions of the calibration boards are not changed during calibration.
And 2, measuring coordinates of the calibration points on the spot and measuring the size of the cross section of the river channel. And (3) calibrating the calibration board erected in the step (1) by using a total station, measuring the coordinate of the calibration board in a three-dimensional space, and measuring the attitude parameter of the camera by using an auxiliary measuring instrument to obtain the pixel coordinate of the calibration point. And simultaneously, measuring the size of the cross section of the river channel.
And step 3, calibrating the camera. In order to determine the corresponding relation between the coordinates of the pixel points and the three-dimensional coordinates, the parameters inside and outside the camera model are determined, and a DLT linear transformation method is introduced. And (3) obtaining the pixel coordinates and the space coordinates of the calibration point after the steps 1 and 2, and establishing the relationship between the two-dimensional pixel point and the three-dimensional coordinate point through DLT linear transformation. And resolving a DLT linear equation according to the world coordinates of the site calibration point and the corresponding world coordinates, and finding an optimal solution to obtain the internal and external parameters calibrated by the camera.
(xi p,yi p)＝f(xi 3D,yi 3D,zi 3D)
In the above formula, the first and second carbon atoms are,
Figure BDA0002554624130000051
is the coordinate of a two-dimensional pixel point,
Figure BDA0002554624130000052
are three-dimensional coordinates.
Specifically, if the coordinates of the pixel points in the image are recorded as (u, v), the coordinates of the corresponding points in the three-dimensional space are recorded as (x, y, z), and the undetermined parameters calibrated by the camera are recorded as li. According to transformation formulae
Figure BDA0002554624130000053
And solving the equation through a preset calibration point to complete the transformation from the three-dimensional coordinate point of the image sequence to the coordinate of the image pixel point.
And 4, shooting a river video and generating a space-time image. The river is shot through the erected camera, a speed measuring line (shown in figure 2) is drawn in a video after the shooting is finished, and a corresponding space-time image is generated through the speed measuring line.
A speed measurement line as shown in fig. 2 is drawn in a captured river video, and characteristics such as water surface ripples and ripples can be considered to move with water flow under the condition that disturbances such as wind are ignored. The resulting change in image intensity may reflect the surface flow velocity of the channel. The basic principle of the STIV flow measurement is that a space-time image is generated by plotting the change of the brightness of the pixel points on the flow velocity line with time on an image (as shown in fig. 3) through a flow velocity line (as shown in fig. 2) plotted along the flow direction in a video shot in an inclined manner, wherein the horizontal axis represents the length of the flow velocity line, and the vertical axis represents the time axis, so that the generated space-time image is called as the STIV image.
And 5, preprocessing the spatio-temporal image. The generated image texture effect is not obvious, and the texture features of the image are enhanced through histogram equalization. And then, carrying out edge identification on the equalized image through a Canny operator, identifying a closed polygon outline, and identifying the obtained sample as p ═ { p ═ piIn which p isiFor each polygon end point coordinate set
Figure BDA0002554624130000061
The specific method comprises the following steps:
step 5.1, the image collected by the camera is three-channel RGB, and for convenient processing, the image is converted into gray scale processing in a calculation mode
Figure BDA0002554624130000062
And 5.2, histogram equalization, wherein linear characteristics presented by the image texture are considered to be reflected on the degree change of the gray level. The histogram equalization method can be adopted to enhance local contrast characteristics under the condition of not influencing the overall contrast. The pixel point with gray scale i has niA probability of
Figure BDA0002554624130000063
n is the total number of pixel points. Let s and r be the gray levels of the images before and after transformation, and perform normalization processing, and introduce the transformation function T as follows
Figure BDA0002554624130000064
The equalized pixel values may be calculated from the histogram before the transformation. See figure 4 for an STIV image after histogram equalization
And 5.3, Canny edge identification, wherein in the step 4, in order to describe the angle characteristics of the texture, a Canny operator is introduced for edge identification. Firstly, filtering the image, wherein the Canny operator mainly works by depending on the brightness derivative of the image, and the unfiltered image has large noise and can cause large influence on the performance of the operator. The processing is usually performed by Gaussian filtering
Figure BDA0002554624130000065
Secondly, solving the gradient direction and the gradient amplitude by using a Sobel Operator (a convolution operation is carried out on the pixel matrix and the Operator matrix of the picture
Figure BDA0002554624130000071
Where F is the pixel matrix of the image, and the recognition result is shown in fig. 5, that is, the sample obtained by recognition is p ═ piIn which p isiFor each polygon end point coordinate set
Figure BDA0002554624130000072
The identified closed polygons serve as a basis for subsequent maximum sequence density estimation.
And 6, locally sealing the slope of the interval. In order to obtain the slope of the speed measurement line of the spatio-temporal image, ridge regression analysis is performed on the edge point rows constituting the closed polygon. Finding the linear slope sequence W ═ W of regression1,w2...wn}. The overall processing steps are shown in FIG. 7, and the speed measuring line preprocessing process comprises the following specific steps
Step 6.1 is that in the STIV image, the proportion of noise is smaller than the proportion of the correct result, and therefore the sample obtained by edge recognition is p ═ { p ═ pi}，
Figure BDA0002554624130000073
For each piPerforming ridge regression analysis, i.e.
Figure BDA0002554624130000074
Wherein Y is { Y ═ Y1,y2..yn},X＝[x1，x2，...，xn]T w＝[w1,w2....,wn]T。
Step 6.2, the formula of the step 6.1 is iteratively solved by gradient descent, namely
Figure BDA0002554624130000075
Wherein eta is a step length parameter. Solved wiThe tangent value of the angle representing this piece of texture can be considered.
And 7, dip angle of the space-time pattern texture. On the basis of step 6, the optimal subinterval of the obtained sequence W is estimated, and the average value of the sequence in the optimal subinterval is taken as the estimated value of the flow velocity, and the specific steps are as follows:
w by point column regression analysis of each contour of step 6iSlope sequence of composition w ═ (w)1,w2,..,wn). Let ε be the allowable error, center value of window be
Figure BDA0002554624130000076
Counting intervals by sliding windows
Figure BDA0002554624130000077
The number of inner samples. Considering the interval covering the maximum number of samples N
Figure BDA0002554624130000078
Using the mean value in this interval
Figure BDA0002554624130000079
As an estimate of the tangent of the texture angle, wherein
Figure BDA00025546241300000710
See fig. 6 for a detailed maximum sequence density estimate.
And 8, river surface flow velocity. According to the slope estimation value calculated in the previous step, the average speed of the speed measurement lines can be calculated by combining the movement distance of the river surface ripple texture pixel points in a fixed time and the frame rate of the camera, and the surface flow velocity of the river can be obtained by carrying out weighted average on a plurality of speed measurement lines; the method comprises the following specific steps:
in the case of an STIV image,
Figure BDA0002554624130000081
the maximum sequence density estimated value calculated in the step 7
Figure BDA0002554624130000082
The calculated flow velocity assumes that the moving characteristic quantity moves along the velocity measurement line within the time T by the distance S in the three-dimensional space, and meanwhile, the pixel pointIf i pixels move in the k frame, the velocity of the velocity measurement line is
Figure BDA0002554624130000083
Wherein fps is the frame rate of the shooting camera and the unit is frame/s,
Figure BDA0002554624130000084
tangent of the slope of the STIV spatio-temporal image, PxIs the spatial distance represented by the pixel.
For the velocity v of different velocity linesi(i ═ 1,2.. n), use
Figure BDA0002554624130000085
To calculate the average flow velocity of the river surface, whereiFor the weight coefficients, they are selected according to the relevant specifications GB 50179-1993 and the like.

Claims (7)

1. A river channel flow velocity measurement method based on edge identification and maximum sequence density estimation is characterized by comprising the following steps:
step 1, erecting a camera and setting a site calibration point: arranging cameras and calibration points on two sides of a river bank to be detected;
step 2, measuring coordinates of the calibration points on the spot and measuring the size of the cross section of the river channel: calibrating the calibration point erected in the step 1 by using a total station to obtain the world coordinate of the calibration point, measuring the attitude parameter of the camera by using an auxiliary measuring instrument to obtain the pixel coordinate of the calibration point, and measuring the size of the cross section of the river channel by using the total station;
step 3, calibrating the camera: obtaining the pixel coordinates and the space coordinates of the calibration points after the step 2, establishing the relationship between two-dimensional pixel points and three-dimensional coordinate points through DLT linear transformation, analyzing a DLT linear equation according to the world coordinates of the site calibration points and the corresponding world coordinates, and finding the optimal solution to obtain the internal and external parameters calibrated by the camera;
(xi p,yi p)＝f(xi 3D,yi 3D,zi 3D)
in the above formula, the first and second carbon atoms are,
Figure FDA0003462061430000011
is the coordinate of a two-dimensional pixel point,
Figure FDA0003462061430000012
is a three-dimensional coordinate;
step 4, shooting a river video and generating a space-time image: shooting a river through an erected camera, drawing a speed measuring line in a video after the shooting is finished, and generating a corresponding space-time image through the speed measuring line;
step 5, space-time image preprocessing: the texture effect of the generated spatio-temporal image is not obvious, the texture features of the spatio-temporal image are enhanced through histogram equalization, then the edge of the spatio-temporal image after equalization is identified through a Canny operator, a closed polygon outline is identified, and the identified sample is p ═ { p ═ piIn which p isiFor each polygon end point coordinate set
Figure FDA0003462061430000013
Step 6, slope of the local closed area: in order to obtain the slope of a speed measurement line of a spatio-temporal image, ridge regression analysis is performed on edge point rows constituting a closed polygon, and a linear slope sequence W of regression is obtained as { W }1,w2...wn}；
Step 7, the inclination angle of the space-time pattern texture: on the basis of the step 6, carrying out estimation of an optimal subinterval on the obtained sequence W, and taking a sequence average value in the optimal subinterval as an estimation value of the flow speed;
step 8, river surface flow rate: according to the slope estimation value calculated in the previous step, the average speed of the speed measurement lines can be calculated by combining the movement distance of the river surface ripple texture pixel points in a fixed time and the frame rate of the camera, and the surface flow velocity of the river can be obtained by carrying out weighted average on a plurality of speed measurement lines;
in the step 4, the specific method for generating the spatio-temporal image comprises the following steps:
the method comprises the steps of drawing the change of brightness of pixel points on a flow velocity line along the time on an image through the flow velocity line drawn along the flowing direction in an obliquely shot video, wherein the horizontal axis represents the length of the flow velocity line, and the vertical axis represents the time axis, so that the generated image is called a space-time image.
2. The method of measuring a flow rate of a river of claim 1, wherein: in the step 2, a camera is installed at a proper position of the river channel, at least 6 ground calibration boards are prepared as calibration points, and calibration hanging boards are arranged on two sides of the river channel, so that the space coordinates of the calibration boards are measured by using a total station under the condition that the positions of the calibration boards are not changed during calibration.
3. The method of measuring a flow rate of a river of claim 2, wherein: in the step 3, in order to determine the corresponding relationship between the coordinates of the pixel points and the three-dimensional coordinates and determine the parameters inside and outside the camera model, a DLT linear transformation method is introduced, the coordinates of the pixel points in the image shot by the camera are recorded as (u, v), the coordinates of the corresponding points in the three-dimensional space are recorded as (x, y, z), and the undetermined parameters calibrated by the camera are recorded as liThe transformation formula is as follows:
Figure FDA0003462061430000021
and solving the equation through a preset calibration point to complete the transformation from the three-dimensional coordinate point of the image sequence to the coordinate of the image pixel point.
4. The method of measuring a flow rate of a river of claim 2, wherein: in the step 5, the specific steps are as follows:
step 5.1, the image collected by the camera is three-channel RGB, and for convenient processing, the image is converted into gray processing, and the calculation mode is as follows:
Figure FDA0003462061430000022
step 5.2, histogram equalization, wherein the pixel point with the gray level of i has niA probability of
Figure FDA0003462061430000023
n is the total number of pixel points, s and r are the gray levels of the image before and after transformation, normalization processing is carried out, and a transformation function T is introduced as follows
Figure FDA0003462061430000024
The equalized pixel values can be calculated from the histogram before transformation;
step 5.3, Canny edge identification, firstly, filtering the image, and processing by adopting Gaussian filtering; secondly, solving the gradient direction and the gradient amplitude by using a Sobel operator, wherein the specific method is to convolute the pixel matrix and the operator matrix of the picture:
Figure FDA0003462061430000025
where F is the pixel matrix of the image and the identified closed polygon is used as the basis for subsequent maximum sequence density estimation.
5. The method of claim 4, wherein: the specific steps in the step 6 are as follows:
step 6.1, for each p identified in step 5iPerforming ridge regression analysis, i.e.
Figure FDA0003462061430000031
Wherein Y is { Y ═ Y1,y2..yn},X＝[x1，x2，...，xn]T w＝[w1,w2....,wn]T；
Step 6.2, the formula in the step 6.1 is iteratively solved by adopting gradient descent, namely
Figure FDA0003462061430000032
In the formula, eta is a step length parameter, and w obtained by solvingiThe tangent value of the angle representing the texture of the block can be considered, where
Figure FDA0003462061430000033
6. The method of claim 5, wherein: the specific steps in the step 7 are as follows:
w by point column regression analysis of each contour of step 6iSlope sequence of composition w ═ (w)1,w2,..,wn) Let ε be the allowable error, let the center value of the window be
Figure FDA0003462061430000034
Counting intervals by sliding windows
Figure FDA0003462061430000035
The number of internal samples, considering the interval covering the maximum number of samples N
Figure FDA0003462061430000036
Using the mean value in this interval
Figure FDA0003462061430000037
As an estimate of the tangent of the texture angle, wherein
Figure FDA0003462061430000038
7. The method of claim 6, wherein: the specific steps in the step 8 are as follows:
in the spatio-temporal image, the spatial-temporal image,
Figure FDA0003462061430000039
the maximum sequence density estimated value calculated in the step 7
Figure FDA00034620614300000310
The calculated flow velocity is that the distance of the moving characteristic quantity moving along the velocity measurement line in the time T is S under the assumption of three-dimensional space, and meanwhile, the pixel point moves by i pixels in the k frame, so that the flow velocity of the velocity measurement line is S
Figure FDA00034620614300000311
Where fps is the frame rate of the camera,
Figure FDA00034620614300000312
is the tangent of the slope of the spatio-temporal image, PxIs the spatial distance represented by the pixel; for the velocity v of different velocity linesi1,2.. n, using
Figure FDA00034620614300000313
To calculate the average flow velocity of the river surface, whereiFor the weight coefficients, they were chosen according to the corresponding specification GB 50179-1993.
CN202010588422.0A 2020-06-24 2020-06-24 River channel flow velocity measurement method based on edge identification and maximum sequence density estimation Active CN111798386B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010588422.0A CN111798386B (en) 2020-06-24 2020-06-24 River channel flow velocity measurement method based on edge identification and maximum sequence density estimation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010588422.0A CN111798386B (en) 2020-06-24 2020-06-24 River channel flow velocity measurement method based on edge identification and maximum sequence density estimation

Publications (2)

Publication Number Publication Date
CN111798386A CN111798386A (en) 2020-10-20
CN111798386B true CN111798386B (en) 2022-03-22

Family

ID=72803054

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010588422.0A Active CN111798386B (en) 2020-06-24 2020-06-24 River channel flow velocity measurement method based on edge identification and maximum sequence density estimation

Country Status (1)

Country Link
CN (1) CN111798386B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11885653B2 (en) 2021-09-24 2024-01-30 Hydro Radar, LLC Flow and level monitoring fluid system including nadir-facing and angle flow sensors with MIMO phase radar sensors

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112258570B (en) * 2020-11-17 2022-12-20 中国科学院空天信息创新研究院 Method and system for measuring and calculating water surface width of river full water area
CN112862898B (en) * 2021-02-05 2022-03-08 慧目(重庆)科技有限公司 Flow velocity measuring method based on computer vision
CN112861856B (en) * 2021-02-05 2022-05-27 慧目(重庆)科技有限公司 Drainage monitoring method based on computer vision and water body monitoring method
CN113012195B (en) * 2021-03-04 2024-02-06 西安电子科技大学 Method for river surface flow velocity estimation based on optical flow calculation and readable storage medium
CN113419075B (en) * 2021-05-31 2022-08-16 广州忘平信息科技有限公司 Ship speed measuring method, system, device and medium based on binocular vision
CN113343923B (en) * 2021-07-01 2024-12-03 江苏舆图信息科技有限公司 A real-time recognition method for drainage status of river outlets based on video images
CN114397476B (en) * 2021-11-15 2022-10-14 河海大学 Effective identification and correction method of flow velocity for spatiotemporal image velocimetry in frequency domain
CN115601664A (en) * 2022-01-28 2023-01-13 华北水利水电大学（Cn） AI unmanned aerial vehicle vision technology-based river surface flow velocity detection method and device
CN114754833A (en) * 2022-04-06 2022-07-15 武汉大水云科技有限公司 A portable video flow measurement device and video flow measurement method thereof
CN114481980B (en) * 2022-04-18 2022-06-24 水利部交通运输部国家能源局南京水利科学研究院 Data processing control method and system suitable for hydraulic engineering
CN115761563B (en) * 2022-09-27 2025-08-15 武汉大水云科技有限公司 River surface flow velocity calculation method and system based on optical flow measurement and calculation
CN115456916B (en) * 2022-11-11 2023-03-24 武汉大水云科技有限公司 Surface flow velocity calculation method based on image quality
CN115641363B (en) * 2022-12-23 2023-03-14 武汉大水云科技有限公司 River surface flow velocity calculation method based on depth matching
CN116466104B (en) * 2023-04-06 2025-09-16 昆明理工大学 Video current measurement method and device based on LK tracking and gray statistics feature method

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR950019642A (en) * 1993-12-24 1995-07-24 김은영 Method and device for simultaneous measurement of flow velocity and depth in rivers
US5952583A (en) * 1998-06-10 1999-09-14 Chang Min Tech Co., Ltd. Method of measuring a river horizontal average flow velocity
CN102564508A (en) * 2011-12-14 2012-07-11 河海大学 Method for implementing online tests of stream flow based on video images
CN103792533A (en) * 2014-01-21 2014-05-14 北京艾力泰尔信息技术有限公司 River cross section multipoint flow-measuring method based on fixed points
CN104777327A (en) * 2015-03-17 2015-07-15 河海大学 Time-space image speed measuring system and method based on auxiliary laser calibration
CN109977829A (en) * 2019-03-18 2019-07-05 江河瑞通（北京）技术有限公司 Current velocity testing method and device based on image recognition analysis technology
KR101996992B1 (en) * 2018-11-13 2019-07-08 주식회사 하이드로셈 Apparatus and Method for Measuring Flow Velocity of River using Optical Flow Image Processing
CN111062978A (en) * 2019-11-27 2020-04-24 武汉大学 Texture recognition method for spatio-temporal image flow measurement based on frequency domain filtering technology

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106097389B (en) * 2016-06-07 2018-08-21 浙江工业大学 River flow velocity measuring method based on image pattern recognition

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR950019642A (en) * 1993-12-24 1995-07-24 김은영 Method and device for simultaneous measurement of flow velocity and depth in rivers
US5952583A (en) * 1998-06-10 1999-09-14 Chang Min Tech Co., Ltd. Method of measuring a river horizontal average flow velocity
CN1281149A (en) * 1998-06-10 2001-01-24 昌民技术株式会社 Method of measuring overage flow rate of river water
CN102564508A (en) * 2011-12-14 2012-07-11 河海大学 Method for implementing online tests of stream flow based on video images
CN103792533A (en) * 2014-01-21 2014-05-14 北京艾力泰尔信息技术有限公司 River cross section multipoint flow-measuring method based on fixed points
CN104777327A (en) * 2015-03-17 2015-07-15 河海大学 Time-space image speed measuring system and method based on auxiliary laser calibration
KR101996992B1 (en) * 2018-11-13 2019-07-08 주식회사 하이드로셈 Apparatus and Method for Measuring Flow Velocity of River using Optical Flow Image Processing
CN109977829A (en) * 2019-03-18 2019-07-05 江河瑞通（北京）技术有限公司 Current velocity testing method and device based on image recognition analysis technology
CN111062978A (en) * 2019-11-27 2020-04-24 武汉大学 Texture recognition method for spatio-temporal image flow measurement based on frequency domain filtering technology

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Design and evaluation of an FFT-based space-time image velocimetry (STIV) for time-averaged velocity measurement;Zhang Z.等;《2019 14th IEEE International Conference on Electronic Measurement & Instruments (ICEMI)》;20200527;第503-514页 *
基于机器视觉的河水表面流速估计;冯全 等;《农业工程学报》;20181031;第34卷(第19期);第140-146页 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11885653B2 (en) 2021-09-24 2024-01-30 Hydro Radar, LLC Flow and level monitoring fluid system including nadir-facing and angle flow sensors with MIMO phase radar sensors

Also Published As

Publication number Publication date
CN111798386A (en) 2020-10-20

Similar Documents

Publication Publication Date Title
CN111798386B (en) River channel flow velocity measurement method based on edge identification and maximum sequence density estimation
KR101996992B1 (en) Apparatus and Method for Measuring Flow Velocity of River using Optical Flow Image Processing
KR101978351B1 (en) System and Method for Measuring Real Time Water Discharge based on CCTV Image
CN102564508B (en) Method for implementing online tests of stream flow based on video images
Dal Sasso et al. Exploring the optimal experimental setup for surface flow velocity measurements using PTV
CN104021676B (en) Vehicle location based on vehicle dynamic video features and vehicle speed measurement method
CN103996171B (en) Fluid motion vector estimation method based on space-time image
CN111914695B (en) Tidal bore monitoring method based on machine vision
CN116047112A (en) Method, device and storage medium for measuring fluid surface velocity
CN119470967B (en) Water area surface flow velocity area monitoring method based on intelligent visual analysis
CN116381650B (en) A Simulation and Testing Method for LiDAR Point Cloud Position and Intensity
CN109579774B (en) An Antenna Downtilt Measurement Method Based on Deep Instance Segmentation Network
CN112261719B (en) Area positioning method combining SLAM technology with deep learning
CN115326026B (en) Method and device for acquiring hydraulic characteristics based on non-contact measurement-hydrodynamic fusion assimilation
CN104408746A (en) Passenger flow statistical system based on depth information
CN109900331A (en) River discharge online test method based on depth camera signature tracking
JP4574657B2 (en) Method and apparatus for measuring water flow
CN117690303A (en) A noise early warning system, device and early warning method based on traffic data collection
CN113792452A (en) A method of video inversion of rainfall intensity based on raindrop velocity
Chodura et al. Evaluation of flight parameters in UAV-based 3D reconstruction for rooftop infrastructure assessment
Massó et al. Enhancing LSPIV accuracy in low-speed flows and heterogeneous seeding conditions using image gradient
Zhang et al. Computer vision-based real-time monitoring for swivel construction of bridges: from laboratory study to a pilot application
Kim et al. Measurement of debris flow velocity in flume using normal image by space-time image velocimetry incorporated with machine learning
CN118135465B (en) Nighttime raindrop falling speed measurement method and system based on monitoring near-infrared video
CN112819817A (en) River flow velocity estimation method based on graph calculation

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
TR01 Transfer of patent right

Effective date of registration: 20220712

Address after: 430070 15a, Zhongke Kaiwu building, No. 19, LUOJIASHAN Road, Hongshan District, Wuhan City, Hubei Province

Patentee after: Wuhan dashuiyun Technology Co.,Ltd.

Address before: 430072 Hubei Province, Wuhan city Wuchang District of Wuhan University Luojiashan

Patentee before: WUHAN University

Patentee before: Wuhan dashuiyun Technology Co.,Ltd.

TR01 Transfer of patent right