A method for reconstructing three-dimensional underwater flow fields under complex lighting conditions

Abstract

This invention discloses a method and apparatus for reconstructing three-dimensional underwater flow fields under complex lighting conditions. The method includes: calculating the initial flow velocity; acquiring a sequence of binocular polarized images, solving the total intensity image, performing multi-scale filtering and adaptive background subtraction, and outputting a particle-enhanced image; combining local dynamic thresholding and gradient extremum search to perform particle identification and overlap decoupling to extract sub-pixel coordinates; correcting polarization distortion using calibration parameters, and establishing corresponding points based on epipolar lines and polarization features, then correcting the particle's three-dimensional spatial coordinates through layered medium refraction; constructing a polarization-space fusion cost function to perform inter-frame particle correlation to solve for the three-dimensional velocity vector; and finally, performing post-processing that incorporates hydrodynamic physical constraints such as geometric, extremum, and local consistency constraints based on the initial flow velocity. This method effectively eliminates interference from complex underwater lighting fluctuations and improves the accuracy of flow field reconstruction.

Description

Technical Field

[0001] This application relates to the fields of artificial intelligence, computer vision, photoelectric detection and fluid measurement technology, and specifically to a method and apparatus for reconstructing three-dimensional underwater flow fields under complex lighting conditions. Background Technology

[0002] Accurate reconstruction of three-dimensional underwater flow fields is a core foundation for underwater fluid mechanics mechanism research, precision fluid device development, underwater equipment performance testing, and operating condition optimization. It directly determines the accuracy of flow field characteristic analysis, the reliability of fluid structure design, and the validity of underwater test data. Non-contact optical measurement technology, with its advantages of no flow field disturbance, large measurement range, and fast dynamic response, has now replaced traditional contact measurement methods and become the mainstream technology for underwater flow field parameter characterization and three-dimensional flow field reconstruction, providing core technical support for refined flow field research under complex underwater operating conditions.

[0003] Current underwater optical flow field measurements mainly include two core technologies: Particle Image Velocimetry (PIV) and Particle Tracking Velocimetry (PTV). PIV, based on a grid averaging algorithm, solves for global flow field velocities, quickly acquiring macroscopic distribution characteristics of the flow field and is suitable for overall measurement of large-scale, strongly disturbed flow fields. However, this technology is limited by the grid averaging mechanism, unable to identify microscale flow details, has low spatial resolution, and requires stringent requirements for the distribution density and uniformity of tracer particles, making it difficult to adapt to actual underwater conditions with low particle concentrations, resulting in significant application limitations. In contrast, PTV technology inverts flow field parameters by tracking the spatiotemporal trajectory of individual tracer particles, avoiding the shortcomings of PIV technology such as grid averaging ambiguity, loss of microscopic features, and poor particle concentration adaptability. It has advantages such as high measurement accuracy, adaptability to low particle concentration conditions, and the ability to resolve microscopic flow characteristics, making it more suitable for underwater three-dimensional refined flow field reconstruction.

[0004] The real underwater environment is characterized by non-uniform and time-varying illumination. Illumination intensity and distribution dynamically change with water depth, suspended solids concentration, and current disturbances, creating complex, non-equilibrium illumination conditions that severely impact the particle imaging quality and 3D flow field reconstruction accuracy of PTV technology. Complex dynamic illumination leads to significant spatiotemporal differences in the brightness and contrast of tracer particle images, resulting in noticeable deviations in the imaging characteristics of the same particle at different spatial locations and measurement times. Existing PTV particle recognition algorithms often use fixed light intensity thresholds for particle segmentation and extraction, exhibiting poor adaptability to different lighting environments: low-illumination areas are prone to missing real particles due to insufficient particle grayscale features, while high-illumination areas are susceptible to false particle identification and contour misjudgment due to overexposure and edge distortion. This prevents sub-pixel-level high-precision particle localization, fundamentally reducing the quality of the original flow field reconstruction data and leading to problems such as broken particle trajectories, matching errors, and flow velocity calculation deviations. Currently, there are no effective solutions in the field to address underwater dynamic illumination interference and adapt to complex illumination measurement scenarios, hindering the high-precision engineering application of 3D underwater PTV flow field reconstruction technology. Summary of the Invention

[0005] The main objective of this application is to provide a method and apparatus for three-dimensional underwater flow field reconstruction under complex lighting conditions. This application aims to solve the technical problems in existing underwater PTV flow field reconstruction technology, such as weak illumination robustness due to non-uniform underwater dynamic illumination and time-varying fluctuations, low particle identification and positioning accuracy, difficulty in effectively segmenting overlapping particles, and poor three-dimensional reconstruction accuracy and temporal trajectory correlation errors caused by multi-medium refraction and polarization distortion.

[0006] To achieve the above objectives, the first aspect of this application provides a method for reconstructing a three-dimensional underwater flow field under complex lighting conditions, comprising: An underwater measurement system including a binocular polarization camera and a laser displacement sensor was constructed. The real physical distance provided by the laser displacement sensor was used as a constraint for joint calibration to obtain calibration parameters. A single tracer particle was released in the flow field. Based on the calibration parameters and the particle motion image sequence, the inter-frame actual physical displacement of the tracer particle was calculated to obtain the initial flow velocity of the flow field. A sequence of binocular polarization images of underwater tracer particles is acquired, the intensity image is calculated, and multi-scale filtering and adaptive background subtraction are performed to separate the tracer particles from the backscattered background and obtain the particle-enhanced image. Using the particle-enhanced image as input, particle recognition and overlap decoupling are performed by combining local dynamic thresholding and gradient extremum search algorithms to extract sub-pixel-level center coordinates; the polarization distortion of the sensor system is corrected using the calibration parameters to obtain the corrected polarization features; epipolar constraints are constructed based on the calibration parameters, and corresponding points are established by combining the sub-pixel-level center coordinates, grayscale features, and the corrected polarization features; the three-dimensional spatial coordinates of the particles are solved using triangulation and layered medium refraction correction. Based on the three-dimensional spatial coordinates of temporally adjacent frames and the corrected polarization features, a polarization-space fusion matching cost function is constructed to perform inter-frame particle correlation and solve for the three-dimensional velocity vector of the particles. The three-dimensional velocity vector is combined with the initial flow velocity of the flow field, and three-dimensional flow field post-processing based on hydrodynamic physical constraints is performed to obtain the final three-dimensional underwater flow field reconstruction result.

[0007] Optionally, the actual physical object distance provided by the laser displacement sensor is used as a constraint for joint calibration to obtain calibration parameters, including: Set the extrinsic translation vector of the camera under different spatial attitudes, extract the axial calculated distance from the camera's optical center to the calibration plate plane, and obtain the actual axial physical distance measured by the laser displacement sensor; A joint optimization objective function with distance constraints is constructed, wherein the objective function is a weighted sum of a reprojection error term and a physical scale penalty term; wherein the reprojection error term is the sum of squared deviations between the actual extracted corner coordinates and the coordinates calculated by the nonlinear projection model, and the physical scale penalty term is the sum of squared deviations between the calculated axial distance and the actual axial physical distance multiplied by a weight penalty coefficient; The joint optimization objective function is solved using a nonlinear iterative algorithm to obtain the high-precision camera intrinsic parameters, extrinsic parameters, and distortion coefficients, which serve as the calibration parameters.

[0008] Optionally, the inter-frame actual physical displacement of the tracer particles is calculated based on the calibration parameters and the particle motion image sequence to obtain the initial flow velocity of the flow field, including: To locate the local region of interest where a single tracer particle is located, assuming that the light spot intensity of the tracer particle follows a two-dimensional symmetric Gaussian distribution, a mathematical fitting model is constructed that includes the peak amplitude of light intensity, the coordinates of the centroid to be solved, the standard deviation, and the background noise. The mathematical fitting model is linearized by taking the logarithm of both sides, and then converted into a quadratic polynomial form with respect to spatial coordinates. The coordinates and grayscale data of all pixels in the local region of interest are solved by the least squares method, and the sub-pixel level centroid coordinates of the tracer particles are obtained by inverse solution of the polynomial fitting coefficients. The average pixel displacement of particles in multiple adjacent frames is statistically analyzed, and combined with the pixel-to-physical scale mapping relationship in the calibration parameters, it is converted into the actual physical displacement between frames. The initial flow velocity of the flow field is then calculated by combining the video frame time interval.

[0009] Optionally, the calculation of the intensity image and the application of multi-scale filtering and adaptive background subtraction to obtain the particle-enhanced image include: The original polarized mosaic image is obtained by using a split-plane polarization camera. The intensity images of the four polarization channels are obtained by using gradient-guided interpolation. The total intensity image is obtained by summing the intensity images of the four polarization channels and taking half of the sum. A Gaussian filter bank with multiple scales is constructed, and corresponding weights are assigned to Gaussian kernels of different scales. Multi-scale weighted convolution response calculation is performed on the total intensity image to suppress high-frequency detection noise while compensating for the difference in particle light intensity distribution caused by depth changes, thus obtaining a denoised image. A local feature sliding window larger than a preset multiple of the particle's maximum diameter is constructed. The local grayscale expectation of the pixels in the denoised image within the window is calculated as the background estimate. The background is subtracted by the difference between the denoised image and the local grayscale expectation. Extract the minimum and peak gray values ​​of the non-zero pixels in the image after background subtraction, construct a linear mapping function to extend the particle gray-scale response to full scale, expand the brightness gradient between the particle center and edge, and output the particle-enhanced image.

[0010] Optionally, the step of combining local dynamic thresholding and gradient extremum search algorithms for particle recognition and overlap decoupling, and extracting sub-pixel-level center coordinates, includes: A local sliding window is established based on the particle-enhanced image. The mean and variance of pixel gray levels are calculated for each window. A pixel-level dynamic segmentation threshold is constructed by combining the water noise adaptive adjustment coefficient. The threshold is equal to the local gray level mean plus the product of the adjustment coefficient and the local gray level variance. Based on the threshold, pixel-by-pixel binarization is performed to generate particle candidate regions. For the particle candidate region, the first-order gray-level partial derivatives in the horizontal and vertical directions of the gray-level image are solved, and the gradient magnitude and gradient direction angle of each pixel are calculated using the differential operator to construct a gray-level gradient vector field covering the entire connected domain. Local grayscale peak traversal and discrimination are performed within each candidate connected region. If there are two or more local peaks within the connected region and they meet the physical constraints, the gradient direction angle sequence is extracted pixel by pixel along the reference path connecting the two peaks, and the gradient direction change rate of adjacent pixels is calculated. If the change rate exceeds the preset gradient turning threshold, the pixel is determined to be a grayscale gradient inflection point. The line connecting the consecutive inflection points is used as the segmentation boundary to complete the decoupling segmentation of overlapping particles, and two-dimensional Gaussian fitting is performed on the decoupled single-particle region to extract the sub-pixel-level center coordinates.

[0011] Optionally, the polarization state distortion of the sensor system is corrected using the calibration parameters to obtain the corrected polarization characteristics, including: A polarization refraction physical model based on Stokes vectors and Mueller matrices is established. The Stokes vectors of the incident and outgoing light are physically correlated through a calibrated four-by-four Mueller matrix. The diagonal elements of the Mueller matrix characterize the degree of change of the linear polarization component of the refracted light and satisfy the properties of total light intensity conservation and circular polarization state invariance. For the water-glass interface and the glass-air interface in the light propagation path, the corresponding Mueller matrices are determined respectively, and a cascaded positive transformation relationship from the incident polarization state to the sensor observation polarization state is constructed. That is, the observation polarization state is equal to the product of the glass-air interface Mueller matrix and the water-glass interface Mueller matrix with respect to the initial polarization state. The distorted polarization data acquired by the sensor is inversely compensated by using the inverse transformation of the cascaded matrix. That is, the observed polarization state is multiplied by the inverse of the Mueller matrix of each interface in turn to obtain the true Stokes vector after compensation. The linear polarization degree and polarization angle are then extracted as the corrected polarization features.

[0012] Optionally, the step of establishing corresponding points by combining the sub-pixel-level center coordinates, grayscale features, and the corrected polarization features, and solving for the three-dimensional spatial coordinates of the particles using triangulation and layered medium refraction correction, includes: Based on the calibration parameters, a basic matrix is ​​constructed, and the normalized distance from the candidate particle point to the epipolar line in the right image is calculated. When the distance is less than the epipolar line search bandwidth threshold, it is used as a candidate matching point. The horizontal and vertical search intervals are calculated by combining the effective range of disparity to narrow the search space. Calculate the gray-level normalized cross-correlation coefficients of the left and right local windows; calculate the Stokes parameters from the intensity images of the four polarization directions, convert the polarization angles into periodic invariant features to construct polarization feature vectors, and calculate the polarization feature differences between the left and right particles; construct a matching cost function by weighted summation of the gray-level cross-correlation coefficients and polarization feature differences, select the particle pair with the minimum cost, and perform bidirectional consistency judgment to establish corresponding points; After solving the initial three-dimensional coordinates based on the principle of triangulation, a layered medium refraction correction model is established for the refraction of the water tank wall and the water body. The direction of the incident direction after refraction through the interface is calculated according to the vector form of the law of refraction. The true back projection ray in the water body is solved, and the midpoint of the line connecting the shortest distance between the left and right refracted rays is found to obtain the three-dimensional spatial coordinates after refraction correction.

