Large-view-field two-dimensional flow field splicing calibration method for particle image velocity measurement

Abstract

The invention relates to the technical field of image processing, solves the technical problems that vectors at splicing positions are abnormal during flow field splicing calculation among multiple cameras in large-view-field measurement and the deviation between the flow velocity and the actual flow velocity is calibrated, and particularly relates to a large-view-field two-dimensional flow field splicing calibration method for particle image velocity measurement. Comprising the steps of data acquisition, calibration parameter calculation, calibration parameter correction, homography matrix calculation, PIV displacement field calculation, displacement field splicing, scale transformation and offset calibration, velocity field calculation and the like. According to the method, the problem that vectors at splicing positions are abnormal during flow field splicing calculation among multiple cameras in PIV-2D2C large-view-field measurement is solved, meanwhile, when an image plane is calibrated to be not completely parallel to a laser plane, the deviation between the flow velocity and the actual flow velocity is calculated, and the PIV measurement precision is improved. Moreover, the method is high in universality, can complete the splicing calibration only through the additional shooting of one calibration image, is simple in operation, and is wide in application prospect.

Description

Technical Field

[0001] The invention relates to the technical field of image processing, and in particular to a large-field-of-view two-dimensional flow field stitching calibration method for particle image velocimetry. Background Art

[0002] In the field of image processing technology, Particle Image Velocimetry (PIV) is an important technology for measuring fluid motion. It is widely used in fluid velocity field measurement due to its full flow field, high precision and no interference with fluid motion. Tracer particles are spread in the area to be measured, and then analyzed using a high-speed camera after illuminating the tracer particles with a light sheet.

[0003] With the continuous advancement of technology and the rapid development of the field, the demand for high precision, large field of view, and high temporal and spatial resolution is also gradually increasing. The large field of view can effectively capture and analyze macroscopic flow structures. This technology can not only simultaneously observe large-scale flow phenomena (such as vortices and laminar flows), but also systematically identify aerodynamic phenomena such as airflow separation and reattachment. This overall observation capability enables researchers to quantify flow patterns and gain a deep understanding of the characteristics and behavior of fluid flows.

[0004] In order to ensure a sufficiently large measurement field of view, it is usually necessary to use multiple cameras to synchronize and simultaneously shoot the flow field area to be measured. For example, the Chinese invention patent with the publication number CN116205791A and the invention name is a method for stitching original particle images for large-field PIV measurement. The processing method adopted in this application is that each camera calculates the velocity vector field based on the original particle image taken, and then stitches the two sets of flow fields according to the reference point coordinates through scale calibration and reference points. This method requires that the two cameras are completely parallel to the laser plane and there is no rotation between the two cameras, which is almost difficult to achieve when the cameras are actually arranged. When the image planes of the two cameras appear in the above-mentioned non-parallel scene, obvious data discontinuity or abnormal data will appear at the data splicing.

[0005] Another method is to use Zhang’s calibration method to calibrate the internal and external parameters between camera groups. The position of the secondary camera image in the main camera image is calculated based on the external parameters. After the stitched image is obtained, the velocity vector field is calculated based on the stitched particle image. However, due to the influence of lighting, shooting angle, etc., there will be obvious brightness changes in the stitching seams of the stitching image; in addition, since the laser surface always has a certain thickness, and PIV generally measures penetrable fluids, there may be splitting between particles at the stitching seams. When the images in the stitching seam area are used for cross-correlation calculations, there is a high probability that abnormal values ​​will be calculated. In addition, at least 10 sets of calibration plate images freely placed in a public area are required for the calibration of camera internal and external parameters, which is difficult to do in some measurement scenarios.

[0006] When the camera image plane and the laser plane are not completely parallel, there will be a deviation between the calculated flow velocity and the actual flow velocity. Assuming that a uniform flow field is measured, there is an angle between the image plane and the laser plane. The speed calculated at the near light surface is too large, and the speed calculated at the far light surface is too small. However, the existing technology has not calibrated this deviation. Summary of the invention

[0007] In view of the shortcomings of the prior art, the present invention provides a large-field-of-view two-dimensional flow field stitching calibration method for particle image velocimetry, which solves the technical problems of anomalies in vectors at the stitching point when calculating the flow field stitching between multiple cameras in PIV-2D2C large-field-of-view measurement, and calibrating the deviation between the flow velocity and the actual flow velocity.

[0008] In order to solve the above technical problems, the present invention provides the following technical solutions: a large-field-of-view two-dimensional flow field stitching calibration method for particle image velocimetry, the method comprising the following steps:

[0009] S1, respectively obtaining a calibration image and a particle image located in the field of view of at least two cameras by a two-dimensional PIV stitching measurement device;

[0010] S2. Calculate the calibration parameters used to describe the conversion relationship between the two-dimensional pixel points in the calibration image and the three-dimensional points in the world coordinate system. The calibration parameters include the coefficients [l1, l2, ..., l 12 ], and the coefficients of the second-order polynomial model [a1, a2, …, a 10 ] and [b1,b2,…,b 10 ];

[0011] S3, calibrating the calibration parameters to obtain a coordinate mapping model for each camera to describe the conversion relationship between the spatial three-dimensional coordinates in the particle image and the distorted pixel coordinates (u, v);

