A water area surface flow field reconstruction method and system based on time sequence trajectory fitting

Abstract

The application discloses a water area surface flow field reconstruction method and system based on time sequence trajectory fitting, and belongs to the technical field of non-contact visual velocity measurement. The method comprises the following steps: acquiring continuous images of a water area and pre-processing, and constructing an effective water surface mask; initializing a grid feature point in the effective water surface mask to obtain an initial feature set; constructing an image pyramid, performing multi-scale optical flow tracking, removing abnormal trajectories based on three-level reliability constraints, and dynamically supplementing points to ensure global feature coverage; performing polynomial fitting on effective trajectories, deriving a six-dimensional time sequence flow velocity vector, generating a flow velocity field through network aggregation, interpolation encryption and color gradient; for marine and lake water areas, high and low frequency wave velocity decoupling is completed, and the relative pressure field of the water surface is solved by combining the incompressible N-S equation and the finite element method. The application does not need artificial tracers, is suitable for weak texture and strong fluctuation water surface, has excellent velocity measurement precision and time sequence stability, and is suitable for surface flow field observation of rivers, channels, lakes and near seas.

Description

Technical Field

[0001] This invention relates to the field of non-contact visual velocimetry technology, and in particular to a method and system for reconstructing the surface flow field of water bodies based on time-series trajectory fitting. Background Technology

[0002] Existing visual velocimetry methods for water surfaces (such as LSPIV, STIV, OFV, etc.) generally employ inter-frame instantaneous difference / optical flow calculations to determine velocity, which suffers from the following core defects: reliance on instantaneous matching, making them prone to inaccuracies and feature loss in environments with weak textures, strong waves, and water surface reflections; lack of long-term temporal constraints, resulting in high velocity noise and poor temporal consistency; uneven spatial distribution of feature points, with local clustering and numerous blank areas; lack of a complete trajectory reliability control mechanism, leading to easy propagation of tracking errors; failure to introduce multi-scale coarse-to-fine propagation and dynamic point replenishment, resulting in poor stability across the entire area; ability to only invert flow velocity, unable to achieve closed-loop analysis of higher-order physical quantities such as flow velocity → pressure; lack of a standardized and reproducible process, making engineering reproduction and field deployment difficult; failure to publicly disclose temporal velocity vectors and explicit expressions of flow velocity vector matrices; lack of velocity interpolation encryption and color gradient visualization processes; and failure to decouple high and low frequencies of advection / waves + net drift in ocean wave environments, resulting in physically imprecise pressure calculations. Summary of the Invention

[0003] Based on this, the purpose of this invention is to provide a method and system for reconstructing the surface flow field of water bodies based on temporal trajectory fitting. This innovative technical approach replaces frame difference with temporal trajectory fitting, forming four major modules and a multi-layer constraint system: original and innovative input preprocessing and gridded feature point initialization, multi-scale optical flow tracking, trajectory reliability constraints, and trajectory fitting and velocity solution.

[0004] To achieve the above objectives, the present invention adopts the following technical solution: This invention provides a method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting. It comprises four core stages: image preprocessing and initialization of gridded feature points, multi-scale optical flow tracing, trajectory reliability constraints, and trajectory fitting and velocity calculation. It strictly distinguishes the physical levels of image detection, kinematic analysis, scene decoupling, and dynamic analysis, and explicitly states that high- and low-frequency wave velocity decoupling is not enabled for natural rivers and artificial channels, but must be enabled for lakes and oceans. The specific steps are as follows: S1. Obtain and sort the continuous image sequence of the water area to be tested, preprocess each frame of the image, and then construct an effective water surface mask to obtain the preprocessed image sequence; S2. Divide the flow measurement area effectively covered by the water surface mask into... A uniform grid is used to construct a locally effective mask for each grid. Feature points are detected in the locally effective mask. A globally unique index is assigned to each feature point. Constraint rules for single-grid feature points are set to constrain the number of features in a single grid to eliminate blind spots and obtain an initial set of spatially uniformly distributed feature points. S3. Construct a maximum of [number] layers from the preprocessed image sequence. The image pyramid is used; the global active point set (Active) is initialized using an initial feature point set. Initialization involves clearing the Active container, traversing all feature points in the initial feature point set, storing all feature points and their indexes in Active, and marking the initial state as updated. The Active set manages all feature points uniformly throughout the process. Forward / backward Lucas-Kanade optical flow tracking is performed from coarse to fine. The global active point set (Active) is maintained frame by frame, while the global indexes of the feature points remain unchanged. For continuous... The low-density grid of the frame is dynamically filled with points; after the single-level processing is completed, the trajectory is propagated to a finer scale as a motion prior until the finest scale is traversed, and finally the full-resolution trajectory set is output. S4. Perform three-level joint verification on the full-resolution trajectory set, and combine the trajectory length threshold to complete the abnormal trajectory removal and abnormal segment correction, and output the effective trajectory set with feature number index; S5. Convert the frame number of the effective trajectory into physical time, perform polynomial fitting on the feature points and differentiate to obtain the pixel velocity, combine it with the scale coefficient of the field calibration to convert it into physical velocity, calculate the velocity amplitude and direction angle, and construct a six-dimensional temporal velocity vector; aggregate the velocity data according to the three-level rules of trajectory level, frame level, and grid level to generate a gridded temporal velocity vector matrix; use spatial interpolation algorithm to complete the discrete velocity densification, combine color gradient mapping and vector arrow superposition to generate standard velocity field cloud map and vector map, complete the kinematic analysis of the whole process, and obtain the total velocity field of temporal observation; for lakes and ocean waters, execute steps S6 and S7 in sequence; for natural rivers and artificial channels, execute step S7. S6. The total velocity observed in the time series is decomposed into a composite component of low-frequency plane advection velocity, wave oscillation velocity and Stokes drift velocity. Time-frequency filtering or dynamic mode decomposition is used to separate the high-frequency and low-frequency plane advection velocities, and only the low-frequency plane advection velocity is retained as the decoupled velocity field. S7. Based on the time-series observed total velocity field / decoupled flow velocity field, divergence correction is performed on the velocity field; the Poisson pressure equation is derived based on the incompressible Navier-Stokes equations, and the computational domain is discretized using the finite element method; three types of boundary conditions are applied sequentially: free water surface, solid wall, and inflow / outflow open boundary, and the linear equations are solved to obtain the relative pressure field of the entire water surface; the relative pressure field of the entire water surface is interpolated, densified, and color-gradientd, and a standard pressure cloud map is output to complete the analysis of the entire dynamic process.

[0005] As a further improvement to the above-described solution of the present invention, step S1 specifically includes: S11. Acquire a continuous image sequence of the water area to be measured. And sorted by filename completion time. Indicates the first t Original frame image, 0≤ t ≤ T -1, T Represents the total number of frames in the image sequence; S12. For each frame image Perform ROI extraction → perspective correction → grayscale conversion → CLAHE contrast enhancement → Gaussian smoothing; for scenes with extremely strong reflections and missing textures on the water surface, a high-boost filter can be selectively added for edge enhancement processing; S13. Construct an effective water surface mask: In the formula, A region of interest mask for the water surface, used to delineate the target flow measurement range; A manually labeled non-water surface exclusion mask used to remove invalid areas; symbol For set intersection operation, the symbol is... This is the negation operation; S14. Finally, the preprocessed image sequence is obtained. , Representing the t The image after frame preprocessing.

[0006] As a further improvement to the above-described solution of the present invention, step S2 specifically includes: S21. Divide the current measurement area covered by the effective water surface mask into... Uniform grid, a single grid is denoted as , where: row index Column index For each grid Constructing a locally effective mask ; S22. In the local effective mask Internal use Shi Tomagi Corner detection algorithms detect feature points. Shi Tomagi Corner detection algorithm configuration parameters include quality coefficient Minimum spacing between feature points Boundary shielding parameters , The value ranges from 0.01 to 0.1. Values ​​range from 5 to 20 pixels; each detected feature point is assigned a globally unique and unchanging index number. Set single-grid feature point constraint rules: The initial number of feature points within a single grid must satisfy: , The actual number of feature points in the grid at row r and column c. The maximum number of feature points per grid cell is set; a minimum of two feature points per grid cell is enforced to eliminate blind spots in observation. S23. Remove coordinates located on the effective water surface mask. All feature points other than those in the original set are combined to form an initial dynamic feature point set that is spatially uniformly distributed. :

[0007] In the formula: Representing the The pixel coordinates of each feature point This represents the total number of initial feature points across the entire region.

[0008] As a further improvement to the above-described solution of the present invention, step S3 specifically includes: S31. The preprocessed image sequence The maximum number of layers to build is The image pyramid, and the image pyramid downsampling rule are:

[0009] In the formula, Representing the The first layer of the pyramid Frame image; Representing the The first layer of the pyramid Frame image; hierarchical index satisfy ,in For the finest scale, This is the coarsest scale; The downsampling function scales down the width and height of the previous image proportionally to obtain the next image. S32. From the coarsest scale of the pyramids We begin by using a coarse-to-fine strategy to traverse layer by layer towards finer scales: (1) Forward LK optical flow tracking At the current level Internal forward LK optical flow tracing is performed, by the first Frame Image To the Frame Image The formula for estimating the motion of feature points is:

[0010] in, For the first Layer Frame number Coordinates of feature points The first digit obtained by forward tracking Layer Frame number Predict the coordinates of each feature point; (2) Reverse LK optical flow backtracking verification At the current level Internally, a reverse LK optical flow backtracking verification is performed on the forward tracking results, by the first... Frame Image To the Frame Image Backtracking, the formula is:

[0011] in, The first obtained by reverse tracing Layer t, frame t Backtracking coordinates of each feature point; calculate and error Only retain those with errors less than the forward-backward consistency threshold. The reliable trajectory is obtained at the current level. Frame number Coordinates of reliable feature points ; The value ranges from 1 to 3 pixels; (3) Cross-scale trajectory transfer The reliable feature point trajectory results filtered at the current level are upsampled and passed to the next level as motion priors, as shown in the formula:

