Particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraint

Abstract

The present invention relates to a particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraints, comprising the following steps: performing multi-plane calibration to obtain a camera calibration function relationship; calculating the three-dimensional line of sight corresponding to the image pixel coordinates of each camera-captured image; acquiring a particle image, calculating the three-dimensional line of sight corresponding to the particle pixel coordinates in the particle image, and performing particle matching and triangulation reconstruction on each particle in parallel to obtain the particle's preliminary three-dimensional spatial position; performing nearest neighbor matching on the preliminary three-dimensional spatial position of each particle in parallel, and performing line of sight interpolation and line of sight dithering to generate the particle's precise three-dimensional spatial position; constructing a four-dimensional particle trajectory based on the precise three-dimensional spatial position of each particle; calculating the particle's velocity vector and acceleration vector, and obtaining a three-dimensional flow field through interpolation. Compared with the prior art, the present invention can compensate for refractive interfaces, improve measurement accuracy, and enhance the computational efficiency of 4D PTV technology based on parallel processing.

Description

Technical Field

[0001] The present invention belongs to the field of particle image velocimetry and relates to a four-dimensional Lagrangian particle trajectory velocimetry method, in particular to a particle trajectory velocimetry method based on multi-plane calibration and line-of-sight constraints. Background Art

[0002] Four-Dimensional Lagrangian Particle Tracking Velocimetry (4D PTV) is a non-contact, non-interference method for quantitatively measuring spatial flow fields. 4D PTV utilizes uniformly dispersed, tiny tracer particles to track the motion of the fluid. The tracer particle motion is recorded through continuous illumination and photography. Using algorithms such as three-dimensional particle matching, iterative particle reconstruction from images, trajectory initialization, Wiener filter prediction, and trajectory connection, 4D PTV extracts particle trajectory information in a time series (the fourth dimension) from continuous tracer particle images, providing an approximation of the velocity of a point in the fluid space. Four-dimensional Lagrangian particle trajectory velocimetry combines PTV technology with a triangulated reconstruction algorithm. It uses multiple perspectives (at least three cameras) to simultaneously capture and record the illuminated particle field, and reconstructs the three-dimensional spatial position of the tracer particles through three-dimensional particle matching and triangulated reconstruction algorithms. The tracer particles are then reprojected onto the image for iterative particle reconstruction to improve the three-dimensional spatial position accuracy of the particles and the image particle utilization rate, thereby being able to process particle images with high particle concentrations. The trajectory of the initial four frames of particle three-dimensional spatial positions is then initialized to obtain the initial four-dimensional Lagrangian particle trajectory. The possible positions of all particle trajectories in the next frame are then predicted based on the initial four-dimensional Lagrangian particle trajectory based on the Wiener filter, and iterative particle reconstruction is performed with the captured particle image to improve the three-dimensional spatial position accuracy of the particles and the image particle utilization rate. The particle trajectory is then optimally fitted to obtain a particle trajectory with higher accuracy. Finally, the first-order and second-order differentials of each spatial dimension of the particle trajectory are taken in the time dimension to obtain the velocity and acceleration of the particle at the three-dimensional spatial position.

[0003] Compared to traditional 2D PIV, 3D stereo PIV, and 3D tomographic PIV technologies, 4D PTV offers higher temporal and spatial resolution, enabling efficient quantitative measurement of instantaneous 3D flow velocity and acceleration, and full-field quantitative reconstruction of 3D flow fields. 4D Lagrangian particle trajectory velocimetry offers advantages such as ease of operation, high measurement accuracy, and greater information capture. However, due to the presence of the refractive interface between the air-plexiglass interface and the fluid during actual measurement, camera-captured images are severely distorted, leading to high calibration errors and very low computational efficiency. Consequently, higher-precision refraction compensation algorithms and more demanding camera angles are required.

[0004] In recent years, research on 4D Lagrangian particle trajectory velocimetry (PTV) has focused on improving the reconstruction accuracy and image utilization of 4D PTV techniques. Numerous high-precision reconstruction algorithms, including triangulation, iterative particle reconstruction (IPR), shake the box (STB), and Wiener filter prediction, have been proposed, improving the quality of 3D particle position reconstruction and the accuracy of 4D particle trajectories. These algorithms perform well when tested with synthetic particle images under ideal conditions. However, when applied to actual measurements, the camera must be perfectly aligned with the fluid-glass-air interface to avoid severe image distortion caused by the air-plexiglass interface and the fluid refraction interface. However, the complex geometric layout of the measurement object during actual fluid measurements makes it difficult to meet the stringent requirements of existing 4D PTV techniques. Furthermore, these algorithms, running on a single CPU core, are inefficient for reconstructing 4D particle trajectories over long time series. Scholars at home and abroad have conducted research on multi-camera calibration compensation to address the severe image distortion caused by the air-plexiglass interface-fluid refraction interface. However, the camera simulations obtained in these studies are more complex and difficult to integrate and apply in 4DPTV technology.

[0005] In summary, the severe distortion and distortion of particle images caused by the air-plexiglass interface-fluid refractive interface in fluid measurement seriously hinders the large-scale application of 4D PTV technology. Summary of the Invention

[0006] The purpose of the present invention is to overcome the defects of the above-mentioned prior art and provide a particle trajectory velocity measurement method based on multi-plane calibration and line of sight constraint, which effectively improves the measurement accuracy and calculation efficiency through multi-plane calibration, line of sight storage and refractive index compensation.

[0007] The purpose of the present invention can be achieved by the following technical solutions:

[0008] A particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraint includes the following steps:

[0009] 1) Calibrate multiple cameras based on multi-plane calibration images to obtain camera calibration function relationships;

