Backflow surface velocity field algorithm

Abstract

The invention relates to the technical field of hydraulic engineering, in particular to a backflow surface velocity field algorithm. Through a particle image velocity measurement technology and in combination with an advanced image processing algorithm, automatic identification and tracking of a complex flow state of a Three Gorges change backwater area are realized. Through key technologies of image preprocessing, feature point matching, flow field interpolation and the like, visual display of a flow velocity vector diagram, the flow velocity and the flow velocity direction can be provided, so that powerful technical support is provided for reservoir operation management, ship navigation safety guidance, environment monitoring and scientific research; besides, the provided backflow surface velocity field algorithm not only can improve the flow state identification capability of the Three Gorges Reservoir, but also is expected to be popularized to other large-scale water conservancy projects and applied to ship navigation safety guidance, and has wide application prospects and commercial values.

Description

Technical Field

[0001] The invention relates to the technical field of water conservancy engineering, in particular to a backflow surface velocity field algorithm. Background Art

[0002] As one of the largest artificial reservoirs in the world, the complex hydrodynamic conditions of the Three Gorges Reservoir have had a profound impact on shipping safety, ecological protection, water resources management, etc. In particular, the flow pattern in the variable backwater area has always been a difficult and hot issue in the study of water conservancy engineering and fluid mechanics due to its large dynamic changes, many influencing factors and high prediction difficulty.

[0003] In the variable backwater area of ​​the Three Gorges Reservoir, the flow velocity is unevenly distributed and the vortex structure is complex. These characteristics pose a threat to the safety of ship navigation and also affect the water quality and ecosystem health of the reservoir area. Therefore, accurately identifying and analyzing these complex flow states is of great significance for optimizing waterway management, formulating scientific scheduling strategies, and protecting the aquatic ecological environment.

[0004] Traditional flow pattern identification methods, such as field observations and physical model tests, are often time-consuming, labor-intensive, and costly, and have difficulty capturing the dynamic changes of the flow field.

[0005] Based on this, the present invention provides a reflow surface velocity field algorithm to solve the above-mentioned technical problems. Summary of the invention

[0006] The purpose of the present invention is to provide a backflow surface velocity field algorithm, which realizes automatic recognition and tracking of the complex flow state of the Three Gorges variable backwater area through particle image velocimetry (PIV) technology combined with advanced image processing algorithms; through key technologies such as image preprocessing, feature point matching, flow field interpolation, etc., it can provide a flow velocity vector diagram, a direct display of the flow velocity size and direction, thereby providing strong technical support for reservoir operation management, ship navigation safety guidance, environmental monitoring and scientific research; in addition, the implementation of this technical solution can not only enhance the flow state recognition capability of the Three Gorges Reservoir, but also is expected to be extended to other large-scale water conservancy projects and applied to ship navigation safety guidance, and has broad application prospects and commercial value.

[0007] To achieve the above object, the present invention provides the following technical solutions:

[0008] The present invention provides a reflow surface velocity field algorithm, comprising the following steps:

[0009] S1. Initialize and set the path: clear MATLAB environment variables, close warning messages, get the current working directory and build the main path;

[0010] S2, image reading: use the imread function to read the image from the specified folder path;

[0011] S3, image preprocessing: enhance the read image through grayscale conversion, median filtering, image subtraction and histogram equalization;

[0012] S4, PIV flow field calculation: PIV algorithm based on FFT is used to calculate the velocity component of the flow field;

[0013] S5, analysis of processing PIV flow field calculation: remove outliers exceeding 3 times the standard deviation and fill in NaN values. At the same time, use a median filter to smooth the velocity vector field;

[0014] S6. Visualize the flow field results: Use the quiver and pcolor functions to display the velocity vector and velocity magnitude of the flow field;

[0015] S7, SURF feature point detection and matching: Track the correspondence between images through SURF feature point matching.

[0016] The present invention is further configured as follows: in step S4, the specific process is as follows:

[0017] S41, initialization parameters: output prompt information, create a cell array s for storing parameter settings for PIV analysis;

[0018] S42, perform PIV analysis: input the pre-processed image and the set parameters, and perform PIV calculation. The function calculates and returns the x and y corresponding to the x-coordinate and y-coordinate of the flow field, and the corresponding velocity components u and v;

[0019] S43, filling the missing values ​​NaN in the velocity components u and v;

[0020] S44, output results: print information, and the PIV flow field calculation is completed.