[0012] in: For the first Layer Frame number Coordinates of a reliable feature point; For the first Layer Frame number The initial coordinates of each feature point are predicted; coefficient 2 is the inverse scaling factor for downsampling the image pyramid by 2 times, corresponding to the mapping relationship from coarse-scale coordinates to fine-scale coordinates; S33. Count the number of active feature points in each frame grid by grid, for continuous... If a frame is determined to be a low-density mesh, a dynamic point-filling operation is performed, as follows: (1) Active point statistics: After optical flow tracking is completed for each frame, statistics are performed grid by grid. Frame grid The number of active feature points within the range is denoted as ; (2) Low-density mesh determination criteria: Set the sparseness determination ratio coefficient If a single grid is continuous Frames satisfy inequalities If the grid is low-density, then feature points need to be added. (3) Point replacement objective constraint: Reapply the low-density grid Shi Tomagi Corner detection algorithms detect feature points and supplement them to... ,in Fill in the proportional coefficient for dynamically added points; (4) Constraint on the number of feature points added in a single operation: The maximum number of feature points added in a single operation is [number missing]. And satisfy To avoid local clustering of feature points; (5) New feature point management: Assign a brand-new globally unique number index to the newly detected feature points, and directly include them in the optical flow tracking queue of all subsequent frames to ensure uniform coverage of features across the entire domain; S34. Maintain the global active point set Active frame by frame, distinguishing between three states: inflow, outflow, and update; once a feature point is assigned an index, it is never modified, repeated, or reset throughout the process; details are as follows: (1) Initialization phase After the pyramid is constructed, initialize the global active point set. Load the initial feature point set Assign a unique global index to each feature point. ; (2) Coordinate with forward LK optical flow, reverse LK optical flow, and reliability screening process At each level When performing forward LK optical flow, reverse LK optical flow, and reliability screening: ① Tracking is effective and the error is less than the forward and backward consistency judgment threshold. Feature points: Update location information, mark status as "updated", retain number, and continue to stay. gather; ① Tracking failure, error exceeding the forward and backward consistency judgment threshold Feature point: The state is marked as "moved out", from Remove from the set, terminating the trajectory of that point; (3) Coordinating cross-scale trajectory transfer Coarse-scale reliable feature points When upsampling is passed to a finer scale: ① The existing feature point global numbering is used directly without renumbering; ② The set inherits the batch of feature points synchronously, and the state continues; (4) Dynamic point replenishment in conjunction with low-density grid For continuous When the low-density grid of the frame is used to supplement new feature points: ① The status of newly added feature points is marked as "inflow"; ② Assign a brand new, unique global ID and add it. Gather and participate in subsequent tracking; S35. After single-level processing, the trajectory is propagated to finer scales as motion priors until the finest scale is reached, finally outputting a full-resolution trajectory set. .

[0013] As a further improvement to the above-mentioned scheme of the present invention, in step S4, the three-level joint verification includes forward and backward consistency constraints, prediction guidance constraints, and trajectory continuity constraints, specifically: S41. Level 1: Forward and Backward Consistency Constraints Forward and backward errors The calculation formula is:

[0014] in, European distance; Representing the Layer Frame number The coordinates of reliable feature points are taken from the global active point set. ; For the first Layer Frame number The coordinates of each feature point are traced back by reverse optical flow. when If the condition is met, the feature point is determined to be a reliable point and is directly retained; otherwise, it proceeds to the second-level constraint verification. The threshold for determining forward and backward consistency has a value range of 1 to 3 pixels. S42. Level Two: Predictive Guiding Constraints Using cross-scale predicted location to correct trajectory points, calculate prediction guidance bias. :

[0015] in, Based on reliable feature point coordinates The first obtained by forward LK optical flow tracing Layer Frame number Predict the coordinates of each feature point; The first result obtained for cross-scale transfer Layer Frame number Predict the coordinates of each feature point; If the correction simultaneously satisfies and If the feature point is positive, it is marked as a recovery point and retained; otherwise, the feature point is directly discarded. The threshold for predicting guidance deviation is set to a value ranging from 2 to 5 pixels. S43. Level 3: Trajectory Continuity Constraints Perform polynomial fitting on the complete trajectory and calculate the root mean square error of the fit. :

[0016] In the formula: For the trajectory number k Frame 1 i Measured feature coordinates of the feature points; For the trajectory number k Frame 1 i Fitted coordinates of feature points This represents the total number of valid frames for the trajectory. like If the trajectory is continuous and reliable, it is considered a continuous and reliable trajectory; otherwise, it is considered an abnormal trajectory; fitting deviation threshold. Value range: 0.5~2 pixels; S44. Set minimum trajectory length ,Require Frame; effective trajectory length less than If a trajectory is deemed abnormal, the entire segment is removed. The effective length of the trajectory is statistically analyzed based on reliable feature points. The generated continuous trajectory segments; S45. For abnormal trajectories with excessive fitting deviation or insufficient length, local re-examination should be performed on the abnormal segments first, specifically as follows: ① Locating the abnormal frame range: based on , , Lock the range of abnormal frames within the trajectory; ② Local Restart Multi-Scale Optical Flow Tracking: For the image region corresponding to the abnormal frame, the original algorithm process is reused, namely, reconstructing the local image pyramid and relying on the global active point set. Re-execute forward + reverse LK optical flow tracing with the original global number of the feature; ③ Re-examination and screening: Recalculate the error; if the new location meets the requirements... , , If all thresholds are used, then new points will be used. The original abnormal point was replaced, and the abnormal segment repair was completed. ④ Determine the flow split: If the constraints are still not met after re-checking, proceed to the timing interpolation stage; S46. For abnormal segments that still do not meet all constraints after local re-inspection, time-series interpolation is used to supplement the points, specifically: ① Select the adjacent valid and reliable points before and after the abnormal segment as the interpolation reference; ② Use a temporal interpolation algorithm to calculate the supplementary coordinates of the abnormal frame; ③ Perform a three-level joint verification on the points generated by interpolation. If the verification passes, they are included in the trajectory; if the verification fails, the abnormal segment is determined to be a permanent anomaly. S47. Trajectories that still fail to satisfy all constraints after re-inspection and interpolation are eliminated entirely, including any of the following cases: ① The number of abnormal segments is large and the span is large. Even after multiple re-detection and interpolation, all threshold constraints are still not met. ②The total effective length of the corrected trajectory is still less than the minimum number of frames. ; ③ The fitting error after interpolation continues to exceed the limit. Trajectory continuity cannot be guaranteed; S48. Final output: Set of valid trajectories with feature numbers All points within this trajectory set originate from reliable feature points at each level. .

[0017] As a further improvement to the above-described solution of the present invention, step S5 specifically includes: S51. Physical Time Conversion: Extract the frame index number of the valid trajectory. Read the original frame rate of the video , frame index sequence number Converted to physical time:

[0018] in, Representing the k The physical time corresponding to the frame; S52. Polynomial Fitting Rule: For feature point pixel coordinates , Perform third-order polynomial fitting respectively, and the fitting formula is as follows:

[0019]

[0020] In the formula: , for Directional polynomial coefficients, for Directional polynomial coefficients, Physical time; S53. Physical velocity conversion: The pixel velocity is obtained by taking the first derivative of the fitted curve, and then converted into the actual physical velocity by combining it with the scale coefficient s determined on-site.

[0021]

[0022] in, The unit is ; represent The physical velocity component in the X direction at time t, represent The physical velocity component in the Y direction at any given time, in m / s; S54. Definition of a six-dimensional temporal velocity vector: Constructing a complete six-dimensional temporal velocity vector containing spatiotemporal coordinates, velocity, velocity amplitude, and velocity direction angle:

[0023] Where: Resultant velocity The flow velocity direction angle is calculated using the four-quadrant arctangent function. ; u i , v i These represent the X-axis physical velocity component and the Y-axis physical velocity component, respectively. U i Represents the magnitude of the resultant velocity; S55. Three-level data aggregation system: ① Trajectory-level output: Collects all six-dimensional temporal velocity vectors of single feature points to form a trajectory-level dataset. ; ② Frame-by-frame output: Statistically calculate the number of valid feature samples, average velocity, median velocity, and mean velocity component for each frame; ③ Grid-level output: Press Grid-aggregated data to construct a gridded time-series velocity vector matrix The point velocity at the trajectory level is spatially grouped by grid to form a temporal velocity sequence for each grid. Representing the Line 1 Column grid in The aggregate value of the velocity vector at time t, Represented as The matrix; generating a grid velocity dataset. ;in: Representing the r Line 1 The horizontal average velocity component of the column grid, Representing the Line 1 The vertical average velocity component of the column grid. Representing the Line 1 The sum velocity of the grid; S56. Flow Field Visualization Process: Adopted Inverse distance weighting The system uses either radial basis function or kriging interpolation to refine the grid of discrete velocity points; it maps velocity amplitudes to color gradient fields, overlays vector arrows, and generates smooth velocity cloud maps and vector maps; it also exports trajectory history files, statistical tables, visualization images, and videos.

[0024] As a further improvement to the above-mentioned solution of the present invention, step S6 specifically includes the following steps: S61. The total temporal velocity is decomposed into a low-frequency advection component and a high-frequency wave composite component using time-frequency filtering or dynamic mode decomposition:

[0025] In the formula: This refers to low-frequency planar advection velocity. This is the composite component of wave oscillation velocity and wave-induced Stokes drift velocity; this step only performs one decoupling of high and low frequencies, and does not further decompose the wave components. The time-frequency filtering method is as follows: [The text abruptly ends here, likely due to an incomplete sentence or a missing section.] Filtering is performed with a cutoff frequency set to 0.05Hz~1.2Hz to filter out wave oscillation components above the cutoff frequency and retain the plane advection velocity below the cutoff frequency, which is the low-frequency plane advection velocity. The dynamic mode decomposition method is as follows: The total velocity of time-series observations is... Construct a snapshot matrix, extract the dominant dynamic mode through singular value decomposition, and separate the low-frequency dominant mode and high-frequency oscillation mode according to frequency characteristics. The frequency range of the low-frequency mode is consistent with that of the filtering method. S63. Set the cross-scale velocity deviation threshold. The unit is m / s, and the value is selected and configured according to the scenario. It is used to verify the validity of the decoupling results. The low-frequency plane advection velocity obtained by high and low frequency decoupling is... It is necessary to perform a deviation check between the coarse-scale / fine-scale velocity obtained from the original multi-scale trajectory transfer: If cross-scale velocity deviation This indicates the low-frequency plane advection velocity. Consistent with the overall movement trend, decoupling is effective, where... Output the original flow rate for multi-scale optical flow; like This indicates that the decoupling algorithm failed due to noise / strong wave interference, so this step was skipped, and the total velocity was observed using time series. This serves as a low-frequency planar advection velocity to prevent erroneous data from being included in subsequent pressure field calculations. After decoupling, only the pure low-frequency plane advection velocity field is retained as the decoupled velocity field and input into the pressure field solution stage.

