Binocular vision three-dimensional reconstruction method based on timing consistency constraint

Abstract

This invention discloses a binocular vision-based 3D reconstruction method based on temporal consistency constraints, belonging to the interdisciplinary fields of computer vision and precision measurement. The method acquires intrinsic and extrinsic parameters through binocular camera arrangement and Zhang's calibration method to eliminate lens distortion; it simultaneously acquires dynamic video and performs grayscale normalization and background suppression. Sub-pixel-level depth recovery employs a gradient difference cost function, initial phase difference estimation, and parabolic fitting to improve pixel accuracy; it improves semi-global matching by fusing grayscale and gradient differences, combining 8-direction cost aggregation and multi-scale guided filtering to enhance accuracy in weakly textured regions. Temporal consistency constraints suppress jumps through disparity difference detection, Kalman filtering, and median filtering to ensure the continuity of the 3D trajectory. Finally, a dynamic point cloud is generated, which, verified by a laser displacement meter, shows reduced displacement and dominant frequency errors, making it suitable for vibration structure health monitoring and dynamic deformation analysis.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of computer vision and precision measurement, specifically relating to a binocular vision 3D reconstruction method based on temporal consistency constraints. Background Technology

[0002] In vibration monitoring of engineering structures, binocular cameras have gradually become an important measurement method due to their non-contact and full-field acquisition capabilities. The traditional binocular measurement process includes: camera calibration, feature point identification and matching, parallax calculation, and point cloud reconstruction. Through the principle of parallax, three-dimensional information of the structure in space can be obtained.

[0003] The most common implementation scheme currently available is a feature-matching-based binocular vision method. This method extracts natural textures or artificial targets as feature points from the left and right camera images, obtains depth information through matching and disparity calculation, and further reconstructs the point cloud to achieve 3D vibration measurement. While this type of method achieves good results under laboratory conditions, it still has shortcomings in practical engineering applications.

[0004] (1) When the surface texture is insufficient or the target is missing, the matching is prone to failure and the depth reconstruction accuracy decreases;

[0005] (2) Parallax calculation is limited to the pixel level. Digital errors accumulate gradually during the calculation, resulting in insufficient resolution in the depth direction.

[0006] (3) Frame-by-frame independent calculation lacks time constraints, and the point cloud trajectory is not stable enough in dynamic measurement.

[0007] Therefore, the existing solution closest to this invention is a binocular 3D measurement method based on feature matching and disparity calculation. Based on this, this invention proposes a binocular vision 3D reconstruction method based on temporal consistency constraints. Summary of the Invention

[0008] To overcome the shortcomings and deficiencies of the existing technology, the present invention adopts the following technical solution:

[0009] The binocular vision 3D reconstruction method based on temporal consistency constraints includes the following steps:

[0010] Step 1: Measurement system setup and parameter calibration;

[0011] Step 2: Dynamic Video Acquisition and Image Preprocessing

[0012] Step 3, Subpixel-level Depth Restoration: For each frame output from Step 2, perform the following sub-steps:

[0013] Step 4: Improve the semi-global matching (SGM) disparity calculation:

[0014] Step 5: Timing Consistency Constraint Correction:

[0015] Step 6: 3D point cloud reconstruction and measurement accuracy verification.

[0016] Preferably, step one is as follows:

[0017] S11. Binocular Camera Setup: Two identical industrial cameras are selected and fixed in front of the structure under test using a tripod. The camera angles are adjusted so that the object being measured is completely within the overlapping area of ​​the two cameras' fields of view. The lens direction is perpendicular to the plane of the structure under test, and the baseline length B is set. S12. Camera Intrinsic and Extrinsic Parameter Calibration: The Zhang calibration method is used. After collecting calibration data, the pixel coordinates of the corner points of the checkerboard calibration board and the world coordinates are extracted and matched using the Harris corner detection algorithm. The camera intrinsic and extrinsic parameters are solved, and distortion correction is performed on the subsequently acquired images based on the calibration parameters.

[0018] Preferably, the baseline length B is set according to the measurement distance; preferably, the structure to be measured is a vibrating thin plate.

[0019] Preferably, step two is as follows:

[0020] S21. Synchronous video acquisition: The binocular camera is triggered by the synchronization signal to synchronously acquire video of the dynamic process of the structure under test, ensuring that the time of each frame of the left and right cameras is strictly aligned; S22. Image preprocessing: The acquired images are subjected to grayscale normalization processing, the grayscale values ​​are mapped to the set range, and the frame difference method is used to calculate the grayscale difference between adjacent frames to filter static background areas, and the corrected left and right image sequences are output.