[0021] The present invention is further configured as follows: in step S42, the piv_FFTmulti function based on the Fast Fourier Transform FFT cross-correlation algorithm is called to perform PIV calculation.

[0022] 4. The reflow surface velocity field algorithm according to claim 2, characterized in that, in step S43, the missing values ​​NaN in the velocity components u and v are filled using 4th nearest neighbor interpolation through the inpaint_nans function.

[0023] The present invention is further configured as follows: in step S5, the specific process is as follows:

[0024] S51, calculate standard deviation: use the std function to calculate the standard deviation of the velocity vector fields u and v respectively, and store them respectively;

[0025] S52, remove outliers;

[0026] S53, fill NaN value: fill NaN value by interpolation method;

[0027] S54, smoothing: using the median filter medfilt2 to smooth the velocity vector field;

[0028] S55, calculate the average speed: use the median filter medfilt2 again to calculate the average speeds u_mean and v_mean.

[0029] S56. Calculate the pulsating speed: subtract the average speed from the instantaneous speed to obtain the pulsating speeds u_fluc and v_fluc.

[0030] S57, output optimization progress: output processing progress information;

[0031] S58. Visualize flow field results: Use the figure function to create a new graphics window, and use the quiver function to draw the velocity vectors u and v at the flow field position (x, y). The vectors are represented in red, and the contents of the current graphics window are saved as an image file in JPEG format.

[0032] The present invention is further configured as follows: in step S52, the outlier removal is set to be a speed value exceeding 3 times the standard deviation.

[0033] The present invention is further configured as follows: in step S7, the specific process is as follows:

[0034] S71. Image scaling: Use the imresize function to reduce the original image by half and store them separately;

[0035] S72. Feature point detection: Use the detectSURFFeatures function to detect SURF features in the scaled image;

[0036] S73, feature description: use the extractFeatures function to calculate the description vector of the detected feature points;

[0037] S74, feature point matching: match feature points between images using the matchFeatures function;

[0038] S75, feature point position adjustment and visualization: Use the Location attribute to extract the position coordinates of the matching feature points in the two images and store them in Lct1 and Lct2 respectively. For each pair of matching points, calculate their differences in x and y coordinates, square them and add them, and then take the square root to get the Euclidean distance between the matching point pairs;

[0039] S76. Draw matching points: Use imshow and plot functions to draw matching points on the Fimg1 image, loop through the plot function to connect each pair of matching points, add a legend, and set the background of the graphics window to white;

[0040] S77. Velocity field interpolation: Use the griddata function to interpolate the velocity of the entire flow field based on the coordinates of the matching point pairs and the displacement between them, and then use the medfilt2 function to perform median filtering on the interpolated velocity components to smooth the image.

[0041] The present invention is further configured as follows: after step S77, it also includes visualizing the flow field: creating a new graphics window and displaying the background image img1, using the pcolor function to color the flow field according to the size of the velocity vector uvm_fmv, and removing the boundary lines between the color blocks.

[0042] Compared with the prior art, the present invention has the following beneficial effects:

[0043] (1) The present invention adopts advanced image preprocessing technology: It adopts advanced image preprocessing technology such as median filtering and histogram equalization to enhance the surface texture of the image and improve the accuracy of subsequent PIV flow field calculation, which is more effective than the simple filter used in the prior art.

[0044] (2) The present invention provides an accurate flow field vector calculation method: a customized piv_FFTmulti function is used to calculate the flow field vector, and sub-pixel level correlation analysis is combined to provide a higher-precision flow velocity vector diagram, which has significantly improved the spatial resolution and velocity accuracy compared to traditional PIV technology.

[0045] (3) The present invention has real-time and dynamic monitoring capabilities: it can process and analyze flow field data in real time, has good adaptability to dynamically changing flow states, and improves the delays and static analysis limitations that may exist in the prior art.

[0046] (4) The present invention optimizes outlier processing: by setting a threshold and using the inpaint_nans function, outliers in flow field data are effectively identified and processed, data cleaning is optimized, and the reliability of flow field analysis is improved.

[0047] (5) The present invention provides a comprehensive feature point matching technology: it integrates the SURF feature point matching technology, improves the robustness of flow state recognition, especially in complex and dynamically changing flow fields, and is more advanced than the single feature point matching method that may be used in the prior art.

[0048] (6) The present invention has user-friendly visual output: it provides intuitive flow field vector diagrams and velocity distribution diagrams, so that non-professional users can also easily understand the flow analysis results.