[0010] 2) calculating the three-dimensional spatial line of sight corresponding to the image pixel coordinates of each camera-captured image according to the camera calibration function relationship, and storing the result in a line of sight library;

[0011] 3) acquiring a particle image, calculating the three-dimensional spatial line of sight corresponding to the particle pixel coordinates in the particle image based on the camera calibration function relationship, and performing particle matching and triangulation reconstruction on each particle in parallel to obtain a preliminary three-dimensional spatial position of the particle;

[0012] 4) performing nearest neighbor matching on the preliminary 3D spatial position of each particle with all 3D spatial sight lines in the sight line library in parallel, and performing sight line interpolation and sight line jittering to obtain the optimal reprojected particle image, thereby generating the precise 3D spatial position of the particle;

[0013] 5) Construct a four-dimensional particle trajectory based on the precise three-dimensional spatial position of each particle and calculate the velocity vector and acceleration vector of each particle at each moment;

[0014] 6) The velocity vector and acceleration vector of each particle are interpolated onto the orthogonal Euler grid to obtain the three-dimensional flow field.

[0015] Furthermore, in step 1), each camera photographs a plurality of calibration planes at different distances to obtain the multi-plane calibration image.

[0016] Furthermore, in step 2), the sightline library is obtained by the following steps:

[0017] Based on the camera calibration function relationship, the image pixel coordinates of each camera-captured image are mapped to the world coordinates of multiple calibration planes, and multiple intersection points of a single pixel and multiple planes are obtained. The multiple intersection points form a three-dimensional line of sight, and the direction unit vector of the three-dimensional line of sight is calculated. The image pixel coordinates and their corresponding direction unit vectors of the line of sight and their intersection points with each calibration plane are stored to form the line of sight library.

[0018] Furthermore, in step 3) and step 4), parallel processing of each particle is implemented based on OpenMP.

[0019] Furthermore, in step 3), the pixel coordinates of the particles in the particle image are determined by:

[0020] The particle pixel coordinates are obtained by searching the particle pixel peak and Gaussian fitting the particle image position center on the particle image actually captured.

[0021] Furthermore, in step 3), the particle matching and triangulation reconstruction specifically include:

[0022] 301) based on the three-dimensional space sight line and direction unit vector of each camera corresponding to each particle, obtain the particle image pixel position projected on each camera by the particle three-dimensional space position;

[0023] 302) mapping each of the particle image pixel positions to the world coordinates of a plurality of calibration planes, obtaining a plurality of intersections between the line of sight corresponding to the particle image pixel position and the plurality of calibration planes, and connecting the plurality of intersections to obtain a plurality of three-dimensional space lines of sight;

[0024] 303) Perform triangulation reconstruction based on each three-dimensional spatial line of sight to obtain the preliminary three-dimensional spatial position of the particle.

[0025] Furthermore, in step 4), the sight line interpolation and sight line jitter specifically include:

[0026] 401) performing three-dimensional space interpolation based on the plurality of sight lines obtained by the nearest neighbor matching to obtain a particle pixel coordinate, generating a particle pixel distribution on a new blank image based on the particle pixel coordinate as the center point of the particle image, and obtaining a reprojected particle image;

[0027] 402) obtaining the three-dimensional spatial sight line of each camera corresponding to the particle pixel coordinates;

[0028] 403) performing spatial translation and jittering on the particle pixel coordinates corresponding to the three-dimensional spatial line of sight of each camera on the reprojected particle image, subtracting the reprojected particle image generated by each jitter from the original particle image to obtain a residual particle image distribution, and taking the reprojected particle image with the smallest total particle brightness on the residual particle image as the optimal reprojected particle image;

[0029] 404) extracting the particle pixel coordinates on the best reprojected particle image, and calculating the particle pixel coordinates according to the camera calibration function relationship, performing line of sight triangulation reconstruction on the three-dimensional line of sight of each camera to obtain the precise three-dimensional spatial position of the particle.

[0030] Furthermore, in step 5), the construction of the four-dimensional particle trajectory specifically includes:

[0031] 501) Based on the precise three-dimensional spatial positions of the particles obtained in the previous frames, starting from the precise three-dimensional spatial positions of all particles in the first frame, all possible particle trajectories of the next frame are sequentially connected to form a particle trajectory tree network;

[0032] 502) traversing each branch of the particle trajectory tree network, and finding the branch with the smallest acceleration change as the optimal particle initial trajectory;

[0033] 503) Based on the initial trajectory of the particle in the previous frames, an N-order polynomial fitting is used to predict the possible position of the particle in the next frame;

[0034] 504) performing nearest neighbor matching on the possible positions and all three-dimensional sight lines in the sight line library, and performing sight line interpolation and sight line reprojection to obtain an optimal reprojected particle image;

[0035] 505) Based on the camera calibration function relationship, the three-dimensional spatial line of sight corresponding to the particle pixel coordinates in the optimal reprojected particle image is calculated, and triangulated reconstruction is performed to obtain the accurate three-dimensional spatial position of the particle in the next frame;

[0036] 506) Connecting the precise three-dimensional spatial position of the next frame of particles into the initial trajectory of the particles;

[0037] 507) Repeat steps 503)-506) until all particles at all times are connected.

[0038] Furthermore, in step 5), the calculation of the velocity vector and acceleration vector of each particle at each moment specifically includes:

[0039] According to the four-dimensional particle trajectory, the velocity vector of the particle at the nth moment is calculated based on the three-dimensional spatial positions of every three adjacent particles on the time series trajectory;

[0040]

[0041] in, is the velocity vector of the particle at the nth moment, X(n-1) and X(n+1) are the three-dimensional spatial positions of the particle at the n-1th and n+1th moments respectively, and Δt is the time interval between each two moments;