[0012] S4. Calculate the homography matrix H from the pixel coordinates in any slave camera of at least two cameras to the pixel coordinates in the master camera, and the homography matrix H from any pixel coordinate in the master camera to the Z=0 plane based on the coordinate mapping model. w ;

[0013] S5, calculate the pixel displacement vector field of the particle image taken by each camera, the pixel displacement vector field includes the original coordinates g of n grid points in the main camera grid n (u n ,v n ), and the displacement Δg of the corresponding point is obtained by cross-correlation n (Δu n ,Δv n );

[0014] S6. In the pixel displacement vector field, determine the stitching data g of the slave camera in the master camera according to the homography matrix H. all_1n (u n ,v n ) and displacement Δg all_1n (Δu n ,Δv n ) complete the splicing;

[0015] S7, splice data g all_1n (u n ,v n ) and displacement Δg all_1n (Δu n ,Δv n ) scale transformation and offset calibration to obtain the physical coordinates G all_1n (x n ,y n ) and displacement ΔG all_1n (Δx n ,Δy n );

[0016] S8, use displacement ΔG all_1n (Δx n ,Δy n ) divided by the image acquisition interval T to obtain the velocity field.

[0017] Furthermore, in step S1, the specific process includes the following steps:

[0018] S11, arrange a two-dimensional PIV stitching measurement device, including a water tank, multiple cameras, a laser, a synchronous controller, a PC terminal, and a single-layer or double-layer calibration plate;

[0019] S12, turn on the laser, place a single-layer or double-layer calibration plate on the laser plane, make the calibration plate and the laser plane basically overlap, then turn off the laser, fine-tune the measurement field of view and depth of field of multiple cameras, and the cameras synchronously collect and save calibration images;

[0020] S13, remove the calibration plate from the water tank, spread the tracer particles that can be illuminated by the laser in the water tank, let the fluid in the water tank circulate for a period of time, turn on the laser, and synchronously collect and save the particle image with the camera.

[0021] Furthermore, in step S2, the specific process includes the following steps:

[0022] S21, based on the projection matrix P in the pinhole camera model, establish a linear transformation to convert the three-dimensional spatial coordinates (X, Y, Z) into undistorted pixel coordinates (u und ,v und )’s coordinate projection model, the expression of the coordinate projection model is:

[0023]

[0024]

[0025] Where P is the projection matrix; s is the scaling factor; [l1,l2,…,l 12 ] are the coefficients in the projection matrix P;

[0026] S22, obtain the detection pixel coordinates of all calibration points in the calibration plate (u i ',v i ') and the corresponding three-dimensional spatial coordinates (X i ,Y i ,Z i ), and bring it into the coordinate projection model to solve the coefficients [l1,l2,…,l 12 ]’s initial value;

[0027] S23, establishing a second-order polynomial model for describing the deviation of camera distortion, the expression of the second-order polynomial model is:

[0028]

[0029]

[0030] Where (u, v) is the distorted pixel coordinate; (Δu, Δv) is the distorted pixel offset; [a1, a2, …, a 10 ], [b1,b2,…,b 10 ] are the coefficients of the second-order polynomial model;

[0031] S24, the three-dimensional coordinates (X i ,Y i ,Z i ) and the projection matrix P are brought into the coordinate projection model to obtain the undistorted pixel coordinates (u iund ,v iund ), and use the detected pixel coordinates (u i ',v i ') minus the undistorted pixel coordinates (u iund ,v iund ) to obtain the distorted pixel offset (Δu i ,Δv i );

[0032] S25, offset the distorted pixel (Δu i ,Δv i ) is brought into the second-order polynomial model and the coefficients [a1, a2, …, a 10 ] and [b1,b2,…,b10 ]’s initial value;

[0033] S26, the three-dimensional coordinates (X i ,Y i ,Z i ) is brought into the coordinate projection model and the second-order polynomial model to obtain the calculated pixel coordinates (u i ,v i );

[0034] S27, construct an optimization function to calculate the pixel coordinates (u i ,v i ) and the detected pixel coordinates (u i ',v i ') The optimization equation is established with the minimum deviation as the goal, that is:

[0035]

[0036] Where m is the number of pixels i;

[0037] S28, use the Levenberg-Marquardt iterative optimization algorithm to optimize the coefficients [l1,l2,…,l 12 ]、[a1,a2,…,a 10 ] and [b1,b2,…,b 10 ] is optimized until the results converge and the calibration parameters are obtained.

[0038] Furthermore, in step S3, the specific process includes the following steps:

[0039] S31, calculating the maximum common area of ​​any two left and right cameras;

[0040] S32, mapping the particle images taken by the two cameras at the same time to the Z=0 plane respectively according to the calibration parameters, and obtaining a mapped particle image of the same size as the original image by interpolation of the maximum common area;

[0041] S33, the mapping particle images generated by the left and right cameras are divided into sub-areas and then cross-correlation calculation is performed to obtain the offset Δd of each sub-area of ​​the left and right cameras. i ;

[0042] S34, according to the offset Δd i Calculate the distance Δz between each vector and the corresponding point of each sub-area in the laser plane and the calibration plate plane i , distance Δz i The calculation formula is:

[0043]

[0044] In the formula, is the pixel u along the horizontal direction of camera n ni Partial derivatives with respect to X, Y, and Z; is the number of pixels v along the vertical direction of camera n ni Partial derivatives with respect to X, Y, and Z;