[0026] As a further improvement to the above-mentioned solution of the present invention, a cross-scale velocity deviation threshold is provided. The possible values ​​are as follows: The random error of conventional water surface video speed measurement is generally ≤0.05~0.10 m / s; Under strong wave and turbulent flow conditions, the overall error can reach 0.10~0.20 m / s; The ocean current field is inverted at multiple scales with an inherent bias of approximately 0.05–0.6 m / s.

[0027] As a further improvement to the above-described solution of the present invention, step S7 specifically includes: S71. Incompressible Fluid Pre-Constraint: Perform divergence correction on the input time-series observed total velocity field / decoupled flow velocity field to strictly satisfy the fluid incompressibility condition: ;in, Represents the gradient / divergence differential operator. Represents the velocity vector; S72. Derivation of Poisson's pressure equation: The Poisson pressure equation is derived from the incompressible Navier-Stokes equations.

[0028] In the formula: This represents the relative pressure field over the water surface, with local atmospheric pressure as the reference zero point. The density of the water is a fixed value. ; For the Laplace operator; S73. Finite Element Mesh Rules: The computational domain is discretized using the finite element method, and the finite element mesh is compared with the front-end velocity statistics. The grid remains completely consistent; S74. Three types of mandatory boundary conditions: ① Dynamic boundary of free water surface: ; ② Bank / Hydraulic Structure Solid Boundary: ; ③ Open boundaries for river / water inflow and outflow: Directly specify the relative pressure at the boundary. Or pressure normal derivative ; in: The normal vector is the unit vector outside the boundary, and the gravitational acceleration has a fixed value. ; S75. Equation Solving and Post-processing: Solving the linear equation system yields the global relative pressure distribution. Interpolation and color gradient mapping are performed on the pressure field to output a standard water surface relative pressure cloud map.

[0029] This invention also provides a water surface flow field reconstruction system based on time-series trajectory fitting, which is used to implement the water surface flow field reconstruction method based on time-series trajectory fitting as described above, including: The image preprocessing module is used to acquire and sort a continuous image sequence of the water area to be tested, preprocess each frame of the image, construct an effective water surface mask, and output the preprocessed image sequence. The gridded feature point initialization module is used to divide the current measurement area covered by the effective water surface mask into... A uniform grid is used to construct a locally effective mask for each grid. Feature points are extracted from the locally effective mask using a feature detection algorithm. A globally unique index is assigned to each feature point. Constraint rules are set for feature points within a single grid to constrain the number of features in a single grid in order to eliminate blind spots in observation, resulting in an initial set of spatially uniformly distributed feature points. A multi-scale optical flow tracking and dynamic point-filling module is used for preprocessing image sequences. The maximum number of layers to build is The image pyramid is used; the global active point set (Active) is initialized using an initial feature point set. Initialization involves clearing the Active container, traversing all feature points in the initial feature point set, storing all feature points and their indexes in Active, and marking the initial state as updated. The Active set manages all feature points uniformly throughout the process. Forward / backward Lucas-Kanade optical flow tracking is performed from coarse to fine. The global active point set (Active) is maintained frame by frame, while the global indexes of the feature points remain unchanged. For continuous... The low-density grid of the frame is dynamically filled with points; after the single-level processing is completed, the trajectory is propagated to a finer level as a motion prior until the finest scale is traversed, and finally the full-resolution trajectory set is output. The trajectory reliability constraint module performs three-level joint verification on the full-resolution trajectory set, and combines trajectory length thresholds to remove abnormal trajectories and correct abnormal segments, outputting a set of valid trajectories with feature number indexes.

[0030] The trajectory fitting and temporal velocity vector solution module converts the frame number of the effective trajectory into physical time, performs polynomial fitting on the feature points and calculates the derivative to obtain the pixel velocity, combines the scale coefficient of the field calibration to convert it into physical velocity, calculates the velocity amplitude and direction angle, and constructs a six-dimensional temporal velocity vector. It aggregates velocity data according to three levels of rules: trajectory level, frame level, and grid level, and generates a gridded temporal velocity vector matrix. It uses a spatial interpolation algorithm to complete the discrete velocity densification, and combines color gradient mapping and vector arrow superposition to generate standard velocity field cloud map and vector map, completes the kinematic analysis of the whole process, and obtains the total velocity field of temporal observation. The wave velocity decoupling module decomposes the total velocity from time-series observations into a composite component of low-frequency plane advection velocity, wave oscillation velocity, and Stokes drift velocity. It employs time-frequency filtering or dynamic mode decomposition to separate high- and low-frequency plane advection velocities, retaining only the low-frequency plane advection velocity as the decoupled velocity field.

[0031] The pressure field solution module is used to perform divergence correction on the velocity field based on the time-series observed total velocity field / decoupled flow velocity field; derive the Poisson pressure equation based on the incompressible Navier-Stokes equations; and discretize the computational domain using the finite element method; apply three types of boundary conditions in sequence—free water surface, solid wall, and inflow / outflow open boundary—and solve the linear equations to obtain the relative pressure field of the entire water surface; perform interpolation refinement and color gradient processing on the relative pressure field of the entire water surface, output a standard pressure contour map, and complete the analysis of the entire dynamic process.

[0032] Compared with the prior art, the present invention has the following beneficial effects: This invention employs a temporal trajectory fitting velocity measurement framework, abandoning traditional frame difference velocity measurement. It uses continuous multi-frame trajectories, polynomial fitting, and first-order differential derivation to obtain smooth velocity, fundamentally improving stability and accuracy. It features low noise, high stability, and extremely high robustness to weak textures. This invention utilizes four standardized original modules: input preprocessing and meshed feature point initialization, multi-scale optical flow tracking, trajectory reliability constraints, and trajectory fitting and velocity solution, forming a complete new framework distinct from all existing visual flow measurement methods.

[0033] This invention employs a multi-scale optical flow tracking mechanism that uses a coarse-to-fine propagation approach. Coarse-scale trajectory upsampling guides fine-scale tracking, significantly improving robustness in weakly textured environments and increasing the success rate of tracking complex water surfaces. The invention utilizes a gridded uniform feature structure combined with a dynamic point-filling mechanism to force uniformity in the feature space and automatically fills in points in long-term low-density areas, ensuring stable coverage across the entire watershed. This invention employs a three-level trajectory reliability constraint: forward and backward consistency, prediction guidance, and trajectory continuity, comprehensively guaranteeing trajectory reliability at the frame, scale, and trajectory levels.

[0034] This invention assigns a globally unique index number to each feature point, tracking the temporal appearance, movement, inflow, and outflow of feature points to construct a stable and traceable dynamic feature set. This invention provides a complete six-dimensional temporal velocity vector with spatiotemporal coordinates and a gridded vector matrix, realizing the mathematical expression of discrete point arrival fields. For ocean and lake waters, this invention decouples high- and low-frequency wave velocities, decomposing the velocity into a low-frequency advection component and a composite component of wave oscillation and wave-induced net drift, without secondary decoupling, ensuring the physical accuracy of the pressure solution. Finally, using Navier-Stokes equations as the physical governing equations and the finite element method as the numerical method, the Poisson equation is constructed and strict boundaries are applied to completely solve the relative pressure field of the water surface.

[0035] This invention requires no artificial tracer, is highly adaptable to weak textures and strong undulations in water surfaces, has high velocity measurement accuracy, good temporal stability, and is reproducible throughout the entire process. It is suitable for surface flow field observation in rivers, canals, lakes, and nearshore oceans. Attached Figure Description

[0036] Picture 1 A flowchart of a method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting, provided in an embodiment of the present invention; Picture 2 A process diagram illustrating a method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting, provided in an embodiment of the present invention; Picture 3 This is a visualization of the velocity field in an embodiment of the present invention; Picture 4 This is a visualization of the surface relative pressure field in an embodiment of the present invention. Detailed Implementation

[0037] To facilitate understanding of the present invention, a more comprehensive description will be given below with reference to specific embodiments. However, the present invention can be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided to provide a more thorough and complete understanding of the disclosure of the present invention.

[0038] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used herein in the description of the invention is for the purpose of describing particular embodiments only and is not intended to limit the invention.

[0039] Reference Picture 1-Picture 2This embodiment proposes a method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting. It includes four core stages: image preprocessing and initialization of gridded feature points, multi-scale optical flow tracing, trajectory reliability constraints, and trajectory fitting and velocity calculation. It strictly distinguishes the physical levels of image detection, kinematic analysis, scene decoupling, and dynamic analysis, and explicitly states that high- and low-frequency wave velocity decoupling is not enabled for natural rivers and artificial channels, but it is mandatory for lakes and oceans. The specific steps are as follows: S1. Image Preprocessing: Obtaining a continuous image sequence of the water area to be tested. The system performs sorting and preprocessing on each frame of the image to construct an effective water surface mask. The preprocessed image sequence is obtained. Specifically: S11. Preprocessing steps for image sequences Acquire continuous image sequences of the water body to be tested , Indicates the first t Original frame image, 0≤ t ≤ T -1, T This represents the total number of frames in the image sequence. All images are sorted chronologically according to their filenames to ensure temporal continuity.

[0040] S12. Standardized Image Preprocessing Workflow For each frame of the original image after sorting The processing follows a fixed sequence: ROI extraction → perspective correction → grayscale conversion → Claeshe contrast enhancement → Gaussian smoothing. For special scenes with extremely strong reflections and missing textures on the water surface, a high-lift filter can be selectively added for edge enhancement.

[0041] S13. Rules for constructing effective water surface masks The effective water surface mask calculation formula is defined as follows:

[0042] In the formula: A region of interest mask for the water surface, used to delineate the target flow measurement range; This is a manually marked non-water surface exclusion mask used to remove invalid areas such as shorelines, floating debris, and obstructions; symbol For set intersection operation, the symbol is... This is an inversion operation. All subsequent feature detection, optical flow tracking, and velocity calculations in this method are performed only on the effective water surface mask. Execute within the specified scope.

[0043] S14. Output of Preprocessing Results After all the above processing steps, a standardized preprocessed image sequence is obtained. ,in Representing the The image after frame preprocessing.

[0044] S2. Mesh feature point initialization: Divide the current measurement area covered by the effective water surface mask into... Uniform mesh, constructing a locally effective mask Feature points are extracted from the local effective mask using a feature detection algorithm, and each feature point is assigned a globally unique index. By setting constraint rules for feature points within a single grid, the number of features in a single grid is constrained to eliminate blind spots in observation, resulting in an initial set of feature points that are spatially uniformly distributed. Specifically: S21. Mesh Generation The effective water surface mask-covered flow measurement area is divided into: R OK 、C A uniform grid of columns, with a single grid denoted as . row index Column index For each grid Constructing a locally effective mask .

[0045] S22. Feature Detection Algorithm Local effective masking Internal use Shi Tomagi Corner detection algorithms detect feature points. Shi Tomagi Corner detection algorithm configuration parameters include quality coefficient Minimum spacing between feature points Boundary shielding parameters ; The value ranges from 0.01 to 0.1. Values ​​range from 5 to 20 pixels.

[0046] S23. Hard constraints on feature point numbering and quantity Assign a globally unique and unchanging index to each detected feature point. Define the constraint rules for feature points within a single grid: The initial number of feature points within a single grid must satisfy the following: , This represents the actual number of feature points in the grid at row r and column c. This is the maximum number of feature points per grid; it forces that the number of feature points per grid be no less than 2 to avoid areas without observation data and eliminate blind spots in observation.

[0047] S24. Feature point set output Elimination coordinates located at the effective water surface mask All feature points other than those in the original dataset are combined to form an initial dynamic feature point set that is spatially uniformly distributed.

[0048] In the formula For the first The pixel coordinates of each feature point This represents the total number of initial feature points across the entire region.

[0049] S3. Multi-scale optical flow tracking and dynamic point filling of low-density meshes: For preprocessed image sequences The maximum number of layers to build is The image pyramid, the image pyramid downsampling rule is... Using an initial feature point set Initialize the global active point set Active by clearing the Active container and iterating through the initial feature point set. All feature points are processed, and all feature points and their indexes are stored in Active, with the initial state marked as updated. Active sets manage all feature points uniformly throughout the process. Forward / backward Lucas-Kanade optical flow tracking is performed from coarse to fine, with coarse-scale trajectory tracking... pass Sampling is used as the initial value for fine scale; the global active point set (Active) is maintained frame by frame, while the global index of the feature points remains unchanged; for continuous... The low-density grid is dynamically padded with points; after single-level processing, the trajectory is propagated to finer scales as motion priors until the finest scale is reached, and finally, a full-resolution trajectory set is output. Specifically, this includes: S31. The preprocessed image sequence The maximum number of layers to build is The image pyramid, and the image pyramid downsampling rule are:

[0050] In the formula, Representing the The first layer of the pyramid Frame image; Representing the The first layer of the pyramid Frame image; hierarchical index satisfy ,in This is the finest scale (full resolution original image). This is the coarsest scale (the highest level of the pyramid). This is a downsampling function that proportionally reduces the width and height of the previous image to obtain the next image.

[0051] The constructed image pyramid provides a framework for the multi-scale motion trajectory of this method.

[0052] S32. From the coarsest scale of the pyramids We begin by using a coarse-to-fine strategy to traverse layer by layer towards finer scales: (1) Forward LK optical flow tracking At the current level Internal forward LK optical flow tracing is performed, by the first Frame Image To the Frame Image The formula for estimating the motion of feature points is:

[0053] in, For the first Layer Frame number Coordinates of feature points The first digit obtained by forward tracking Layer Frame number Predict the coordinates of each feature point; (2) Reverse LK optical flow backtracking verification At the current level Internally, a reverse LK optical flow backtracking verification is performed on the forward tracking results, by the first... Frame Image To the Frame Image Backtracking, the formula is:

[0054] in, The first obtained by reverse tracing Layer t, frame t Backtracking coordinates of each feature point; calculate and error Only retain those with errors less than the forward-backward consistency threshold. The reliable trajectory is obtained at the current level. Frame number Coordinates of reliable feature points ; The value ranges from 1 to 3 pixels.

[0055] (3) Cross-scale trajectory transfer (coarse to fine) The reliable feature point trajectory results filtered at the current level are upsampled and passed to the next level (a finer scale) as motion priors, as follows:

[0056] in: For the first Layer (coarse scale) Frame number Coordinates of a reliable feature point; For the first Layer (fine scale) The predicted initial coordinates of each feature point; coefficient 2 is the inverse scaling factor of the image pyramid downsampling by 2 times, corresponding to the mapping relationship from coarse-scale coordinates to fine-scale coordinates.

[0057] S33. Count the number of active feature points in each frame grid by grid, for continuous... Frames identified as low-density meshes undergo dynamic patching. The complete determination rules, parameter constraints, and execution flow are as follows: (1) Active point statistics: After optical flow tracking is completed for each frame, statistics are performed grid by grid. Frame grid The number of active feature points within the range is denoted as ; (2) Low-density mesh determination criteria: Set the sparseness determination ratio coefficient If a single grid is continuous Frames satisfy inequalities If the grid is low-density, then feature points need to be added. (3) Point replacement objective constraint: Reapply the low-density grid Shi Tomagi Corner detection algorithms detect feature points and supplement them to... ,in Fill in the proportional coefficient for dynamically added points; (4) Constraint on the number of feature points added in a single operation: The maximum number of feature points added in a single operation is [number missing]. And satisfy To avoid local clustering of feature points; (5) New feature point management: Assign a brand-new globally unique number index to the newly detected feature points and directly include them in the optical flow tracking queue of all subsequent frames to ensure uniform coverage of features across the entire domain.

[0058] S34. Maintain the global active point set Active frame by frame, distinguishing between three states: inflow, outflow, and update; once all feature points are assigned an index, they are never modified, repeated, or reset throughout the process; the entire rule is a global top-level constraint, as follows: (1) Initialization phase After the pyramid is constructed, initialize the global active point set. Load the initial feature point set Assign a unique global number to each feature point This number remains unchanged throughout the entire process, frame, and scale.

[0059] (2) Coordinate with forward LK optical flow, reverse LK optical flow, and reliability threshold screening process At each level When performing forward LK optical flow, reverse LK optical flow, and reliability screening: ① Tracking is effective and the error is less than the forward and backward consistency judgment threshold. Feature points: Update location information, mark status as "updated", retain number, and continue to stay. gather; ② Tracking failure, error exceeding the forward and backward consistency judgment threshold Feature point: The state is marked as "moved out", from Remove from the set and terminate the trajectory of that point.

[0060] (3) Coordinating cross-scale trajectory transfer Coarse-scale reliable feature points When upsampling is passed to a finer scale: ① The existing feature point global numbering is used directly without renumbering; ② The set inherits the batch of feature points synchronously, and the state continues.

[0061] (4) Dynamic point replenishment in conjunction with low-density grid For continuous When the low-density grid of the frame is used to supplement new feature points: ① The status of newly added feature points is marked as "inflow"; ② Assign a brand new, unique global ID and add it. The group will participate in subsequent tracking.

[0062] S35. After single-level processing, the trajectory is propagated to finer scales as motion priors until the finest scale is reached, finally outputting a full-resolution trajectory set. .

[0063] S4. Trajectory Reliability Constraints: For full-resolution trajectory sets The system performs a three-level joint verification of forward and backward consistency constraints, prediction guidance constraints, and trajectory continuity constraints. Simultaneously, it incorporates trajectory length thresholds to remove abnormal trajectories and correct abnormal segments, outputting a set of valid trajectories with feature-indexed numbers. Specifically, this includes: S41. Level 1: Forward and Backward Consistency Constraints Forward and backward errors The calculation formula is:

[0064] in, European distance; Representing the Layer (coarse scale) Frame number The coordinates of reliable feature points are taken from the global active point set. The global numbering of feature points remains unchanged; For the first Layer Frame number The coordinates of each feature point are traced back by reverse optical flow. when If the condition is met, the feature point is determined to be a reliable point and is directly retained; otherwise, it proceeds to the second-level constraint verification. The threshold for determining forward and backward consistency is 1 to 3 pixels.

[0065] S42. Level Two: Predictive Guiding Constraints Using cross-scale predicted location to correct trajectory points, calculate prediction guidance bias. :

[0066] in, For reliable feature points The first obtained by forward LK optical flow tracing Layer Frame number Predicted coordinates of each feature point; The first result obtained for cross-scale transfer Layer Frame number Predicted coordinates of each feature point; If the correction simultaneously satisfies and If the feature point is positive, it is marked as a recovery point and retained; otherwise, the feature point is directly discarded. The threshold for predicting guidance deviation is set to a value ranging from 2 to 5 pixels.

[0067] S43. Level 3: Trajectory Continuity Constraints Perform polynomial fitting on the complete trajectory and calculate the root mean square error of the fit. :