[0042] Based on the velocity vector of every two adjacent particles on the time series trajectory, the acceleration of the previous particle is calculated;

[0043]

[0044] in, and are the velocity vectors of the particle at the nth and n+1th moments, respectively. The acceleration vector of the particle at the nth moment.

[0045] Furthermore, in step 6), obtaining the three-dimensional flow field specifically includes:

[0046] Based on scattered point interpolation, the velocity vector and acceleration vector of each particle at each moment are interpolated to the orthogonal Euler grid nodes to obtain the velocity vector of the grid nodes and construct the Euler velocity field;

[0047] The Euler velocity field is corrected for error vectors to obtain the three-dimensional flow field.

[0048] The present invention can be applied to particle reconstruction and velocity field reconstruction of four-dimensional Lagrangian particle trajectory velocimetry (4D PIV) technology or equipment, provides a basis for the widespread application of 4D PTV technology in flow field reconstruction, and is of great significance for accelerated experimental fluid research.

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

[0050] 1. The present invention replaces the traditional camera model calibration with multi-plane calibration and stores the line of sight of each image pixel, which can compensate for the refractive interface existing in the actual measurement process and greatly improve the measurement accuracy.

[0051] 2. This invention decomposes time-consuming processes such as 3D particle matching, triangulation reconstruction, and line of sight jittering into multi-threaded parallel computing, greatly improving computing efficiency;

[0052] 3. The present invention establishes four-dimensional particle trajectories by first constructing the initial particle trajectories of the previous several frames, and then predicting, comparing, and accurately locating the possible positions of the particles in the next frame. This effectively improves the overall calculation accuracy while ensuring calculation efficiency.

[0053] 4. This invention constructs a three-dimensional flow field of particles through operations such as particle matching, triangulation reconstruction, and line-of-sight dithering. This can significantly improve the measurement accuracy of 4D PTV at the air-plexiglass interface-fluid refractive interface. It also greatly improves the computational efficiency of 4DPTV technology, reducing the tedious data processing work required by users when using 4D PTV technology-related hardware equipment, thereby saving a considerable amount of time and accelerating experimental fluid dynamics research. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 This is a flow chart of the particle trajectory velocity measurement method of the present invention;

[0055] Figure 2 A flowchart of stereo matching and triangulation reconstruction according to the present invention;

[0056] Figure 3 Flowchart of jittering sight line and triangulation reconstruction of the present invention;

[0057] Figure 4 A flow chart showing the connection of the initial trajectory of the present invention;

[0058] Figure 5 : The calibration images in the embodiment of the present invention, where (a) is the calibration image at Z = +15 mm, and (b) is the calibration image at Z = -15 mm;

[0059] Figure 6 Schematic diagram of the line of sight corresponding to each pixel in the calibration process according to an embodiment of the present invention;

[0060] Figure 7 Schematic diagram of particle image and distribution Gaussian fitting in an embodiment of the present invention;

[0061] Figure 8 Schematic diagram of the stereo matching process in an embodiment of the present invention;

[0062] Figure 9 is the particle image stereo matching result in an embodiment of the present invention;

[0063] Figure 10 Schematic diagram of the line of sight triangulation reconstruction process in an embodiment of the present invention;

[0064] Figure 11 This is a schematic diagram of line of sight interpolation in an embodiment of the present invention;

[0065] Figure 12 Schematic diagram of line of sight translation jitter in an embodiment of the present invention;

[0066] Figure 13 This is a schematic diagram of trajectory initialization in an embodiment of the present invention;

[0067] Figure 14 Schematic diagram of long trajectory and speed calculation in an embodiment of the present invention;

[0068] Figure 15 Schematic diagram of particle point interpolation to an orthogonal Euler grid in an embodiment of the present invention;

[0069] Figure 16 This is the three-dimensional velocity field result in an embodiment of the present invention. DETAILED DESCRIPTION

[0070] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments. This embodiment is implemented based on the technical solution of the present invention, and provides a detailed implementation method and specific operation process, but the protection scope of the present invention is not limited to the following embodiments.

[0071] Reference Figure 1 As shown, this embodiment provides a particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraint, which is used to implement four-dimensional Lagrangian particle trajectory velocity measurement, including the following steps:

[0072] S1. Calibrate multiple cameras based on the multi-plane calibration image to obtain the camera calibration function relationship.

[0073] 11) For the standard checkerboard calibration plane in the multi-camera three-dimensional complex flow field measurement system, each camera captures multiple calibration planes at different distances to obtain a multi-plane calibration image.

[0074] In this embodiment, for the standard checkerboard calibration plane in the multi-camera three-dimensional complex flow field measurement system, each camera shoots two calibration planes at different distances to obtain a dual-plane calibration image, such as Figure 5 Shown are two calibration images taken by camera 1 at Z=+15 mm and Z=-15 mm, respectively. The distance between the dots is ΔX=5 mm, and ΔY=5 mm.

[0075] 12) Perform camera calibration based on the multi-plane calibration image obtained in step 11) to obtain a calibration function relationship between the calibration image or particle image obtained by the camera and the plane physical three-dimensional space. The calibration function polynomial is:

[0076] u=a0+a1X+a2X 2 +a3X 3 +a4Y+a5Y 2 +a6Y 3 +a7XY+a8X 2 Y+a9XY 2

[0077] v=b0+b1X+b2X 2 +b3X 3 +b4Y+b5Y 2 +b6Y 3 +b7XY+b8X 2 Y+b9XY 2

[0078] X=c0+c1u+c2u 2 +c3u 3 +c4v+c5v 2 +c6v 3 +c7uv+c8u 2 v+c9uv 2