[0045] S35, refitting a new Z=0 plane using all intersection points on the maximum common area, using the distance Δz of the non-zero disparity vector i Move and rotate the coordinate system so that the new Z=0 plane coincides with the actual laser plane;

[0046] S36, calculating the rotation and translation matrix of the corresponding points in the new Z=0 plane and the actual laser plane, and performing perspective transformation on the spatial coordinates of the calibration points in the calibration plate to obtain the spatial coordinates of the calibration points when the laser plane is taken as the Z=0 plane;

[0047] S37, re-execute step S2 according to the spatial coordinates of the calibration points to obtain the corrected calibration parameters;

[0048] S38, bringing the corrected calibration parameters into the coordinate projection model and the second-order polynomial model to obtain a coordinate mapping model.

[0049] Further, in step S31, the maximum common area is the largest inscribed rectangle in the common area of ​​the fields of view of the two cameras.

[0050] Furthermore, in step S4, the specific process includes the following steps:

[0051] S41, for each pair of corresponding points p from the pixel coordinates in the slave camera to the pixel coordinates in the master camera 1i (u 1i ,v 1i ) and p 2i (u 2i ,v 2i ), the constructed equation is:

[0052]

[0053] Where [h1,h2,…,h9] is the coefficient of the homography matrix H;

[0054] S42. Generate a 2n×9 matrix A through n>4 corresponding points, and solve the coefficients [h1,h2,…,h9] according to the least squares method to obtain the homography matrix H. The expression of the homography matrix H is:

[0055]

[0056] S43, divide the homography matrix H by h9 to normalize and obtain the homography matrix H w .

[0057] Furthermore, in step S6, the specific process includes the following steps:

[0058] S61, the original coordinates g of the point to be measured from the camera 2n (u 2n ,v 2n ) plus the displacement Δg 2n (Δu 2n ,Δv 2n ), the coordinates of the slave camera in the main camera are obtained by transforming it to the main camera coordinate system through the homography matrix H. 1n ′(u 1n ′,v 1n ′) and displacement Δg 1n ′(Δu 1n ′,Δv 1n ′);

[0059] S62, according to the original coordinates g in the main camera 1n (u 1n ,v 1n ) and the coordinates g transformed from the camera 1n ′(u 1n ′,v 1n ') determine the maximum area;

[0060] S63, expand the grid points divided in the main camera, and convert the coordinates g 1n ′(u 1n ′,v 1n ′) is interpolated into the grid of the main camera through triangulation, and finally the spliced ​​data g is obtained all_1n (u n ,v n ) and displacement Δg all_1n (Δu n ,Δv n ).

[0061] Furthermore, in step S61, the conversion formula is:

[0062]

[0063]

[0064] Where s is the scaling factor and ps is the scaling factor from the camera to the main camera.

[0065] Furthermore, in step S7, the scale transformation and offset calibration formula is:

[0066]

[0067]

[0068] In the formula, H w It is the homography matrix from any pixel coordinate in the main camera to the Z=0 plane.

[0069] By means of the above technical solution, the present invention provides a large-field-of-view two-dimensional flow field stitching calibration method for particle image velocimetry, which has at least the following beneficial effects:

[0070] 1. The present invention solves the problem of abnormal vectors at the splicing point when calculating the flow field splicing between multiple cameras in PIV-2D2C large field of view measurement. At the same time, when the calibration image plane and the laser plane are not completely parallel, the deviation between the calculated flow velocity and the actual flow velocity is improved, thereby improving the measurement accuracy of PIV.

[0071] 2. The method proposed in the present invention has strong versatility and only requires taking one additional calibration image to complete the stitching calibration. It is simple to operate and has broad application prospects. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] The drawings described herein are used to provide a further understanding of the present application and constitute a part of the present application. The illustrative embodiments of the present application and their descriptions are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

[0073] Figure 1 It is a flow chart of the large field of view two-dimensional flow field stitching calibration method of the present invention;

[0074] Figure 2 It is a structural schematic diagram of the two-dimensional PIV splicing measurement device in the present invention;

[0075] Figure 3 A schematic diagram of the placement of the calibration plate in the two-dimensional PIV splicing measurement device of the present invention;

[0076] Figure 4 It is a schematic diagram of the displacement deviation introduced when the calibration plate and the laser plane do not coincide with each other in the present invention;

[0077] Figure 5 A schematic diagram for determining the maximum common area in the present invention;

[0078] Figure 6 A schematic diagram of dividing the mapped particle image according to sub-areas in the present invention;

[0079] Figure 7 A schematic diagram of a calibration image collected by a two-dimensional PIV stitching measurement device in the present invention;

[0080] Figure 8 is a schematic diagram of a particle image collected by a two-dimensional PIV stitching measurement device in the present invention;

[0081] Fig. 9is a schematic diagram of images to be stitched in the left camera of the present invention;

[0082] Fig.10 is a schematic diagram of images to be stitched in the right camera of the present invention;

[0083] Fig.11 It is a schematic diagram of the flow field after the images to be stitched in the left and right cameras in the present invention are stitched. DETAILED DESCRIPTION

[0084] In order to make the above-mentioned purposes, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific implementation methods, so that the implementation process of how the present application uses technical means to solve technical problems and achieve technical effects can be fully understood and implemented accordingly.