[0049] (7) The present invention has excellent environmental adaptability and broad application prospects: The present invention is specially optimized for the complex flow state of the Three Gorges variable backwater area, has strong environmental adaptability, and can be extended to other similar hydrological environments. DETAILED DESCRIPTION

[0050] The following will be combined with the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0051] Glossary: ​​PIV (Particle Image Velocimetry): By scattering tracer particles in the flow field and using a pulsed laser to illuminate the cross-sectional area of ​​the flow field, the particle image is recorded by an imaging system, and then image processing techniques such as autocorrelation or cross-correlation are used to analyze the image to obtain the flow field velocity distribution.

[0052] LSPIV (Large Scale Particle Image Velocimetry): Improves PIV technology and applies it to larger-scale flow velocity measurements under field conditions.

[0053] SURF (Speeded Up Robust Features): is an interest point detection and descriptor algorithm for image matching. It is based on the scale space representation of the image and can detect stable feature points at different scales and directions.

[0054] Example:

[0055] This embodiment provides a reflow surface velocity field algorithm, and an implementation process thereof is as follows:

[0056] 1. Initialize and set the path:

[0057] Clear the MATLAB environment variables, turn off warning messages, get the current working directory and build the home path.

[0058] 2. Image reading:

[0059] Use the imread function to read two images from the specified folder path, denoted as img1 and img2.

[0060] 3. Image preprocessing (surface texture enhancement calculation): The surface texture of the image is significantly enhanced through grayscale conversion, median filtering, image subtraction and histogram equalization. This helps to improve the accuracy of subsequent PIV flow field calculations because the flow features in the image become more obvious and easier to identify.

[0061] Among them, 31) Grayscale conversion: Use rgb2gray function to convert the read color image into a grayscale image and store them in variables img_gray1 and img_gray2 respectively.

[0062] 32) Median filtering: Apply the median filter medfilt2 with a filter size of [200, 200] to the grayscale images img_gray1 and img_gray2 to obtain the filtered image I1.

[0063] 33) Image Subtraction: Use the imsubtract function to subtract the filtered image I1 from the original grayscale image img_gray1 to get the difference image I2, which highlights the difference between the two images.

[0064] 34) Histogram equalization: Apply histogram equalization histeq to the difference image I2 to enhance the contrast of the image and obtain the final texture enhanced images img_piv1 and img_piv2. Output the message "The image has completed the texture enhancement processing procedure and the output calculation picture effect".

[0065] 4. PIV flow field calculation: Using the FFT-based PIV algorithm, the velocity component of the flow field can be calculated quickly and accurately. This method is more efficient than the traditional correlation algorithm, especially when processing large-scale data, which improves the accuracy of PIV flow field calculation.

[0066] 41) Initialize parameters: Output the prompt message "Image is being processed for PIV flow field calculation..." and create a cell array s to store the parameter settings for PIV analysis.

[0067] Int.area 1, Step size 1, Subpix.finder, etc. are the parameter names of the PIV algorithm, which correspond to:

[0068] Int.area 1: The window size of the first layer correlation operation, set to 256.

[0069] Step size 1: The step size of the first layer correlation operation, set to 128.

[0070] Subpix.finder: Whether to perform sub-pixel level correlation calculations, set to 1 (yes).

[0071] Mask and ROI: used to set the mask and region of interest (ROI) respectively. They are empty here, indicating that they are not used.

[0072] Nr.ofpasses: The number of iterations of the correlation operation, set to 3.

[0073] Int.area 2 and Int.area 3: The window sizes for subsequent iterations, set to 128 and 64 respectively.

[0074] Window deformation: Window deformation option, set to '*linear', indicating the use of linear interpolation.

[0075] Repeated Correlation: Whether to repeat the correlation calculation, set to 1 (yes).

[0076] DisableAutocorrelation: Whether to disable autocorrelation calculation, set to 1 (yes).

[0077] 42) Perform PIV analysis: pass in two preprocessed images img_piv1 and img_piv2, as well as the parameter s set above, and call the piv_FFTmulti function based on the fast Fourier transform (FFT) cross-correlation algorithm to perform PIV calculations. The function returns the x and y coordinates of the flow field, as well as the corresponding velocity components u and v.

[0078] 43) Use the inpaint_nans function to fill the missing NaN values ​​in the velocity components u and v using 4th nearest neighbor interpolation.

[0079] 44) Output result: The printed information “Image PIV flow field calculation is completed, output flow field vector diagram” indicates that the PIV flow field calculation is completed.