[0013] Optionally, based on the three-dimensional spatial coordinates of temporally adjacent frames and the corrected polarization features, a polarization-spatial fusion matching cost function is constructed to perform inter-frame particle correlation, and the three-dimensional velocity vector of the particles is solved, including: For effective particles and candidate particles in two adjacent frames, a three-dimensional spatial distance cost term and an inter-frame polarization degree feature difference cost term are defined respectively. Dimensionless weight coefficients are introduced, and the two types of cost terms are normalized by the maximum spatial distance and the maximum polarization degree difference in the candidate matching pair set to eliminate the dimension difference. Construct a particle matching cost function that integrates polarization and space, which is equal to the weighted sum of the normalized spatial distance cost and the normalized polarization degree difference cost; Traverse all cross-frame candidate particle matching combinations, select the particle pair corresponding to the minimum fusion matching cost as the optimal matching result to complete the temporal association, and calculate the temporal three-dimensional displacement and instantaneous three-dimensional velocity vector of the particles in combination with the inter-frame time interval of the camera.

[0014] Optionally, the final three-dimensional underwater flow field reconstruction result, obtained by combining the three-dimensional velocity vector with the initial flow velocity and performing three-dimensional flow field post-processing based on hydrodynamic physical constraints, includes: Establish geometric effective domain constraints based on finite space: Define the effective measurement area of ​​the water tank and the set of wall planes, and eliminate velocity vectors that exceed the boundary or violate the physical rule that the wall is impenetrable; Establish physical extreme value constraints for velocity: Define the reference velocity of the flow field based on the initial flow velocity and the average inter-frame physical displacement, construct the physical upper limit of velocity in combination with the three-dimensional positioning error, and eliminate abnormal vectors whose velocity magnitude exceeds the reasonable range of physics; Establish local consistency constraints for velocity length: Use a robust statistical method based on the neighborhood median and the median absolute deviation to calculate the normalized ratio of the target vector velocity length residual to the neighborhood length residual scale, and identify outlier vectors in length. Establish a local consistency constraint for velocity direction based on the median value of the unit spherical direction: Define the angular distance between two unit direction vectors, select the direction vector with the smallest angular distance from the neighborhood direction set as the local principal direction, calculate the normalized ratio of the angular residual between the target velocity direction and the local principal direction to the fluctuation scale of the neighborhood direction, and determine the outlier vector of the direction.

[0015] A second aspect of this application provides a three-dimensional underwater flow field reconstruction device for complex lighting conditions, comprising: The initial flow velocity calculation module is used to build an underwater measurement system including a binocular polarization camera and a laser displacement sensor. It uses the real physical distance provided by the laser displacement sensor as a constraint to perform joint calibration and obtain calibration parameters. A single tracer particle is released in the flow field, and the inter-frame actual physical displacement of the tracer particle is calculated based on the calibration parameters and the particle motion image sequence to obtain the initial flow velocity of the flow field. The particle image enhancement module is used to acquire a sequence of binocular polarization images of underwater tracer particles, calculate the intensity image and perform multi-scale filtering and adaptive background subtraction, separate the tracer particles from the backscattered background, and obtain the particle-enhanced image. The three-dimensional spatial coordinate solving module is used to take the particle-enhanced image as input, combine the local dynamic threshold and gradient extremum search algorithm to perform particle recognition and overlap decoupling, and extract sub-pixel-level center coordinates; use the calibration parameters to correct the polarization distortion of the sensor system to obtain the corrected polarization features; construct epipolar constraints based on the calibration parameters, establish corresponding points by combining the sub-pixel-level center coordinates, grayscale features and the corrected polarization features, and solve the three-dimensional spatial coordinates of the particles using triangulation and layered medium refraction correction. The output module is used to construct a polarization-space fusion matching cost function based on the three-dimensional spatial coordinates of temporally adjacent frames and the corrected polarization features to perform inter-frame particle correlation and solve for the three-dimensional velocity vector of the particles; the three-dimensional velocity vector is combined with the initial flow velocity of the flow field to perform three-dimensional flow field post-processing based on hydrodynamic physical constraints to obtain the final three-dimensional underwater flow field reconstruction result.

[0016] This application provides a method and apparatus for reconstructing three-dimensional underwater flow fields under complex lighting conditions. The method includes constructing an underwater measurement system comprising a binocular polarization camera and a laser displacement sensor; performing joint calibration using the actual physical distance provided by the laser displacement sensor as a constraint to obtain calibration parameters; releasing a single tracer particle in the flow field; calculating the inter-frame actual physical displacement of the tracer particle based on the calibration parameters and the particle motion image sequence to obtain the initial flow velocity; acquiring the binocular polarization image sequence of the underwater tracer particle; calculating the intensity image and performing multi-scale filtering and adaptive background subtraction to separate the tracer particle from the backscattered background to obtain a particle-enhanced image; and using the particle-enhanced image as input, combining local dynamic thresholding and gradient extremum search. The algorithm performs particle recognition and overlap decoupling to extract sub-pixel-level center coordinates; it corrects the polarization distortion of the sensor system using calibration parameters to obtain the corrected polarization features; it constructs epipolar constraints based on the calibration parameters, and establishes corresponding points by combining sub-pixel-level center coordinates, grayscale features, and corrected polarization features; it then uses triangulation and layered medium refraction to correct and solve for the particle's three-dimensional spatial coordinates; based on the three-dimensional spatial coordinates of temporally adjacent frames and the corrected polarization features, it constructs a polarization-space fusion matching cost function to perform inter-frame particle association and solve for the particle's three-dimensional velocity vector; finally, it combines the three-dimensional velocity vector with the initial flow velocity of the flow field and performs three-dimensional flow field post-processing based on hydrodynamic physical constraints to obtain the final three-dimensional underwater flow field reconstruction result. Compared to existing technologies, this application, based on the polarization optical imaging mechanism, captures the inherent polarization characteristics of tracer particles, overcoming the inherent limitations of traditional PTV technology that relies on light intensity information for particle identification and positioning. It can effectively suppress complex lighting interference caused by non-uniform underwater illumination, time-varying fluctuations, and suspended impurities in the water. It solves problems such as missed particle detection, overexposure distortion, and false particle misjudgment in traditional fixed light intensity threshold algorithms, improves the accuracy of tracer particle identification and sub-pixel positioning stability under complex lighting conditions, and makes up for the technical shortcomings of traditional PTV technology in adapting to dynamic underwater lighting scenarios. It can provide high-precision and high-reliability flow field data support for underwater hydrodynamic mechanism research, precision fluid device development, and underwater equipment performance testing and optimization. Attached Figure Description

[0017] Figure 1 This is a schematic flowchart of an embodiment of the three-dimensional underwater flow field reconstruction method for complex lighting conditions provided in this application.

[0018] The realization of the purpose, functional features and advantages of this application will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are for illustrative purposes only and are not intended to limit the scope of this application. Where there is no conflict, the embodiments and features described herein can be combined with each other.

[0020] It should be noted that, in the description of the embodiments of this application, in order to make the expression more concise, avoid redundancy and highlight key improvements, the word "the" is not used in any embodiments. Instead, terms such as "the," "the above," "this," "this," and "corresponding" are used to refer to the aforementioned elements. This is intended to ensure that the specification text corresponding to the claims conforms to the requirements of rigor, standardization, and examination guidance style, while strictly meeting the requirements of specific drafting specifications.

[0021] Before detailing the specific implementation scheme of this application, we first provide a rigorous physical and technical definition of the core technical terms involved in this application. These terms have a stable relationship in subsequent steps: Calibration parameters: refer to the set of parameters including the intrinsic parameter matrix and extrinsic parameter matrix (rotation matrix and translation vector) of the binocular polarization camera, as well as the radial and tangential distortion coefficients of the lens. In this application, these parameters are high-precision mapping parameters obtained by joint optimization under the physical distance constraints provided by the laser displacement sensor.

[0022] Tracer particles: These are tiny spherical physical particles that are released into an underwater flow field and move synchronously with the fluid. These particles not only have a specific geometric scale, but also scatter light with stable polarization physical properties when illuminated.

[0023] Particle-enhanced images refer to images obtained by performing total intensity calculation, multi-scale noise reduction, background spatial degradation field estimation, and global contrast stretching on the original multi-channel polarized mosaic image. These images can significantly amplify the brightness gradient of particles and suppress backscattering noise.

[0024] Polarization characteristics: including the degree of linear polarization (DoLP) and the angle of polarization (AoLP). In this application, these characteristics are corrected in the spatial domain by the Mueller matrix cascade inverse transformation of the "water-glass-air" multi-medium refractive interface to characterize the true surface scattering polarization properties of the tracer particles.

[0025] Corresponding points: refer to pairs of corresponding pixels in the left and right imaging planes of a binocular polarization camera, generated by the projection of the same three-dimensional physical tracer particle in space.

[0026] Three-dimensional velocity vector: refers to the instantaneous velocity of the tracer particle in three-dimensional space, including motion components in the horizontal, vertical and depth directions. It is calculated by combining the optimal particle correlation results between time-series frames and the inter-frame sampling interval.

[0027] Physical constraints refer to the set of fluid dynamics and spatial geometry rules, consisting of the finite space boundary of the experimental container, the fluid kinematic limit, the local continuity topology of the flow velocity, and the incompressible continuity equation of the water body, used to filter out trajectory noise and erroneous associations.

[0028] Traditional underwater particle tracking and velocimetry (PTV) relies heavily on single-dimensional light intensity grayscale information for particle extraction when reconstructing the flow field. However, in actual underwater conditions, the presence of multi-scattering media causes significant attenuation and backscattering of light during propagation. The disturbances of water flow, bubbles, and uneven distribution of the light source create non-uniform and dynamically fluctuating background noise on the image plane (often referred to as the "veil effect"). This complex, non-equilibrium illumination results in significant spatial differences in the brightness and contrast of the tracer particles in the image.

[0029] Using traditional fixed light intensity threshold segmentation leads to the following problems: First, tracer particles in low-light and far-depth-of-field areas have extremely low signal-to-noise ratios and are easily submerged in background noise, resulting in a large number of missed detections of real particles. Second, particles in high-light and near-depth-of-field areas experience feature saturation due to overexposure, causing multiple adjacent particles to overlap and connect in the image. Traditional center localization methods cannot decouple these overlapping particles, resulting in severe localization errors. Third, the layered refraction of the tank wall and the water medium causes nonlinear geometric distortion in binocular imaging, rendering traditional epipolar constraints ineffective. These problems fundamentally limit the accuracy of matching and lead to numerous trajectory breaks and incorrect matches during time-series tracking.

[0030] To overcome the aforementioned technical deficiencies, this application introduces the physical characteristics of polarization imaging into the PTV technology system. Polarization information, as a physical dimension independent of light intensity, is naturally immune to external light intensity fluctuations, providing a robust physical constraint for underwater flow field reconstruction.

[0031] refer to Figure 1 The first embodiment of this application provides a method for reconstructing a three-dimensional underwater flow field under complex lighting conditions, thereby solving the technical problems mentioned in the background art. This method can be executed by a processor located on a terminal or server, and the execution process of this method can be as follows: Step S101: Set up the system and perform joint calibration to obtain the initial flow rate.

[0032] An underwater measurement system was constructed, comprising a binocular polarization camera and a laser displacement sensor (LDS). The actual physical distance provided by the LDS was used as a constraint for joint calibration to obtain calibration parameters. A single tracer particle was released into the flow field, and the inter-frame actual physical displacement of the tracer particle was calculated based on the calibration parameters and the particle motion image sequence to obtain the initial flow velocity.

[0033] In this embodiment, the process of obtaining calibration parameters by using the actual physical object distance provided by the laser displacement sensor as a constraint for joint calibration can be as follows: Define the extrinsic translation vector of the camera under different spatial attitudes, extract the axial calculated distance representing the distance from the camera's optical center to the calibration plate plane, and obtain the actual axial physical distance measured by the laser displacement sensor; A joint optimization objective function with distance constraints is constructed. The objective function is a weighted sum of the reprojection error term and the physical scale penalty term. The reprojection error term is the sum of squared deviations between the actual extracted corner coordinates and the coordinates calculated by the nonlinear projection model. The physical scale penalty term is the sum of squared deviations between the calculated axial distance and the actual axial physical distance multiplied by the weight penalty coefficient. The joint optimization objective function is solved using a nonlinear iterative algorithm to obtain high-precision camera intrinsic and extrinsic parameters, as well as distortion coefficients, which serve as calibration parameters.