[0085] In order to solve the technical problems of abnormal vectors at the splicing point when calculating the flow field between multiple cameras in PIV-2D2C large field of view measurement, and calibrating the deviation between the flow velocity and the actual flow velocity, please refer to Figure 1-Figure 11 This embodiment proposes a large-field-of-view two-dimensional flow field stitching calibration method for particle image velocimetry, which solves the problem of abnormal vectors at the stitching point when the flow field stitching calculation between multiple cameras in PIV-2D2C large-field-of-view measurement. At the same time, when the calibration image plane and the laser plane are not completely parallel, the deviation between the calculated flow velocity and the actual flow velocity is improved, thereby improving the measurement accuracy of PIV. The method includes the following steps:

[0086] S1. Obtain a calibration image and a particle image located in the field of view of at least two cameras respectively through a two-dimensional PIV stitching measurement device. The specific process includes the following steps:

[0087] S11. Arrange a two-dimensional PIV splicing measurement device, including a water tank, multiple cameras, lasers, a synchronous controller, a PC, and a single-layer or double-layer calibration plate. Figure 2 As shown, a schematic diagram of using two cameras to arrange a two-dimensional PIV stitching measurement device is given. The camera used in this embodiment is a high-speed camera, and the laser is a continuous pulse laser. The water tank is a circulating flow device, and the fluid in the water tank can flow at a set speed. There is no need to deliberately control the parallelism of the camera, and an angle between the camera image plane and the laser plane is allowed. Multiple cameras are synchronized through a synchronous controller and then connected to the PC to ensure that images are collected at the same time. Then the water flow in the long water tank is photographed, and high-quality image data is obtained.

[0088] S12, turn on the laser, place the single-layer or double-layer calibration plate on the laser plane, make the calibration plate and the laser plane basically overlap, then turn off the laser, fine-tune the measurement field of view and depth of field of multiple cameras, and ensure that the two cameras to be spliced ​​can capture the calibration plate. Figure 3The camera synchronously collects and saves the calibration image. The collected calibration image is shown in Figure 7 shown.

[0089] S13, remove the calibration plate from the water tank, spread the tracer particles that can be illuminated by the laser in the water tank, let the fluid in the water tank circulate for a period of time, turn on the laser, and synchronously collect and save the particle image with the camera. The collected particle image is as follows: Figure 8 As shown, the particle images only show two frames.

[0090] S2. Calculate the calibration parameters used to describe the conversion relationship between the two-dimensional pixel points in the calibration image and the three-dimensional points in the world coordinate system. The calibration parameters include the coefficients [l1, l2, ..., l 12 ], and the coefficients of the second-order polynomial model [a1, a2, …, a 10 ] and [b1,b2,…,b 10 The specific process includes the following steps:

[0091] S21, based on the projection matrix P in the pinhole camera model, establish a linear transformation to convert the three-dimensional spatial coordinates (X, Y, Z) into undistorted pixel coordinates (u und ,v und )’s coordinate projection model, the expression of the coordinate projection model is:

[0092]

[0093]

[0094] Where P is the projection matrix; s is the scaling factor; [l1,l2,…,l 12 ] are the coefficients in the projection matrix P;

[0095] S22, obtain the detection pixel coordinates (u) of all calibration points in the calibration plate by means of dot center recognition and other methods. i ',v i ') and the corresponding three-dimensional spatial coordinates (X i ,Y i ,Z i ), and bring it into the coordinate projection model to solve the coefficients [l1,l2,…,l 12 ]’s initial value;

[0096] S23. Since the camera lens inevitably has distortion, a second-order polynomial is used here to describe the deviation caused by the distortion. A second-order polynomial model is established to describe the deviation of the camera distortion. The expression of the second-order polynomial model is:

[0097]

[0098]

[0099] Where (u, v) is the distorted pixel coordinate; (Δu, Δv) is the distorted pixel offset; [a1, a2, …, a 10 ]、[b1,b2,…,b 10 ] are the coefficients of the second-order polynomial model;

[0100] S24, the three-dimensional coordinates (X i ,Y i ,Z i ) and the projection matrix P are brought into the coordinate projection model to obtain the undistorted pixel coordinates (u iund ,v iund ), and use the detected pixel coordinates (u i ',v i ') minus the undistorted pixel coordinates (u iund ,v iund ) to obtain the distorted pixel offset (Δu i ,Δv i );

[0101] S25, offset the distorted pixel (Δu i ,Δv i ) is brought into the second-order polynomial model and the coefficients [a1, a2, …, a 10 ] and [b1,b2,…,b 10 ]’s initial value;

[0102] S26, the three-dimensional coordinates (X i ,Y i ,Z i ) is brought into the coordinate projection model and the second-order polynomial model to obtain the calculated pixel coordinates (u i ,v i );

[0103] S27, construct an optimization function to calculate the pixel coordinates (u i ,v i ) and the detected pixel coordinates (u i ',v i ') The optimization equation is established with the minimum deviation as the goal, that is:

[0104]

[0105] Where m is the number of pixels i;

[0106] S28, use the Levenberg-Marquardt iterative optimization algorithm to optimize the coefficients [l1,l2,…,l 12 ]、[a1,a2,…,a10 ] and [b1,b2,…,b 10 ] is optimized until the results converge and the calibration parameters are obtained.