[0080] 5. Process the analysis results of PIV flow field calculation: In this operation, outliers exceeding 3 times the standard deviation are eliminated and NaN values ​​are filled, which helps to remove inaccurate data caused by noise or erroneous tracking, ensures the reliability of flow field data, and thus improves the quality of analysis results; at the same time, a median filter is used in this operation to smooth the velocity vector field, reduce noise, and make the velocity field data more stable and continuous, which is crucial for subsequent flow field analysis and interpretation.

[0081] 51) Calculate the standard deviation: Use the std function to calculate the standard deviation of the velocity vector fields u and v, respectively, and store them in the variables ustd and vstd, which are used to measure the statistics of the fluctuation size of the velocity vector field.

[0082] 52) Eliminate outliers: u(abs(u)>3*ustd)=nan and v(abs(v)>3*vstd)=nan Set velocity values ​​exceeding 3 times the standard deviation to NaN (Not a Number), which are considered to be outliers caused by noise or erroneous particle tracking.

[0083] 53) Fill NaN values: u = inpaint_nans(u,4) and v = inpaint_nans(v,4) use interpolation to fill NaN values. The number 4 here means using 4 nearest neighbor interpolation.

[0084] 54) Smoothing: Use the median filter medfilt2 with a filter size of [3,3] to smooth the velocity vector field to reduce noise and obtain the instantaneous velocities u and v.

[0085] 55) Calculate the average speed: Use the median filter medfilt2 with a filter size of [10,10] again to calculate the average speeds u_mean and v_mean.

[0086] 56) Calculate the pulsating speed: Subtract the average speed from the instantaneous speed to obtain the pulsating speed u_fluc and v_fluc.

[0087] 57) Output optimization progress: Output processing progress information "Image PIV flow field result optimization, processing..." and "Image PIV flow field result optimization completed, output effect diagram.

[0088] 58) Visualize the flow field results: Use the figure function to create a new graphics window, and use the quiver function to draw the velocity vectors u and v at the flow field position (x, y). The vectors are represented in red. print(gcf,[mainPath,'\result output\img_piv1_flow velocity vector.jpg'],'-djpeg','-r600') saves the contents of the current graphics window as an image file in JPEG format.

[0089] 6. Visualize the flow field results: The velocity vector and velocity magnitude of the flow field are intuitively displayed through the quiver and pcolor functions, making complex flow field data easier to understand and analyze, which is conducive to presenting and communicating research results.

[0090] 61) Calculate the magnitude of the velocity vector: Use uvm = sqrt(u.^2 + v.^2), that is, the formula Calculate the magnitude of the velocity vector for each pixel.

[0091] 62) Visualize the flow field: figure; image(img1); hold on to create a new graphics window and display the background image img1. Use the pcolor function to color the flow field according to the size of the velocity vector uvm and store it in the variable h1. Use the EdgeColor property to set h1.EdgeColor = 'none' to set the edge color of h1 to none, thereby removing the boundary lines between the color blocks and making the graphics look smooth and clear. Use the colormap function to set the color map to a gradient color from blue to red 'jet' (jet is a variable commonly used to represent speed or other continuous changes). Use the caxis function to set the range of the color axis so that the color map covers a range from 0 to 1.5 times the median value of uvm. h = colormap bar displays a color bar to explain the velocity magnitude represented by different colors; s00 = streamslice (x, y, u_ROI, v_ROI, 10, 'linear', 'arrows') draws streamlines on the flow field to display the direction and magnitude of the velocity vector, where the 'linear' and 'arrows' parameters specify the style of the streamlines; set (s00, 'Color', 'k', 'LineWidth', 0.25) sets the color of the streamlines to black and the line width to 0.25; set (gcf, 'color', 'w') sets the background color of the graphics window to white; and gives the graphics window a text title of "Backflow velocity magnitude and flow velocity vector diagram".

[0092] 63) Save visualization results: Save the contents of the current graphics window as an image file in JPEG format.

[0093] 64) Pause and close the graphics window: Use pause(15) to pause the execution for 15 seconds so that the user can view the results, and use close('all') to close all graphics windows.

[0094] 7. SURF (Speeded Up Robust Features) feature point detection and matching: Tracking the correspondence between images through SURF feature point matching is conducive to understanding the geometric transformation and movement between images, while enhancing the accuracy of feature point matching.

[0095] 71) Image scaling: To reduce the amount of computation and simulate the image features at different viewpoints or scales to a certain extent, the imresize function is used to shrink the original images img_piv1 and img_piv2 by half and store them in variables Fimg1 and Fimg2 respectively.