[0079] Y=d0+d1u+d2u 2 +d3u 3 +d4v+d5v 2 +d6v 3 +d7uv+d8u 2 v+d9uv 2

[0080] Where u and v are the pixel coordinates of the calibration image or particle image in the width and height directions, respectively; X and Y are the spatial coordinates in the X and Y axis directions on the calibration plane, respectively; a0-a9 are the coefficients of the physical space coordinates to the width direction of the calibration image or particle image; b0-b9 are the coefficients of the physical space coordinates to the height direction of the calibration image or particle image; c0-c9 are the coefficients of the calibration image or particle image pixel coordinates to the physical space coordinates in the X direction; and d0-d9 are the coefficients of the calibration image or particle image pixel coordinates to the physical space coordinates in the Y direction.

[0081] In this embodiment, there are 4 cameras and 2 calibration planes, and a total of 320 parameters are obtained.

[0082] S2. Calculate the three-dimensional spatial line of sight corresponding to the image pixel coordinates of each camera-captured image based on the camera calibration function relationship obtained in step 1), and store the calculated line of sight as a line of sight library.

[0083] 21) Mapping the image pixel coordinates of each camera-captured image to the world coordinates of multiple calibration planes using the camera calibration function relationship obtained in step 12), obtaining multiple intersection points between a single pixel and multiple planes, and connecting these intersection points to form a three-dimensional line of sight.

[0084] In this embodiment, for the standard checkerboard calibration plane in the multi-camera three-dimensional complex flow field measurement system, each camera shoots two calibration planes at different distances to obtain a dual-plane calibration image, such as Figure 6 The figure shows a schematic diagram of the line of sight corresponding to each pixel in the calibration process of this embodiment. The line of sight of each pixel will bend when passing through the refractive interface. The intersection of the line of sight of each pixel and the two calibration planes is calculated through the camera calibration function relationship obtained in step 12) to form a fixed line of sight.

[0085] 22) Calculate the direction unit vector of the sight line obtained in step 21), and store all image pixel coordinates and their corresponding direction unit vectors and their intersections with the two calibration planes to form the sight line library. The unit vector is calculated as follows:

[0086]

[0087] in, is the sight unit vector corresponding to j pixels, j is a natural number 1, 2, 3..., P 1,j and P 2,j are the intersection points of the sight lines corresponding to j pixels and any two calibration planes in the multiplane, ||P 1,j -P 2,j ||2 represents the geometric distance between the line of sight and the intersection of any two calibration planes in the multiplane.

[0088] S3. Obtain a particle image and calculate the three-dimensional spatial line of sight corresponding to the particle pixel coordinates in the particle image based on the camera calibration function relationship. Use a stereo matching algorithm to perform particle matching and triangulation reconstruction on each particle in parallel to obtain the particle's preliminary three-dimensional spatial position (i.e., the three-dimensional spatial point closest to the line of sight).

[0089] In this embodiment, parallel processing of particles is achieved based on OpenMP parallel programming, and stereo matching and triangulation reconstruction of each particle are completed by one OpenMP thread.

[0090] Take the four cameras of this embodiment as an example, Figure 2 As shown, step S3 specifically includes the following steps:

[0091] 31) Perform particle pixel peak search and Gaussian fitting of the particle image position center on the particle image actually captured.

[0092] like Figure 7 The image shown is the particle image distribution actually obtained in this embodiment, and Gaussian fitting is performed on it to obtain the particle pixel coordinates.

[0093] 32) The particle pixel coordinates of all camera particle images are calibrated according to the camera calibration function relationship obtained in step 1), and the three-dimensional spatial line of sight of each camera corresponding to the particle pixel coordinates is calculated and stored.

[0094] 33) Based on the line of sight and direction unit vector of the first camera corresponding to the single particle p1 obtained in step 32), three-dimensional stereo matching is performed, the line of sight of particle p1 corresponding to camera 1 is projected onto camera 2, and the epipolar line L12 of camera 1 to camera 2 is obtained. The particle image pixel position p12 within a certain area near the epipolar line L12 is found and stored.

[0095] 34) Project the line of sight of camera 1 corresponding to the single particle p1 onto camera 3 to obtain the epipolar line L13 of camera 1 to camera 3. Then project the line of sight of camera 2 corresponding to the particle p12 near the epipolar line L12 in step 33) onto camera 3 to obtain the epipolar line L23 of camera 2 to camera 3. Find the particle image pixel position p123 within a certain area near the intersection of these two epipolar lines.

[0096] In this embodiment, Figure 8The figure shows a detailed process diagram of particle matching. Based on the line of sight of camera 1 corresponding to the single particle p1 obtained in step 32), it is projected onto the second camera to obtain the epipolar line L12 of camera 1 to camera 2. The particle image pixel position p12 in a certain area near the epipolar line L12 is searched and stored. The line of sight of camera 1 corresponding to particle p1 is then projected onto camera 3 to obtain the epipolar line L13 of camera 1 to camera 3. The line of sight of camera 2 corresponding to the particle in the particle image pixel position p12 is then projected onto camera 3 to obtain the epipolar line L23 of camera 2 to camera 3. The particle image pixel position p123 in a certain area near the intersection of these two epipolar lines is searched.

[0097] 35) Project the line of sight of camera 1 for a single particle p1 onto camera 4 to obtain the epipolar line L14 from camera 1 to camera 4. Next, project the line of sight of camera 2 for particle p12 near the epipolar line from step 33) onto camera 4 to obtain the epipolar line L24 from camera 2 to camera 4. Finally, project the line of sight of camera 3 for particle p123 near the epipolar line intersection from step 34) onto camera 4 to obtain the epipolar line L34 from camera 2 to camera 4. Find the particle image pixel positions within a certain area near the intersection of these three epipolar lines. The particle image pixel position p1234 closest to the intersection is the particle image pixel position projected onto this image. This process continues for more cameras.