[0107] In this embodiment, the coordinate mapping model is a coordinate projection model and a second-order polynomial model. By constructing a coordinate mapping model and optimizing the parameters in the model, the conversion relationship between the two-dimensional pixel points in the calibration image and the three-dimensional points in the world coordinate system is obtained, and the origin of the world coordinate system is made to coincide with the origin of the calibration plate, and the plane where the calibration plate is located is the laser plane of the Z=0 plane.

[0108] S3. Correct the calibration parameters to obtain the coordinate mapping model used by each camera to describe the conversion relationship between the spatial three-dimensional coordinates in the particle image and the distorted pixel coordinates (u, v). Since the laser surface and the calibration plate cannot be completely overlapped during calibration, there are Figure 3 The displacement deviation shown in the figure needs to be corrected for the calibration parameters obtained in the previous step. The specific process includes the following steps:

[0109] S31, calculate the maximum common area of ​​any two left and right cameras. Figure 5 As shown, in this embodiment, according to the coordinate transformation model, the spatial coordinates X and Y corresponding to the four vertices of the images in the left and right cameras in the Z=0 plane are calculated respectively, and the common area of ​​the two camera fields of view is obtained through 8 points. By traversing all possible rectangular positions in the area, it is checked whether each rectangle is completely in any irregular area, and the largest inscribed rectangle is found by judging the area of ​​the rectangle, which is the maximum common area of ​​the two cameras.

[0110] S32, mapping the particle images taken by the two cameras at the same time to the Z=0 plane respectively according to the calibration parameters, and obtaining a mapped particle image of the same size as the original image by interpolation of the maximum common area;

[0111] S33, the mapping particle images generated by the left and right cameras are divided into sub-areas and then cross-correlation calculation is performed to obtain the offset Δd of each sub-area of ​​the left and right cameras. i As a further explanation, theoretically, if the calibration plate coincides with the laser plane and the surface laser source only illuminates the tracer particles in the laser plane, the positions of the particles in the particle images mapped by the maximum common area of ​​the two cameras should be exactly the same, that is, the deviation vector is 0.

[0112] like Figure 6 As shown, the left camera image is divided into sub-areas according to a certain grid size, and the sub-area in the left camera image (i.e., the interpretation window in the figure) is slid in the vicinity of the corresponding position in the right camera (i.e., the search window in the figure) with a size of one pixel and the mutual correlation coefficient is calculated, and the distance between the position with the largest mutual correlation coefficient and the search window is the offset.

[0113] Among them, the calculation formula of the mutual correlation coefficient ZNCC can be expressed as:

[0114]

[0115] in, is the mean of sub-region A, that is:

[0116]

[0117] is the mean of sub-region B, that is:

[0118]

[0119] Where N is the number of sub-regions; A(i,j) and B(i,j) are pixels (i,j) in sub-regions A and B respectively.

[0120] S34, according to the offset Δd i Calculate the distance Δz between each vector and the corresponding point of each sub-area in the laser plane and the calibration plate plane i , distance Δz i The calculation formula is:

[0121]

[0122] In the formula, is the pixel u along the horizontal direction of camera n ni Partial derivatives with respect to X, Y, and Z; is the number of pixels v along the vertical direction of camera n ni Partial derivatives with respect to X, Y, and Z;

[0123] S35, refit a new Z=0 plane (the plane where the laser is located) using all the intersection points on the maximum common area, using the distance Δz of the non-zero parallax vector i Move and rotate the coordinate system so that the new Z=0 plane coincides with the actual laser plane;

[0124] S36, calculate the rotation and translation matrix of the corresponding points in the new Z=0 plane and the actual laser plane, and perform perspective transformation on the spatial coordinates of the calibration points in the calibration plate to obtain the spatial coordinates of the calibration points when the laser plane is taken as the Z=0 plane.

[0125] More specifically, after calculating the rotation and translation matrix, the perspective transformation of the spatial coordinates can be achieved, as shown in the following equation, x i ,y i ,z i It is expressed as the spatial coordinates of the calibration points, n is the number of calibration points, a0, a1, a2 are the plane rotation parameters, then:

[0126]

[0127] Let z i =Δz i , according to the above equations, the plane rotation parameters a0, a1, a2 are solved, and the rotation and translation matrix and the Z-direction offset can be expressed as:

[0128]

[0129]

[0130]

[0131] The spatial coordinates of the calibration points after perspective transformation can be expressed as:

[0132]

[0133] S37. Re-execute step S2 according to the spatial coordinates of the calibration points to obtain the corrected calibration parameters. In this embodiment, the calibration parameters are corrected by bringing the spatial coordinates of the calibration points into the calibration parameter calculation process given in step S2, which will not be described in detail here.

[0134] S38, bringing the corrected calibration parameters into the coordinate projection model and the second-order polynomial model to obtain a coordinate mapping model. Here, the corrected calibration parameters are brought into the coordinate projection model and the second-order polynomial model respectively, and a coordinate mapping model for each camera to describe the conversion relationship between the spatial three-dimensional coordinates in the particle image and the distorted pixel coordinates (u, v) can be obtained, thereby eliminating the defect that when the camera image plane is not completely parallel to the laser plane, the calculated flow velocity will deviate from the actual flow velocity.

[0135] S4. Calculate the homography matrix H from the pixel coordinates in any slave camera of at least two cameras to the pixel coordinates in the master camera, and the homography matrix H from any pixel coordinate in the master camera to the Z=0 plane based on the coordinate mapping model. w .