[0096] 72) Feature point detection: Use the detectSURFFeatures function to detect SURF features in the two scaled images, and return a point object containing the location and scale of the feature point, which is stored in the variables points1 and points2 respectively. NumOctaves, NumScaleLevels and MetricThreshold are the parameters for detecting feature points, where:

[0097] 'NumOctaves',4: specifies the number of octaves used in the difference Gaussian scale space. An octave is a level in the image scale space. Increasing the number of octaves can make feature detection more robust, but it will increase the amount of computation.

[0098] 'NumScaleLevels',6: specifies the number of scale levels in each octave. More scale levels can improve the sensitivity of feature detection, but also increase the amount of computation.

[0099] 'MetricThreshold',1000: Sets the threshold for filtering feature points. This threshold is based on the determinant value of the Hessian matrix of the feature point. A higher threshold means that only very significant feature points will be detected.

[0100] 73) Feature description: Use the extractFeatures function to calculate the description vector of the detected feature points.

[0101] 74) Feature point matching: matchFeatures function is used to match feature points between two images, and the matched point pairs are stored in variables matched_pts1 and matched_pts2. Location attribute is used to extract the position coordinates of the matching feature points in the two images and store them in Lct1 and Lct2 respectively. For each pair of matching points, the difference between their x and y coordinates is calculated, squared and added, and then the square root is taken to get the Euclidean distance between the matching point pairs. Find function is used to find matching point pairs with a distance greater than 30 pixels and store them in [df,~]. Matching point pairs with a distance greater than 30 pixels are deleted to reduce the impact of false matches. Then a new graphics window is created and showMatchedFeatures function is used to display the filtered matching feature points for intuitive inspection of the matching results. Legend function is used to add a legend to the graphics window to distinguish the matching points of the two images.

[0102] 75) Feature point position adjustment and visualization: The location coordinates of the matching feature points in the two images are extracted using the Location attribute and stored in Lct1 and Lct2 respectively. For each pair of matching points, the difference in their x and y coordinates is calculated, squared and added, and then the square root is taken to get the Euclidean distance between the matching point pairs. If the image is reduced by half by imresize before processing, the following three lines of code will enlarge the coordinates and distances back to the original size:

[0103] Lct1 = Lct1 * 2;

[0104] Lct2 = Lct2*2;

[0105] d = d * 2;

[0106] 76) Draw matching points: Use imshow and plot functions to draw matching points on the Fimg1 image. Red circles represent Lct1 and blue crosses represent Lct2. Loop through the plot function to connect each pair of matching points with black lines. Add a legend and set the background of the graphics window to white.

[0107] 77) Velocity field interpolation: Use the griddata function to interpolate the velocity components u_fmv and v_fmv of the entire flow field according to the coordinates of the matching point pairs and the displacement between them, and then use the medfilt2 function to perform median filtering on the interpolated velocity components, smooth the image, and estimate the continuous velocity field from discrete matching points, thereby improving the resolution and coverage of the velocity field data. Use the inpolygon function to define a polygonal region (mask) defined by the coordinates specified by myMask{1,1} and myMask{1,2}, which is used to identify points in a specific area of ​​the flow field, calculate the magnitude of the velocity vector, and set the magnitude of the velocity vector in the mask area to NaN; similarly, u_fmv_ROI and v_fmv_ROI are the x and y components of the velocity vector, respectively, and they are also set to NaN in the mask area.

[0108] 78) Visualize the flow field: Create a new graphics window and display the background image img1. Use the pcolor function to color the flow field according to the magnitude of the velocity vector uvm_fmv and remove the boundary lines between the color blocks. Set the color map to 'jet'. Set the range of the color axis so that the color map covers a range from 0 to 1.5 times the median of uvm_fmv. Display a color bar to explain the velocity magnitude represented by different colors. Use the streamslice function to draw streamlines on the flow field to show the direction and magnitude of the velocity vector. Set the streamline color to black and the line width to 0.25. Set the background color of the graphics window to white.

[0109] In the description of this specification, the description with reference to the terms "one embodiment", "example", "specific example", etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representation of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in any one or more embodiments or examples in a suitable manner.

[0110] The preferred embodiments of the present invention disclosed above are only used to help illustrate the present invention. The preferred embodiments do not describe all the details in detail, nor do they limit the invention to the specific implementation methods described. Obviously, many modifications and changes can be made according to the content of this specification. This specification selects and specifically describes these embodiments in order to better explain the principles and practical applications of the present invention, so that those skilled in the art can understand and use the present invention well. The present invention is limited only by the claims and their full scope and equivalents.