[0068] In the formula: For the trajectory number k Frame 1 iMeasured feature coordinates of feature points (tracing back to reliable feature points selected at each level) ); For the trajectory number k Frame 1 i Fitted coordinates of feature points This represents the total number of valid frames for the trajectory. like If the trajectory is continuous and reliable, it is considered a continuous and reliable trajectory; otherwise, it is considered an abnormal trajectory; fitting deviation threshold. Value range: 0.5~2 pixels; S44. Set minimum trajectory length ,Require Frame; effective trajectory length less than If a trajectory is deemed abnormal, the entire segment is removed. The effective length of the trajectory is statistically analyzed based on reliable feature points. The generated continuous trajectory segments; S45. For abnormal trajectories with excessive fitting deviation or insufficient length, local re-examination should be performed on the abnormal segments first, specifically as follows: ① Locating the abnormal frame range: based on , , Lock the range of abnormal frames within the trajectory; ② Local Restart Multi-Scale Optical Flow Tracking: For the image region corresponding to the abnormal frame, the original algorithm process is reused, namely, reconstructing the local image pyramid and relying on the global active point set. Re-execute forward + reverse LK optical flow tracing with the original global number of the feature; ③ Re-examination and screening: Recalculate the error; if the new location meets the requirements... , , If all thresholds are used, then new points will be used. The original abnormal point was replaced, and the abnormal segment repair was completed. ④ Determine the flow split: If the constraints are still not met after re-checking, proceed to the timing interpolation stage; S46. For abnormal segments that still do not meet all constraints after local re-inspection, time-series interpolation is used to supplement the points, specifically: ① Select adjacent valid and reliable points before and after the abnormal segment (all of which have been verified). () as the interpolation reference; ② Use a temporal interpolation algorithm (linear interpolation / polynomial interpolation, consistent with the trajectory fitting algorithm mentioned above) to calculate the supplementary coordinates of the abnormal frame; ③ Perform a three-level joint verification on the points generated by interpolation. If the verification passes, they are included in the trajectory; if the verification fails, the abnormal segment is determined to be a permanent anomaly. S47. Trajectories that still fail to satisfy all constraints after re-inspection and interpolation are eliminated entirely, including any of the following cases: ① The number of abnormal segments is large and the span is large. Even after multiple re-detection and interpolation, all threshold constraints are still not met. ②The total effective length of the corrected trajectory is still less than the minimum number of frames. ; ③ The fitting error after interpolation continues to exceed the limit. Trajectory continuity cannot be guaranteed; S48. Final output: Set of valid trajectories with feature numbers All points within this trajectory set originate from reliable feature points at each level. The global numbering of feature points remains unchanged throughout the entire process.

[0069] S5. Trajectory Fitting and Temporal Velocity Vector Calculation: The trajectory frame number is converted into physical time. Polynomial fitting is performed on the feature points and the derivative is obtained to obtain the pixel velocity. Combined with the scale coefficient calibrated on site, it is converted into physical velocity. The velocity amplitude and direction angle are calculated to construct a six-dimensional temporal velocity vector. Velocity data is aggregated according to three levels of rules: trajectory level, frame level, and grid level, to generate a gridded temporal velocity vector matrix. Spatial interpolation algorithm is used to complete the discrete velocity densification. Combined with color gradient mapping and vector arrow superposition, standard velocity field cloud map and vector map are generated to complete the kinematic analysis and obtain the total velocity field of temporal observation. For lakes and ocean waters, steps S6 and S7 are executed in sequence to decouple high and low frequency wave velocities. For natural rivers and artificial channels, step S7 is executed, and the total velocity field of temporal observation is used to enter the next stage.

[0070] Step S5 specifically includes: S51. Physical Time Conversion: Extract the frame index number of the valid trajectory. Read the original frame rate of the video , frame index sequence number Converted to physical time:

[0071] in, Representing the k The physical time corresponding to the frame.

[0072] S52. Polynomial Fitting Rules For feature point pixel coordinates , Perform third-order polynomial fitting respectively, and the fitting formula is as follows:

[0073]

[0074] In the formula: , for Directional polynomial coefficients, for Directional polynomial coefficients, For physical time.

[0075] S53. Physical Flow Rate Conversion The pixel velocity is obtained by taking the first derivative of the fitted curve, and then converted into physical velocity by combining it with the scale coefficient s determined on-site.

[0076]

[0077] in, The unit is ; represent The physical velocity component in the X direction at time t, represent The physical velocity component in the Y direction at any given time, in m / s.

[0078] S54. Definition of Six-Dimensional Temporal Velocity Vector Construct a complete six-dimensional temporal velocity vector containing spatiotemporal coordinates, flow velocity, velocity amplitude, and direction angle:

[0079] Where: Resultant velocity amplitude The flow velocity direction angle is calculated using the atan2 function. ; u i , v i These represent the X-axis physical velocity component and the Y-axis physical velocity component, respectively. U i This represents the magnitude of the resultant velocity.

[0080] S55. Three-level data aggregation system ① Trajectory-level output: Collects the temporal velocity vectors of all single feature points to form a trajectory-level dataset. ; ② Frame-by-frame output: Statistically calculate the number of valid feature samples, average velocity, median velocity, and mean velocity component for each frame; ③ Grid-level output: Press Grid-aggregated data to construct a gridded time-series velocity vector matrix The point velocity at the trajectory level is spatially grouped by grid to form a temporal velocity sequence for each grid. Representing the Line 1 Column grid in The aggregate value of the velocity vector at time t, Represented as The matrix; generating a grid velocity dataset. ;in: Representing the r Line 1 The horizontal average velocity component of the column grid, Representing the Line 1 The vertical average velocity component of the column grid. Representing the Line 1 The sum velocity of the grid; S56. Flow Field Visualization Process: Adopted Inverse distance weighting The system uses either radial basis function or kriging interpolation to refine the grid of discrete velocity points; it maps velocity amplitudes to color gradient fields, overlays vector arrows, and generates smooth velocity cloud maps and vector maps; it also exports trajectory history files, statistical tables, visualization images, and videos.

[0081] S6. Scene-Differentiated High- and Low-Frequency Wave Velocity Decoupling. This step is selectively performed based on the type of water body: natural rivers and artificial channels can skip this step and proceed to step S7 using the total velocity field from time-series observations; lakes and oceans must undergo high- and low-frequency wave velocity decoupling before proceeding to step S7. Specifically, high- and low-frequency wave velocity decoupling involves decomposing the total velocity from time-series observations into a composite component of low-frequency plane advection velocity, wave oscillation velocity, and Stokes drift velocity. Time-frequency filtering or Dynamic Mode Decomposition (DMD) is used to separate the high- and low-frequency plane advection velocities, retaining only the low-frequency plane advection velocity to obtain the decoupled velocity field, which is then used for subsequent pressure field calculations.

[0082] S61. The total velocity observed in the time series is decomposed into a low-frequency advection component and a high-frequency wave composite component using time-frequency filtering or dynamic mode decomposition:

[0083] In the formula: This represents the low-frequency plane advection velocity (low-frequency advection component, used for pressure field calculation). This is a composite component of wave oscillation velocity and wave-induced Stokes drift velocity (high-frequency component, not involved in pressure calculation); this method only performs one decoupling of high and low frequencies, and does not perform a second decomposition of wave components.

[0084] ① Time-frequency filtering methods (such as Butterworth low-pass filtering) are used to measure the total velocity of time-series observations. Filtering is performed with a cutoff frequency set to 0.05Hz~1.2Hz to filter out wave oscillation components above the cutoff frequency and retain the plane advection velocity below the cutoff frequency, i.e., low-frequency plane advection velocity. ② Dynamic Mode Decomposition (DMD) is used to analyze the total velocity of time-series observations. A snapshot matrix is ​​constructed, and the dominant dynamic modes are extracted through singular value decomposition. The low-frequency dominant mode (low-frequency advection component) and the high-frequency oscillation mode (wave component) are separated according to frequency characteristics. The frequency range of the low-frequency mode is consistent with that of the filtering method. Choose one of the two methods to execute, with time-frequency filtering being preferred, and DMD being used in strong noise scenarios to improve decoupling robustness.

[0085] S63. Set the cross-scale velocity deviation threshold. The unit is m / s, and the value is selected and configured according to the scenario. It is used to verify the validity of the decoupling results. The low-frequency plane advection velocity obtained by high and low frequency decoupling is... It is necessary to perform a deviation check between the coarse-scale / fine-scale velocity obtained from the original multi-scale trajectory transfer: If cross-scale velocity deviation This indicates that the low-frequency plane advection velocity is consistent with the overall motion trend, and the decoupling is effective. In the formula... Output the original flow rate for multi-scale optical flow; like This indicates that the decoupling algorithm failed due to noise / strong wave interference, so this step should be skipped and the original speed should be used. This serves as a low-frequency planar advection velocity to prevent erroneous data from being included in subsequent pressure field calculations. Cross-scale velocity deviation threshold The possible values ​​are as follows: The random error of conventional water surface video speed measurement is generally ≤0.05~0.10 m / s; Under strong wave and turbulent flow conditions, the overall error can reach 0.10~0.20 m / s; Multi-scale inversion of ocean current fields has an inherent bias of approximately 0.05–0.6 m / s; After decoupling, only the pure low-frequency plane advection velocity field is retained as input to the pressure field solution stage.

[0086] S7. Dynamics of Relative Pressure Field at Water Surface: Based on the time-series observed total velocity field / decoupled flow velocity field, first perform divergence correction on the velocity field to satisfy the requirements of incompressible fluid. The physical constraints are defined; the Poisson pressure equation is derived based on the incompressible Navier-Stokes equations, and the computational domain is discretized using the finite element method; three types of boundary conditions are applied sequentially: free water surface, solid wall, and inflow / outflow open boundary, and the linear equations are solved to obtain the relative pressure field of the entire water surface; the pressure field is interpolated, refined, and color-gradientd to output a standard pressure contour map, completing the analysis of the entire dynamic process. Specifically, this includes: S71. Incompressible Fluid Pre-Constraint: Perform divergence correction on the input time-series observed total velocity field / decoupled flow velocity field to strictly satisfy the fluid incompressibility condition: ;in, Represents the gradient / divergence differential operator. Represents the velocity vector; S72. Derivation of Poisson's pressure equation: The Poisson pressure equation is derived from the incompressible Navier-Stokes equations.

[0087] In the formula: The relative pressure field of the water surface is represented, with local atmospheric pressure as the reference zero point; the water density is a fixed value. ; It is a velocity vector; For the Laplace operator; S73. Finite Element Mesh Rules: The computational domain is discretized using the finite element method, and the finite element mesh is compared with the front-end velocity statistics. The grid remains completely consistent; S74. Three types of mandatory boundary conditions (applied uniformly across the entire domain): ① Dynamic boundary of free water surface: ; ② Bank / Hydraulic Structure Solid Boundary: ; ③ Open boundaries of river / waterway inflow and outflow: Directly define the relative pressure field of the water surface at the boundary. Or pressure normal derivative ; in: The normal vector is the unit vector outside the boundary, and the gravitational acceleration has a fixed value. ; S75. Equation Solving and Post-processing: Solving the linear equation system yields the global relative pressure distribution. Interpolation and color gradient mapping are performed on the pressure field to output a standard water surface relative pressure cloud map.