[0136] In this embodiment, N spatial points P are selected from the medium spacing according to the maximum common area between two adjacent cameras. n (X n ,Y n ,0), respectively, and the calibration parameters calculated above are introduced to obtain the spatial point P through the coordinate projection model. n (X n ,Y n ,0) The coordinates of the projected pixel point p in the main camera 1n (u 1n ,v 1n) and the projected pixel point in the slave camera, thereby calculating the homography matrix H from the point in the slave camera to the point in the main camera. Then calculate the projected pixel point p in the main camera in the previous step 1n (u 1n ,v 1n ) to the spatial point P in the Z = 0 plane (laser plane) n (X n ,Y n ,0) w The specific process includes the following steps:

[0137] S41, for each pair of corresponding points p from the pixel coordinates in the slave camera to the pixel coordinates in the master camera 1i (u 1i ,v 1i ) and p 2i (u 2i ,v 2i ), the constructed equation is:

[0138]

[0139] Where [h1,h2,…,h9] is the coefficient of the homography matrix H;

[0140] S42. Generate a 2n×9 matrix A through n>4 corresponding points, and solve the coefficients [h1,h2,…,h9] by the least squares method to obtain the homography matrix H. The least squares method solves overdetermined equations. Common software or algorithm libraries (matlab, opencv, eigen, etc.) have functions that can be directly calculated, so I will not go into details here. The expression of the homography matrix H is:

[0141]

[0142] S43, divide the homography matrix H by h9 to normalize and obtain the homography matrix H w .

[0143] S5. Use the classic PIV single-channel cross-correlation method to calculate the pixel displacement vector field of the particle image taken by each camera. The pixel displacement vector field includes the original coordinates g of the n grid points in the main camera grid. n (u n ,v n ), and the displacement Δg of the corresponding point is obtained by cross-correlation n (Δu n ,Δv n ). Here, the classic PIV single-channel cross-correlation calculation is used. This method is a conventional existing method and can be directly adopted by those skilled in the art, so it will not be described in detail here.

[0144] S6. In the pixel displacement vector field, determine the stitching data g of the slave camera in the master camera according to the homography matrix H. all_1n (u n ,v n ) and displacement Δg all_1n (Δu n ,Δv n ) to complete the splicing. The specific process includes the following steps:

[0145] S61, the original coordinates g of the point to be measured from the camera 2n (u 2n ,v 2n ) plus the displacement Δg 2n (Δu 2n ,Δv 2n ), the coordinates of the slave camera in the main camera are obtained by transforming it to the main camera coordinate system through the homography matrix H. 1n ′(u 1n ′,v 1n ′) and displacement Δg 1n ′(Δu 1n ′,Δv 1n ′);

[0146] The conversion formula is:

[0147]

[0148]

[0149] Where s is the scaling factor and ps is the scaling factor from the camera to the main camera.

[0150] S62, according to the original coordinates g in the main camera 1n (u 1n ,v 1n ) and the coordinates g transformed from the camera 1n ′(u 1n ′,v 1n ′) to determine the maximum area. This embodiment calculates the minimum value u of all coordinates in the horizontal direction min =min(u 11 ~u 1n ,u 11 ′,u 1n ′), the minimum value v in the vertical direction min =min(v 11 ~v 1n ,v 11 ′,v 1n ′);

[0151] The maximum value u in the horizontal direction max =max(u 11 ~u1n ,u 11 ′,u 1n ′), the maximum value v in the vertical direction max =max(v 11 ~v 1n ,v 11 ′,v 1n ′), where u min ~u max is the horizontal range of the maximum area, v min ~v max is the vertical range of the maximum area, and the maximum area is finally determined.

[0152] S63, expand the grid points divided in the main camera, and convert the coordinates g 1n ′(u 1n ′,v 1n ′) is interpolated into the grid of the main camera through triangulation, and finally the spliced ​​data g is obtained all_1n (u n ,v n ) and displacement Δg all_1n (Δu n ,Δv n ). In this embodiment, assuming that the original sub-area half-width is N, the sub-area center distance is p, and the horizontal coordinate of the leftmost sub-area of ​​the main camera is ul, if ul-pN>u min , then continue to expand to the left until ul-pN<u min Stop when the horizontal coordinate of the right sub-area is ur, if ur+p+N<u max , then continue to expand to the right until ur+p+N>u max The same applies to the vertical direction v.

[0153] In order to further understand the splicing proposed in this embodiment, the coordinates g 1n ′(u 1n ′,v 1n ′) displacement data Δg 1n ′(Δu 1n ′,Δv 1n ′) Interpolate to the grid closest to the main camera through triangulation to find the coordinates g 1n ′(u 1n ′,v 1n ′) are the three coordinates closest to the grid to be interpolated. Assuming that the coordinates of the three endpoints are a1(x1,y1), a2(x2,y2), and a3(x3,y3), and z1, z2, and z3 are the displacements of the points, the displacement z of the grid point is:

[0154]

[0155] S7, splice data g all_1n (u n ,v n ) and displacement Δg all_1n (Δu n ,Δv n ) scale transformation and offset calibration to obtain the physical coordinates G all_1n (x n ,y n ) and displacement ΔG all_1n (Δx n ,Δy n );