[0021] Preferably, step three is as follows:

[0022] S31. Local Gray-Level Difference: Set a matching window, calculate the horizontal and vertical gradients of corresponding pixels in the left and right images, and construct a gradient difference cost function to quantify the gradient similarity of corresponding pixels in the left and right images; S32. Phase Difference Estimation: Perform Fourier transform on the gray-level signals within the matching window, extract the phase of the left and right images and calculate the phase difference to initially estimate the sub-pixel offset; S33. Interpolation Fitting Correction: Using the integer disparity d corresponding to the initial phase difference estimate as the center, select three integer disparity points d-1, d, and d+1, substitute the cost function values ​​of the three into the parabolic equation to solve for the coefficients, calculate the optimal sub-pixel offset, and obtain the sub-pixel disparity; S34. Depth Calculation and Stability Control: Based on the triangulation principle, substitute the sub-pixel disparity into the depth formula to calculate the depth Z, and then apply a moving average filter to the sub-pixel disparity of the same pixel in adjacent frames to remove abnormal jump values ​​and output a stable depth value;

[0023] Preferably, in step S31, when calculating the horizontal and vertical gradients, if a pixel is located at the edge of the image, one-sided difference is used to avoid exceeding the boundary.

[0024] Preferably, step four is as follows: S41, Basic cost construction: fuse gray-level difference and gradient difference to construct a weighted cost function; S42, Multi-path cost aggregation: perform cost aggregation in multiple directions, and calculate the cumulative cost after aggregation along each direction using the cumulative cost formula; S43, Multi-scale guided filtering and edge constraint: downsample the corrected image to at least two scales, use the original image as a guide map to divide the region into regions, apply different window filtering, upsample and weighted fuse, and then add a cross-edge penalty term through the edge indicator function to obtain the final aggregated cost; S44, Disparity optimization and hole filling: mark the wrong matching points by checking the left and right consistency, and use bilinear interpolation combined with the temporal information of the previous frame to compensate for the wrong matching points and occluded area holes to generate an optimized disparity map.

[0025] Preferably, step five is as follows: S51, disparity difference detection: calculate the disparity difference between adjacent frames, set a difference threshold T, mark pixels with a disparity difference greater than T as inconsistent regions, and those less than or equal to T as consistent regions; S52, temporal smoothing correction: apply multi-frame temporal median filtering to small-scale inconsistent regions, construct a Kalman filtering model for large-scale inconsistent regions, and output the filtered and corrected disparity; S53, trajectory consistency protection: substitute the corrected disparity into the trigonometric formula to obtain spatial points, calculate the distance between spatial points in adjacent frames, determine trajectory abrupt changes, construct a linear prediction model based on multi-frame trajectories, replace abnormal spatial points with predicted values, and output a stable three-dimensional trajectory.

[0026] Preferably, in step S53, the condition for determining a sudden change in trajectory is: calculating the spatial distance between adjacent frames. ,like The average displacement over the past 5 frames was determined to be a sudden change in trajectory; where ΔP is the distance between spatial points in adjacent frames. Let P(t) be the Euclidean norm; P(t) be the coordinates of a point in three-dimensional space at frame t; P(t-1) be the coordinates of a point in three-dimensional space at frame (t-1); preferably, the linear prediction model is constructed based on the trajectory of the first 5 frames, and the model formula is: Replace the outlier P(t) with the predicted value to output a stable three-dimensional trajectory. ;in, Let be the predicted 3D spatial coordinates of the t-th frame, a be the slope of the prediction model, and b be the intercept of the prediction model.

[0027] Preferably, step six is ​​as follows:

[0028] S61. 3D point cloud generation: Based on the stable parallax output in step five and the camera parameters calibrated in step one, triangulation is performed on each frame of image to generate a dynamic 3D point cloud sequence of the structure under test; S62. Vibration time history extraction and verification: Key measurement points are selected in the point cloud, and the displacement time history curves of their vibration directions are extracted. A laser displacement meter is used as a benchmark to compare the camera measurement and laser measurement results to verify the measurement accuracy.

[0029] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are:

[0030] 1. This invention achieves sub-pixel parallax accuracy of 0.01~0.1 pixels by using gradient difference cost function + initial phase difference estimation + parabolic fitting. It combines multi-scale guided filtering and edge constraint optimization to match weak texture regions, and Zhang's calibration to eliminate lens distortion. Finally, it compares and verifies the reduction of displacement error and main frequency error with laser displacement meter, thus meeting the high-precision requirements of vibration structure deformation analysis.

[0031] 2. This invention employs a three-level temporal constraint mechanism of disparity difference detection + Kalman filtering + median filtering to effectively suppress disparity jumps of >0.5 pixels in dynamic vibration scenarios; it corrects trajectory mutation points by constructing a linear prediction model based on the trajectory of the first 5 frames, ensuring that the distance change between spatial points in adjacent frames is ≤ the threshold, thereby achieving continuous and stable output of three-dimensional trajectory, which is suitable for health monitoring of dynamic structures such as bridges and thin plates. Attached Figure Description

[0032] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0033] Figure 1 A flowchart of the binocular vision 3D reconstruction method based on temporal consistency constraints of the present invention is shown;

[0034] Figure 2 A flowchart of the improved semi-global matching (SGM) disparity calculation method of the present invention is shown. Detailed Implementation

[0035] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0036] Furthermore, the described features, structures, or characteristics can be combined in any suitable manner in one or more exemplary embodiments. Numerous specific details are provided in the following description to give a full understanding of exemplary embodiments of this disclosure. However, those skilled in the art will recognize that the technical solutions of this disclosure can be practiced with one or more of the specific details omitted, or other methods, components, steps, etc., can be employed. In other instances, well-known structures, methods, implementations, or operations are not shown or described in detail to avoid obscuring various aspects of this disclosure.

[0037] Example 1:

[0038] See Figure 1 As shown in the figure, the binocular vision 3D reconstruction method based on temporal consistency constraints in this embodiment has the following process:

[0039] Step 1: Measurement system setup and parameter calibration (basic preparation stage).

[0040] S11, Binocular camera setup.

[0041] Two industrial cameras of the same model (such as HIKMV-CS028-10UM, focal length 12mm, frame rate 100 frames / second, 8-bit depth) were selected and fixed in front of the structure to be tested (such as a vibrating thin plate) using a tripod.

[0042] Adjust the camera angle to ensure that the object being measured is completely within the overlapping area of ​​the two cameras' fields of view, with the lens direction perpendicular to the plane of the structure to be measured (such as the plane of a thin plate). The baseline length B is set according to the measurement distance (example: 150.2 mm).

[0043] S12. Camera internal and external parameter calibration (Zhang's calibration method)

[0044] Collect calibration data: Place the checkerboard calibration board (5mm×5mm, 11×9 grids) in front of the camera and take 60 sets of calibration images at different distances (0.5~2m) and different angles (-30°~30°).

[0045] Corner extraction and matching: The Harris corner detection algorithm is used to extract the checkerboard corners of each calibration image, and the pixel coordinates and world coordinates of the corners are matched.

[0046] Parameter calculation: The camera's intrinsic parameters (focal length f, principal point coordinates, radial / tangential distortion coefficients) and extrinsic parameters (rotation matrix R, translation vector T) are solved using Zhang's calibration algorithm; Distortion correction: Based on the calibration parameters, all subsequently acquired images are preprocessed to eliminate the influence of lens radial / tangential distortion.

[0047] Step 2: Dynamic Video Acquisition and Image Preprocessing (Data Input Stage)

[0048] S21, Synchronous video capture.

[0049] Trigger control: The binocular camera is triggered by a synchronization signal to synchronously acquire video of the dynamic process of the structure under test (example: video of the vibration of a thin plate excited by a small hammer, 10 seconds long, 1000 frames in total), ensuring that the time of each frame of the left and right cameras is strictly aligned (no timing misalignment).

[0050] S22, Image preprocessing.

[0051] Gray-level normalization: Maps the image gray-level values ​​from [0,255] to the [0,1] range to eliminate the interference of light intensity differences on the gray-level signal;

[0052] Background suppression: The frame difference method is used to calculate the gray-level difference between adjacent frames, filter static background areas, preserve the dynamic texture features of the structure under test, and output the corrected left and right image sequences. , (t=1~1000, t is the frame number).

[0053] Step 3: Subpixel-level depth recovery (high-precision depth estimation stage).

[0054] For each frame output from step two, calculate the subpixel depth using the following sub-steps:

[0055] S31, Local gray-level difference (cost function construction).

[0056] Set a matching window (e.g., 5×5 or 7×7) and calculate the horizontal gradient of corresponding pixels in the left and right images. Vertical gradient ;

[0057] Construct the gradient difference cost function: ;in, The gradient difference cost is used to quantify the gradient similarity of corresponding pixels in the left and right images; the smaller the value, the higher the matching degree. d is an integer disparity candidate value (unit: pixels), representing the horizontal offset of the right image pixel relative to the left image pixel. These are the image pixel coordinates (unit: pixels), corresponding to the row and column positions in the image plane; The horizontal grayscale gradient (unit: grayscale value / pixel) at pixel (x,y) in the left image reflects the rate of change of grayscale in the horizontal direction; The horizontal grayscale gradient (unit: grayscale value / pixel) at pixel (xd,y) in the right image is a candidate matching pair with the pixel (x,y) in the left image; The vertical grayscale gradient (unit: grayscale value / pixel) at pixel (x,y) in the left image reflects the rate of change of grayscale in the vertical direction. The vertical grayscale gradient (unit: grayscale value / pixel) at pixel (xd,y) in the right image is a candidate matching pair with pixel (x,y) in the left image.

[0058] The horizontal gradient in this method Vertical gradient The calculation is performed using the simple difference method, and the specific formula is as follows, directly related to the grayscale value I:

[0059] 1. Horizontal gradient Calculation (with left figure) (For example).

[0060] The horizontal gradient describes the grayscale change of a pixel in the horizontal direction (x-axis direction), and is calculated using the following formula: ;

[0061] Meaning: The horizontal gradient of pixel (x, y) in the left image is equal to the grayscale value of the pixel (x+1, y) to its right. Subtract the gray value of the left pixel (x-1, y). ;

[0062] If a pixel is located at the edge of the image (e.g., x=0 or x=image width-1), then one-sided difference is used (e.g., ...). ), to avoid crossing the boundary.