[0098] In this embodiment, the line of sight of particle p1 corresponding to camera 1 is projected onto camera 4 to obtain the epipolar line L14 from camera 1 to camera 4. The line of sight of the particle at particle image pixel position p12 corresponding to camera 2 is then projected onto camera 4 to obtain the epipolar line L24 from camera 2 to camera 4. The line of sight of the particle at particle image pixel position p123 corresponding to camera 3 is then projected onto camera 4 to obtain the epipolar line L34 from camera 3 to camera 4. The nearest particle image pixel position p1234 near the intersection of these three epipolar lines is then found. Figure 9 FIG. 1 is a schematic diagram of epipolar lines on a particle image during the stereo matching process in this embodiment.

[0099] 36) The particle pixel positions p1, p12, p123, and p1234 on each camera particle image obtained in step 35) are first mapped to the world coordinates of multiple calibration planes using the camera calibration function relationship obtained in step 12), and multiple intersection points of the lines of sight corresponding to the particle pixel positions and the multiple planes are obtained. These intersection points are connected to form four three-dimensional spatial lines of sight.

[0100] Calculate the direction unit vectors of the four lines of sight, and store the direction unit vectors of the four lines of sight and the intersection points of the four lines of sight with multiple calibration planes. The unit vector of each line of sight is calculated as follows:

[0101]

[0102] in, is the sight unit vector of the particle corresponding to the i-th camera, i is a natural number, and its value is 1, 2, 3, 4, P 1,i and P 2,i are the intersection points of the line of sight of the i-th camera corresponding to a particle and any two calibration planes in the multiplane, ||P 1,i -P 2,i ||2 represents the geometric distance between the line of sight of the particle corresponding to the i-th camera and the intersection of any two calibration planes in the multiplane.

[0103] Based on the four obtained line-of-sight vectors, triangulation reconstruction is performed. Due to the errors in the measurement process, the triangulation reconstruction does not necessarily intersect at one point. The three-dimensional space intersection point closest to the four lines of sight is obtained, and the three-dimensional space position of the particle is preliminarily obtained (that is, the three-dimensional space point closest to all lines of sight). The calculation process of the three-dimensional space intersection point closest to the four lines of sight is as follows:

[0104]

[0105] Where X is the three-dimensional spatial position of the particle, I is the third-order identity matrix, and T represents the transpose of the vector.

[0106] In this embodiment, the particle pixel positions p1, p12, p123, p1234 are adjusted according to Figure 10 As shown in the figure, the triangulated reconstruction of the four lines of sight is performed. Due to the errors in the measurement process, the triangulated reconstruction does not necessarily intersect at one point. The three-dimensional spatial intersection point closest to the four lines of sight is obtained, and the three-dimensional spatial position of the particle is preliminarily obtained.

[0107] S4. Perform nearest neighbor matching on the preliminary three-dimensional spatial position of each particle obtained in step S3 and all three-dimensional spatial sight lines in the sight line library in step S2 in parallel, and perform sight line interpolation and sight line jitter to obtain the optimal reprojected particle image, thereby generating the precise three-dimensional spatial position of the particle.

[0108] In this embodiment, the parallel processing of particles is implemented based on OpenMP parallel programming, and the line of sight jitter of each particle is completed by an OpenMP thread.

[0109] like Figure 3 As shown, step S4 specifically includes the following steps:

[0110] 41) Perform nearest neighbor matching on the three-dimensional spatial position of the particle obtained in step S3 and the corresponding sight lines of all pixel coordinates of all calibration images or particle images obtained and stored in step S2 to obtain the four nearest sight lines (sight lines corresponding to integer pixel coordinates) from the particle space coordinates, perform three-dimensional spatial interpolation of sight lines based on the four nearest sight lines, and then obtain the particle pixel coordinates of the four interpolated integer pixel coordinates of sight lines, generate a particle pixel distribution on a new blank image based on the interpolated particle pixel coordinates as the center point of the particle image, obtain a reprojected particle image, and then obtain the sight line and direction unit vector corresponding to this particle pixel coordinate for each camera.

[0111] like Figure 11 The figure shows a schematic diagram of line of sight interpolation in this embodiment. The four lines of sight closest to the three-dimensional spatial position of the particle (lines of sight corresponding to integer pixel coordinates) are obtained, and the three-dimensional spatial interpolation of the lines of sight is performed based on the four lines of sight, thereby obtaining the particle pixel coordinates corresponding to the four integer pixel coordinates interpolated on the reprojected particle image.

[0112] 42) The pixel coordinates of the particles in the reprojected particle image obtained in step 41) are spatially translated and jittered corresponding to the line of sight of each camera. The jittering process is an iterative process. The line of sight is first translated significantly in three spatial dimensions, and then translated slightly in three spatial dimensions. During the jittering process, the pixel coordinates of the particles in the reprojected particle image will also change. The residual particle image distribution is obtained by subtracting the reprojected particle image from the original particle image. The optimal reprojected particle image (i.e., the image closest to the original particle image) is obtained when the sum of the particle brightness on the residual particle image is minimized.

[0113] like Figure 12 The figure shows a schematic diagram of the line of sight translation jitter in this embodiment. The particle pixel coordinates will also change during the jitter process. The residual image distribution is obtained by subtracting the reprojected particle image from the original particle image. The optimal reprojected particle image (i.e., the image closest to the original particle image) is obtained when the sum of the particle brightness on the residual particle image is minimized.