[0088] Next, the publicly available CRR river dataset from the U.S. Geological Survey (USGS) was used to verify the velocity measurement accuracy, flow field reconstruction, and relative pressure field solution capabilities of the method of this invention in natural and complex waters. Quantitative verification was completed based on ADCP measured data.

[0089] (1) Experimental conditions and equipment 1) Observation area: The natural channel of the Colorado River in the United States, with a width of about 49.38m and a significant velocity gradient in the water flow, which is a typical natural river condition.

[0090] 2) Acquisition equipment: The drone takes photos from the front, without perspective distortion, and no additional perspective correction is required.

[0091] 3) Video parameters: resolution 1920×1080, frame rate 23.98fps, total frames 767.

[0092] 4) Physical calibration: Scale factor s = 0.018 m / pixel.

[0093] 5) Verification benchmark: ADCP near-surface velocity data (hydrological standard measurement equipment).

[0094] (2) Algorithm parameter configuration 1) Preprocessing: grayscale conversion, CLAHE contrast enhancement, 5×5 Gaussian smoothing to construct an effective water surface mask.

[0095] 2) Gridded feature initialization: Grid R×C=20×25 grids, maximum number of features per grid N cell =2, each grid has at least 2 feature points, using Shi Tomagi Corner detection algorithm.

[0096] 3) Multi-scale optical flow tracking: number of pyramid layers L max =4, LK window ω=30×30, forward and backward consistency threshold =2 pixels, prediction bias threshold =4 pixels.

[0097] 4) Trajectory reliability constraint: minimum trajectory length =5 frames, cubic polynomial fitting m=3, fitting deviation threshold =1 pixel.

[0098] 5) Pressure calculation: water density The acceleration due to gravity is g = 9.8 m / s². 2 The solution is based on the incompressible Navier-Stokes equations and the finite element method.

[0099] (3) Processing flow 1) Perform ROI extraction, grayscale conversion, CLAHE enhancement, and Gaussian smoothing on a continuous image sequence to generate an effective water surface mask; 2) Initialize uniformly meshed feature points within the mask to ensure uniform spatial coverage; 3) Construct a 4-layer image pyramid, perform multi-scale LK optical flow tracking from coarse to fine, and perform dynamic point filling on continuous low-density grids; 4) After three levels of constraints—forward and backward consistency, prediction guidance, and trajectory continuity—abnormal trajectories are eliminated to obtain a set of valid trajectories; 5) Perform cubic polynomial fitting and first-order differentiation on the effective trajectory to convert it into a physical time-series velocity vector; 6) The meshed temporal velocity vector matrix is ​​obtained by aggregation, and the flow field diagram is generated by IDW interpolation and densification. The velocity field visualization effect is shown in the figure below. Picture 3 As shown; 7) Divergence correction is applied to the velocity field to satisfy the fluid incompressibility constraint. The relative pressure field at the water surface is solved based on the Poisson equation and the finite element method. The visualization of the relative pressure field at the surface is shown in the figure below. Picture 4 As shown.

[0100] (4) Experimental results 1) The method of this invention has an average relative error of -1.76% relative to the ADCP reference, an average absolute error (MAE) of 0.085 m / s, and a root mean square error (RMSE) of 0.109 m / s. 2) It can accurately reproduce the natural flow velocity distribution of the riverbank with low flow velocity in the middle and high flow velocity in the middle. The curve is smooth without abnormal jumps. The stability of the weak textured riverbank area is significantly better than the traditional STIV and OFV methods. 3) Output gridded temporal velocity field and relative pressure field cloud map of water surface, fully realizing the reconstruction of the physical field of the entire process from temporal trajectory → velocity → pressure.

[0101] This embodiment demonstrates that the method of the present invention has high velocity measurement accuracy, good temporal stability, and strong adaptability to weak texture under natural and complex river conditions. It can complete the reconstruction of standardized and reproducible surface flow field and relative pressure field of water bodies, meeting the actual observation needs of river engineering.

[0102] Based on the aforementioned method for reconstructing the surface flow field of water bodies using time-series trajectory fitting, this embodiment also proposes a system for reconstructing the surface flow field of water bodies using time-series trajectory fitting. The system consists of eight mutually cooperating functional modules, fully covering the entire algorithm process, kinematics, scene decoupling, and dynamics. The functions of each module are as follows: Image preprocessing module: Acquires continuous image sequences of the water area to be tested. The images are then sorted and preprocessed for each frame to construct an effective water surface mask. Output the preprocessed image sequence .

[0103] Mesh feature initialization module: Divides the current measurement area covered by the effective water surface mask into... Uniform grid, constructing a locally effective mask for each grid. Feature points are detected within a locally effective mask, and each feature point is assigned a globally unique index. By setting single-grid feature point constraint rules to constrain the number of features in a single grid to eliminate blind spots in observation, an initial set of spatially uniformly distributed feature points is obtained. .

[0104] Multi-scale optical flow tracking and dynamic point filling module: for preprocessed image sequences The maximum number of layers to build is The image pyramid, the image pyramid downsampling rule is... Perform forward / backward Lucas-Kanade optical flow tracking from coarse to fine, with coarse-scale trajectory. pass Sampling is used as the initial value for fine scale; the global active point set (Active) is maintained frame by frame, while the global index of the feature points remains unchanged; for continuous... The low-density grid is dynamically padded with points; after single-level processing, the trajectory is propagated to finer scales as motion priors until the finest scale is reached, and finally, a full-resolution trajectory set is output. .

[0105] Trajectory reliability constraint module: for full-resolution trajectory sets The system performs a three-level joint verification of forward and backward consistency constraints, prediction guidance constraints, and trajectory continuity constraints. Simultaneously, it incorporates trajectory length thresholds to remove abnormal trajectories and correct abnormal segments, outputting a set of valid trajectories with feature-indexed numbers. .

[0106] The trajectory fitting and temporal velocity vector solution module converts the frame number of the effective trajectory into physical time, performs polynomial fitting on the feature points and calculates the derivative to obtain the pixel velocity, combines the scale coefficient calibrated on site to convert it into physical velocity, calculates the velocity amplitude and direction angle, and constructs a six-dimensional temporal velocity vector; aggregates velocity data according to three levels of rules: trajectory level, frame level, and grid level, and generates a gridded temporal velocity vector matrix; uses a spatial interpolation algorithm to complete the discrete velocity densification, and combines color gradient mapping and vector arrow overlay to generate standard velocity field cloud map and vector map, completing the analysis of the entire kinematic process.

[0107] The flow field visualization module is responsible for using IDW / RBF / Kriging interpolation to complete flow velocity field encryption, color gradient mapping, and vector arrow overlay, generating flow velocity cloud maps and vector maps, and exporting trajectory files, statistical reports, images, and video results.

[0108] Wave velocity decoupling module: A dedicated velocity decoupling module for lakes or oceans, used to decompose the total velocity of time-series observations into a composite component of low-frequency plane advection velocity, wave oscillation velocity and Stokes drift velocity. It uses time-frequency filtering or dynamic mode decomposition to separate high- and low-frequency plane advection velocities, retaining only the low-frequency plane advection velocity as the decoupled velocity field for subsequent pressure field solution.

[0109] Pressure field solution module: Based on the time-series observed total velocity field / decoupled flow velocity field, it first performs divergence correction on the velocity field to satisfy the requirements of incompressible fluid. The physical constraints are determined; the Poisson pressure equation is derived based on the incompressible Navier-Stokes equations, and the computational domain is discretized using the finite element method; three types of boundary conditions are applied sequentially: free water surface, solid wall, and inflow / outflow open boundary, and the linear equations are solved to obtain the relative pressure field of the water surface in the whole field; the relative pressure field of the water surface in the whole field is interpolated, densified, and color-gradientd to output a standard pressure cloud map, thus completing the analysis of the entire dynamic process.

[0110] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0111] The embodiments described above are merely illustrative of several implementations of the present invention, and while the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the invention patent. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these all fall within the protection scope of the present invention. Therefore, the protection scope of this invention patent should be determined by the appended claims.

Claims

1. A method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting, characterized in that, It includes: S1. Obtain and sort the continuous image sequence of the water area to be tested, preprocess each frame of the image, and then construct an effective water surface mask to obtain the preprocessed image sequence; S2. Divide the flow measurement area effectively covered by the water surface mask into... A uniform grid is used to construct a locally effective mask for each grid. Feature points are detected in the locally effective mask. A globally unique index is assigned to each feature point. Constraint rules for single-grid feature points are set to constrain the number of features in a single grid to eliminate blind spots and obtain an initial set of spatially uniformly distributed feature points. S3. Construct a maximum of [number] layers from the preprocessed image sequence. The image pyramid is constructed; the global active point set (Active) is initialized using an initial feature point set. The initialization method is to clear the Active container, traverse all feature points in the initial feature point set, store all feature points and their indexes in Active, and mark the initial state as updated. The Active set manages all feature points uniformly throughout the process; forward / backward Lucas-Kanade optical flow tracking is performed from coarse to fine; the global active point set (Active) is maintained frame by frame, and the global indexes of the feature points remain unchanged. For continuous Dynamic point filling in low-density grids within frames; After single-level processing is completed, the trajectory is propagated to finer scales as motion priors until the finest scale is reached, and finally a full-resolution trajectory set is output. S4. Perform three-level joint verification on the full-resolution trajectory set, and combine the trajectory length threshold to complete the abnormal trajectory removal and abnormal segment correction, and output the effective trajectory set with feature number index; S5. Convert the frame number of the effective trajectory into physical time, perform polynomial fitting on the feature points and differentiate to obtain the pixel velocity, combine it with the scale coefficient of the field calibration to convert it into physical velocity, calculate the velocity amplitude and direction angle, and construct a six-dimensional temporal velocity vector; aggregate the velocity data according to the three-level rules of trajectory level, frame level, and grid level to generate a gridded temporal velocity vector matrix; use spatial interpolation algorithm to complete the discrete velocity densification, combine color gradient mapping and vector arrow superposition to generate standard velocity field cloud map and vector map, complete the kinematic analysis of the whole process, and obtain the total velocity field of temporal observation; for lakes and ocean waters, execute steps S6 and S7 in sequence. For natural rivers and artificial channels, proceed to step S7; S6. The total velocity observed in the time series is decomposed into a composite component of low-frequency plane advection velocity, wave oscillation velocity and Stokes drift velocity. Time-frequency filtering or dynamic mode decomposition is used to separate the high-frequency and low-frequency plane advection velocities, and only the low-frequency plane advection velocity is retained as the decoupled velocity field. S7. Based on the time-series observed total velocity field / decoupled flow velocity field, divergence correction is performed on the velocity field; the Poisson pressure equation is derived based on the incompressible Navier-Stokes equations, and the computational domain is discretized using the finite element method; three types of boundary conditions are applied sequentially: free water surface, solid wall, and inflow / outflow open boundary, and the linear equations are solved to obtain the relative pressure field of the entire water surface; the relative pressure field of the entire water surface is interpolated, densified, and color-gradientd, and a standard pressure cloud map is output to complete the analysis of the entire dynamic process.