[0034] In this embodiment, the process of calculating the inter-frame actual physical displacement of the tracer particles based on calibration parameters and particle motion image sequences to obtain the initial flow velocity of the flow field can be as follows: To locate the local region of interest where a single tracer particle is located, assuming that the light spot intensity of the tracer particle follows a two-dimensional symmetric Gaussian distribution, a mathematical fitting model is constructed that includes the peak amplitude of light intensity, the coordinates of the centroid to be solved, the standard deviation, and the background noise. Linearize the mathematical fitting model by taking the logarithm of both sides, and convert it into a quadratic polynomial form with respect to spatial coordinates; The least squares method is used to solve the overdetermined system of equations for the coordinates and grayscale data of all pixels in the local region of interest. The subpixel-level centroid coordinates of the tracer particles are obtained by inverse solving the polynomial fitting coefficients. The average pixel displacement of particles in multiple adjacent frames is statistically analyzed, and the pixel-to-physical scale mapping relationship in the calibration parameters is combined to convert it into the actual physical displacement between frames. The initial flow velocity of the flow field is then calculated by combining the inter-frame time interval of the video.

[0035] Specifically, in an optional embodiment of this application, the sub-process of using the actual physical object distance provided by the laser displacement sensor as a constraint for joint calibration to obtain the calibration parameters is as follows: Set the camera on the checkerboard calibration board. The extrinsic translation vectors under each spatial pose are: ,in Characterize the axial calculated distance from the camera's optical center to the calibration plate plane obtained from the solution; set the first The true axial physical distance obtained by LDS measurement under each attitude and outlier removal is .

[0036] The processor can construct a joint optimization objective function that includes LDS distance constraints: (1) Among them, the first term on the right side of formula (1) is the reprojection error term, which ensures the fitting accuracy of the two-dimensional pixel plane, and the second term on the right side is the LDS physical scale penalty term, which forces the solution depth of the three-dimensional space to converge to the absolute physical truth provided by LDS. The total number of valid calibration board images acquired; This represents the total number of valid corner points on a single calibration board image. For the first The first image The actual sub-pixel coordinates of each corner point are extracted; It is a nonlinear projection function that incorporates the radial and tangential distortion models of the camera lens; The intrinsic parameter matrix of the camera to be optimized; The first The camera rotation matrix and translation vector corresponding to each image; For the first The three-dimensional physical coordinates of each corner point in the calibration plate's world coordinate system; This is the weighting penalty coefficient for the distance constraint, and its value is determined by the measurement accuracy and confidence level of the LDS sensor.

[0037] The processor then uses the Levenberg-Marquardt algorithm on the objective function. Iterative solutions are performed to obtain highly accurate optimized intrinsic parameters. External reference Including the distortion coefficient, the sensor system is calibrated.

[0038] Furthermore, a single tracer particle is released in the flow field to extract the initial flow velocity, and the sub-process is as follows: In the image sequence split according to the shooting frame rate, the local region of interest (ROI) where a single tracer particle is located is identified, and the gray-level distribution of each pixel within the ROI is obtained. Assuming that the light spot intensity of the tracer particle follows a two-dimensional symmetric Gaussian distribution, the following mathematical fitting model is constructed: (2) In the formula, Image coordinates The pixel grayscale value at that location, This is the amplitude parameter of the peak particle light intensity. Let be the sub-pixel centroid coordinates of the particle to be solved. These are the standard deviations of the Gaussian distribution in the horizontal and vertical directions (characterizing the particle imaging diameter), respectively. This represents background ambient noise.

[0039] (3) In the formula, These are the polynomial fitting coefficients. The overdetermined equations for the coordinates and grayscale data of all pixels within the ROI region are solved using the least squares method. Finally, the sub-pixel-level centroid coordinates of the tracer particles are obtained through the inverse solution of the fitting coefficients. (4) The processor can sequentially extract the subpixel centroid coordinates of the individual tracer particle in each frame of the image sequence, calculate the average pixel displacement of the particle centroids between adjacent frames, and convert it into actual physical displacement between frames by combining the mapping scale factor (unit: mm / pixel) from the two-dimensional pixel plane to the three-dimensional world in the calibration parameters. Finally, it calculates the initial flow velocity of the underwater flow field by combining the time interval between adjacent frames of the video. This initial flow velocity is the global physical reference for the velocity-hydrodynamic constraints in this application.

[0040] Step S102: Image enhancement with multi-scale and adaptive background subtraction.

[0041] A sequence of binocular polarization images of underwater tracer particles is acquired, intensity images are calculated, and multi-scale filtering and adaptive background subtraction are performed to separate the tracer particles from the backscattered background, resulting in a particle-enhanced image.

[0042] The process of calculating the intensity image and performing multi-scale filtering and adaptive background subtraction to obtain the particle-enhanced image can be described as follows: The original polarized mosaic image is obtained by using a split-plane polarization camera. The intensity images of the four polarization channels are obtained by using gradient-guided interpolation. The total intensity image is obtained by summing the intensity images of the four polarization channels and taking half of the sum. A Gaussian filter bank with multiple scales is constructed, and corresponding weights are assigned to Gaussian kernels of different scales. Multi-scale weighted convolution response calculation is performed on the total intensity image to suppress high-frequency detection noise while compensating for the difference in particle light intensity distribution caused by depth changes, thus obtaining a denoised image.

[0043] Construct a local feature sliding window that is larger than a preset multiple of the particle's maximum diameter, calculate the local gray level expectation of the pixels in the denoised image within the window as the background estimate, and subtract the denoised image from the local gray level expectation to achieve background subtraction. Extract the minimum and peak gray values ​​of the non-zero pixels in the image after background subtraction, construct a linear mapping function to extend the particle gray-scale response to full scale, expand the brightness gradient between the particle center and edge, and output a particle-enhanced image.

[0044] Specifically, in an optional embodiment of this application, the process of calculating the intensity image from the polarized mosaic image is as follows: First, the processor can obtain the original polarized mosaic image using a focal plane polarization camera, and then extract the polarized mosaic image using gradient-guided interpolation. Intensity images of four polarization channels Finally, the intensity image is calculated. : (5) Furthermore, considering that underwater tracer particles are distributed at different depths along with the fluid, their projected diameter on the image plane exhibits scale uncertainty. To achieve adaptive noise reduction for all particles, the processor can construct a system containing... A filter bank of several scales. For the intensity image of tracer particles. The processor can calculate the multi-scale weighted response using the following formula. : (6) In the formula, The scaling factor is Two-dimensional Gaussian kernel, The contribution weights for the corresponding scales. The value range is determined by the equivalent pixel diameter of the particle in the near-field and far-field views. Through multi-scale convolution, this application can effectively compensate for the differences in particle light intensity distribution caused by underwater defocusing and depth changes while suppressing high-frequency detection noise, and prevent feature loss or false contour generation caused by single-scale filtering.

[0045] Furthermore, the processor can be constructed with a size of A local feature sliding window, where the window size The particle diameter must be greater than three times the maximum particle diameter to avoid suppressing the effective signal, and the local grayscale expectation of the pixels within the window must be calculated. : (7) Then, the processor can perform adaptive background subtraction using background subtraction to separate particles from low-frequency backscattering: (8) The process dynamically tracks and subtracts the low-frequency background field composed of uneven illumination and stray light, making the tracer particles stand out from the degraded background as isolated, bright patches.

[0046] Next, the processor can process the image after background subtraction. Perform global grayscale statistics and extract the minimum grayscale value of non-zero pixels. With peak grayscale .

[0047] The processor can then construct the following linear mapping function to generate the final preprocessed image. : (9) This step expands the brightness gradient between the particle center and the edge by forcibly extending the particle grayscale response to full scale (0-255), thereby improving the robustness of centroid localization.

[0048] Step S103: Particle recognition overlap decoupling, polarization distortion correction and three-dimensional spatial positioning.

[0049] This step takes the particle-enhanced image as input, combines local dynamic thresholding and gradient extremum search algorithms to perform particle recognition and overlap decoupling, and extracts sub-pixel-level center coordinates; uses the calibration parameters to correct the polarization distortion of the sensor system, and obtains the corrected polarization features; constructs epipolar constraints based on the calibration parameters, and establishes corresponding points by combining the sub-pixel-level center coordinates, grayscale features and the corrected polarization features; and uses triangulation and layered medium refraction correction to solve for the three-dimensional spatial coordinates of the particles.

[0050] In this embodiment, the processor can combine local dynamic thresholding and gradient extremum search algorithms to perform particle recognition and overlap decoupling. The process of extracting sub-pixel level center coordinates can be as follows: A local sliding window is established based on particle-enhanced images. The mean and variance of pixel gray levels are statistically analyzed window by window. A pixel-level dynamic segmentation threshold is constructed by combining the water noise adaptive adjustment coefficient. The threshold is equal to the local gray level mean plus the product of the adjustment coefficient and the local gray level variance. Based on the threshold, pixel-by-pixel binarization is performed to generate particle candidate regions. For the particle candidate region, the first-order gray-level partial derivatives in the horizontal and vertical directions of the gray-level image are solved. The gradient magnitude and gradient direction angle of each pixel are calculated using the differential operator to construct a gray-level gradient vector field covering the entire connected domain. Local gray-level peak traversal and discrimination are performed within each candidate connected region. If there are two or more local peaks within the connected region and they meet the physical constraints, the gradient direction angle sequence is extracted pixel by pixel along the baseline path connecting the two peaks, and the gradient direction change rate of adjacent pixels is calculated. If the change rate exceeds the preset gradient turning threshold, the pixel is determined to be a gray-level gradient inflection point. The line connecting the continuous inflection points is used as the segmentation boundary to complete the decoupling segmentation of overlapping particles, and two-dimensional Gaussian fitting is performed on the decoupled single-particle region to extract sub-pixel level center coordinates.

[0051] In this embodiment, the process of correcting the polarization state distortion of the sensor system using calibration parameters to obtain the corrected polarization characteristics can be as follows: A polarization refraction physical model based on Stokes vectors and Mueller matrices is established. The Stokes vectors of the incident and outgoing light are physically correlated through a calibrated four-by-four Mueller matrix. The diagonal elements of the Mueller matrix characterize the degree of change of the linear polarization component of the refracted light and satisfy the properties of total light intensity conservation and circular polarization state invariance. For the water-glass interface and the glass-air interface in the light propagation path, the corresponding Mueller matrices are determined respectively, and a cascaded positive transformation relationship from the incident polarization state to the sensor observation polarization state is constructed. That is, the observation polarization state is equal to the product of the glass-air interface Mueller matrix and the water-glass interface Mueller matrix with respect to the initial polarization state. Inverse compensation is performed on the distorted polarization data acquired by the sensor using the inverse transformation of the cascaded matrix. That is, the inverse matrix of each interface Mueller matrix is ​​multiplied by the observed polarization state in turn to obtain the true Stokes vector after compensation. The degree of linear polarization and the polarization angle are then extracted as the corrected polarization features.

[0052] In this embodiment, the process of establishing corresponding points by combining sub-pixel-level center coordinates, grayscale features, and corrected polarization features, and then using triangulation and layered medium refraction correction to solve for the three-dimensional spatial coordinates of the particles can be as follows: Based on the calibration parameters, a basic matrix is ​​constructed, and the normalized distance from the candidate particle point to the epipolar line in the right image is calculated. When the distance is less than the epipolar line search bandwidth threshold, it is used as a candidate matching point. The horizontal and vertical search intervals are calculated by combining the effective range of disparity to narrow the search space.

[0053] Calculate the gray-level normalized cross-correlation coefficients of the left and right local windows; calculate the Stokes parameters from the intensity images of the four polarization directions, convert the polarization angles into periodic invariant features to construct polarization feature vectors, and calculate the polarization feature differences between the left and right particles; construct a matching cost function by weighted summation of the gray-level cross-correlation coefficients and polarization feature differences, select the particle pair with the minimum cost, and perform a two-way consistency determination to establish corresponding points.

[0054] After solving the initial three-dimensional coordinates based on the principle of triangulation, a layered medium refraction correction model is established for the refraction of the tank wall and water. The direction of the incident direction after refraction through the interface is calculated according to the vector form of the law of refraction. The true back projection ray in the water is solved, and the midpoint of the line connecting the shortest distance between the left and right refracted rays is found to obtain the three-dimensional spatial coordinates after refraction correction.