[0114] 43) Based on step 42), the particle pixel coordinates are obtained based on the particle pixel peak search and the particle image position center Gaussian fitting according to the optimal reprojected particle image. The particle pixel coordinates are used to calibrate the camera calibration function relationship obtained according to step 12) to calculate the particle pixel coordinates corresponding to each camera's three-dimensional spatial line of sight, and perform line of sight triangulation reconstruction to obtain the precise three-dimensional spatial position of the particle (i.e., the three-dimensional spatial intersection of each line of sight).

[0115] S5. Repeat steps S3 and S4 four times to obtain the three-dimensional spatial position of the particles in the first four frames, connect the initial particle trajectories, and find the trajectory with the smallest acceleration change among all possible trajectories as the optimal initial trajectory, such as Figure 4 As shown, specifically:

[0116] 51) Based on the three-dimensional spatial positions of the particles in the first four frames, starting from the three-dimensional spatial positions of all particles in the first frame, all possible particle trajectories of the next frame are sequentially connected to form a particle trajectory tree network.

[0117] 52) Traverse each branch of the particle trajectory tree network and calculate three velocity vectors based on the two adjacent particle points on each branch. The specific process is to calculate the velocity vector of the particle at the nth moment based on the three-dimensional spatial position of each two adjacent particles on the initial trajectory (nth moment, n+1th moment):

[0118]

[0119] in, is the velocity vector of the particle at the nth moment, X(n+1) is the three-dimensional spatial position where the particle may appear at the n+1th moment, and Δt is the time interval between every two moments.

[0120] According to the three velocity values, an acceleration value is calculated for each two adjacent velocities:

[0121]

[0122] in, and are the nth and n+1th velocity vectors respectively, The acceleration vector of the particle at the nth moment.

[0123] 53) Traverse each branch of the particle trajectory tree network and find the branch with the smallest acceleration change as the optimal particle initial trajectory.

[0124] In this embodiment, based on the three-dimensional spatial positions of the particles in the first four frames, starting from the three-dimensional spatial position p1 of the particle in the first frame, the optimal initial particle trajectory with the smallest acceleration change is as follows: Figure 13 The first four particles are shown.

[0125] S6. Based on the optimal initial particle trajectory obtained in the first four frames, the possible position of the particle in the next frame is predicted based on N-order polynomial fitting, and the predicted position is projected onto the particle image of each camera to obtain a fifth frame of reprojected image. Compared with the fifth frame of measured image, the particle's three-dimensional spatial position is jittered corresponding to the line of sight of each camera to obtain the optimal reprojected image. The specific steps include:

[0126] 61) Based on the obtained optimal initial trajectory of the first four frames, the three-dimensional spatial position of the particle in the next frame is predicted based on N-order polynomial fitting.

[0127] In this embodiment, based on the particle trajectory track1 obtained by connecting the three-dimensional spatial positions of the particles in the first four frames, the three-dimensional spatial position p5 of the particle that may appear in the next frame is predicted based on the N-order polynomial fitting, as shown in FIG. Figure 13 As shown by the fifth particle, the predicted particles do not necessarily coincide with the actual measured particles.

[0128] 62) Perform nearest neighbor matching on the three-dimensional spatial position of the particle obtained in step 61) and the corresponding lines of sight of all image pixel coordinates obtained in step 2), and perform line of sight interpolation and reprojection according to the direction unit vector of the line of sight corresponding to the nearest neighbor pixel coordinate, and project it onto the particle image of each camera to obtain the fifth frame of reprojected image, subtract the reprojected image from the original image to obtain the residual image distribution, jitter the line of sight so that the total brightness of the particles on the residual image is minimized, and finally obtain the best reprojected image (i.e., the image closest to the original image).

[0129] In this embodiment, the predicted particles are Figure 12 The figure shows a schematic diagram of the line of sight translation jitter in this embodiment. The particle pixel coordinates will also change during the jitter process. The residual image distribution is obtained by subtracting the reprojected image from the original image. The optimal reprojected image (i.e., the image closest to the original image) is obtained when the sum of the particle brightness on the residual image is minimized.

[0130] 63) Based on step 62), the particle pixel coordinates are obtained based on the particle pixel peak search and the Gaussian fitting of the particle image position circle center according to the optimal reprojection image. The particle pixel coordinates are used to calibrate the camera calibration function relationship obtained according to step 1) to calculate the particle pixel coordinates and perform line of sight triangulation reconstruction corresponding to each camera's three-dimensional line of sight to obtain the accurate three-dimensional space position of the particle (i.e., the three-dimensional space point closest to each of the four lines of sight).

[0131] 64) Connect the precise three-dimensional spatial position of the particle obtained in step 63) to the existing initial trajectory to form a long-time series of four-dimensional Lagrangian particle trajectories. If the captured particle images have more moments, repeat step S6 until the particles at all moments are connected to the trajectory to construct a four-dimensional particle trajectory.

[0132] S7. Based on the four-dimensional particle trajectory, calculate the velocity and acceleration of each particle point at each moment on the time series trajectory, specifically including:

[0133] 71) According to the four-dimensional Lagrangian particle trajectory obtained in step 64), based on the three-dimensional spatial positions of every three adjacent particles on the time series trajectory (n-1th moment, nth moment, n+1th moment), the velocity of the particle at the nth moment is calculated.

[0134] Figure 14This is a schematic diagram of the long trajectory and velocity calculation in this embodiment. According to the four-dimensional Lagrangian particle trajectory, based on the three-dimensional spatial positions of the three particles on the time series trajectory (the 49th moment, the 50th moment, and the 51st moment), the velocity of the particle at the 50th moment is calculated. Based on the three-dimensional spatial positions of the three particles (the 50th moment, the 51st moment, and the 52nd moment) on the time series trajectory, the velocity of the particle at the 51st moment is