[0063] Similarly, the horizontal gradient of the candidate matching pixel (xd, y) in the right figure is as follows: .

[0064] 2. Vertical gradient Calculation (with left figure) (For example).

[0065] The vertical gradient describes the grayscale change of a pixel in the vertical direction (y-axis direction), and its calculation formula is: ;

[0066] Meaning: The vertical gradient of pixel (x, y) in the left image is equal to the gray value of the pixel (x, y+1) below it. Subtract the gray value of the upper pixel (x, y-1) If a pixel is located at the edge of the image (e.g., y=0 or y=image height-1), then one-sided difference is used (e.g., ...). Similarly, the vertical gradient of the candidate matching pixel (xd, y) in the right figure is as follows: .

[0067] in, The grayscale value at pixel (x,y) in the left image (range [0,1] after preprocessing); The grayscale value at pixel (xd,y) in the right image (range [0,1] after preprocessing).

[0068] S32, Phase difference estimation (initial judgment of sub-pixel displacement).

[0069] Perform a Fourier transform on the grayscale signal within the matching window to extract the phase of the left image. Phase with the right figure Calculate the phase difference ;

[0070] Preliminary estimate of subpixel offset: The value range is [-0.5, 0.5] pixels.

[0071] S33, Interpolation Fitting Correction (Subpixel Parallax Calculation).

[0072] Determine integer disparity candidates: Using the integer disparity d corresponding to the initial phase difference estimate as the center, select three integer disparity points: d-1, d, and d+1.

[0073] Parabolic fitting: Substitute the cost function values ​​C(d-1), C(d), and C(d+1) of the three disparity points into the parabolic equation. Solve for the coefficients a, b, and c;

[0074] Calculate the optimal sub-pixel offset: Ultimately, subpixel parallax (Precision 0.01~0.1 pixels).

[0075] S34, Depth Calculation and Stability Control.