[0055] Specifically, in an optional embodiment of this application, the process of particle recognition and overlap decoupling is as follows: In order to decouple the overlapping and adhered regions of multiple particles, the processor can first construct a physical model for underwater tracer particle optical imaging.

[0056] Define the pixel coordinates of the image plane as The ideal noiseless imaging grayscale distribution of a single tracer particle satisfies a two-dimensional Gaussian distribution model: (10) In the formula, The pixel coordinates of the ideal center of the particle. For particle imaging peak grayscale amplitude, The particle imaging dispersion radius is denoted as .

[0057] The processor can combine the complex underwater imaging environment with the slowly varying scattering background grayscale component and the imaging random noise component to construct a particle imaging grayscale model adapted to the actual underwater scene: (11) In the formula, The background grayscale is the slowly varying scattering of light in the water body. This refers to Gaussian random noise in underwater imaging. When two or more tracer particles overlap in their line-of-sight or are densely adhered in a plane, the image grayscale exhibits a superposition of multiple Gaussian functions, with no clear grayscale valley boundaries in the image space. Traditional fixed-threshold segmentation and single-peak detection methods cannot identify multi-particle targets in the superimposed state, ultimately leading to problems such as multi-particle connected component fusion, missed detection, and false detection. This invention designs an adaptive segmentation and overlapping particle decoupling strategy based on this physical mechanism.

[0058] To extract particles, the processor can construct an adaptive pixel-level segmentation threshold using local neighborhood grayscale statistics. This threshold is based on the pixel-level characteristics within the image. Establish a size centered on . A local sliding window is used. The window size is determined based on the maximum imaging particle size of the tracer particles under experimental conditions to ensure that the window range can completely cover the single-particle imaging area and the surrounding local background area.

[0059] The processor can calculate the gray-level statistical features of local windows, statistically analyzing the mean and variance of pixel gray levels window by window to characterize the local background gray-level distribution characteristics. The specific calculation formula is as follows: (12) (13) Next, the processor can construct an adaptive segmentation threshold model. Combining prior features of underwater imaging signal-to-noise ratio, an adaptive adjustment coefficient for water noise is introduced to construct a pixel-level dynamic segmentation threshold, achieving adaptive adaptation to water environments with different turbidity levels. The threshold model is as follows: (14) In the formula, The adaptive adjustment coefficient for water noise is determined through multiple calibration experiments with different water turbidity levels. It is used to dynamically compensate for the impact of underwater environmental noise fluctuations on the segmentation effect.

[0060] Specifically, the processor can binarize the particle-enhanced image pixel by pixel based on the dynamic threshold to generate particle candidate regions.

[0061] Based on the above adaptive threshold, the original grayscale image is binarized pixel by pixel to separate the particle target region and the background region, generating a binary image of the particle candidate region. The specific determination rules for the binary image are as follows: (15) This step replaces the traditional global fixed threshold with a local dynamic threshold, which can effectively suppress underwater slowly varying background interference and weak random noise, fully preserve the edge detail features of low-contrast particles, realize the complete extraction of candidate connected regions of tracer particles across the entire scale, and avoid the technical problems of missed detection of small particles and broken regions of weak particles.

[0062] Next, the processor removes background pixels from the candidate connected components of the particles in the binary image.

[0063] Specifically, based on the assumption of continuous differentiability of image gray levels, the processor can construct a high-precision gray-level gradient vector field across the entire connected domain. This provides a core feature for differentiating single particles from overlapping and adhered particles. The first-order gray-level partial derivatives of this connected domain in the horizontal and vertical directions are then calculated. (16) The 3×3 Sobel operator is used to perform numerical differentiation operations on discrete pixels, suppressing gradient noise introduced by discrete image sampling, and calculating the gradient magnitude and gradient direction angle pixel by pixel: (17) (18) Finally, a gray-level gradient vector field covering the entire connected domain is generated. Among them, the single-particle connected domain conforms to the characteristics of a single-peak Gaussian gray-level distribution, and the gradient vector converges unidirectionally from the edge of the region to the center of the particle peak, with no gradient reversal inflection point in the entire domain; the overlapping and sticky multi-particle connected domain is a multi-peak gray-level superposition field, with significant gradient direction jump inflection points and gray-level valley regions between adjacent peak intervals. The above-mentioned differentiated topological features are the core discrimination criteria for multi-particle decoupling in this invention.

[0064] Local grayscale peak values ​​are traversed and identified within each candidate connected component. If the grayscale value of the current pixel within the connected component is greater than that of all its 8 neighboring pixels, then the point is recorded as a local grayscale peak value. If two or more local grayscale peak values ​​are detected within a single connected component, and the spatial distance between these two peak values ​​is greater than the minimum imaging particle size of a single particle, then the connected component is determined to be a multi-particle overlapping and adhesion region.

[0065] Next, the processor can perform the steps of multi-peak topology determination and decoupling of overlapping particle gradient inflection points.

[0066] Specifically, based on the gray-level gradient vector field of connected components, the processor can accurately distinguish and decouple connected components by detecting local gray-level peak topology and determining gradient inflection point thresholds. The above steps can specifically include three sub-steps: peak detection, connected component attribute determination, and gradient inflection point segmentation. Sub-step 1: Local grayscale peak detection. Perform 8-neighborhood peak detection on all valid pixels within each candidate connected component. If the grayscale value of the current pixel is greater than the grayscale values ​​of all its 8 neighboring pixels, then the pixel is determined to be a local grayscale peak. The peak detection criteria are as follows: (19) Sub-step 2: Connected Component Particle Attribute Determination. The number of valid local peaks within a single connected component is counted, and type determination is performed based on prior particle physics constraints. If only one valid local grayscale peak exists within a connected component, it is determined to be a single-particle connected component. If two or more local peaks exist within a connected component, and both conditions are met—the peak grayscale difference being less than a preset signal-to-noise ratio threshold and the peak spatial spacing being greater than the minimum imaging particle size of a single particle—the connected component is determined to be a multi-particle overlapping and adhering region, thus initiating the gradient inflection point decoupling process.

[0067] Sub-step 3: Gradient inflection point detection and decoupling of overlapping particles. For multi-peak connected components, select any two effective local peaks, use the line connecting the two points as the baseline detection path, and extract the gradient direction angle sequence pixel by pixel along the baseline path. Calculate the gradient direction change rate of adjacent pixels: (20) The gradient direction change rate and the preset gradient turning threshold are used. Comparison is performed; among which, the gradient turning threshold is used. Determined through particle imaging dispersion radius calibration experiments. If a pixel within the path satisfies... If the pixel is identified as a grayscale gradient inflection point, then the grayscale gradient inflection point corresponds to the grayscale valley separation position in multi-particle superposition imaging. Using the line connecting the continuous inflection points as the region segmentation boundary, the multi-particle connected domain in the fused state is accurately split into several independent single-particle sub-connected domains, completing the decoupling and segmentation processing of underwater overlapping and adhering tracer particles.

[0068] If the rate of change exceeds a preset gradient turning threshold, the pixel is identified as a grayscale gradient inflection point (corresponding to a grayscale valley value of multiple particles). The line connecting these continuous gradient inflection points is used as the segmentation boundary to complete the decoupling segmentation of overlapping particles. Finally, the above two-dimensional Gaussian fitting is performed on the decoupled single-particle region to extract high-precision, unbiased sub-pixel-level center coordinates.

[0069] The purpose of polarization distortion correction in sensor systems is to restore the true polarization state of underwater targets. During underwater detection, light rays are refracted sequentially through multiple layers of media (water-glass-air), causing distortion in the polarization state information acquired by the sensor. To eliminate this degradation effect, we first analyze the propagation law of light waves at the interfaces of these multiple media, establishing a physical model of polarization refraction based on the Stokes vector and the Mueller matrix, revealing the intrinsic mechanism of the deflection of linear polarization components at the refraction interfaces. Then, for the water-glass and glass-air interfaces, we determine calibrated 4×4 Mueller matrices describing the polarization state deflection, and construct a cascaded forward transformation relationship from the incident polarization state to the sensor-observed polarization state. Finally, we use the inverse transformation of the cascaded matrix to perform inverse compensation on the distorted polarization data acquired by the sensor, obtaining the compensated Stokes vector, thereby eliminating the polarization distortion caused by multi-layer media refraction and providing accurate polarization state information input for subsequent steps.

[0070] Specifically, firstly, the processor can establish a physical model of polarization refraction based on the Stokes vector and the Mueller matrix. In the polarization transformation of an optical system, the effect on the polarization state of light waves can be regarded as a mathematical transformation of the Stokes vector. Assume that the incident light reflected from the underwater target is linearly polarized, and that the multiple media (water, glass, air) through which the refraction passes are all non-magnetic and isotropic. Let the Stokes vector of the incident light be... The Stokes vector of the emitted light after refraction through the medium is Both are obtained through a calibrated 4×4 Mueller matrix. Establish physical connections: (twenty one) During the matrix transformation process described above, since the refraction process satisfies the law of conservation of total light intensity, the elements in the first row and first column of the Mueller matrix satisfy... This leads to the derivation of the Mueller matrix. of Matrix diagonal elements and The linear polarization components of refracted light are altered to varying degrees in the horizontal and vertical directions, the extent of which depends on the physical properties of the refracting surface. Since the refraction process does not change the inherent properties of circular polarization, the matrix elements... By constructing this model, the degradation characteristics of light wave polarization states in underwater multilayer media can be accurately described.

[0071] Next, the processor can construct a cascaded positive transformation relationship of polarization states for the two interfaces of "water-glass" and "glass-air".

[0072] Let the Mueller matrices of the "water-glass" interface and the "glass-air" interface in the light propagation path be denoted as follows: and The initial polarization state of the target reflected light on the Poincaré sphere undergoes a first polarization change upon passing through the water-glass interface and a second polarization change upon passing through the glass-air interface. Therefore, the complete polarization state change process of the target reflected light upon reaching the sensor can be expressed as: (twenty two) Finally, the processor can use the inverse transformation of the cascaded matrix to perform inverse compensation on the distorted polarization data.

[0073] To obtain true polarization information unaffected by multilayer media interference, the distorted polarization state captured by the sensor is analyzed. Perform the inverse transformation of the concatenated matrix to calculate the compensated Stokes vector. : (twenty three) Through the aforementioned inverse compensation calculation, the polarization state shift trajectory of the target characteristic light on the Poincaré sphere is restored in reverse. This correction step effectively removes the composite polarization distortion caused by the complex waterproof layer and air medium, providing accurate polarization state information for subsequent steps.

[0074] Furthermore, polar constraints are constructed based on calibration parameters to establish corresponding left and right points. The specific process is as follows: First, based on the obtained intrinsic and extrinsic parameters and distortion correction parameters of the stereo camera, a stereo epipolar geometry constraint is constructed. Let the intrinsic parameter matrices of the left and right cameras be respectively... , The rotation matrix and translation vector of the right camera relative to the left camera are respectively , For the first image identified in the left image The sub-pixel centroid of a particle, its homogeneous pixel coordinates are represented as follows: (twenty four) In the formula, Let be the homogeneous pixel coordinates of the subpixel centroid of the i-th particle in the left image; This represents the horizontal pixel coordinates of the particle in the left image; This represents the vertical pixel coordinates of the particle in the left image.

[0075] The homogeneous pixel coordinates of the candidate particles in the right image are represented as follows: (25) In the formula, Let be the homogeneous pixel coordinates of the j-th candidate particle in the right image; This represents the horizontal pixel coordinates of the candidate particle in the right image; The vertical pixel coordinates of the candidate particle in the right image.

[0076] The projection rays of the same spatial particle point on the normalized imaging planes of the left and right cameras are as follows: (26) In the formula, The direction of the projection ray on the normalized imaging plane of the left camera; This is the inverse of the intrinsic parameter matrix of the left camera; These are the homogeneous pixel coordinates of the left image.

[0077] (27) In the formula, The direction of the projection ray on the normalized imaging plane of the right camera; This is the inverse of the intrinsic parameter matrix of the right camera; These are the homogeneous pixel coordinates of the right image.

[0078] For the same spatial point, the optical centers of the left and right cameras, the projection ray of the left camera, and the projection ray of the right camera are coplanar, thus satisfying the epipolar geometric constraints: (28) In the formula, This is the transpose of the normalized projection ray direction of the right camera; Let be the antisymmetric matrix constructed from the translation vector t; This is the rotation matrix of the right camera relative to the left camera; The normalized projection ray direction for the left camera.

[0079] (29) In the formula, For the translation vector Constructed antisymmetric matrix; This represents the component of the translation vector in the X direction; This represents the component of the translation vector in the Y direction; Let Z be the component of the translation vector in the Z direction.

[0080] in, For the translation vector Construct an antisymmetric matrix.