[0156] The scale transformation and offset calibration formula is:

[0157]

[0158]

[0159] In the formula, H w It is the homography matrix from any pixel coordinate in the main camera to the Z=0 plane (laser plane).

[0160] S8, use displacement ΔG all_1n (Δx n ,Δy n ) divided by the image acquisition interval T to obtain the velocity field. Fig. 9 , Fig.10 and Fig.11 As shown in the figure, by integrating the algorithm and importing it into the PIV analysis software Rflow for analysis, the analysis results show that the two images to be stitched displayed in the left camera and the right camera have significant effects, and there is no splitting phenomenon between the particles in the stitching seam area in the image. Fig.11 As shown in the figure, even when the camera image plane is not completely parallel to the laser plane, there is no obvious data discontinuity or abnormal data at the data splicing point. At the same time, the deviation between the calculated flow velocity and the actual flow velocity is eliminated, and the active calibration of the deviation is realized. Therefore, this method has strong versatility, and only needs to take an additional calibration image to complete the splicing calibration. It is simple to operate and has broad application prospects.

[0161] Those skilled in the art can understand that all or part of the steps in the above-mentioned embodiment method can be completed by instructing the relevant hardware through a program, so the present application can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present application can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program codes.

[0162] The above implementation methods have been described in detail. Specific examples are used herein to illustrate the principles and implementation methods of the present invention. The description of the above embodiments is only used to help understand the method of the present invention and its core idea. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as a limitation on the present invention.

Claims

1. A large-field-of-view two-dimensional flow field stitching calibration method for particle image velocimetry, characterized in that: The method comprises the following steps: S1, respectively obtaining a calibration image and a particle image located in the field of view of at least two cameras by a two-dimensional PIV stitching measurement device; S2. Calculate the calibration parameters used to describe the conversion relationship between the two-dimensional pixel points in the calibration image and the three-dimensional points in the world coordinate system. The calibration parameters include the coefficients [l1, l2, ..., l 12 ], and the coefficients of the second-order polynomial model [a1, a2, …, a 10 ] and [b1,b2,…,b 10 ]; S3, calibrating the calibration parameters to obtain a coordinate mapping model for each camera to describe the conversion relationship between the spatial three-dimensional coordinates in the particle image and the distorted pixel coordinates (u, v); S4. Calculate the homography matrix H from the pixel coordinates in any slave camera of at least two cameras to the pixel coordinates in the master camera, and the homography matrix H from any pixel coordinate in the master camera to the Z=0 plane based on the coordinate mapping model. w ; S5, calculate the pixel displacement vector field of the particle image taken by each camera, the pixel displacement vector field includes the original coordinates g of n grid points in the main camera grid n (u n ,v n ), and the cross-correlation gives the displacement Δg of the corresponding point n (Δu n ,Δv n ); S6. In the pixel displacement vector field, determine the stitching data g of the slave camera in the master camera according to the homography matrix H. all_1n (u n ,v n ) and displacement Δg all_1n (Δu n ,Δv n ) complete the splicing; S7, splice data g all_1n (u n ,v n ) and displacement Δg all_1n (Δu n ,Δv n ) scale transformation and offset calibration to obtain the physical coordinates G all_1n (x n ,y n ) and displacement ΔG all_1n (Δx n ,Δy n ); S8, use displacement ΔG all_1n (Δx n ,Δy n ) divided by the image acquisition interval T to obtain the velocity field.

2. The large-field-of-view two-dimensional flow field stitching calibration method according to claim 1, characterized in that: In step S1, the specific process includes the following steps: S11, arrange a two-dimensional PIV stitching measurement device, including a water tank, multiple cameras, a laser, a synchronous controller, a PC terminal, and a single-layer or double-layer calibration plate; S12, turn on the laser, place a single-layer or double-layer calibration plate on the laser plane, make the calibration plate and the laser plane basically overlap, then turn off the laser, fine-tune the measurement field of view and depth of field of multiple cameras, and the cameras synchronously collect and save calibration images; S13, remove the calibration plate from the water tank, spread the tracer particles that can be illuminated by the laser in the water tank, let the fluid in the water tank circulate for a period of time, turn on the laser, and synchronously collect and save the particle image with the camera.