2. The method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting according to claim 1, characterized in that, Step S1 specifically includes: S11. Acquire a continuous image sequence of the water area to be measured. And sorted by filename completion time. Indicates the first t Original frame image, 0≤ t ≤ T -1, T Represents the total number of frames in the image sequence; S12. For each frame of the original image Perform ROI extraction → perspective correction → grayscale conversion → CLAHE contrast enhancement → Gaussian smoothing; for scenes with extremely strong reflections and missing textures on the water surface, a high-boost filter can be selectively added for edge enhancement processing; S13. Construct an effective water surface mask: In the formula, A region of interest mask for the water surface, used to delineate the target flow measurement range; A manually labeled non-water surface exclusion mask used to remove invalid areas; symbol For set intersection operation, the symbol is... This is the negation operation; S14. Finally, the preprocessed image sequence is obtained. , Representing the t The image after frame preprocessing.

3. The method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting according to claim 1, characterized in that, Step S2 specifically includes: S21. Divide the current measurement area covered by the effective water surface mask into... Uniform grid, a single grid is denoted as , where: row index Column index For each grid Constructing a locally effective mask ; S22. In the local effective mask Internal use Shi Tomasi Corner detection algorithms detect feature points. Shi Tomasi Corner detection algorithm configuration parameters include quality coefficient Minimum spacing between feature points Boundary shielding parameters , The value ranges from 0.01 to 0.

1. Values ​​range from 5 to 20 pixels; each detected feature point is assigned a globally unique and unchanging index number. Set single-grid feature point constraint rules: The initial number of feature points within a single grid must satisfy: , The actual number of feature points in the grid at row r and column c. The maximum number of feature points per grid cell is set; a minimum of two feature points per grid cell is enforced to eliminate blind spots in observation. S23. Remove coordinates located on the effective water surface mask. All feature points other than those in the original set are combined to form an initial dynamic feature point set that is spatially uniformly distributed. : In the formula: Representing the The pixel coordinates of each feature point This represents the total number of initial feature points across the entire region.

4. The method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting according to claim 1, characterized in that, Step S3 specifically includes: S31. The preprocessed image sequence The maximum number of layers to build is The image pyramid, and the image pyramid downsampling rule are: In the formula, Representing the The first layer of the pyramid Frame image; Representing the The first layer of the pyramid Frame image; hierarchical index satisfy ,in For the finest scale, This is the coarsest scale; The downsampling function scales down the width and height of the previous image proportionally to obtain the next image. S32. From the coarsest scale of the pyramids We begin by using a coarse-to-fine strategy to traverse layer by layer towards finer scales: (1) Forward LK optical flow tracking At the current level Internal forward LK optical flow tracing is performed, by the first Frame Image To the Frame Image The formula for estimating the motion of feature points is: in, For the first Layer Frame number Coordinates of feature points The first digit obtained by forward tracking Layer Frame number Predict the coordinates of each feature point; (2) Reverse LK optical flow backtracking verification At the current level Internally, a reverse LK optical flow backtracking verification is performed on the forward tracking results, by the first... Frame Image To the Frame Image Backtracking, the formula is: in, The first obtained by reverse tracing Layer t, frame t Backtracking coordinates of each feature point; calculate and error Only retain those with errors less than the forward-backward consistency threshold. The reliable trajectory is obtained at the current level. Frame number Coordinates of reliable feature points ; (3) Cross-scale trajectory transfer The reliable feature point trajectory results filtered at the current level are upsampled and passed to the next level as motion priors, as shown in the formula: in: For the first Layer Frame number Coordinates of a reliable feature point; For the first Layer Frame number The initial coordinates of each feature point are predicted; coefficient 2 is the inverse scaling factor for downsampling the image pyramid by 2 times, corresponding to the mapping relationship from coarse-scale coordinates to fine-scale coordinates; S33. Count the number of active feature points in each frame grid by grid, for continuous... If a frame is determined to be a low-density mesh, a dynamic point-filling operation is performed, as follows: (1) Active point statistics: After optical flow tracking is completed for each frame, statistics are performed grid by grid. Frame grid The number of active feature points within the range is denoted as ; (2) Low-density mesh determination criteria: Set the sparseness determination ratio coefficient If a single grid is continuous Frames satisfy inequalities If the grid is low-density, then feature points need to be added. (3) Point replacement objective constraint: Reapply the low-density grid Shi Tomasi Corner detection algorithms detect feature points and supplement them to... ,in Fill in the proportional coefficient for dynamically added points; (4) Constraint on the number of feature points added in a single operation: The maximum number of feature points added in a single operation is [number missing]. And satisfy To avoid local clustering of feature points; (5) New feature point management: Assign a brand-new globally unique number index to the newly detected feature points, and directly include them in the optical flow tracking queue of all subsequent frames to ensure uniform coverage of features across the entire domain; S34. Maintain the global active point set Active frame by frame, distinguishing between three states: inflow, outflow, and update; once a feature point is assigned an index, it is never modified, repeated, or reset throughout the process; details are as follows: (1) Initialization phase After the pyramid is constructed, initialize the global active point set. Load the initial feature point set Assign a unique global index to each feature point. ; (2) Coordinate with forward LK optical flow, reverse LK optical flow, and reliability screening process At each level When performing forward LK optical flow, reverse LK optical flow, and reliability screening: ① Tracking is effective and the error is less than the forward and backward consistency judgment threshold. Feature points: Update location information, mark status as "updated", retain the ID, and continue to remain. gather; ② Tracking failure, error exceeding the forward and backward consistency judgment threshold Feature point: The state is marked as "moved out", from Remove from the set, terminating the trajectory of that point; (3) Coordinating cross-scale trajectory transfer Coarse-scale reliable feature points When upsampling is passed to a finer scale: ① The existing feature point global numbering is used directly without renumbering; ② The set inherits the batch of feature points synchronously, and the state continues; (4) Dynamic point replenishment in conjunction with low-density grid For continuous When the low-density grid of the frame is used to supplement new feature points: ① The status of newly added feature points is marked as "inflow"; ② Assign a brand new, unique global ID and add it. Gather and participate in subsequent tracking; S35. After single-level processing, the trajectory is propagated to finer scales as motion priors until the finest scale is reached, finally outputting a full-resolution trajectory set. .

5. The method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting according to claim 1, characterized in that, In step S4, the three-level joint verification includes forward and backward consistency constraints, prediction guidance constraints, and trajectory continuity constraints, specifically: S41. Level 1: Forward and Backward Consistency Constraints Forward and backward errors The calculation formula is: in, European distance; Representing the Layer Frame number The coordinates of reliable feature points are taken from the global active point set. ; For the first Layer Frame number The coordinates of each feature point are traced back by reverse optical flow. when If the condition is met, the feature point is determined to be a reliable point and is directly retained; otherwise, it proceeds to the second-level constraint verification. The threshold for determining forward and backward consistency has a value range of 1 to 3 pixels. S42. Level Two: Predictive Guiding Constraints Using cross-scale predicted location to correct trajectory points, calculate prediction guidance bias. : in, Based on reliable feature point coordinates The first obtained by forward LK optical flow tracing Layer Frame number Predict the coordinates of each feature point; The first result obtained for cross-scale transfer Layer Frame number Predict the coordinates of each feature point; If the correction simultaneously satisfies and If the feature point is positive, it is marked as a recovery point and retained; otherwise, the feature point is directly discarded. The threshold for predicting guidance deviation is set to a value ranging from 2 to 5 pixels. S43. Level 3: Trajectory Continuity Constraints Perform polynomial fitting on the complete trajectory and calculate the root mean square error of the fit. : In the formula: For the trajectory number k Frame 1 i Measured feature coordinates of the feature points; For the trajectory number k Frame 1 i Fitted coordinates of feature points This represents the total number of valid frames for the trajectory. like If the trajectory is continuous and reliable, it is considered a continuous and reliable trajectory; otherwise, it is considered an abnormal trajectory; fitting deviation threshold. Value range: 0.5~2 pixels; S44. Set minimum trajectory length ,Require Frame; effective trajectory length less than If a trajectory is deemed abnormal, the entire segment is removed. The effective length of the trajectory is statistically analyzed based on reliable feature points. The generated continuous trajectory segments; S45. For abnormal trajectories with excessive fitting deviation or insufficient length, local re-examination should be performed on the abnormal segments first, specifically as follows: ① Locating the abnormal frame range: based on , , Lock the range of abnormal frames within the trajectory; ② Local Restart Multi-Scale Optical Flow Tracking: For the image region corresponding to the abnormal frame, the original algorithm process is reused, namely, reconstructing the local image pyramid and relying on the global active point set. Re-execute forward + reverse LK optical flow tracing with the original global number of the feature; ③ Re-examination and screening: Recalculate the error; if the new location meets the requirements... , , If all thresholds are used, then new points will be used. The original abnormal point was replaced, and the abnormal segment repair was completed. ④ Determine the flow split: If the constraints are still not met after re-checking, proceed to the timing interpolation stage; S46. For abnormal segments that still do not meet all constraints after local re-inspection, time-series interpolation is used to supplement the points, specifically: ① Select the adjacent valid and reliable points before and after the abnormal segment as the interpolation reference; ② Use a temporal interpolation algorithm to calculate the supplementary coordinates of the abnormal frame; ③ Perform a three-level joint verification on the points generated by interpolation. If the verification passes, they are included in the trajectory; if the verification fails, the abnormal segment is determined to be a permanent anomaly. S47. Trajectories that still fail to satisfy all constraints after re-inspection and interpolation are eliminated entirely, including any of the following cases: ① The number of abnormal segments is large and the span is large. Even after multiple re-detection and interpolation, all threshold constraints are still not met. ②The total effective length of the corrected trajectory is still less than the minimum number of frames. ; ③ The fitting error after interpolation continues to exceed the limit. Trajectory continuity cannot be guaranteed; S48. Final output: Set of valid trajectories with feature numbers All points within this trajectory set originate from reliable feature points at each level. .