[0081] Converting the normalized coordinates back to pixel coordinates yields the fundamental matrix. : (30) In the formula, The fundamental matrix of the binocular system; This is the transpose of the inverse matrix of the right camera's intrinsic parameters; It is a translation vector antisymmetric matrix; This is the binocular extrinsic rotation matrix; This is the inverse matrix of the intrinsic parameters of the left camera.

[0082] Particles with the same name in the left and right images should satisfy the following epipolar constraints: (31) In the formula, This is the transpose of the homogeneous coordinates of the j-th candidate particle in the right image; The basic matrix; Let be the homogeneous coordinates of the i-th particle in the left image.

[0083] For particle points in the left image Its corresponding polar line in the right image is: (32) In the formula, The polar line is mapped from the i-th particle in the left image to the right image. The basic matrix; Let i be the homogeneous coordinates of the i-th particle in the left image; The horizontal coordinate coefficient of the i-th right epipolar line; The vertical coordinate coefficient of the i-th right epipolar line; is the constant term in the polar equation.

[0084] Then the candidate matching points in the right image should satisfy: (33) In the formula, The horizontal coordinate coefficient of the i-th right epipolar line; The vertical coordinate coefficient of the i-th right epipolar line; This is a constant term for the i-th right epipolar line; Let J be the horizontal coordinate of the j-th candidate point in the right image; Let be the vertical coordinate of the j-th candidate point in the right image.

[0085] Considering calibration errors, particle centroid positioning errors, and underwater refraction disturbances, actual matching does not require candidate points to strictly fall on the epipolar line. Therefore, it is necessary to calculate the normalized distance from the candidate particle point to the epipolar line: (34) In the formula, The epipolar distance between the i-th particle in the left image and the j-th candidate point in the right image; For epipolar-constrained algebraic residuals; These are the horizontal coordinate coefficients of the polar lines; This represents the polar coordinate coefficient.

[0086] When the following formula is satisfied, the first image in the right image will be... The particle is the first in the left image. Candidate matching points for 1 particle: (35) In the formula, The epipolar distance between candidate matching points; This represents the polar search bandwidth threshold.

[0087] in, The epipolar search bandwidth threshold is determined by the dual-target weight projection error, particle sub-pixel localization error, and experimental medium refraction error, and its value ranges from [value missing]. .

[0088] Finally, robust search interval constraints based on disparity range are established. After stereo correction, the epipolar lines of the binocular images are corrected to horizontal lines, and the vertical coordinates of corresponding points on the left and right sides are approximately equal, with disparity existing only in the horizontal direction. According to the principle of binocular triangulation, the relationship between disparity d and depth Z is as follows: (36) In the formula, The horizontal parallax of corresponding points on the left and right; The horizontal coordinates of the left image; The horizontal coordinates of the right image; The equivalent focal length in the horizontal direction; The binocular baseline length; The particle point depth.

[0089] Where B is the baseline length. It is the horizontal equivalent focal length.

[0090] For the effective depth range of the flow field The effective range of parallax can be deduced as follows: (37) In the formula, This is the lower limit of parallax. This represents the upper limit of parallax. To effectively measure the upper limit of depth; This is the lower limit for effective measurement depth; The equivalent focal length in the horizontal direction; This represents the binocular baseline length, i.e., the range of disparity values. .

[0091] Combining the relationship between parallax and coordinates The theoretical constraint interval for the horizontal coordinates of the corresponding point in the right image can be obtained: (38) In the formula, The horizontal coordinates of the corresponding point in the right image; The horizontal coordinates of the particles in the left image; This represents the upper limit of parallax. This represents the lower limit of parallax.

[0092] Considering factors such as camera calibration error, stereo correction residual, and particle positioning noise, a horizontal tolerance is introduced. and vertical tolerance This expands the search range to robust constraints. Specifically, the horizontal search interval and vertical coordinate constraints are as follows: (39) (40) Let be the sub-pixel coordinates of the i-th particle in the left image. The coordinates of the candidate particles in the right image are given. For parallax direction tolerance, The vertical tolerance can be set according to the actual noise level (usually 1~3 pixels).

[0093] In this embodiment, the process of three-dimensional spatial mapping is as follows: First, the processor can construct grayscale-polarization joint similarity matching constraints. Let the particles in the left image... The local neighborhood window is Candidate particles in the right image The local neighborhood window is The corresponding grayscale images are respectively and The gray-level normalized cross-correlation coefficient is defined as: (41) in, and These are the average grayscale values ​​of the left and right local windows, respectively. To prevent the use of tiny positive numbers with a denominator of zero, the Stokes parameters are calculated from the intensity images of the four polarization directions for polarization images: (42) in, The image light intensity values ​​are for the four polarization directions: 0°, 45°, 90°, and 135°.

[0094] This yields the degree of linear polarization (DoLP) and the angle of polarization (AoLP): (43) Considering the polarization angle have Periodicity makes direct comparisons prone to angular jump errors. Therefore, the polarization angle is converted into a periodically invariant characteristic. (44) Then the first i The particle in the left image and the first j The polarization feature vectors of the candidate particles in the right figure are represented as follows: (45) in, This represents the normalized intensity characteristics. The polarization characteristic difference is defined as: (46) This is the polarization feature weight matrix.

[0095] The processor can then construct the particle matching cost function: (47) In the formula, The combined matching cost between the i-th particle in the left image and the j-th candidate particle in the right image; The weight of the epipolar distance term; The polar distance; The normalized scale for epipolar distance; Weights for grayscale-related items; The grayscale correlation coefficient between the left and right local windows; Weights for polarization feature terms; Differences in polarization characteristics; This is the normalization scale for polarization characteristics. , This is the polarization characteristic scale normalization parameter.

[0096] For each particle in the left image In the candidate set The particle in the right image that minimizes the matching cost is selected as the initial matching target: (48) To avoid false matches, the processor can further perform bidirectional consistency checks (left-to-right optimal and right-to-left optimal): (49) At the same time, a uniqueness criterion is introduced, assuming the optimal and suboptimal costs are... and If the following conditions are met: Then the matching pair is considered valid.

[0097] For cases with a large number of particles, the problem is structured as a global allocation problem: (50) In the formula, This is the optimal global matching matrix; The matching matrix to be optimized; This is the matching indicator for the i-th left image particle and the j-th right image particle; This represents the number of candidate particles in the left image. The number of candidate particles in the right image; This represents the cost of matching the corresponding candidate.

[0098] The constraints are Furthermore, the sum of each row and column is no greater than 1. Through this joint determination, the reliability of three-dimensional localization of tracer particles under underwater scattering conditions can be improved.

[0099] The process of solving for three-dimensional coordinates based on the principle of triangulation is as follows: Establish a 3D coordinate system with the optical center of the left camera as the origin, and let the projection matrices of the left and right cameras be respectively: (51) The spatial three-dimensional homogeneous coordinates of the matching pair Satisfy projection relationship: (52) The homogeneous equation system is obtained by eliminating the scale factor through cross product. ,in: (53) In the formula, This is the transpose of the linear trigonometric measurement coefficient matrix; The horizontal coordinates of the left image; The vertical coordinates of the left image; The horizontal coordinates of the right image; The vertical coordinates of the right image; , This is the transpose of each row vector of the left camera projection matrix; This is the transpose of the row vectors of the right camera projection matrix.

[0100] Through singular value decomposition Solving the least squares problem yields the initial three-dimensional coordinates: (54) In the formula, The initial three-dimensional coordinates are obtained from linear triangulation. These are the components of the initial 3D point along the three coordinate axes; These are the four components of a homogeneous coordinate vector.

[0101] Optimize by constructing a model that minimizes the reprojection error: (55) In the formula, The optimized 3D homogeneous coordinates; Let be the covariance matrix of the observation error of the k-th camera; Let k be the projection matrix of the k-th camera; Let K be the coordinates of the observed pixel in the k-th camera; These are the camera indices, representing the left and right cameras respectively.

[0102] A layered medium refraction correction model is established to address the refraction caused by the tank wall and the water. Based on Fermat's principle, Snell's law can be derived by taking the extreme value of the total optical path. (56) In the formula, The refractive index of air; The refractive index of glass; The refractive index of water; The angle of incidence on the air side; The angle of refraction of the glass layer; The angle of refraction within the water.

[0103] In vector form, the incident direction via normal vector Interface Reflection (set up , ): (57) In the formula, The unit direction vector of the refracted light ray; The unit direction vector of the incident ray before refraction; It is the cosine of the angle between the incident direction and the interface normal; The ratio of the refractive indices of adjacent media; Let be the unit normal vector of the refractive interface.

[0104] For the For each camera, the air-side incident direction corresponding to each pixel is: (58) In the formula, Let be the direction of the unit incident ray from the k-th camera on the air side; This is the transpose of the rotation matrix of the k-th camera; It is the inverse matrix of the k-th camera intrinsic parameter matrix; Let be the coordinates of the k-th camera observation pixel.

[0105] The point where light intersects with the air-glass interface is ,in: (59) In the formula, Let be the parameters of the intersection point between the k-th ray and the air-glass interface; The unit normal vector of the air-glass interface; Let the k-th camera be the optical center; For the air-glass interface plane parameters; This refers to the direction of the incident ray on the air side.

[0106] After two refractions, the true back-projected rays in the water body are obtained. .

[0107] Find the midpoint of the line connecting the shortest distances of the two refracted rays, let , , , Solve the system of linear equations: (60) In the formula, The coefficient of the first term in the system of linear equations; These are the coefficients of the cross terms in a system of linear equations; The coefficient of the second term in the system of linear equations; The optimal distance parameter for the refracted light from the left camera; The optimal distance parameter for the refracted light from the right camera; The directions of refracted light rays from the left and right cameras in the water; The direction of the common perpendicular line of the left and right light rays.

[0108] The three-dimensional coordinates after refraction correction are: (61) In the formula, A point in three-dimensional space obtained by refracting light rays from the left and right; The spatial position of the refracted light from the left camera at the optimal parameters; This represents the spatial position of the refracted light from the right camera at the optimal parameters.

[0109] If combined with LDS ranging, the depth residual can be fitted using a second-order polynomial: (62) In the formula, This represents the depth residual or depth correction amount between the binocular reconstructed depth and the actual LDS depth. The true depth measured by the laser displacement sensor; The estimated depth obtained from binocular reconstruction; To deeply correct the coefficients of the second-order polynomial.

[0110] After obtaining the parameters using least squares, the final depth coordinates are corrected as follows: The final result is a high-precision set of three-dimensional particle coordinates in the underwater space, which is used for subsequent velocity field reconstruction.

[0111] Calculate the air-side rays corresponding to the pixels of the left and right cameras. After two refractions ("air-glass" and "glass-water"), obtain the true back-projected rays in the water. Solve for the midpoint of the line connecting the shortest distances of the two refracted rays to obtain the high-precision three-dimensional spatial coordinates after refraction correction.

[0112] Step S104: Temporal particle correlation, velocity solution and physical constraint flow field post-processing.

[0113] Based on the three-dimensional spatial coordinates of adjacent frames in time and the corrected polarization features, a polarization-space fusion matching cost function is constructed to perform inter-frame particle correlation and solve for the three-dimensional velocity vector of the particles. The three-dimensional velocity vector is combined with the initial flow velocity of the flow field, and three-dimensional flow field post-processing based on hydrodynamic physical constraints is performed to obtain the final three-dimensional underwater flow field reconstruction result.

[0114] In this embodiment, based on the three-dimensional spatial coordinates of temporally adjacent frames and the corrected polarization features, a polarization-spatial fusion matching cost function is constructed to perform inter-frame particle correlation. The process of solving the three-dimensional velocity vector of the particles can be as follows: For effective particles and candidate particles in two adjacent frames, a three-dimensional spatial distance cost term and an inter-frame polarization degree feature difference cost term are defined respectively. Dimensionless weight coefficients are introduced, and the two types of cost terms are normalized by the maximum spatial distance and the maximum polarization degree difference in the candidate matching pair set to eliminate the dimension difference. Construct a particle matching cost function that integrates polarization and space, which is equal to the weighted sum of the normalized spatial distance cost and the normalized polarization degree difference cost; Traverse all cross-frame candidate particle matching combinations, select the particle pair corresponding to the minimum fusion matching cost as the optimal matching result to complete the temporal association, and calculate the temporal three-dimensional displacement and instantaneous three-dimensional velocity vector of the particles in combination with the inter-frame time interval of the camera.