[0076] Depth Conversion: Based on the principle of triangulation, subpixel parallax is substituted into the formula to calculate depth Z (along the camera's optical axis): Where Z is the depth of the point to be measured along the optical axis of the camera (unit: mm, determined according to the actual measurement distance), that is, the vertical distance from the point to be measured to the camera plane in three-dimensional space; f is the camera focal length (unit: mm), an intrinsic parameter obtained by Zhang's calibration algorithm in step S12, with an example value of 12.05 mm; B is the baseline length of the binocular camera (unit: mm), an extrinsic parameter obtained by calibration in step S12 (distance between the optical centers of the two cameras), with an example value of 150.2 mm; The subpixel corrected disparity (unit: pixels) is obtained by parabolic fitting in step S33, with an accuracy range of 0.01~0.1 pixels.

[0077] Temporal filtering: subpixel disparity of the same pixel in adjacent frames , A moving average filter is used to remove abnormal jump values ​​and output a stable depth value. .

[0078] Step 4: Improve the semi-global matching (SGM) disparity calculation (disparity optimization stage)

[0079] See Figure 2 As shown, the improved semi-global matching (SGM) disparity calculation is as follows:

[0080] Based on the subpixel parallax in step 3, further optimize the parallax accuracy in areas with weak texture / low contrast:

[0081] S41. Basic cost construction: Fusing grayscale difference and gradient difference to construct a weighted cost function: ;in, The basic matching cost is a comprehensive cost that combines grayscale difference and gradient difference; the smaller the value, the higher the matching reliability. This is a weighting factor (unitless) used to balance the contributions of gray-scale difference and gradient difference, with a value range of 0.3 to 0.5, and an example value of 0.4. The gray-level similarity is quantized by the absolute difference (unitless) of the gray levels of corresponding pixels in the left and right images; other parameters (d, x, y, ...) , , The meaning is the same as step S31.

[0082] S42, Multi-path cost aggregation.

[0083] Cost aggregation is performed in eight directions (left→right, right→left, up→down, down→up, and four diagonals). The aggregation formula is as follows: ;in, This represents the cumulative cost of pixel (x,y) under disparity d after aggregation along direction r; r represents the cost aggregation direction (8 in total: left→right, right→left, up→down, down→up, top left→bottom right, bottom right→top left, top right→bottom left, bottom right→bottom left, bottom left→top right). , This represents the step size (unitless) in the direction r, such as left → right. , Up → Down , ; This represents other disparity candidate values ​​(in pixels), used to iterate through all possible disparities to find the minimum cumulative cost; This represents a small penalty parameter (unitless), used to penalize cases where the disparity difference between adjacent pixels is 1. This represents a large penalty parameter (unitless), used to penalize cases where the disparity difference between adjacent pixels exceeds 1 (suppressing mismatches across edges). Indicates an indicator function (unitless); Indicates when The value is 1 if the condition is met, otherwise it is 0. Indicates when The parameter is set to 1 if the condition is met, and 0 otherwise; the other parameters ( The meanings of x, y, and d are the same as before.

[0084] S43, Multi-scale Guided Filtering and Edge Constraints.

[0085] Multi-scale filtering: The corrected image is downsampled to two scales, 1 / 2 and 1 / 4. Using the original image as a guide image, low-contrast areas are filtered using a 9×9 window (to enhance texture), and edge areas are filtered using a 3×3 window (to preserve details). After upsampling, weighted fusion is performed to obtain the aggregate cost. Edge Preservation: Edge pixels are extracted using the Canny operator, and a cross-edge penalty term is added to the cost update. : ;in, This represents the final aggregation cost after adding edge constraints; This represents the aggregation cost after multi-scale guided filtering fusion; This represents the edge penalty parameter (unitless), used to enhance the matching reliability of edge regions; The edge indicator function (unitless) takes the value 1 when pixel (x,y) is an edge pixel (detected by the Canny operator) and 0 when it is a non-edge pixel; the other parameters (d, x, y) have the same meaning as before.

[0086] S44, Parallax Optimization and Hole Filling.

[0087] Left-right consistency check: parallax of the left image Parallax compared to the right figure ,like Marked as an incorrect match point;

[0088] Hole handling: For holes in mismatched points and occluded areas, bilinear interpolation combined with temporal information from the previous frame is used for compensation to generate a complete optimized disparity map. .

[0089] Step 5: Timing Consistency Constraint Correction (Dynamic Trajectory Stabilization Stage).

[0090] For dynamic measurement scenarios, ensure continuous and stable 3D trajectory:

[0091] S51, disparity difference detection.

[0092] Calculate the disparity difference between adjacent frames: ;

[0093] Threshold determination: Set a difference threshold T = 0.5 pixels. Pixels with Δd(t) > T are marked as inconsistent regions. The pixels retain t-frame parallax (marked as consistent regions).

[0094] S52, Timing Smoothing Correction.

[0095] Small inconsistency areas (≤3 inconsistent pixels within a 3×3 window): Use 3-frame temporal median filtering: ;

[0096] in, This represents the time-corrected parallax (in pixels), used to eliminate jumps in small, inconsistent areas. The median function (unitless) is used to calculate the median value of the input sequence and suppress outliers; t is the frame number (unitless), indicating the current video frame being processed. This represents the optimized disparity of frame t (in pixels, output from step 4.4). , These represent the optimized disparity (in pixels) for the (t-1)th and (t-2)th frames, respectively.

[0097] Large-scale inconsistent regions (e.g., areas of intense vibration): Construct a Kalman filter model: State equations: Where d(t) is the disparity state quantity of frame t (unit: pixels), representing the estimated value of the true disparity; d(t-1) is the disparity state quantity of frame t-1 (unit: pixels); u(t-1) is the disparity change rate of frame t-1 (unit: pixels / frame), reflecting the trend of disparity change over time; w(t-1) is the process noise (unit: pixels), which follows a Gaussian distribution and is used to describe the uncertainty of the model.

[0098] Observation equation: ; (v is the observation noise); where z(t) is the disparity observation value of the t-th frame (unit: pixels); Let v(t) be the optimized disparity of the t-th frame (in pixels, i.e., the observation data); v(t) be the observation noise (in pixels), which follows a Gaussian distribution and is used to describe the observation error (such as matching noise).

[0099] Output filtering and parallax correction As the final correction for parallax .

[0100] S53, Trajectory Consistency Protection.

[0101] 3D coordinate mapping: Substituting into the trigonometric formulas, we obtain the spatial points. ;

[0102] Mutation detection: Calculate the spatial distance between adjacent frames ,like The average displacement of the past 5 frames is used to determine if there is a sudden change in the trajectory; where ΔP is the distance between spatial points in adjacent frames (unit: mm), used to determine whether there is a sudden change in the trajectory. is the Euclidean norm (unitless), used to calculate the straight-line distance between two points in three-dimensional space; P(t) is the coordinate of a point in three-dimensional space in frame t (unit: mm), i.e. ,Depend on Substituting into the trigonometric formula, we get: P(t-1) is the coordinate of a three-dimensional point in the (t-1)th frame (unit: mm).

[0103] Prediction compensation: Constructing a linear prediction model based on the trajectory of the first 5 frames Replace the outlier P(t) with the predicted value to output a stable three-dimensional trajectory. .in, denoted as (t) is the predicted 3D spatial coordinates of the t-th frame (in mm), used to correct trajectory abrupt changes; a is the slope of the prediction model (in mm / frame), obtained by fitting the trajectory of the first 5 frames, reflecting the motion trend of the spatial points; b is the intercept of the prediction model (in mm), obtained by fitting the trajectory of the first 5 frames; the meanings of the remaining parameters (t) are the same as before.

[0104] Step 6: 3D point cloud reconstruction and measurement accuracy verification (result output stage).

[0105] S61, 3D point cloud generation.

[0106] Based on the stable parallax output of step five Using the camera parameters calibrated in step one, triangulation is performed on each frame of the image to generate a dynamic three-dimensional point cloud sequence of the structure under test.

[0107] S62. Vibration time history extraction and verification.

[0108] Extract the trajectory of the target point: Select key measuring points (such as two points at the bottom of the thin plate) in the point cloud and extract their displacement time history curves in the Z direction (vibration direction);

[0109] Accuracy Comparison: Using a laser displacement meter (such as Keyence LK-H150, with an accuracy of 0.1μm) as a benchmark, the results of camera measurement and laser measurement are compared:

[0110] Clock frequency consistency: The clock frequency difference between the two is ≤1% (example: both are 12.5Hz);

[0111] Displacement error: The difference in displacement between the two is ≤3%, verifying the subpixel-level measurement accuracy.

[0112] The beneficial effects of this embodiment are as follows: high-precision dynamic 3D reconstruction is achieved by improving parallax calculation and temporal consistency constraints. Subpixel-level depth recovery combined with multi-scale guided filtering improves the matching reliability of weak texture regions; temporal filtering and Kalman filtering suppress parallax jumps and ensure dynamic trajectory stability. Binocular system calibration and distortion correction eliminate lens errors, and synchronous video acquisition ensures temporal alignment, meeting the requirements of high-precision dynamic measurement and suitable for vibration structure health monitoring and deformation analysis.

[0113] All formulas in this invention are dimensionless and calculated numerically. The preset parameters in the formulas can be set by those skilled in the art according to the actual situation.

[0114] The weighting coefficients of this invention are used to measure the degree of influence of different factors or variables on a certain outcome or decision. The weighting coefficient is defined as the numerical value assigned to each factor when comparing and evaluating multiple factors, reflecting their importance or priority. These weighting coefficients can be determined according to specific circumstances and needs, and are usually jointly formulated and confirmed by professionals or relevant stakeholders. By reasonably setting the weighting coefficients, programs or systems can be helped to make decisions or predictions more accurately.

[0115] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

[0116] The preferred embodiments of the present invention disclosed above are merely illustrative of the invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the invention to specific implementations. Clearly, many modifications and variations can be made based on the content of this specification. This specification selects and specifically describes these embodiments to better explain the principles and practical applications of the invention, thereby enabling those skilled in the art to better understand and utilize the invention. The invention is limited only by the claims and their full scope and equivalents.

Claims

1. A binocular vision three-dimensional reconstruction method based on timing consistency constraints, characterized in that, The method comprises the following steps: Step one, measurement system setup and parameter calibration; Step two, dynamic video acquisition and image preprocessing; Step three, sub-pixel level depth recovery: based on the image sequence preprocessed in step two, sub-pixel parallax is obtained through gradient calculation, phase difference estimation and interpolation fitting, and stable depth value is obtained by combining sliding average filtering, to provide initial depth reference for subsequent improved semi-global matching (SGM) parallax calculation; Step four, improved semi-global matching (SGM) disparity calculation: based on the sub-pixel disparity obtained in step three, a weighted cost function is constructed and the disparity is optimized through multi-path aggregation, multi-scale guided filtering and edge constraint, while filling the occlusion area holes to generate an optimized disparity map; specifically: the weighted cost function ; wherein, is a basic matching cost, which is a comprehensive cost of the gray difference and the gradient difference, and the smaller the value is, the higher the matching reliability is; is a weight factor, which is used to balance the contribution of the gray difference and the gradient difference; is an absolute gray difference of corresponding pixels of left and right images, which quantifies the gray similarity; is a gradient difference cost function; the disparity is optimized through multi-direction cost aggregation, multi-scale guided filtering and edge constraint, the multi-direction cost aggregation includes cost accumulation along eight directions of left→right, right→left, up→down, down→up and four diagonal directions, while filling the occlusion area holes to generate an optimized disparity map; Step five, time sequence consistency constraint correction: based on the optimized parallax map in step four, frame inconsistency area is located by parallax difference detection, and abnormal parallax and spatial points are corrected by combining time sequence filtering and trajectory consistency protection, to output stable parallax and three-dimensional trajectory; Step six, three-dimensional point cloud reconstruction and measurement accuracy verification: combined with the calibration parameters in step one and the stable parallax in step five, dynamic three-dimensional point cloud sequence of the structure to be measured is generated by triangulation, and the displacement time course of key measuring points is extracted and compared with the laser displacement meter as a reference device, to verify the measurement accuracy of the method and complete the whole three-dimensional reconstruction process.

2. The binocular vision three-dimensional reconstruction method based on timing consistency constraints according to claim 1, characterized in that, Step one is specifically as follows: S11, binocular camera arrangement: two industrial cameras of the same type are selected, which are fixed in front of the structure to be measured by a tripod, the camera viewing angle is adjusted so that the measurement object is completely in the overlapping area of the viewing fields of the two cameras, the lens direction is perpendicular to the plane of the structure to be measured, and the baseline length B is set; S12, camera internal and external parameter calibration: Zhang's calibration method is adopted, the pixel coordinates and world coordinates of the corner points of the checkerboard calibration plate are extracted and matched by the Harris corner point detection algorithm after collecting the calibration data, the camera internal and external parameters are solved, and the subsequent collected images are corrected for distortion based on the calibration parameters; 3. The binocular vision 3D reconstruction method based on timing consistency constraint according to claim 1, characterized in that, The baseline length B is set according to the measurement distance; and the structure to be measured is a vibrating thin plate. Step two is specifically as follows:

4. The binocular vision three-dimensional reconstruction method based on timing consistency constraints according to claim 1, characterized in that, S21, synchronous video acquisition: the binocular cameras are triggered by a synchronization signal to synchronously collect the video of the dynamic process of the structure to be measured, so that the time of each frame of images of the left and right cameras is strictly aligned; S22, image preprocessing: the collected images are subjected to gray scale normalization processing, the gray scale value is mapped to a set interval, and then the frame difference method is used to calculate the gray scale difference between adjacent frames to filter the static background area, and the corrected left and right image sequences are output. Step three is specifically as follows: S31, local gray difference: a matching window is set, the horizontal gradient and the vertical gradient of the corresponding pixels of the left and right images are calculated, a gradient difference cost function is constructed to quantify the gradient similarity of the corresponding pixels of the left and right images; S32, phase difference estimation: the gray scale signals in the matching window are subjected to Fourier transform, the phase of the left image and the phase of the right image are extracted, and the phase difference is calculated to preliminarily estimate the sub-pixel offset; S33, interpolation fitting correction: taking the integer parallax d corresponding to the preliminary estimation of the phase difference as the center, selecting three integer parallax points of d-1, d and d+1, substituting the cost function values of the three points into a parabolic equation to solve the coefficients, calculating the optimal sub-pixel offset, and obtaining the sub-pixel parallax; S34, depth calculation and stability control: based on the principle of triangulation, the sub-pixel parallax is substituted into the depth formula to calculate the depth Z, and the sliding average filtering is used on the sub-pixel parallax of the same pixel in adjacent frames to remove abnormal jump values, and the stable depth value is output; In the step S31, when calculating the horizontal gradient and the vertical gradient, if the pixel is located at the edge of the image, one-sided difference is used to avoid border crossing.

5. The binocular vision 3D reconstruction method based on timing consistency constraint according to claim 1, characterized in that, The step four is specifically as follows: S41, base cost construction: fusing the gray difference and the gradient difference to construct a weighted cost function; S42, multi-path cost aggregation: using multi-direction to aggregate the cost, and calculating the accumulated cost after aggregation along each direction through an accumulated cost formula; S43, multi-scale guided filtering and edge constraint: down-sampling the corrected image into at least two scales, using different window filtering on the original image as a guide image to be regionally up-sampled and fused, and then increasing a cross-edge penalty term through an edge indicator function to obtain the final aggregated cost; S44, disparity optimization and hole filling: marking the error matching points through left-right consistency checking, and using bilinear interpolation combined with the time sequence information of the previous frame to compensate the error matching points and the occlusion area holes to generate an optimized disparity map.

6. The binocular vision 3D reconstruction method based on timing consistency constraint according to claim 1, characterized in that, The step five is specifically as follows: S51, disparity difference detection: calculating the disparity difference of adjacent frames, setting a difference threshold T, marking the pixels with disparity difference greater than T as inconsistent regions, and marking the pixels with disparity difference less than or equal to T as consistent regions; S52, time sequence smoothing correction: using multi-frame time sequence median filtering on small-range inconsistent regions, and constructing a Kalman filtering model on large-range inconsistent regions to output the filtered and corrected disparity; S53, trajectory consistency protection: substituting the corrected disparity into a triangular formula to obtain a space point, calculating the distance between the space points of adjacent frames, determining a linear prediction model based on multiple frames of trajectory after trajectory mutation, replacing the abnormal space points with the predicted values, and outputting stable three-dimensional trajectories.

7. The binocular vision three-dimensional reconstruction method based on timing consistency constraints according to claim 6, characterized in that, In step S53, the condition for determining trajectory mutation is that the distance between adjacent frame spatial points is calculated , and if is greater than a set value, it is determined that the trajectory is mutated; wherein is the distance between adjacent frame spatial points; is the Euclidean norm; is the three-dimensional spatial point coordinate of the t-th frame; P(t-1) is the three-dimensional spatial point coordinate of the t-1-th frame; the linear prediction model is constructed based on the previous 5 frames of trajectory, and the model formula is , the abnormal P(t) is replaced by the predicted value, and a stable three-dimensional trajectory is output; wherein is the predicted three-dimensional spatial point coordinate of the t-th frame, a is the slope of the prediction model, and b is the intercept of the prediction model.

8. The binocular vision three-dimensional reconstruction method based on timing consistency constraints according to claim 1, characterized in that, The step six is specifically as follows: S61, three-dimensional point cloud generation: based on the stable disparity output by the step five and the camera parameters calibrated by the step one, performing triangular measurement on each frame of image to generate a dynamic three-dimensional point cloud sequence of the structure to be measured; S62, vibration time course extraction and verification: selecting key measuring points in the point cloud, extracting the displacement time course curve of the vibration direction, and comparing the camera measurement result with the laser measurement result by using a laser displacement meter as a reference to verify the measurement accuracy.

Images