[0135] 72) Based on the particle point velocity obtained in step 71), the acceleration of the intermediate particle is calculated based on the velocities of every two adjacent particles on the time series trajectory.

[0136] The acceleration of the particle at the 50th moment in this embodiment is

[0137] S8. Based on the velocity and acceleration of each particle, interpolate them onto the orthogonal Euler grid to obtain the three-dimensional flow field, specifically including:

[0138] 81) The velocity vectors of the three-dimensional scattered particles obtained in step 71) and step 72) are interpolated based on scattered points. and the acceleration vector Interpolate to the orthogonal Euler grid node to obtain the velocity vector of the grid node Note: The subscript p represents three-dimensional scattered particles, and the subscript m represents orthogonal grid nodes.

[0139] Figure 15 The velocity vector of the particle point in the embodiment of the present invention is shown as follows: and the acceleration vector Interpolation to an orthogonal Euler grid diagram, using the nearest neighbor interpolation algorithm to obtain the velocity vector on the orthogonal Euler grid node

[0140] 82) Perform error vector correction on the Euler velocity field obtained in step 81).

[0141] In this embodiment, each velocity vector is traversed and compared with surrounding velocity vectors. If there is a large difference, the median of the surrounding velocity vectors is used for replacement and correction.

[0142] 83) Obtain the three-dimensional flow field, output and save the calculated three-dimensional particle field and three-dimensional velocity field results, and end the calculation.

[0143] The four-dimensional particle trajectory and the three-dimensional velocity field results at a certain moment in this embodiment are as follows: Figure 16 shown.

[0144] The performance improvement of using CPU single core and OpenMP multi-core parallel in this embodiment is shown in Table 1.

[0145] Table 1 Performance comparison of CPU single-core and OpenMP multi-core parallelization

[0146]

[0147] If the above method is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, or the part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes several instructions for enabling a computer device (which can be a personal computer, server, or network device, etc.) to execute all or part of the steps of the method described in each embodiment of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk, and other media that can store program code.

[0148] In another embodiment, an electronic device is provided, comprising one or more processors, a memory, and one or more programs stored in the memory, wherein the one or more programs include instructions for executing the particle trajectory measurement method based on multi-plane calibration and line-of-sight constraints as described above.

[0149] The present invention also provides a corresponding Figure 1The invention relates to a particle trajectory speed measuring device with multi-plane calibration and sight-line constraint, comprising a calibration module, a sight-line library generation module, a preliminary position acquisition module, a precise position acquisition module, a particle trajectory acquisition module and a three-dimensional flow field construction module, wherein the calibration module is used to calibrate multiple cameras based on multi-plane calibration images to obtain a camera calibration function relationship; the sight-line library generation module is used to calculate the three-dimensional space sight line corresponding to the image pixel coordinates of each camera-taken image according to the camera calibration function relationship, and store it as a sight-line library; the preliminary position acquisition module is used to obtain a particle image, calculate the three-dimensional space sight line corresponding to the particle pixel coordinates in the particle image based on the camera calibration function relationship, and parallelly calculate the particle trajectory speed measuring device with multi-plane calibration and sight-line constraint, comprising a calibration module, a sight-line library generation module, a preliminary position acquisition module, a precise position acquisition module, a particle trajectory acquisition module and a three-dimensional flow field construction module, wherein the calibration module is used to calibrate multiple cameras based on a multi-plane calibration image to obtain a camera calibration function relationship; the sight-line library generation module is used to calculate the three-dimensional space sight line corresponding to the image pixel coordinates of each camera-taken image according to the camera calibration function relationship, and store it as a sight-line library; the preliminary position acquisition module is used to obtain a particle image, calculate the three-dimensional space sight line corresponding to the particle pixel coordinates in the particle image based on the camera calibration function relationship, and parallelly calculate the particle trajectory speed measuring device with multi-plane calibration and sight-line constraint, Particle matching and triangulation reconstruction are performed on each particle to obtain the preliminary three-dimensional spatial position of the particle; the precise position acquisition module is used to perform nearest neighbor matching on the preliminary three-dimensional spatial position of each particle with all three-dimensional spatial sight lines in the sight line library in parallel, and perform sight line interpolation and sight line jitter to obtain the optimal reprojected particle image, and then generate the precise three-dimensional spatial position of the particle; the particle trajectory acquisition module is used to construct a four-dimensional particle trajectory based on the precise three-dimensional spatial position of each particle, and calculate the velocity vector and acceleration vector of each particle at each moment; the three-dimensional flow field construction module is used to interpolate the velocity vector and acceleration vector of each particle onto the orthogonal Euler grid to obtain the three-dimensional flow field.

[0150] The above describes in detail the preferred embodiments of the present invention. It should be understood that those skilled in the art can make numerous modifications and variations based on the concepts of the present invention without inventive effort. Therefore, any technical solutions that can be derived by those skilled in the art through logical analysis, reasoning, or limited experimentation based on the concepts of the present invention and the prior art should be within the scope of protection defined by the claims.

Claims