3. The large-field-of-view two-dimensional flow field stitching calibration method according to claim 1, characterized in that: In step S2, the specific process includes the following steps: S21, based on the projection matrix P in the pinhole camera model, establish a linear transformation to convert the three-dimensional spatial coordinates (X, Y, Z) into undistorted pixel coordinates (u und ,v und )’s coordinate projection model, the expression of the coordinate projection model is: Where P is the projection matrix; s is the scaling factor; [l1,l2,…,l 12 ] are the coefficients in the projection matrix P; S22, obtain the detection pixel coordinates of all calibration points in the calibration plate (u i ',v i ') and the corresponding three-dimensional coordinates (X i ,Y i ,Z i ), and bring it into the coordinate projection model to solve the coefficients [l1,l2,…,l 12 ]’s initial value; S23, establishing a second-order polynomial model for describing the deviation of camera distortion, the expression of the second-order polynomial model is: Where (u, v) is the distorted pixel coordinate; (Δu, Δv) is the distorted pixel offset; [a1, a2, …, a 10 ], [b1,b2,…,b 10 ] are the coefficients of the second-order polynomial model; S24, the three-dimensional coordinates (X i ,Y i ,Z i ) and the projection matrix P are brought into the coordinate projection model to obtain the undistorted pixel coordinates (u iund ,v iund ), and use the detected pixel coordinates (u i ',v i ') minus the undistorted pixel coordinates (u iund ,v iund ) to obtain the distorted pixel offset (Δu i ,Δv i ); S25, offset the distorted pixel (Δu i ,Δv i ) is brought into the second-order polynomial model and the coefficients [a1, a2, …, a 10 ] and [b1,b2,…,b 10 ]’s initial value; S26, the three-dimensional coordinates (X i ,Y i ,Z i ) is brought into the coordinate projection model and the second-order polynomial model to obtain the calculated pixel coordinates (u i ,v i ); S27, construct an optimization function to calculate the pixel coordinates (u i ,v i ) and the detected pixel coordinates (u i ',v i ') The optimization equation is established with the minimum deviation as the goal, that is: Where m is the number of pixels i; S28, use the Levenberg-Marquardt iterative optimization algorithm to optimize the coefficients [l1,l2,…,l 12 ]、[a1,a2,…,a 10 ] and [b1,b2,…,b 10 ] is optimized until the results converge and the calibration parameters are obtained.

4. The large-field-of-view two-dimensional flow field stitching calibration method according to claim 3, characterized in that: In step S3, the specific process includes the following steps: S31, calculating the maximum common area of ​​any two left and right cameras; S32, mapping the particle images taken by the two cameras at the same time to the Z=0 plane respectively according to the calibration parameters, and obtaining a mapped particle image of the same size as the original image by interpolation of the maximum common area; S33, the mapping particle images generated by the left and right cameras are divided into sub-areas and then cross-correlation calculation is performed to obtain the offset Δd of each sub-area of ​​the left and right cameras. i ; S34, according to the offset Δd i Calculate the distance Δz between each vector and the corresponding point of each sub-area in the laser plane and the calibration plate plane i , distance Δz i The calculation formula is: In the formula, is the pixel u along the horizontal direction of camera n ni Partial derivatives with respect to X, Y, and Z; is the number of pixels v along the vertical direction of camera n ni Partial derivatives with respect to X, Y, and Z; S35, refitting a new Z=0 plane using all intersection points on the maximum common area, using the distance Δz of the non-zero disparity vector i Move and rotate the coordinate system so that the new Z=0 plane coincides with the actual laser plane; S36, calculating the rotation and translation matrix of the corresponding points in the new Z=0 plane and the actual laser plane, and performing perspective transformation on the spatial coordinates of the calibration points in the calibration plate to obtain the spatial coordinates of the calibration points when the laser plane is taken as the Z=0 plane; S37, re-execute step S2 according to the spatial coordinates of the calibration points to obtain the corrected calibration parameters; S38, bringing the corrected calibration parameters into the coordinate projection model and the second-order polynomial model to obtain a coordinate mapping model.

5. The large-field-of-view two-dimensional flow field stitching calibration method according to claim 4, characterized in that: In step S31, the maximum common area is the largest inscribed rectangle in the common area of ​​the two camera fields of view.

6. The large-field-of-view two-dimensional flow field stitching calibration method according to claim 4, characterized in that: In step S4, the specific process includes the following steps: S41, for each pair of corresponding points p from the pixel coordinates in the slave camera to the pixel coordinates in the master camera 1i (u 1i ,v 1i ) and p 2i (u 2i ,v 2i ), the constructed equation is: Where [h1,h2,…,h9] is the coefficient of the homography matrix H; S42. Generate a 2n×9 matrix A through n>4 corresponding points, and solve the coefficients [h1,h2,…,h9] according to the least squares method to obtain the homography matrix H. The expression of the homography matrix H is: S43, divide the homography matrix H by h9 to normalize and obtain the homography matrix H w .

7. The large-field-of-view two-dimensional flow field stitching calibration method according to claim 6, characterized in that: In step S6, the specific process includes the following steps: S61, the original coordinates g of the point to be measured from the camera 2n (u 2n ,v 2n ) plus the displacement Δg 2n (Δu 2n ,Δv 2n ), the coordinates of the slave camera in the main camera are obtained by transforming it to the main camera coordinate system through the homography matrix H. 1n ′(u 1n ′,v 1n ′) and displacement Δg 1n ′(Δu 1n ′,Δv 1n ′); S62, according to the original coordinates g in the main camera 1n (u 1n ,v 1n ) and the coordinates g transformed from the camera 1n ′(u 1n ′,v 1n ') determine the maximum area; S63, expand the grid points divided in the main camera, and convert the coordinates g 1n ′(u 1n ′,v 1n ′) is interpolated into the grid of the main camera through triangulation, and finally the spliced ​​data g is obtained all_1n (u n ,v n ) and displacement Δg all_1n (Δu n ,Δv n ).

8. The large-field-of-view two-dimensional flow field stitching calibration method according to claim 7, characterized in that: In step S61, the conversion formula is: Where s is the scaling factor and ps is the scaling factor from the camera to the main camera.

9. The large-field-of-view two-dimensional flow field stitching calibration method according to claim 8, characterized in that: In step S7, the scale transformation and offset calibration formula is: In the formula, H w It is the homography matrix from any pixel coordinate in the main camera to the Z=0 plane.

Images