Claims

1. The reflow surface velocity field algorithm is characterized by: The following steps are involved: S1. Initialize and set the path: clear MATLAB environment variables, close warning messages, get the current working directory and build the main path; S2, image reading: use the imread function to read the image from the specified folder path; S3, image preprocessing: enhance the read image through grayscale conversion, median filtering, image subtraction and histogram equalization; S4, PIV flow field calculation: PIV algorithm based on FFT is used to calculate the velocity component of the flow field; S5, analysis of processing PIV flow field calculation: remove outliers exceeding 3 times the standard deviation and fill in NaN values. At the same time, use a median filter to smooth the velocity vector field; S6. Visualize the flow field results: Use the quiver and pcolor functions to display the velocity vector and velocity magnitude of the flow field; S7, SURF feature point detection and matching: Track the correspondence between images through SURF feature point matching.

2. The reflow surface velocity field algorithm according to claim 1, characterized in that: In step S4, the specific process is as follows: S41, initialization parameters: output prompt information, create a cell array s for storing parameter settings for PIV analysis; S42, perform PIV analysis: input the pre-processed image and the set parameters, and perform PIV calculation. The function calculates and returns the x and y corresponding to the x-coordinate and y-coordinate of the flow field, and the corresponding velocity components u and v; S43, filling the missing values ​​NaN in the velocity components u and v; S44, output results: print information, and the PIV flow field calculation is completed.

3. The reflow surface velocity field algorithm according to claim 2, characterized in that: In the step S42, the piv_FFTmulti function based on the Fast Fourier Transform FFT cross-correlation algorithm is called to perform PIV calculation.

4. The backflow surface velocity field algorithm according to claim 2, characterized in that: In step S43, the missing values ​​NaN in the velocity components u and v are filled using the inpaint_nans function using 4th order nearest neighbor interpolation.

5. The reflow surface velocity field algorithm according to claim 1, characterized in that: In step S5, the specific process is as follows: S51, calculate standard deviation: use the std function to calculate the standard deviation of the velocity vector fields u and v respectively, and store them respectively; S52, remove outliers; S53, fill NaN value: fill NaN value by interpolation method; S54, smoothing: using the median filter medfilt2 to smooth the velocity vector field; S55, calculate the average speed: use the median filter medfilt2 again to calculate the average speeds u_mean and v_mean. S56. Calculate the pulsating speed: subtract the average speed from the instantaneous speed to obtain the pulsating speeds u_fluc and v_fluc. S57, output optimization progress: output processing progress information; S58. Visualize flow field results: Use the figure function to create a new graphics window, and use the quiver function to draw the velocity vectors u and v at the flow field position (x, y). The vectors are represented in red, and the contents of the current graphics window are saved as an image file in JPEG format.

6. The backflow surface velocity field algorithm according to claim 5, characterized in that: In the step S52, the outlier removal is set to the speed value exceeding 3 times the standard deviation.

7. The reflow surface velocity field algorithm according to claim 1, characterized in that: In step S7, the specific process is as follows: S71. Image scaling: Use the imresize function to reduce the original image by half and store them separately; S72. Feature point detection: Use the detectSURFFeatures function to detect SURF features in the scaled image; S73, feature description: use the extractFeatures function to calculate the description vector of the detected feature points; S74, feature point matching: match feature points between images using the matchFeatures function; S75, feature point position adjustment and visualization: Use the Location attribute to extract the position coordinates of the matching feature points in the two images and store them in Lct1 and Lct2 respectively. For each pair of matching points, calculate their differences in x and y coordinates, square them and add them, and then take the square root to get the Euclidean distance between the matching point pairs; S76. Draw matching points: Use imshow and plot functions to draw matching points on the Fimg1 image, loop through the plot function to connect each pair of matching points, add a legend, and set the background of the graphics window to white; S77. Velocity field interpolation: Use the griddata function to interpolate the velocity of the entire flow field based on the coordinates of the matching point pairs and the displacement between them, and then use the medfilt2 function to perform median filtering on the interpolated velocity components to smooth the image.

8. The backflow surface velocity field algorithm according to claim 7, characterized in that: After step S77, the process also includes visualizing the flow field: creating a new graphics window and displaying the background image img1, using the pcolor function to color the flow field according to the size of the velocity vector uvm_fmv, and removing the boundary lines between the color blocks.