[0115] Establish geometric effective domain constraints based on finite space: Define the effective measurement area of ​​the water tank and the set of wall planes, and eliminate velocity vectors that exceed the boundary or violate the physical rule that the wall is impenetrable; Establish physical extreme value constraints for velocity: Define the reference velocity of the flow field based on the initial flow velocity and the average inter-frame physical displacement, construct the physical upper limit of velocity by combining the three-dimensional positioning error, and eliminate abnormal vectors whose velocity magnitude exceeds the reasonable range of physics; Establish local consistency constraints for velocity length: Use a robust statistical method based on the neighborhood median and the median absolute deviation to calculate the normalized ratio of the target vector velocity length residual to the neighborhood length residual scale, and identify outlier vectors in length. Establish a local consistency constraint for velocity direction based on the median value of the unit spherical direction: Define the angular distance between two unit direction vectors, select the direction vector with the smallest angular distance from the neighborhood direction set as the local principal direction, calculate the normalized ratio of the angular residual between the target velocity direction and the local principal direction to the fluctuation scale of the neighborhood direction, and determine the outlier vector of the direction.

[0116] The process of processor performing underwater tracer particle time series data acquisition and preprocessing can be as follows: A split-focus plane polarization camera was used to acquire continuous temporal images of underwater tracer particles, and the camera sampling frame rate was set to... The time interval between adjacent frames is The acquired time-series images are sequentially processed with camera calibration, lens distortion correction, and binocular stereo 3D resolution. In the particle 3D coordinate resolution stage, the effective tracer particles are pre-screened based on the inherent optical characteristics of particle polarization, automatically removing invalid noise particles corresponding to backscattered stray light from the water and suspended micro-impurities. Finally, the 3D spatial coordinate parameters and temporal polarization characteristic parameters of the effective tracer particles in each frame are output.

[0117] Definition of the first Within the frame image The three-dimensional coordinate vector of an effective tracer particle is: (63) In the formula: These are the X, Y, and Z spatial coordinates of the particle in a three-dimensional Cartesian coordinate system.

[0118] The polarization parameter corresponding to the effective tracer particles is The range of values ​​is In a stable underwater flow environment, the scattering polarization characteristics of the same tracer particle exhibit temporal continuity and stability, with no abrupt changes in polarization degree parameters between adjacent frames. This can serve as an effective feature constraint for precise inter-frame particle correlation matching.

[0119] Specifically, in an optional embodiment of this application, the sub-process of correlating temporally adjacent frames and solving for the three-dimensional velocity vector is as follows: Regarding the first Effective particles of a frame With the Candidate particles of the frame Define the three-dimensional spatial distance cost term: (64) Define the cost term for inter-frame polarization degree feature difference: (65) Introducing dimensionless weighting coefficients satisfy The maximum spatial distance in the candidate matching pair set is used respectively. Maximum polarization difference Normalization is performed on the two types of cost terms to eliminate dimensional differences, and a polarization-space fusion particle matching cost function is constructed: (66) Iterate through all cross-frame candidate particle matching combinations and select the fusion matching cost. The particle pair corresponding to the minimum value is taken as the optimal matching result, thus completing the temporal particle association.

[0120] The process by which the processor solves for the particle's three-dimensional velocity vector can be as follows: Based on the previously obtained optimal inter-frame matched particle pairs, the temporal three-dimensional displacement and instantaneous velocity of the particles are calculated. (Matched particle pairs) The corresponding inter-frame 3D displacement vector is: (67) Combined with camera frame intervals Solve for the particle's three-dimensional instantaneous velocity vector: (68) In the formula: These are the instantaneous velocity components of the particle along the X, Y, and Z axes of the three-dimensional Cartesian coordinate system.

[0121] Furthermore, in order to eliminate velocity field distortion caused by underwater bubble interference, particle mismatch, or random positioning noise, the solved three-dimensional velocity vector is combined with the initial flow velocity of the flow field, and the following post-processing with seven dimensions of fluid dynamic physical constraints is performed.

[0122] Dimension 1: Define a three-dimensional flow field vector database.

[0123] The calculated three-dimensional flow field vector database is shown in the following formula: (69) in, Indicates the first The three-dimensional spatial position of a tracer particle at the current moment. This indicates the three-dimensional physical displacement of the particle between two adjacent frames. This represents the three-dimensional velocity vector corresponding to the particle. The quality factor is the time interval between two adjacent frames. This is a vector quality evaluation parameter that integrates matching confidence, 3D reconstruction error, reprojection error, and polarization similarity.

[0124] Dimension Two: Establishing geometric effective domain constraints based on finite space.

[0125] The effective measurement area of ​​a horizontally placed rectangular water tank is defined as shown in the following formula: (70) For any velocity vector Let its current spatial location be... The predicted position for the next moment obtained by velocity extrapolation is: If satisfied or This indicates that the starting position or extrapolated predicted endpoint of the velocity vector exceeds the effective measurement area of ​​the water tank. The velocity vector is determined to not meet the experimental geometric boundary constraints, is marked as a geometrically abnormal vector, and is removed.

[0126] Furthermore, the five inner walls of the water tank are constructed as a set of planes, as shown in the following expression: (71) in, For the first The unit normal vector of each wall. These are the characteristic parameters of the corresponding plane.

[0127] Define spatial points To the The signed distances between the walls are as follows: (72) When the particle is in the region adjacent to the wall, that is, when the condition is met... At that time, the normal component of the particle's velocity vector is checked: (73) If the velocity vector's tendency to move along the outer normal of the wall will cause the particle to pass through the boundary of the water tank in the next moment, that is, if the following condition is satisfied: (74) If the velocity vector violates the physical rule that the wall is impenetrable, it will be discarded.

[0128] Dimension 3: Establish velocity physical extreme value constraints between average inter-frame physical displacement and propeller initial flow velocity.

[0129] The experimental scenario is a horizontal cuboid water tank, with the internal fluid flowing underwater, driven by a constant-speed, low-speed propeller. Under ideal conditions, free from particle mismatch, bubble reflection interference, and 3D reconstruction jumps, the fluid velocity must remain within a physically reasonable range constrained by both the propeller's driving capability and the measured average particle motion characteristics. Based on this, a global physical extremum constraint is first applied to the velocity modulus.

[0130] Let the average inter-frame physical displacement obtained through particle motion statistics of adjacent frames be . The corresponding fluid characteristic velocity is Let the initial velocity of the propeller be... To avoid the underestimation of the velocity upper limit caused by the weighted averaging algorithm, the flow field reference velocity is defined as the upper envelope of both velocities, as shown in the following expression: (75) In the formula, The reference velocity for the flow field; The statistical velocity is obtained from the average displacement of the particles; The initial flow velocity of the propeller; This is the average displacement statistic of particles in adjacent frames; The time interval between adjacent frames.

[0131] If only obtain or One of the parameters, then Take the velocity value corresponding to the obtained parameters.

[0132] The 3D positioning error will be transmitted to the velocity calculation result. Let the covariance of the 3D positioning error between two adjacent frames be... and The inter-frame displacement error covariance is: (76) In the formula, This is the composite matrix of displacement covariance or uncertainty between adjacent frames; This is the particle position uncertainty matrix for the current frame; This is the particle position uncertainty matrix for the next frame.

[0133] The corresponding displacement error scale is: (77) In the formula, The scale for displacement uncertainty; This is the trace of the displacement uncertainty matrix.

[0134] Further derivation yields the velocity error scale as follows: (78) In the formula, The velocity uncertainty scale; The scale for displacement uncertainty; The time interval between adjacent frames.

[0135] Based on the above parameters, the physical upper limit of velocity is constructed as follows: (79) In the formula, This is the physical upper limit of the velocity modulus; This is the reference speed amplification factor; The reference velocity for the flow field; This is the uncertainty margin coefficient; This is the scale for velocity uncertainty.

[0136] in, For flow rate safety margin, This is the tolerance factor for error.

[0137] For the A three-dimensional velocity vector, with velocity magnitude as follows: (80) In the formula, Let be the magnitude of the i-th velocity vector; Let i be the i-th three-dimensional velocity vector; Let be the components of the i-th velocity vector in the X, Y, and Z directions.

[0138] The corresponding inter-frame physical displacement magnitude is: (81) In the formula, Let be the inter-frame displacement magnitude of the i-th particle; Let be the inter-frame displacement vector of the i-th particle; Let i be the velocity modulus; The time interval between adjacent frames.

[0139] If satisfied Or equivalence determination conditions If the velocity modulus of the vector significantly exceeds the physically reasonable range of water flow under the current propeller-driven operating conditions, it is identified as a physical extreme anomaly vector and is removed.

[0140] Regarding the lower speed limit, this method does not directly eliminate low-speed velocity vectors. The minimum resolvable speed is defined as follows: (82) In the formula, The threshold for low-speed reliability discrimination; This is the low-speed threshold coefficient; The velocity uncertainty scale; The scale for displacement uncertainty; The time interval between adjacent frames.

[0141] in, This is the low-speed reliability discrimination coefficient. If This indicates that the direction of the velocity vector may be dominated by localization noise and will not participate in subsequent directional outlier detection, but it will not be directly eliminated due to its low speed. This allows us to preserve weak flows, stagnant zones, or low-velocity structures near the wall that are far from the propeller region.

[0142] Dimension 4: Establish speed-length local consistency constraints based on the idea of ​​general outlier detection.

[0143] After completing the global velocity physical extreme value constraint screening, outlier detection of velocity moduli is further carried out based on the continuity characteristics of local fluid flow. Low-speed propeller-driven water flow exhibits spatial continuity; except for the propeller near-field, vortex core region, and shear layer region, the velocity moduli at adjacent spatial locations do not show abrupt changes at single points. Based on this principle, a robust statistical method using the neighborhood median and median absolute deviation is employed to detect outliers in the velocity moduli.

[0144] For the A velocity vector, with its three-dimensional spatial position Choose a radius centered on All neighborhood vectors within the range, or select the nearest one. For each vector, a neighborhood set is constructed, and the selection method is as follows: or (83) In the formula, Let i be the set of spatial neighborhoods of the i-th velocity vector; Let be the three-dimensional spatial positions of the i-th and j-th particles; The neighborhood search radius; This is the K-nearest neighbor supplement set used when the radius neighborhood is insufficient; For The set of K nearest neighbors centered on the center.

[0145] To ensure the accuracy of the statistical results, the target vector itself is removed when calculating the neighborhood statistics, resulting in a pure neighborhood velocity magnitude set: (84) In the formula, Let i be the set of velocity magnitudes in the neighborhood of the i-th particle; Let be the magnitude of the j-th velocity vector in the neighborhood; Let be the neighborhood set of the i-th particle.

[0146] Extracting local robust reference values ​​for neighborhood velocity moduli: (85) In the formula, The median of the neighborhood velocity magnitude; For sets Take the median value.

[0147] Define the velocity length residual of the target vector as: (86) In the formula, The absolute residual of the i-th velocity modulus relative to the median of its neighborhood; Let i be the velocity modulus; It is the median of the neighborhood velocity magnitude.

[0148] The neighborhood length residual scale is defined as follows: (87) In the formula, For the local robustness scale of the velocity modulus residual; The neighborhood velocity modulus; The median of the neighborhood velocity magnitude; This is the base of the velocity noise.

[0149] in, This is the base of the velocity noise, used to avoid numerical problems where the denominator is zero. It is a positive number of the same order of magnitude as the velocity error scale.

[0150] Based on the above parameters, the normalized length residual is constructed as follows: (88) In the formula, The velocity length outlier normalized residual; For the absolute residual of the velocity modulus; For velocity modulus residual robustness scale.

[0151] If satisfied If the velocity length of this vector exhibits an isolated abrupt change compared to the flow trend in its neighborhood, it is marked as a length outlier vector. The threshold is the length outlier threshold. This criterion relies on the median and the absolute deviation of the median to construct a local reference value. Abnormally large or small velocity data will not interfere with the neighborhood statistical results. Compared with traditional mean and variance criteria, it has a stronger ability to resist abnormal interference.

[0152] Step 5: Establish local consistency constraints on velocity direction based on the medoid of the unit spherical direction.

[0153] To avoid interference from the velocity magnitude on the direction determination results, this constraint normalizes the velocity vector before performing a direction consistency check. For constraints satisfying... The effective velocity vector, whose unit direction vector is defined as: (89) In the formula, Let be the unit direction vector of the i-th velocity vector; Let i be the i-th three-dimensional velocity vector; It is the second norm of the velocity vector.

[0154] for The low dynamic velocity vector, because its direction is easily dominated by 3D positioning noise, does not participate in the direction outlier detection, and only retains the discrimination results of its length constraint and geometric constraint.

[0155] The direction of velocity belongs to a three-dimensional unit spherical space. To avoid geometric bias caused by the independent median calculation of the three components, this method uses the unit spherical direction (medoid) as the local principal direction. First, the angular distance between two unit direction vectors is defined: (90) In the formula, The angle between two unit direction vectors; Let be the two unit direction vectors to be compared; This is a cutoff function that restricts the cosine value to the range [-1, 1].

[0156] Among them, the minute correction amount This is used to avoid the problem of the inverse cosine function input value exceeding the domain due to numerical errors.