1. A particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraint, characterized in that: The following steps are involved: 1) Calibrate multiple cameras based on multi-plane calibration images to obtain the camera calibration function relationship; 2) Calculating the three-dimensional line of sight corresponding to the image pixel coordinates of each camera-captured image based on the camera calibration function relationship, and storing the result in a line of sight library; 3) Acquire a particle image, calculate the three-dimensional spatial line of sight corresponding to the particle pixel coordinates in the particle image based on the camera calibration function relationship, and perform particle matching and triangulation reconstruction on each particle in parallel to obtain a preliminary three-dimensional spatial position of the particle; 4) performing nearest neighbor matching on the preliminary 3D spatial position of each particle with all 3D spatial lines of sight in the line of sight library in parallel, and performing line of sight interpolation and line of sight dithering to obtain the optimal reprojected particle image, thereby generating the precise 3D spatial position of the particle; 5) Construct a four-dimensional particle trajectory based on the precise three-dimensional spatial position of each particle and calculate the velocity vector and acceleration vector of each particle at each moment; 6) Interpolate the velocity vector and acceleration vector of each particle onto the orthogonal Euler grid to obtain the three-dimensional flow field; In step 5), the construction of the four-dimensional particle trajectory specifically includes: 501) Based on the precise three-dimensional spatial positions of the particles obtained in the previous frames, starting from the precise three-dimensional spatial positions of all particles in the first frame, all possible particle trajectories of the next frame are sequentially connected to form a particle trajectory tree network; 502) traverse each branch of the particle trajectory tree network and find the branch with the smallest acceleration change as the optimal particle initial trajectory; 503) Based on the initial trajectory of the particle in the previous frames, an N-order polynomial fitting is used to predict the possible position of the particle in the next frame; 504) performing nearest neighbor matching on the possible position and all three-dimensional spatial sight lines in the sight line library, and performing sight line interpolation and sight line reprojection to obtain an optimal reprojected particle image; 505) Based on the camera calibration function relationship, the three-dimensional spatial line of sight corresponding to the particle pixel coordinates in the optimal reprojected particle image is calculated, and triangulated reconstruction is performed to obtain the accurate three-dimensional spatial position of the particle in the next frame; 506) Connecting the precise three-dimensional spatial position of the next frame of particles into the initial trajectory of the particles; 507) Repeat steps 503)-506) until all particles at all moments are connected.

2. The particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraint according to claim 1 is characterized in that: In step 1), each camera captures a plurality of calibration planes at different distances to obtain the multi-plane calibration image.

3. The particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraint according to claim 1 is characterized in that: In step 2), the sightline library is obtained by the following steps: Based on the camera calibration function relationship, the image pixel coordinates of each camera-captured image are mapped to the world coordinates of multiple calibration planes, and multiple intersection points of a single pixel and multiple planes are obtained. The multiple intersection points form a three-dimensional line of sight, and the direction unit vector of the three-dimensional line of sight is calculated. The image pixel coordinates and their corresponding direction unit vectors of the line of sight and their intersection points with each calibration plane are stored to form the line of sight library.

4. The particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraint according to claim 1 is characterized in that: In step 3) and step 4), parallel processing of each particle is implemented based on OpenMP.

5. The particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraint according to claim 1 is characterized in that: In step 3), the pixel coordinates of the particles in the particle image are determined by: The particle pixel coordinates are obtained by searching the particle pixel peak and Gaussian fitting the particle image position center on the particle image actually captured.

6. The particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraint according to claim 1 is characterized in that: In step 3), the particle matching and triangulation reconstruction specifically include: 301) Based on the three-dimensional space line of sight and direction unit vector of each camera corresponding to each particle, obtain the particle image pixel position projected on each camera by the particle three-dimensional space position; 302) Mapping each of the particle image pixel positions to the world coordinates of a plurality of calibration planes, obtaining a plurality of intersections between the lines of sight corresponding to the particle image pixel position and the plurality of calibration planes, and connecting the plurality of intersections to obtain a plurality of three-dimensional lines of sight; 303) Perform triangulation reconstruction based on each 3D spatial line of sight to obtain the preliminary 3D spatial position of the particle.

7. The particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraint according to claim 1 is characterized in that: In step 4), the sight line interpolation and sight line jitter specifically include: 401) performing three-dimensional spatial interpolation based on the plurality of sight lines obtained by the nearest neighbor matching to obtain a particle pixel coordinate, generating a particle pixel distribution on a new blank image based on the particle pixel coordinate as the center point of the particle image, and obtaining a reprojected particle image; 402) Obtaining the three-dimensional spatial line of sight of each camera corresponding to the particle pixel coordinates; 403) On the reprojected particle image, the particle pixel coordinates are spatially translated and jittered corresponding to the three-dimensional spatial line of sight of each camera, and the reprojected particle image generated by each jitter is subtracted from the original particle image to obtain a residual particle image distribution, and the reprojected particle image with the smallest total particle brightness on the residual particle image is taken as the optimal reprojected particle image; 404) Extracting the particle pixel coordinates on the optimal reprojected particle image, and calculating the particle pixel coordinates according to the camera calibration function relationship, performing line of sight triangulation reconstruction on the three-dimensional line of sight of each camera to obtain the precise three-dimensional spatial position of the particle.

8. The particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraint according to claim 1 is characterized in that: In step 5), the calculation of the velocity vector and acceleration vector of each particle at each moment specifically includes: According to the four-dimensional particle trajectory, the three-dimensional spatial position of each three adjacent particles on the time series trajectory is calculated. n The velocity vector of the particle at a moment; in, For the n The velocity vector of the particle at that moment is X ( n -1) and X ( n +1) respectively n -1 and n +1 moment's three-dimensional spatial position of the particle, is the time interval between every two moments; Based on the velocity vector of every two adjacent particles on the time series trajectory, the acceleration of the previous particle is calculated; in, and Respectively n and n +1 moment particle velocity vector, No. n The acceleration vector of the particle at that moment.

9. The particle trajectory velocity measurement method based on multi-plane calibration and line-of-sight constraint according to claim 1, characterized in that: In step 6), obtaining the three-dimensional flow field specifically includes: Based on scattered point interpolation, the velocity vector and acceleration vector of each particle at each moment are interpolated to the orthogonal Euler grid nodes to obtain the velocity vector of the grid nodes and construct the Euler velocity field; The Euler velocity field is corrected for error vectors to obtain the three-dimensional flow field.

Images