6. The method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting according to claim 1, characterized in that, Step S5 specifically includes: S51. Physical Time Conversion: Extract the frame index number of the valid trajectory. Read the original frame rate of the video , frame index sequence number Converted to physical time: in, Representing the k The physical time corresponding to the frame; S52. Polynomial Fitting Rule: For feature point pixel coordinates , Perform third-order polynomial fitting respectively, and the fitting formula is as follows: In the formula: , for Directional polynomial coefficients, for Directional polynomial coefficients, Physical time; S53. Physical velocity conversion: The pixel velocity is obtained by taking the first derivative of the fitted curve, and then converted into the actual physical velocity by combining it with the scale coefficient s determined on-site. in, The unit is ; represent The physical velocity component in the X direction at time t, represent The physical velocity component in the Y direction at any given time, in m / s; S54. Definition of a six-dimensional temporal velocity vector: Constructing a complete six-dimensional temporal velocity vector containing spatiotemporal coordinates, velocity, velocity amplitude, and velocity direction angle: Where: Resultant velocity The flow velocity direction angle is calculated using the four-quadrant arctangent function. ; u i , v i These represent the X-axis physical velocity component and the Y-axis physical velocity component, respectively. U i Represents the magnitude of the resultant velocity; S55. Three-level data aggregation system: ① Trajectory-level output: Collects all six-dimensional temporal velocity vectors of single feature points to form a trajectory-level dataset. ; ② Frame-by-frame output: Statistically calculate the number of valid feature samples, average velocity, median velocity, and mean velocity component for each frame; ③ Grid-level output: Press Grid-aggregated data to construct a gridded time-series velocity vector matrix The point velocity at the trajectory level is spatially grouped by grid to form a temporal velocity sequence for each grid. Representing the Line 1 Column grid in The aggregate value of the velocity vector at time t, Represented as The matrix; generating a grid velocity dataset. ;in: Representing the r Line 1 The horizontal average velocity component of the column grid, Representing the Line 1 The vertical average velocity component of the column grid. Representing the Line 1 The sum velocity of the grid; S56. Flow Field Visualization Process: Adopted Inverse distance weighting The system uses either radial basis function or kriging interpolation to refine the grid of discrete velocity points; it maps velocity amplitudes to color gradient fields, overlays vector arrows, and generates smooth velocity cloud maps and vector maps; it also exports trajectory history files, statistical tables, visualization images, and videos.

7. The method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting according to claim 1, characterized in that, Step S6 is as follows: S61. Use time-frequency filtering or dynamic mode decomposition to calculate the total velocity of time-series observations. Decomposed into a low-frequency advection component and a high-frequency wave component: In the formula: This refers to low-frequency planar advection velocity. This is the composite component of wave oscillation velocity and wave-induced Stokes drift velocity; this step only performs one decoupling of high and low frequencies, and does not further decompose the wave components. The time-frequency filtering method is as follows: [The text abruptly ends here, likely due to an incomplete sentence or a missing section.] Filtering is performed with a cutoff frequency set to 0.05Hz~1.2Hz to filter out wave oscillation components above the cutoff frequency and retain the plane advection velocity below the cutoff frequency, which is the low-frequency plane advection velocity. The dynamic mode decomposition method is as follows: The total velocity of time-series observations is... Construct a snapshot matrix, extract the dominant dynamic mode through singular value decomposition, and separate the low-frequency dominant mode and high-frequency oscillation mode according to frequency characteristics. The frequency range of the low-frequency mode is consistent with that of the filtering method. S63. Set the cross-scale velocity deviation threshold. The unit is m / s, and the value is selected and configured according to the scenario. It is used to verify the validity of the decoupling results. The low-frequency plane advection velocity obtained by high and low frequency decoupling is... It is necessary to perform a deviation check between the coarse-scale / fine-scale velocity obtained from the original multi-scale trajectory transfer: If cross-scale velocity deviation This indicates the low-frequency plane advection velocity. Consistent with the overall movement trend, decoupling is effective, where... Output the original flow rate for multi-scale optical flow; like This indicates that the decoupling algorithm failed due to noise / strong wave interference, so this step was skipped, and the total velocity was observed using time series. This serves as a low-frequency planar advection velocity to prevent erroneous data from being included in subsequent pressure field calculations. After decoupling, only the pure low-frequency plane advection velocity field is retained as the decoupled velocity field and input into the pressure field solution stage.

8. The method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting according to claim 7, characterized in that, The cross-scale velocity deviation threshold The possible values ​​are as follows: The random error of conventional water surface video speed measurement is generally ≤0.05~0.10 m / s; Under strong wave and turbulent flow conditions, the overall error can reach 0.10~0.20 m / s; The ocean current field is inverted at multiple scales with an inherent bias of approximately 0.05–0.6 m / s.

9. The method for reconstructing the surface flow field of water bodies based on time-series trajectory fitting according to claim 1, characterized in that, Step S7 specifically includes: S71. Incompressible Fluid Pre-Constraint: Perform divergence correction on the input time-series observed total velocity field / decoupled flow velocity field to strictly satisfy the fluid incompressibility condition: ;in, Represents the gradient / divergence differential operator. Represents the velocity vector; S72. Derivation of Poisson's pressure equation: The Poisson pressure equation is derived from the incompressible Navier-Stokes equations. In the formula: This represents the relative pressure field over the water surface, with local atmospheric pressure as the reference zero point. The density of the water is a fixed value. ; For the Laplace operator; S73. Finite Element Mesh Rules: The computational domain is discretized using the finite element method, and the finite element mesh is compared with the front-end velocity statistics. The grid remains completely consistent; S74. Three types of mandatory boundary conditions: ① Dynamic boundary of free water surface: ; ② Bank / Hydraulic Structure Solid Boundary: ; ③ Open boundaries for river / water inflow and outflow: Directly specify the relative pressure at the boundary. Or pressure normal derivative ; in: The normal vector is the unit vector outside the boundary, and the gravitational acceleration has a fixed value. ; S75. Equation Solving and Post-processing: Solving the linear equation system yields the global relative pressure distribution. Interpolation and color gradient mapping are performed on the pressure field to output a standard water surface relative pressure cloud map.

10. A system for reconstructing the surface flow field of a body of water based on time-series trajectory fitting, characterized in that, It is used to implement the water surface flow field reconstruction method based on time-series trajectory fitting as described in any one of claims 1-9, comprising: The image preprocessing module is used to acquire and sort a continuous image sequence of the water area to be tested, preprocess each frame of the image, construct an effective water surface mask, and output the preprocessed image sequence. The gridded feature point initialization module is used to divide the current measurement area covered by the effective water surface mask into... A uniform grid is used to construct a locally effective mask for each grid. Feature points are extracted from the locally effective mask using a feature detection algorithm. A globally unique index is assigned to each feature point. Constraint rules are set for feature points within a single grid to constrain the number of features in a single grid in order to eliminate blind spots in observation, resulting in an initial set of spatially uniformly distributed feature points. The multi-scale optical flow tracking and dynamic point filling module is used to construct a maximum number of layers for preprocessed image sequences. The image pyramid is used; the global active point set (Active) is initialized using an initial feature point set. Initialization involves clearing the Active container, traversing all feature points in the initial feature point set, storing all feature points and their indexes in Active, and marking the initial state as updated. The Active set manages all feature points uniformly throughout the process. Forward / backward Lucas-Kanade optical flow tracking is performed from coarse to fine. The global active point set (Active) is maintained frame by frame, while the global indexes of the feature points remain unchanged. For continuous... The low-density grid of the frame is dynamically filled with points; after the single-level processing is completed, the trajectory is propagated to a finer level as a motion prior until the finest scale is traversed, and finally the full-resolution trajectory set is output. The trajectory reliability constraint module is used to perform three-level joint verification on the full-resolution trajectory set, and at the same time, combined with the trajectory length threshold, it completes the abnormal trajectory elimination and abnormal segment correction, and outputs a valid trajectory set with feature number index. The trajectory fitting and temporal velocity vector solution module converts the frame number of the effective trajectory into physical time, performs polynomial fitting on the feature points and calculates the derivative to obtain the pixel velocity, combines the scale coefficient of the field calibration to convert it into physical velocity, calculates the velocity amplitude and direction angle, and constructs a six-dimensional temporal velocity vector. It aggregates velocity data according to three levels of rules: trajectory level, frame level, and grid level, and generates a gridded temporal velocity vector matrix. It uses a spatial interpolation algorithm to complete the discrete velocity densification, and combines color gradient mapping and vector arrow superposition to generate standard velocity field cloud map and vector map, completes the kinematic analysis of the whole process, and obtains the total velocity field of temporal observation. The wave velocity decoupling module is used to decompose the total velocity of time-series observations into a composite component of low-frequency plane advection velocity, wave oscillation velocity, and Stokes drift velocity. It uses time-frequency filtering or dynamic mode decomposition to separate high- and low-frequency plane advection velocities, retaining only the low-frequency plane advection velocity as the decoupled velocity field. The pressure field solution module is used to perform divergence correction on the velocity field based on the time-series observed total velocity field / decoupled flow velocity field; derive the Poisson pressure equation based on the incompressible Navier-Stokes equations; and discretize the computational domain using the finite element method; apply three types of boundary conditions in sequence—free water surface, solid wall, and inflow / outflow open boundary—and solve the linear equations to obtain the relative pressure field of the entire water surface; perform interpolation refinement and color gradient processing on the relative pressure field of the entire water surface, output a standard pressure contour map, and complete the analysis of the entire dynamic process.

Images