[0157] In the neighborhood set After removing the target vector itself, the direction vector with the smallest median angular distance is selected from the remaining neighborhood direction set and determined as the local direction medoid, as shown in the following expression: (91) In the formula, The representative neighbor number is the one that minimizes the median difference in the neighborhood directions; Let be the neighborhood set of the i-th velocity vector; It is the angle between two velocity directions in the neighborhood.

[0158] The corresponding local mainstream direction is The angular residual between the target velocity direction and the local principal direction is defined as: (92) In the formula, Let be the angle between the i-th velocity direction and the local reference direction; Let i be the direction vector of the unit velocity. This is a local reference direction determined by the neighboring representative.

[0159] Define the neighborhood directional fluctuation scale as: (93) In the formula, This is a local robustness scale for the orientation angle residual; The angle between the neighborhood velocity direction and the local reference direction; This is the base of the angular noise.

[0160] in, This is the base of the angular noise, used to avoid numerical problems where the denominator is zero.

[0161] The normalized directional residuals are constructed as follows: (94) In the formula, Normalized residuals for outliers in the direction; This is the deviation angle in the current velocity direction; This is the robustness scale for the orientation angle residual.

[0162] If satisfied If the velocity direction of the vector changes abruptly from the mainstream direction in the neighborhood, it is determined to be an isolated outlier vector.

[0163] The physical logic behind this directional constraint is that the velocity direction of a real continuous water body can change continuously in the vortex core, shear layer, and recirculation zone, but will not exhibit random directional jumps that are completely contrary to the surrounding neighborhood at a single spatial point. Employing a unit spherical directional medoid statistical method can effectively solve the geometric deviation problem in three-dimensional directional statistics of the traditional component median method, making directional outlier identification more closely match the spatial directional distribution characteristics of the velocity vector.

[0164] Dimension Six: Perform adaptive threshold correction for the near-field and high-curvature regions of the propeller.

[0165] The flow field in this experiment was driven by a low-speed propeller, exhibiting significant velocity gradients and abrupt changes in direction in the near-field, wake, vortex core, and shear layer. If a uniform threshold for length and direction outliers is applied across the entire field, the true vortex flow structure may be misidentified as anomaly vectors. To address this issue, an adaptive threshold correction mechanism is introduced for regions with highly variable local flow.

[0166] Let the center position of the propeller be... The near-field region of the propeller is defined as follows: (95) In the formula, This is the set of particle positions within the near-field influence region of the propeller; Let i be the three-dimensional position of the i-th particle; The location of the propeller center or near-field reference center; The radius of influence of the propeller in the near field.

[0167] in, The radius of influence of the propeller in the near field.

[0168] Simultaneously, a local first-order linear model is used to fit the velocity field within each vector neighborhood to estimate the local velocity gradient, as expressed below: (96) In the formula, Let be the local velocity gradient matrix at the i-th position; The difference in speed between neighboring regions; The difference in neighborhood location; Let be the neighborhood set of the i-th particle.

[0169] in, The local velocity gradient matrix has the following specific form: (97) In the formula, This is the local velocity gradient matrix; Let be the partial derivatives of the velocity component u with respect to x, y, and z; Let v be the partial derivative of the velocity component v with respect to x, y, and z. Let w be the partial derivative of the velocity component w with respect to x, y, and z.

[0170] The velocity gradient matrix is ​​solved using the weighted least squares method. (98) In the formula, This is the local velocity gradient matrix obtained by weighted least squares estimation; Spatial distance weights; The difference in speed between neighboring regions; This represents the difference in location within the neighborhood.

[0171] in, The spatial distance weight is determined as follows: (99) In the formula, The spatial distance weight between the i-th point and the j-th neighboring point; Let i and j be the positions of the i-th and j-th particles; This is the distance-weighted smoothing scale.

[0172] Based on the obtained velocity gradient matrix, calculate the local flow field curl: (100) In the formula, Let i be the curl of the velocity field at the i-th position; For the curl operator of the velocity field; The velocity vector has three components; The direction of the three-dimensional spatial coordinates.

[0173] If any of the following conditions are met: or or (101) In the formula, The particles are located in the near-field influence region of the propeller; For local curl modulus; This is the curl threshold; The F norm of the velocity gradient matrix; This is the velocity gradient threshold.

[0174] The spatial point is then determined to be in the propeller's near-field, high-curvature region, or high-velocity gradient region. For these regions, a relaxed adaptive outlier threshold is applied. , (102) In the formula, The adaptive length outlier threshold for the i-th vector; The adaptive directional outlier threshold for the i-th vector; This is the local operating condition adjustment coefficient; The base length outlier threshold; The outlier threshold is based on the direction of the base.

[0175] in For ordinary low-gradient regions, then take .

[0176] The revised length outlier criterion is: (103) In the formula, The velocity length outlier normalized residual; This is an adaptive length outlier threshold.

[0177] The revised outlier criterion for direction is: (104) In the formula, The outlier normalized residual is the velocity direction. This is the adaptive outlier threshold.

[0178] This adaptive thresholding mechanism can accurately remove isolated outlier vectors while fully preserving the real vortex, backflow, and shear flow structures generated by propeller drive, thus avoiding the loss of flow field characteristics caused by excessive post-processing.

[0179] Dimension 7: Establish soft constraints for the consistency of divergence in incompressible fluids.

[0180] The medium used in this experiment is water, which can be approximated as an incompressible fluid under low-speed flow conditions, satisfying the continuity equation for incompressible fluids: (105) In the formula, Let the velocity field divergence be denoted as . Let be the partial derivative of the velocity component u with respect to x; Let v be the partial derivative of the velocity component v with respect to y; Let w be the partial derivative of the velocity component w with respect to z.

[0181] Based on the aforementioned solution of the local velocity gradient matrix The local flow field divergence can be calculated: (106) In the formula, This is the divergence estimate at the i-th position; The trace of the local velocity gradient matrix; Let be the component expansion of the velocity field divergence.

[0182] Considering issues such as sparse particle distribution, incomplete boundary neighborhood data, and local interpolation errors, using divergence consistency as a hard rejection condition could easily lead to excessive rejection of true vectors. Therefore, this method uses divergence consistency as a confidence correction factor to construct a soft constraint. The normalized divergence residual is defined as follows: (107) In the formula, This is a normalization index for divergence anomalies; This is a local divergence estimate; The F norm of the local velocity gradient matrix; To prevent divergence noise bases with a denominator of zero.

[0183] in, To prevent tiny positive numbers with a denominator of zero. If If the vector is not found to be out of the loop, the confidence level of the velocity vector will be reduced. Only when the vector satisfies the criteria for either length outlier or direction outlier will the divergence anomaly be used as auxiliary evidence of anomaly to complete the vector removal.

[0184] This processing method incorporates the prior physical knowledge that water is incompressible, effectively improving the physical rationality of flow field post-processing, and can avoid the problem of erroneous deletion of real flow structures caused by unstable local gradient estimation in the near field of propellers and sparse particle regions.

[0185] If the divergence residual is large, the quality evaluation parameter (quality factor) of the velocity vector is reduced proportionally. Only when the velocity vector falls on the boundary outlier edge in a previous dimension and the divergence residual is abnormally high is it identified as a noise vector and removed. This soft divergence constraint incorporates incompressible fluid dynamics priors while avoiding false deletions caused by hard thresholds.

[0186] By post-processing the physical constraints of the above six dimensions, all noise and outliers in the flow field are filtered out, and finally, a highly robust and accurate three-dimensional underwater flow field reconstruction result is output.

[0187] To achieve the same design concept and solve the same technical problems as the methods described above, a second embodiment of this application provides a three-dimensional underwater flow field reconstruction device for complex lighting conditions. This device may include the following functional modules: The initial flow velocity calculation module is used to build an underwater measurement system that includes a binocular polarization camera and a laser displacement sensor. It uses the real physical distance provided by the laser displacement sensor as a constraint for joint calibration to obtain calibration parameters. A single tracer particle is released in the flow field, and the actual physical displacement between frames of the tracer particle is calculated based on the calibration parameters and the particle motion image sequence to obtain the initial flow velocity of the flow field.

[0188] The particle image enhancement module is used to acquire a sequence of binocular polarization images of underwater tracer particles, calculate the intensity image, perform multi-scale filtering and adaptive background subtraction, separate the tracer particles from the backscattered background, and obtain the particle-enhanced image.

[0189] The three-dimensional spatial coordinate solving module is used to take the particle-enhanced image as input, combine local dynamic thresholding and gradient extremum search algorithms to perform particle recognition and overlap decoupling, and extract sub-pixel-level center coordinates; use the calibration parameters to correct the polarization distortion of the sensor system to obtain the corrected polarization features; construct epipolar constraints based on the calibration parameters, and establish corresponding points by combining the sub-pixel-level center coordinates, grayscale features and the corrected polarization features; and use triangulation and layered medium refraction correction to solve for the particle's three-dimensional spatial coordinates.

[0190] The output module is used to construct a polarization-space fusion matching cost function based on the three-dimensional spatial coordinates of the temporally adjacent frames and the corrected polarization features to perform inter-frame particle correlation and solve for the three-dimensional velocity vector of the particles; the three-dimensional velocity vector is combined with the initial flow velocity of the flow field to perform three-dimensional flow field post-processing based on hydrodynamic physical constraints to obtain the final three-dimensional underwater flow field reconstruction result.

[0191] It is understood that each functional module in this embodiment corresponds one-to-one with each step in the above method embodiment. Therefore, all the technical features, specific details, formula calculations, and technical effects of the method embodiment are present in the corresponding modules of this embodiment, and will not be elaborated upon here.

[0192] The third embodiment of this application provides an electronic device whose process can be as follows: a memory and a processor. The memory stores a computer program, and when the processor executes the computer program, it can implement all the steps of the three-dimensional underwater flow field reconstruction method for complex lighting conditions described in the first embodiment.

[0193] Specifically, the electronic device can take the form of various computing devices, such as servers, high-performance graphics workstations, industrial control computers, or embedded fluid computing terminals.

[0194] The electronic device includes a processor, memory, network interface, and database connected via a system bus. The processor provides computational and control capabilities, supporting the operation of the entire 3D flow field reconstruction system. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system, database, and computer programs. The database stores binocular polarization camera calibration parameters, tracer particle image sequences, polarization feature data, overlapping particle decoupling parameters, and 3D velocity vector data of the flow field. The internal memory provides high-speed, high-frequency cache support for the operation of the computer programs in the non-volatile storage media.

[0195] The network interface is used to establish a communication connection with external data sources (such as high-frequency high-speed polarization camera interfaces and laser rangefinder serial ports) to acquire or import particle images and physical distance data in real time. When the computer program is executed by the processor, it can drive the underlying hardware structures to implement the three-dimensional underwater flow field reconstruction method for complex lighting conditions provided in this application.

[0196] The fourth embodiment of this application provides a computer-readable storage medium on which a computer program is stored. When executed by a processor, the computer program can implement all the steps of the three-dimensional underwater flow field reconstruction method for complex lighting conditions described in the first embodiment.

[0197] The computer-readable storage medium provided in this application embodiment can be non-volatile or volatile. For example, it can include: optical disk, hard disk, magneto-optical disk, solid-state flash memory (SSD), random access memory (RAM), read-only memory (ROM), etc. When the instructions stored in the storage medium are called and executed by the computer's processor, they can efficiently control electronic devices to complete the fine extraction of underwater particles based on polarization physics characteristics, high-precision epipolar geometry positioning under multi-medium refraction, and velocity vector field noise reduction based on six-dimensional physical constraints of hydrodynamics, thereby achieving robust and accurate three-dimensional underwater flow field reconstruction.

[0198] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. This program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments of the methods described above. Furthermore, any references to memory, storage, databases, or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory.

[0199] 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.

[0200] The embodiments described above are merely examples of several implementation methods of this application, 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 this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this patent application should be determined by the appended claims. The above are merely preferred embodiments of this application and do not limit the patent scope of this application. Any equivalent structural or procedural transformations made using the content of this application's specification and drawings, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of this application.

Claims

1. A method for reconstructing three-dimensional underwater flow fields under complex lighting conditions, characterized in that, include: An underwater measurement system including a binocular polarization camera and a laser displacement sensor was constructed. The real physical distance provided by the laser displacement sensor was used as a constraint for joint calibration to obtain calibration parameters. A single tracer particle was released in the flow field. Based on the calibration parameters and the particle motion image sequence, the inter-frame actual physical displacement of the tracer particle was calculated to obtain the initial flow velocity of the flow field. A sequence of binocular polarization images of underwater tracer particles is acquired, the intensity image is calculated, and multi-scale filtering and adaptive background subtraction are performed to separate the tracer particles from the backscattered background and obtain the particle-enhanced image. Using the particle-enhanced image as input, particle recognition and overlap decoupling are performed by combining local dynamic thresholding and gradient extremum search algorithms to extract sub-pixel-level center coordinates; the polarization distortion of the sensor system is corrected using the calibration parameters to obtain the corrected polarization features; epipolar constraints are constructed based on the calibration parameters, and corresponding points are established by combining the sub-pixel-level center coordinates, grayscale features, and the corrected polarization features; the three-dimensional spatial coordinates of the particles are solved using triangulation and layered medium refraction correction. Based on the three-dimensional spatial coordinates of temporally adjacent frames and the corrected polarization features, a polarization-space fusion matching cost function is constructed to perform inter-frame particle correlation and solve for the three-dimensional velocity vector of the particles. The three-dimensional velocity vector is combined with the initial flow velocity of the flow field, and three-dimensional flow field post-processing based on hydrodynamic physical constraints is performed to obtain the final three-dimensional underwater flow field reconstruction result.

2. The method for reconstructing three-dimensional underwater flow fields under complex lighting conditions according to claim 1, characterized in that, Using the actual physical object distance provided by the laser displacement sensor as a constraint, joint calibration is performed to obtain calibration parameters, including: Set the extrinsic translation vector of the camera under different spatial attitudes, extract the axial calculated distance from the camera's optical center to the calibration plate plane, and obtain the actual axial physical distance measured by the laser displacement sensor; A joint optimization objective function with distance constraints is constructed, wherein the objective function is a weighted sum of a reprojection error term and a physical scale penalty term; wherein the reprojection error term is the sum of squared deviations between the actual extracted corner coordinates and the coordinates calculated by the nonlinear projection model, and the physical scale penalty term is the sum of squared deviations between the calculated axial distance and the actual axial physical distance multiplied by a weight penalty coefficient; The joint optimization objective function is solved using a nonlinear iterative algorithm to obtain the high-precision camera intrinsic parameters, extrinsic parameters, and distortion coefficients, which serve as the calibration parameters.

3. The method for reconstructing three-dimensional underwater flow fields under complex lighting conditions according to claim 1, characterized in that, Based on the calibration parameters and the particle motion image sequence, the inter-frame actual physical displacement of the tracer particles is calculated, and the initial flow velocity of the flow field is obtained, including: To locate the local region of interest where a single tracer particle is located, assuming that the light spot intensity of the tracer particle follows a two-dimensional symmetric Gaussian distribution, a mathematical fitting model is constructed that includes the peak amplitude of light intensity, the coordinates of the centroid to be solved, the standard deviation, and the background noise. The mathematical fitting model is linearized by taking the logarithm of both sides, and then converted into a quadratic polynomial form with respect to spatial coordinates. The coordinates and grayscale data of all pixels in the local region of interest are solved by the least squares method, and the sub-pixel level centroid coordinates of the tracer particles are obtained by inverse solution of the polynomial fitting coefficients. The average pixel displacement of particles in multiple adjacent frames is statistically analyzed, and combined with the pixel-to-physical scale mapping relationship in the calibration parameters, it is converted into the actual physical displacement between frames. The initial flow velocity of the flow field is then calculated by combining the video frame time interval.

4. The method for reconstructing three-dimensional underwater flow fields under complex lighting conditions according to claim 1, characterized in that, The intensity image is calculated and subjected to multi-scale filtering and adaptive background subtraction to obtain a particle-enhanced image, including: The original polarized mosaic image is obtained by using a split-plane polarization camera. The intensity images of the four polarization channels are obtained by using gradient-guided interpolation. The total intensity image is obtained by summing the intensity images of the four polarization channels and taking half of the sum. A Gaussian filter bank with multiple scales is constructed, and corresponding weights are assigned to Gaussian kernels of different scales. Multi-scale weighted convolution response calculation is performed on the total intensity image to suppress high-frequency detection noise while compensating for the difference in particle light intensity distribution caused by depth changes, thus obtaining a denoised image. A local feature sliding window larger than a preset multiple of the particle's maximum diameter is constructed. The local grayscale expectation of the pixels in the denoised image within the window is calculated as the background estimate. The background is subtracted by the difference between the denoised image and the local grayscale expectation. Extract the minimum and peak gray values ​​of the non-zero pixels in the image after background subtraction, construct a linear mapping function to extend the particle gray-scale response to full scale, expand the brightness gradient between the particle center and edge, and output the particle-enhanced image.

5. The method for reconstructing three-dimensional underwater flow fields under complex lighting conditions according to claim 1, characterized in that, The process of combining local dynamic thresholding and gradient extremum search algorithms for particle recognition and overlap decoupling, and extracting sub-pixel-level center coordinates, includes: A local sliding window is established based on the particle-enhanced image. The mean and variance of pixel gray levels are calculated for each window. A pixel-level dynamic segmentation threshold is constructed by combining the water noise adaptive adjustment coefficient. The threshold is equal to the local gray level mean plus the product of the adjustment coefficient and the local gray level variance. Based on the threshold, pixel-by-pixel binarization is performed to generate particle candidate regions. For the particle candidate region, the first-order gray-level partial derivatives in the horizontal and vertical directions of the gray-level image are solved, and the gradient magnitude and gradient direction angle of each pixel are calculated using the differential operator to construct a gray-level gradient vector field covering the entire connected domain. Local grayscale peak traversal and discrimination are performed within each candidate connected region. If there are two or more local peaks within the connected region and they meet the physical constraints, the gradient direction angle sequence is extracted pixel by pixel along the reference path connecting the two peaks, and the gradient direction change rate of adjacent pixels is calculated. If the change rate exceeds the preset gradient turning threshold, the pixel is determined to be a grayscale gradient inflection point. The line connecting the consecutive inflection points is used as the segmentation boundary to complete the decoupling segmentation of overlapping particles, and two-dimensional Gaussian fitting is performed on the decoupled single-particle region to extract the sub-pixel-level center coordinates.

6. The method for reconstructing three-dimensional underwater flow fields under complex lighting conditions according to claim 1, characterized in that, The polarization state distortion of the sensor system is corrected using the calibration parameters to obtain the corrected polarization characteristics, including: A polarization refraction physical model based on Stokes vectors and Mueller matrices is established. The Stokes vectors of the incident and outgoing light are physically correlated through a calibrated four-by-four Mueller matrix. The diagonal elements of the Mueller matrix characterize the degree of change of the linear polarization component of the refracted light and satisfy the properties of total light intensity conservation and circular polarization state invariance. For the water-glass interface and the glass-air interface in the light propagation path, the corresponding Mueller matrices are determined respectively, and a cascaded positive transformation relationship from the incident polarization state to the sensor observation polarization state is constructed. That is, the observation polarization state is equal to the product of the glass-air interface Mueller matrix and the water-glass interface Mueller matrix with respect to the initial polarization state. The distorted polarization data acquired by the sensor is inversely compensated by using the inverse transformation of the cascaded matrix. That is, the observed polarization state is multiplied by the inverse of the Mueller matrix of each interface in turn to obtain the true Stokes vector after compensation. The linear polarization degree and polarization angle are then extracted as the corrected polarization features.

7. The method for reconstructing three-dimensional underwater flow fields under complex lighting conditions according to claim 1, characterized in that, The process of establishing corresponding points by combining the sub-pixel-level center coordinates, grayscale features, and the corrected polarization features, and then using triangulation and layered medium refraction correction to solve for the particle's three-dimensional spatial coordinates, includes: Based on the calibration parameters, a basic matrix is ​​constructed, and the normalized distance from the candidate particle point to the epipolar line in the right image is calculated. When the distance is less than the epipolar line search bandwidth threshold, it is used as a candidate matching point. The horizontal and vertical search intervals are calculated by combining the effective range of disparity to narrow the search space. Calculate the gray-level normalized cross-correlation coefficients of the left and right local windows; calculate the Stokes parameters from the intensity images of the four polarization directions, convert the polarization angles into periodic invariant features to construct polarization feature vectors, and calculate the polarization feature differences between the left and right particles; construct a matching cost function by weighted summation of the gray-level cross-correlation coefficients and polarization feature differences, select the particle pair with the minimum cost, and perform bidirectional consistency judgment to establish corresponding points; After solving the initial three-dimensional coordinates based on the principle of triangulation, a layered medium refraction correction model is established for the refraction of the water tank wall and the water body. The direction of the incident direction after refraction through the interface is calculated according to the vector form of the law of refraction. The true back projection ray in the water body is solved, and the midpoint of the line connecting the shortest distance between the left and right refracted rays is found to obtain the three-dimensional spatial coordinates after refraction correction.

8. The method for reconstructing three-dimensional underwater flow fields under complex lighting conditions according to claim 1, characterized in that, Based on the three-dimensional spatial coordinates of temporally adjacent frames and the corrected polarization features, a polarization-spatial fusion matching cost function is constructed to perform inter-frame particle correlation, and the three-dimensional velocity vector of the particles is solved, including: For effective particles and candidate particles in two adjacent frames, a three-dimensional spatial distance cost term and an inter-frame polarization degree feature difference cost term are defined respectively. Dimensionless weight coefficients are introduced, and the two types of cost terms are normalized by the maximum spatial distance and the maximum polarization degree difference in the candidate matching pair set to eliminate the dimension difference. Construct a particle matching cost function that integrates polarization and space, which is equal to the weighted sum of the normalized spatial distance cost and the normalized polarization degree difference cost; Traverse all cross-frame candidate particle matching combinations, select the particle pair corresponding to the minimum fusion matching cost as the optimal matching result to complete the temporal association, and calculate the temporal three-dimensional displacement and instantaneous three-dimensional velocity vector of the particles in combination with the inter-frame time interval of the camera.

9. The method for reconstructing three-dimensional underwater flow fields under complex lighting conditions according to claim 1, characterized in that, The three-dimensional velocity vector is combined with the initial flow velocity of the flow field to perform three-dimensional flow field post-processing based on fluid dynamics physical constraints, including: Establish geometric effective domain constraints based on finite space: Define the effective measurement area of ​​the water tank and the set of wall planes, and eliminate velocity vectors that exceed the boundary or violate the physical rule that the wall is impenetrable; Establish physical extreme value constraints for velocity: Define the reference velocity of the flow field based on the initial flow velocity and the average inter-frame physical displacement, construct the physical upper limit of velocity in combination with the three-dimensional positioning error, and eliminate abnormal vectors whose velocity magnitude exceeds the reasonable range of physics; Establish local consistency constraints for velocity length: Use a robust statistical method based on the neighborhood median and the median absolute deviation to calculate the normalized ratio of the target vector velocity length residual to the neighborhood length residual scale, and identify outlier vectors in length. Establish a local consistency constraint for velocity direction based on the median value of the unit spherical direction: Define the angular distance between two unit direction vectors, select the direction vector with the smallest angular distance from the neighborhood direction set as the local principal direction, calculate the normalized ratio of the angular residual between the target velocity direction and the local principal direction to the fluctuation scale of the neighborhood direction, and determine the outlier vector of the direction.

10. A three-dimensional underwater flow field reconstruction device for complex lighting conditions, characterized in that, include: The initial flow velocity calculation module is used to build an underwater measurement system including a binocular polarization camera and a laser displacement sensor. It uses the real physical distance provided by the laser displacement sensor as a constraint to perform joint calibration and obtain calibration parameters. A single tracer particle is released in the flow field, and the inter-frame actual physical displacement of the tracer particle is calculated based on the calibration parameters and the particle motion image sequence to obtain the initial flow velocity of the flow field. The particle image enhancement module is used to acquire a sequence of binocular polarization images of underwater tracer particles, calculate the intensity image and perform multi-scale filtering and adaptive background subtraction, separate the tracer particles from the backscattered background, and obtain the particle-enhanced image. The three-dimensional spatial coordinate solving module is used to take the particle-enhanced image as input, combine the local dynamic threshold and gradient extremum search algorithm to perform particle recognition and overlap decoupling, and extract sub-pixel-level center coordinates; use the calibration parameters to correct the polarization distortion of the sensor system to obtain the corrected polarization features; construct epipolar constraints based on the calibration parameters, establish corresponding points by combining the sub-pixel-level center coordinates, grayscale features and the corrected polarization features, and solve the three-dimensional spatial coordinates of the particles using triangulation and layered medium refraction correction. The output module is used to construct a polarization-space fusion matching cost function based on the three-dimensional spatial coordinates of temporally adjacent frames and the corrected polarization features to perform inter-frame particle correlation and solve for the three-dimensional velocity vector of the particles; the three-dimensional velocity vector is combined with the initial flow velocity of the flow field to perform three-dimensional flow field post-processing based on hydrodynamic physical constraints to obtain the final three-dimensional underwater flow field reconstruction result.

Images