Spatial rotation axis calibration method and system for 3D level measurement

Abstract

The present invention provides a spatial rotation axis calibration method and system for 3D object level measurement, which relates to the field of measurement technology, including arranging multiple high-resolution industrial cameras to ensure that at least two cameras simultaneously observe at least one characteristic plane of the polyhedron calibration target during the rotation of the rotation axis; driving the rotation axis to rotate one circle at a non-uniform angle interval, the rotation interval is determined according to the golden section ratio, and a multi-view image sequence of the polyhedron calibration target is collected at each angle position; using the multi-view geometry principle and camera array parameters, the two-dimensional image coordinates are reconstructed into high-precision three-dimensional space coordinates; based on the reconstructed three-dimensional space coordinates, a feature point trajectory fitting optimization model considering gravity deformation and temperature drift is constructed, and a dynamically compensated rotation axis coordinate system is established using the optimal rotation axis parameters; in actual 3D object level measurement, according to the current angle of the rotation axis and environmental parameters, the coordinate transformation is corrected in real time using the comprehensive error compensation model.

Description

Technical Field

[0001] The present invention relates to measurement technology, and in particular to a spatial rotation axis calibration method and system for 3D object level measurement. Background Art

[0002] In the measurement of three-dimensional material level and stockpile (such as coal pile, ore, etc.), in order to measure the complete shape of the material and accurately calculate the volume, it is necessary to shoot point clouds from multiple angles (using laser scanning or radar scanning). In some cases, the material can be placed on a large turntable for shooting, and then the point clouds shot from multiple angles can be spliced. When splicing, it is necessary to first determine the axis position and inclination angle of the large turntable, that is, the axis calibration, that is, to determine the spatial linear equation of the axis.

[0003] He Wantao, Xu Yongan, Lv Haidong, etc. used a method based on multi-angle (>=3) spatial coordinates of a same-name point rotating around the axis of rotation to fit the circle equation of the spatial trajectory of the same-name point and solve the linear equation of the axis of rotation. However, if the measurement of the spatial position of the same-name point is inaccurate, the calculated axis of rotation is also inaccurate, and this error is proportional to the measurement error of the spatial point position. Zhang Xiaochen, Ye Hao, Park SY, etc. gave a method based on multiple same-name points, which can average out this measurement error and solve the axis of rotation equation with higher accuracy, but the matching of multiple same-name points is a problem. In some cases, same-name points cannot be matched. For example, we can only measure the equation of a spatial plane (such as plane point cloud fitting), and the plane rotates around the axis of rotation to obtain multiple spatial plane equations.

[0004] In fact, based on multiple spatial plane equations rotating around the axis, the linear equation of the axis can also be solved. The plane equation can be obtained by placing a plane plate on the turntable, obtaining the plane point cloud, and then fitting the point cloud plane. Xu Pei and YEY gave a method for calibrating the axis of rotation using a plane placed vertically on the turntable, but the plane is required to be strictly perpendicular to the turntable plane, which is physically difficult to achieve. Summary of the invention

[0005] The embodiments of the present invention provide a spatial rotation axis calibration method and system for 3D object level measurement, which can solve the problems in the prior art.

[0006] According to a first aspect of the embodiments of the present invention,

[0007] Provides a spatial rotation axis calibration method for 3D level measurement, including:

[0008] A polyhedron calibration target is installed on the rotation axis to be calibrated, wherein the polyhedron calibration target includes at least four feature planes that are not coplanarly distributed, and each feature plane is provided with at least three coded feature points that are not colinearly distributed; multiple high-resolution industrial cameras are arranged to ensure that at least two cameras simultaneously observe at least one feature plane of the polyhedron calibration target during the rotation of the rotation axis; the intrinsic and extrinsic parameters of the camera array are calibrated to establish a global coordinate system of the camera array; the rotation axis is driven to rotate one circle at a non-uniform angle interval, and the rotation interval is determined according to the golden section ratio, and a multi-view image sequence of the polyhedron calibration target is collected at each angle position;

[0009] Adaptive threshold segmentation and sub-pixel edge extraction are performed on the collected image sequence to identify the coded feature points and extract the two-dimensional image coordinates of each feature point; the two-dimensional image coordinates are reconstructed into high-precision three-dimensional space coordinates using the multi-view geometry principle and camera array parameters; based on the reconstructed three-dimensional space coordinates, a feature point trajectory fitting optimization model considering gravity deformation and temperature drift is constructed, and the model takes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the gravity deformation parameters and the temperature drift coefficient as parameters to be optimized; the optimization model is solved using the robust Levenberg-Marquardt algorithm to obtain the optimal rotation axis parameters; using the optimal rotation axis parameters, a dynamically compensated rotation axis coordinate system is established;

[0010] The six-degree-of-freedom posture transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotating axis coordinate system at each angular position is calculated; based on the posture transformation matrix sequence, the actual rotation angle, axial play, radial runout and tilt error of the rotating axis are extracted; the wavelet analysis method is used to perform time-frequency analysis on the extracted actual rotation angle, axial play, radial runout and tilt error of the rotating axis to identify periodic error components and random error components; based on the periodic error components, a rotating axis geometric error model considering the mechanical structure characteristics is established; the random error components are modeled using a neural network algorithm to construct an adaptive error prediction model; the geometric error model and the adaptive error prediction model are combined to form a comprehensive error compensation model; in actual 3D position measurement, the comprehensive error compensation model is used to correct the coordinate transformation in real time according to the current angle of the rotating axis and environmental parameters.

[0011] In an optional embodiment,

[0012] Adaptive threshold segmentation and sub-pixel edge extraction are performed on the acquired image sequence to identify the coded feature points and extract the two-dimensional image coordinates of each feature point; using the multi-view geometry principle and camera array parameters, the two-dimensional image coordinates are reconstructed into high-precision three-dimensional space coordinates, including:

[0013] Performing adaptive threshold segmentation based on multi-scale analysis on the acquired multi-view image sequence, including performing wavelet decomposition on each image to obtain sub-images of different scales, applying the Otsu algorithm to each sub-image to obtain a local optimal threshold, obtaining a global optimal threshold by weighted fusion of thresholds of each scale, and using the global optimal threshold to segment the original image;

[0014] Improved Zernike moment sub-pixel edge extraction is performed on the segmented image, including using the Canny algorithm to obtain the initial edge, constructing a circular window of variable size in the neighborhood of the edge point, calculating the Zernike moment up to the 7th order in the window, fitting the edge curve equation using the least squares method, and solving the equation to obtain the sub-pixel edge coordinates;

[0015] Identify the coded feature points of the two-dimensional code mark designed based on the pseudo-random sequence in the extracted edges, the identification process includes detecting the positioning pattern using Hough transform, determining the two-dimensional code area based on the positioning pattern, extracting the code element using adaptive binarization, applying perspective transform to correct the distortion, decoding and verifying to obtain the unique identification of the feature point, and extracting the two-dimensional image coordinates of each coded feature point;

[0016] Reconstructing the two-dimensional image coordinates into three-dimensional space coordinates by using an improved bundle adjustment algorithm, wherein the improved bundle adjustment algorithm comprises using a Schur complement method to decompose the optimization problem into two sub-problems of camera parameters and feature point coordinates, applying a preconditioned conjugate gradient method to solve the camera parameter sub-problem, and using a sparse Cholesky decomposition and parallel calculation to solve the feature point coordinate sub-problem;

[0017] Establish a distortion correction model, including using the Brown model to describe radial distortion and tangential distortion, and using thermal imager data to establish a dynamic compensation model for the effect of temperature on camera intrinsic parameters;

[0018] Substituting the reconstructed three-dimensional space coordinates into the distortion correction model to obtain corrected high-precision three-dimensional space coordinates;

[0019] The design of the two-dimensional code mark includes generating an M sequence with a length of 63 as a basic code, rearranging the M sequence into a 9x7 matrix, adding positioning patterns at the four corners of the matrix, and using the BCH (63, 51) code to perform error correction encoding on the data.

[0020] In an optional embodiment,

[0021] Reconstructing the two-dimensional image coordinates into three-dimensional space coordinates using an improved bundle adjustment algorithm includes:

[0022] The multi-view 3D reconstruction problem is modeled as a nonlinear least squares problem, with the objective function being the sum of squares of the reprojection error. The SIFT feature matching and PnP algorithm are used to obtain the initial camera pose and feature point 3D coordinate estimation.

[0023] The Schur complement method is used to decompose the optimization problem into camera parameter sub-problems and feature point coordinate sub-problems, and the camera parameter increment is obtained using the Schur complement formula;

[0024] The preconditioned conjugate gradient method is applied to solve the camera parameter subproblem, including:

[0025] Construct a preconditioner matrix, using incomplete Cholesky decomposition as a preconditioner;

[0026] Initialize residual and search direction;

[0027] Iterate the solution until convergence or the maximum number of iterations is reached;

[0028] The sparse Cholesky decomposition and parallel calculation are used for the feature point coordinate sub-problem, including:

[0029] Construct the Hessian matrix of the feature point coordinates;

[0030] Perform symbolic analysis to determine filling patterns and elimination order;

[0031] Use multithreading to parallelize numerical decomposition;

[0032] Use the decomposition results to solve the feature point coordinate increments;

[0033] Combine the camera parameter increment and feature point coordinate increment to update the current estimate;

[0034] Introducing an adaptive damping factor strategy, including:

[0035] Initialize the damping factor;

[0036] Calculate the actual reduction amount for the current iteration;

[0037] Dynamically adjust the damping factor according to the actual reduction amount;

[0038] Determine whether to accept the result of this iteration. If so, update the estimated value. Otherwise, increase the damping factor and solve again.

[0039] In an optional embodiment,

[0040] Based on the reconstructed three-dimensional spatial coordinates, a feature point trajectory fitting optimization model considering gravity deformation and temperature drift is constructed. The model takes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the gravity deformation parameters and the temperature drift coefficient as parameters to be optimized; the robust Levenberg-Marquardt algorithm is used to solve the optimization model to obtain the optimal rotation axis parameters; using the optimal rotation axis parameters, the dynamic compensation rotation axis coordinate system is established, including:

[0041] Based on the reconstructed three-dimensional space coordinates, a feature point trajectory fitting optimization model is constructed. The model takes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the gravity deformation parameter and the temperature drift coefficient as parameters to be optimized.

[0042] Define the ideal position calculation formula of the feature point at time t, which includes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the initial position and angular velocity of the feature point;

[0043] A gravity deformation correction term is introduced, which is a function of the gravity deformation parameter, the gravity acceleration vector and the angular velocity; at the same time, a temperature drift correction term is introduced, which is a function of the temperature drift coefficient, the current temperature and the initial time;

[0044] Based on the ideal position calculation formula, gravity deformation correction term and temperature drift correction term, define the actual position calculation formula of the feature point at time t;

[0045] Construct the objective function, which is the sum of squared errors between the observed positions of all feature points at all times and the calculated actual positions;

[0046] The robust Levenberg-Marquardt algorithm is used to solve the optimization model, including: initializing parameters and damping factors, calculating the current residual and Jacobian matrix, defining weights using the Huber loss function and calculating weighted matrix products, solving incremental equations and updating parameters, dynamically adjusting the damping factor according to changes in the objective function value, and iterating the optimization until convergence or the maximum number of iterations is reached;

[0047] Based on the optimal rotation axis parameters, a dynamically compensated rotation axis coordinate system is established, including: defining the rotation axis direction as the Z axis, calculating the X axis perpendicular to the Z axis and the gravity direction, obtaining the Y axis by cross-producting the X axis and the Z axis, defining the origin considering the influence of gravity deformation and temperature drift, and finally constructing the rotation matrix and translation vector;

[0048] The constructed rotation matrix and translation vector are used to transform the coordinates of any point in the world coordinate system into the dynamic compensation rotation axis coordinate system.

[0049] In an optional embodiment,

[0050] Calculating the six-degree-of-freedom posture transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotation axis coordinate system at each angular position; extracting the actual rotation angle, axial play, radial runout and tilt error of the rotation axis based on the posture transformation matrix sequence includes:

[0051] Calculate the six-degree-of-freedom pose transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotation axis coordinate system at each angular position. The calculation process includes:

[0052] Obtaining the coordinates of the feature points on the polyhedron calibration target in the dynamically compensated rotation axis coordinate system and the target coordinate system;

[0053] Calculate the centroid of two sets of coordinate points and obtain the centroid coordinates;

[0054] Construct the covariance matrix and perform singular value decomposition;

[0055] Calculate the rotation matrix based on the singular value decomposition result, and determine whether correction is needed based on the determinant of the rotation matrix;

[0056] Calculate the translation vector and finally get the six-degree-of-freedom posture transformation matrix;

[0057] Based on the pose transformation matrix sequence, the error parameters of the rotation axis are extracted, including:

[0058] Calculate the actual rotation angle through the rotation matrix in two adjacent posture transformation matrices;

[0059] Compare the components of the translation vectors of two adjacent positions in the Z-axis direction to obtain the axial movement;

[0060] Calculate the displacement of the translation vector on the XY plane to obtain the radial runout;

[0061] Analyze the change of the third column of the rotation matrix with the angle, use the least squares method to fit the conical surface, and solve the direction of the average rotation axis;

[0062] The tilt error is obtained by calculating the angle between the average rotation axis direction and the ideal Z axis.

[0063] In an optional embodiment,

[0064] Based on the periodic error components, a rotating axis geometric error model considering the mechanical structure characteristics is established; the random error components are modeled using a neural network algorithm to construct an adaptive error prediction model; the geometric error model and the adaptive error prediction model are combined to form a comprehensive error compensation model; in actual 3D level measurement, according to the current angle of the rotating axis and environmental parameters, the coordinate transformation is corrected in real time using the comprehensive error compensation model, including:

[0065] A rotating axis geometric error model considering the mechanical structure characteristics is established. The model includes:

[0066] Fourier series are used to describe the periodic error components of axial play, radial runout and tilt errors. The Bayesian information criterion is used to determine the optimal harmonic order to balance model complexity and fitting accuracy.

[0067] Construct an adaptive error prediction model, which includes:

[0068] Use long short-term memory network to model random error components; use current angle, environmental parameters and historical error data as network input; use multi-layer long short-term memory network structure, including dropout layer to improve generalization ability; use mean square error as loss function, and use Adam algorithm for optimization;

[0069] A comprehensive error compensation model is formed, which includes:

[0070] Combine the outputs of the geometric error model and the adaptive error prediction model; generate a comprehensive error estimate taking into account the parameters of current angle, temperature, and humidity;

[0071] Correction of coordinate transformation in real time during actual 3D level measurement. The process includes:

[0072] The ideal coordinate transformation matrix is ​​calculated according to the current angle of the rotation axis; the error compensation matrix is ​​calculated using the comprehensive error compensation model; the ideal coordinate transformation matrix is ​​multiplied by the error compensation matrix to obtain a corrected coordinate transformation matrix; and the coordinate transformation of the measurement point is performed using the corrected coordinate transformation matrix.

[0073] According to a second aspect of the embodiments of the present invention,

[0074] Provides a spatial rotation axis calibration system for 3D level measurement, including:

[0075] The first unit is used to install a polyhedron calibration target on the rotation axis to be calibrated, wherein the polyhedron calibration target includes at least four feature planes that are not coplanarly distributed, and each feature plane is provided with at least three coded feature points that are not colinearly distributed; a plurality of high-resolution industrial cameras are arranged to ensure that at least two cameras simultaneously observe at least one feature plane of the polyhedron calibration target during the rotation of the rotation axis; the intrinsic and extrinsic parameters of the camera array are calibrated to establish a global coordinate system of the camera array; the rotation axis is driven to rotate one circle at a non-uniform angle interval, the rotation interval is determined according to the golden section ratio, and a multi-view image sequence of the polyhedron calibration target is collected at each angle position;

[0076] The second unit is used to perform adaptive threshold segmentation and sub-pixel edge extraction on the acquired image sequence, identify the encoded feature points, and extract the two-dimensional image coordinates of each feature point; reconstruct the two-dimensional image coordinates into high-precision three-dimensional space coordinates using the multi-view geometry principle and camera array parameters; based on the reconstructed three-dimensional space coordinates, construct a feature point trajectory fitting optimization model that considers gravity deformation and temperature drift, and the model takes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the gravity deformation parameters and the temperature drift coefficient as parameters to be optimized; use the robust Levenberg-Marquardt algorithm to solve the optimization model to obtain the optimal rotation axis parameters; use the optimal rotation axis parameters to establish a dynamically compensated rotation axis coordinate system;

[0077] The third unit is used to calculate the six-degree-of-freedom posture transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotating axis coordinate system at each angular position; based on the posture transformation matrix sequence, the actual rotation angle, axial play, radial runout and tilt error of the rotating axis are extracted; the wavelet analysis method is used to perform time-frequency analysis on the extracted actual rotation angle, axial play, radial runout and tilt error of the rotating axis to identify periodic error components and random error components; based on the periodic error components, a rotating axis geometric error model considering the mechanical structure characteristics is established; the random error components are modeled using a neural network algorithm to construct an adaptive error prediction model; the geometric error model and the adaptive error prediction model are combined to form a comprehensive error compensation model; in actual 3D position measurement, the comprehensive error compensation model is used to correct the coordinate transformation in real time according to the current angle of the rotating axis and environmental parameters.

[0078] According to a third aspect of the embodiments of the present invention,

[0079] An electronic device is provided, comprising:

[0080] processor;

[0081] a memory for storing processor-executable instructions;

[0082] The processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0083] According to a fourth aspect of the embodiments of the present invention,

[0084] A computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the aforementioned method is implemented.

[0085] The effects of this application are as follows:

[0086] The calibration accuracy and robustness of the turntable are improved. By acquiring the spatial point cloud data of the plane plate at different angles and performing plane fitting, multiple high-precision spatial plane equations are obtained. Using the constraint relationship between the plane normal vector and the point-to-plane distance, a distance equation group is established, and the parameters of the axis line equation are obtained by solving. This method makes full use of the geometric constraint information of the plane feature, overcomes the influence of point cloud data noise and outliers, and realizes the high-precision calibration of the turntable rotation axis.

[0087] The turntable calibration process is simplified and the hardware cost is reduced. The traditional turntable calibration method usually requires a precise calibration target or a specific calibration object, while the present invention only needs to use a simple flat plate as a calibration tool, without the need to specially design and make a calibration object, thus simplifying the calibration process and reducing the hardware cost and operation difficulty.

[0088] The efficient 3D reconstruction of large and complex objects is achieved. The calibrated rotation axis linear equation parameters are used to seamlessly integrate the 3D data obtained at different angles of the measured object through coordinate transformation and point cloud splicing to obtain a complete 3D model. This method avoids the tedious point cloud registration process and greatly improves the efficiency of 3D reconstruction. It is suitable for 3D measurement of complex objects such as large workpieces and cultural relics, which are difficult to directly obtain complete data.

[0089] The flexibility and adaptability of three-dimensional measurement are improved. The turntable calibration method of the present invention does not rely on specific scanning equipment and measurement environment, and is applicable to various types of non-contact three-dimensional scanners such as structured light and laser. Automatic splicing of point cloud data can be achieved through turntable calibration, and the size, shape and placement posture of the measured object can be flexibly selected, which expands the application scope of three-dimensional measurement technology.

[0090] It provides technical support for high-precision 3D detection and reverse engineering. Through turntable calibration, high-precision linear equations of the rotation axis are obtained. Combined with coordinate transformation and point cloud stitching technology, precise 3D detection and quality control of complex workpieces can be achieved. At the same time, the generated high-quality 3D model provides an important data foundation for reverse engineering design, computer-aided manufacturing and other fields, promoting the intelligent and digital development of the manufacturing industry. BRIEF DESCRIPTION OF THE DRAWINGS

[0091] Figure 1 It is a schematic flow chart of a method for calibrating a spatial rotation axis for 3D level measurement according to an embodiment of the present invention;

[0092] Figure 2 The figure is a schematic diagram of the structure of a spatial rotation axis calibration system for 3D object level measurement according to an embodiment of the present invention. DETAILED DESCRIPTION

[0093] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0094] The technical solution of the present invention is described in detail with specific embodiments below. The following specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described in detail in some embodiments.

[0095] Figure 1 FIG. 1 is a flow chart of a method for calibrating a spatial rotation axis for 3D level measurement according to an embodiment of the present invention. Figure 1 As shown, the method includes:

[0096] S101. Install a polyhedron calibration target on the rotation axis to be calibrated, wherein the polyhedron calibration target includes at least four feature planes that are not coplanarly distributed, and each feature plane is provided with at least three coded feature points that are not colinearly distributed; arrange multiple high-resolution industrial cameras to ensure that at least two cameras simultaneously observe at least one feature plane of the polyhedron calibration target during the rotation of the rotation axis; perform intrinsic and extrinsic calibration of the camera array to establish a global coordinate system of the camera array; drive the rotation axis to rotate one circle at a non-uniform angle interval, the rotation interval is determined according to the golden ratio, and collect a multi-view image sequence of the polyhedron calibration target at each angle position;

[0097] S102. Perform adaptive threshold segmentation and sub-pixel edge extraction on the acquired image sequence, identify the coded feature points, and extract the two-dimensional image coordinates of each feature point; reconstruct the two-dimensional image coordinates into high-precision three-dimensional space coordinates using the multi-view geometry principle and camera array parameters; construct a feature point trajectory fitting optimization model that considers gravity deformation and temperature drift based on the reconstructed three-dimensional space coordinates, and the model takes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the gravity deformation parameters and the temperature drift coefficient as parameters to be optimized; solve the optimization model using the robust Levenberg-Marquardt algorithm to obtain the optimal rotation axis parameters; use the optimal rotation axis parameters to establish a dynamically compensated rotation axis coordinate system;

[0098] S103. Calculate the six-degree-of-freedom posture transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotating axis coordinate system at each angular position; extract the actual rotation angle, axial play, radial runout and tilt error of the rotating axis based on the posture transformation matrix sequence; use the wavelet analysis method to perform time-frequency analysis on the extracted actual rotation angle, axial play, radial runout and tilt error of the rotating axis, and identify periodic error components and random error components; based on the periodic error components, establish a rotating axis geometric error model that takes into account the mechanical structure characteristics; use a neural network algorithm to model the random error components and construct an adaptive error prediction model; combine the geometric error model and the adaptive error prediction model to form a comprehensive error compensation model; in actual 3D position measurement, use the comprehensive error compensation model to correct the coordinate transformation in real time according to the current angle of the rotating axis and environmental parameters.

[0099] In an optional embodiment,

[0100] Adaptive threshold segmentation and sub-pixel edge extraction are performed on the acquired image sequence to identify the coded feature points and extract the two-dimensional image coordinates of each feature point; using the multi-view geometry principle and camera array parameters, the two-dimensional image coordinates are reconstructed into high-precision three-dimensional space coordinates, including:

[0101] Performing adaptive threshold segmentation based on multi-scale analysis on the acquired multi-view image sequence, including performing wavelet decomposition on each image to obtain sub-images of different scales, applying the Otsu algorithm to each sub-image to obtain a local optimal threshold, obtaining a global optimal threshold by weighted fusion of thresholds of each scale, and using the global optimal threshold to segment the original image;

[0102] Improved Zernike moment sub-pixel edge extraction is performed on the segmented image, including using the Canny algorithm to obtain the initial edge, constructing a circular window of variable size in the neighborhood of the edge point, calculating the Zernike moment up to the 7th order in the window, fitting the edge curve equation using the least squares method, and solving the equation to obtain the sub-pixel edge coordinates;

[0103] Identify the coded feature points of the two-dimensional code mark designed based on the pseudo-random sequence in the extracted edges, the identification process includes detecting the positioning pattern using Hough transform, determining the two-dimensional code area based on the positioning pattern, extracting the code element using adaptive binarization, applying perspective transform to correct the distortion, decoding and verifying to obtain the unique identification of the feature point, and extracting the two-dimensional image coordinates of each coded feature point;

[0104] Reconstructing the two-dimensional image coordinates into three-dimensional space coordinates by using an improved bundle adjustment algorithm, wherein the improved bundle adjustment algorithm comprises using a Schur complement method to decompose the optimization problem into two sub-problems of camera parameters and feature point coordinates, applying a preconditioned conjugate gradient method to solve the camera parameter sub-problem, and using a sparse Cholesky decomposition and parallel calculation to solve the feature point coordinate sub-problem;

[0105] Establish a distortion correction model, including using the Brown model to describe radial distortion and tangential distortion, and using thermal imager data to establish a dynamic compensation model for the effect of temperature on camera intrinsic parameters;

[0106] Substituting the reconstructed three-dimensional space coordinates into the distortion correction model to obtain corrected high-precision three-dimensional space coordinates;

[0107] The design of the two-dimensional code mark includes generating an M sequence with a length of 63 as a basic code, rearranging the M sequence into a 9x7 matrix, adding positioning patterns at the four corners of the matrix, and using the BCH (63, 51) code to perform error correction encoding on the data.

[0108] Exemplarily, adaptive threshold segmentation and sub-pixel edge extraction are performed on the acquired image sequence, coded feature points are identified, and the two-dimensional image coordinates of each feature point are extracted; the detailed technical contents of reconstructing the two-dimensional image coordinates into high-precision three-dimensional space coordinates by using multi-view geometry principles and camera array parameters are as follows:

[0109] First, adaptive threshold segmentation based on multi-scale analysis is performed. Discrete wavelet transform is applied to each image, usually using Daubechies wavelet, and 3-5 layers of decomposition are performed to obtain low-frequency and high-frequency sub-images of different scales. The Otsu algorithm is applied to each sub-image to calculate the local optimal threshold. The Otsu algorithm determines the optimal threshold by maximizing the inter-class variance. Then, different weights are assigned to the thresholds of each sub-image according to the importance of different scales. Usually, low-frequency sub-images have higher weights, such as 0.5, and high-frequency sub-images have lower weights, such as 0.1-0.2. The weighted thresholds are summed to obtain the global optimal threshold. Finally, the global threshold is used to perform binary segmentation on the original image.

[0110] Next, the improved Zernike moment sub-pixel edge extraction is performed. First, the Canny algorithm is used to obtain the initial edge. The Canny algorithm includes Gaussian filtering, calculation of gradient amplitude and direction, non-maximum suppression, double threshold detection and edge connection. A circular window of variable size is constructed near the obtained edge point, and the window radius is usually 5-15 pixels. Zernike moments up to the 7th order are calculated within the window. Zernike moments are a set of orthogonal polynomial basis functions. The calculated Zernike moments are fitted into the edge curve equation using the least squares method. Solving the equation obtains the edge coordinates with sub-pixel accuracy, with an accuracy of about 0.1 pixel.

[0111] Then identify the coded feature points of the two-dimensional code mark designed based on the pseudo-random sequence. First, use the Hough transform to detect the positioning pattern. The Hough transform can effectively identify geometric shapes such as straight lines and circles in the image. Determine the two-dimensional code area based on the detected positioning pattern. Apply adaptive binarization to the area to extract code elements. Adaptive binarization can overcome the influence of uneven image brightness. Use perspective transformation to correct the distortion caused by the shooting angle. Decode and verify the unique identification of the feature point, and extract the two-dimensional image coordinates of each coded feature point.

[0112] The improved bundle adjustment algorithm is used to reconstruct the 2D image coordinates into 3D space coordinates. The improved bundle adjustment algorithm first uses the Schur complement method to decompose the optimization problem into two sub-problems: camera parameters and feature point coordinates. The preconditioned conjugate gradient method is applied to solve the camera parameter sub-problem, and the preconditioned matrix is ​​usually selected as a diagonal approximation. The sparse Cholesky decomposition is used for parallel calculation of the feature point coordinate sub-problem, which can significantly improve the computational efficiency. Iterate the optimization until convergence or the maximum number of iterations is reached, such as 100 times.

[0113] Establish a distortion correction model. Use the Brown model to describe radial distortion and tangential distortion, usually considering 2-3 order radial distortion coefficients and 2 tangential distortion coefficients. Use thermal imager data to establish a dynamic compensation model for the effect of temperature on camera internal parameters, usually establishing a linear or quadratic polynomial relationship between temperature changes and changes in internal parameters such as focal length and principal point coordinates.

[0114] Finally, the reconstructed 3D space coordinates are substituted into the distortion correction model to obtain the corrected high-precision 3D space coordinates. The correction process usually uses an iterative method until the coordinate change is less than a preset threshold, such as 0.001mm.

[0115] The design process of the QR code tag is as follows: First, generate an M sequence of length 63 as the basic code. The M sequence is a pseudo-random binary sequence with good autocorrelation. The M sequence is rearranged into a 9x7 matrix, which can form a compact QR code structure. Add positioning patterns at the four corners of the matrix, usually using concentric squares as positioning patterns. Use BCH (63, 51) code to encode the data for error correction. BCH code has strong error correction capabilities and can effectively resist the influence of noise and local damage.

[0116] In practical applications, this method can be used in high-precision three-dimensional measurement systems. For example, in an industrial inspection scenario, 5 industrial cameras are used to form a camera array with a camera resolution of 4096x3000 pixels. 100 coded feature points are arranged on the surface of the object to be measured. After the collected image is processed, the positioning accuracy of the two-dimensional feature points can reach 0.05 pixels. After three-dimensional reconstruction and distortion correction, the final three-dimensional coordinate accuracy can reach 0.02mm. The system can maintain stable operation in the temperature range of -10°C to 40°C, and the measurement error caused by temperature changes does not exceed 0.005mm.

[0117] In an optional embodiment,

[0118] Reconstructing the two-dimensional image coordinates into three-dimensional space coordinates using an improved bundle adjustment algorithm includes:

[0119] The multi-view 3D reconstruction problem is modeled as a nonlinear least squares problem, with the objective function being the sum of squares of the reprojection error. The SIFT feature matching and PnP algorithm are used to obtain the initial camera pose and feature point 3D coordinate estimation.

[0120] The Schur complement method is used to decompose the optimization problem into camera parameter sub-problems and feature point coordinate sub-problems, and the camera parameter increment is obtained using the Schur complement formula;

[0121] The preconditioned conjugate gradient method is applied to solve the camera parameter subproblem, including:

[0122] Construct a preconditioner matrix, using incomplete Cholesky decomposition as a preconditioner;

[0123] Initialize residual and search direction;

[0124] Iterate the solution until convergence or the maximum number of iterations is reached;

[0125] The sparse Cholesky decomposition and parallel calculation are used for the feature point coordinate sub-problem, including:

[0126] Construct the Hessian matrix of the feature point coordinates;

[0127] Perform symbolic analysis to determine filling patterns and elimination order;

[0128] Use multithreading to parallelize numerical decomposition;

[0129] Use the decomposition results to solve the feature point coordinate increments;

[0130] Combine the camera parameter increment and feature point coordinate increment to update the current estimate;

[0131] Introducing an adaptive damping factor strategy, including:

[0132] Initialize the damping factor;

[0133] Calculate the actual reduction amount for the current iteration;

[0134] Dynamically adjust the damping factor according to the actual reduction amount;

[0135] Determine whether to accept the result of this iteration. If so, update the estimated value. Otherwise, increase the damping factor and solve again.

[0136] Exemplarily, the detailed technical content of reconstructing the two-dimensional image coordinates into three-dimensional space coordinates using the improved bundle adjustment algorithm is as follows:

[0137] First, the multi-view 3D reconstruction problem is modeled as a nonlinear least squares problem. The objective function is defined as the sum of squares of the reprojection error, where the reprojection error refers to the difference between the image coordinates calculated by the current estimated camera parameters and the 3D point coordinates and the actual observed image coordinates. To obtain an initial estimate, the SIFT (Scale Invariant Feature Transform) algorithm is used to extract image feature points and perform cross-view matching. Then a PnP (Perspective-n-Point) algorithm, such as EPnP or UPnP, is applied to estimate the initial pose of the camera. Combined with the known camera intrinsic parameters, the initial 3D coordinate estimate of the feature points is obtained by triangulation.

[0138] Next, the Schur complement method is used to decompose the optimization problem into camera parameter subproblems and feature point coordinate subproblems. This decomposition takes advantage of the sparse structure of the problem and can significantly improve computational efficiency. Specifically, the camera parameter increment is derived using the Schur complement formula. The coefficient matrix of this equation is usually called the reduced camera matrix, and its dimension is comparable to the number of camera parameters, which is much smaller than the scale of the original problem.

[0139] For the camera parameter subproblem, the preconditioned conjugate gradient method is applied to solve it. First, the preconditioned matrix is ​​constructed. The commonly used choice is incomplete Cholesky decomposition. The role of the preconditioned matrix is ​​to improve the condition number of the coefficient matrix and accelerate convergence. Initialize the residual vector to the right vector minus the product of the coefficient matrix and the current solution, and set the initial search direction to the product of the preconditioned matrix and the residual. Then enter the iterative solution process: calculate the step size, update the solution vector and residual, and calculate the new search direction. Iterate until the residual norm is less than the preset threshold (such as 1e-6) or the maximum number of iterations (such as 50 times) is reached.

[0140] For the feature point coordinate subproblem, sparse Cholesky decomposition is used for parallel calculation. First, the Hessian matrix of the feature point coordinates is constructed, which has a block diagonal structure. Symbolic analysis is performed to determine the filling mode and the optimal elimination order. This step can significantly reduce the complexity of subsequent numerical calculations. Multi-threading technology is used to perform numerical decomposition in parallel, and each thread is responsible for processing a part of the feature points. Finally, the decomposition results are used to solve the feature point coordinate increments.

[0141] Combine the obtained camera parameter increments and feature point coordinate increments to update the current estimate. In order to improve the robustness and convergence of the algorithm, an adaptive damping factor strategy is introduced. Initialize the damping factor to a small positive value, such as 1e-3. In each iteration, calculate the actual reduction of the objective function and compare it with the expected reduction. If the actual reduction is close to the expected, reduce the damping factor (such as multiplying by 0.5); if the actual reduction is much less than expected, increase the damping factor (such as multiplying by 2). Decide whether to accept the result of this iteration based on the reduction. If accepted, update the estimate; otherwise, increase the damping factor and solve again.

[0142] In practical applications, the algorithm can be used for large-scale 3D reconstruction tasks. For example, in a 3D reconstruction project of a certain urban building complex, 100 high-resolution cameras (8192x6144 pixels) were used to capture about 10,000 images. An average of 5,000 SIFT feature points were extracted from each image, and about 5 million 3D points were obtained after matching. The initial PnP solution and triangulation took about 30 minutes. The improved bundle adjustment algorithm was run on a workstation equipped with a 64-core CPU and 256GB memory, with the maximum number of iterations set to 100 and the convergence threshold set to an average reprojection error of less than 0.5 pixels. The algorithm converged in about 2 hours, and the final 3D reconstruction model had an average reprojection error of 0.3 pixels, which is equivalent to an accuracy of about 1 cm in real space. By comparing with the laser scanning results, the size error of the reconstructed building is within the range of ±2 cm, which meets the needs of most urban planning and cultural relics protection.

[0143] The key improvements of the algorithm are problem decomposition and parallel computing strategies, which enable it to efficiently process large-scale data sets. The adaptive damping factor strategy improves the stability and convergence of the algorithm in complex scenarios. The combination of these technologies makes high-precision large-scale 3D reconstruction possible, providing strong technical support for urban modeling, digital protection of cultural relics and other fields.

[0144] In an optional embodiment,

[0145] Based on the reconstructed three-dimensional spatial coordinates, a feature point trajectory fitting optimization model considering gravity deformation and temperature drift is constructed. The model takes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the gravity deformation parameters and the temperature drift coefficient as parameters to be optimized; the robust Levenberg-Marquardt algorithm is used to solve the optimization model to obtain the optimal rotation axis parameters; using the optimal rotation axis parameters, the dynamic compensation rotation axis coordinate system is established, including:

[0146] Based on the reconstructed three-dimensional space coordinates, a feature point trajectory fitting optimization model is constructed. The model takes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the gravity deformation parameter and the temperature drift coefficient as parameters to be optimized.

[0147] Define the ideal position calculation formula of the feature point at time t, which includes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the initial position and angular velocity of the feature point;

[0148] A gravity deformation correction term is introduced, which is a function of the gravity deformation parameter, the gravity acceleration vector and the angular velocity; at the same time, a temperature drift correction term is introduced, which is a function of the temperature drift coefficient, the current temperature and the initial time;

[0149] Based on the ideal position calculation formula, gravity deformation correction term and temperature drift correction term, define the actual position calculation formula of the feature point at time t;

[0150] Construct the objective function, which is the sum of squared errors between the observed positions of all feature points at all times and the calculated actual positions;

[0151] The robust Levenberg-Marquardt algorithm is used to solve the optimization model, including: initializing parameters and damping factors, calculating the current residual and Jacobian matrix, defining weights using the Huber loss function and calculating weighted matrix products, solving incremental equations and updating parameters, dynamically adjusting the damping factor according to changes in the objective function value, and iterating the optimization until convergence or the maximum number of iterations is reached;

[0152] Based on the optimal rotation axis parameters, a dynamically compensated rotation axis coordinate system is established, including: defining the rotation axis direction as the Z axis, calculating the X axis perpendicular to the Z axis and the gravity direction, obtaining the Y axis by cross-producting the X axis and the Z axis, defining the origin considering the influence of gravity deformation and temperature drift, and finally constructing the rotation matrix and translation vector;

[0153] The constructed rotation matrix and translation vector are used to transform the coordinates of any point in the world coordinate system into the dynamic compensation rotation axis coordinate system.

[0154] Exemplarily, based on the reconstructed three-dimensional space coordinates, a feature point trajectory fitting optimization model considering gravity deformation and temperature drift is constructed, and the detailed technical contents of establishing a dynamically compensated rotating axis coordinate system are as follows:

[0155] First, a feature point trajectory fitting optimization model is constructed. The model contains multiple parameters to be optimized: the rotation axis direction vector (a unit vector, represented by three components), the three-dimensional coordinates of a point on the rotation axis (used to determine the position of the axis), the gravity deformation parameter (describing the degree of influence of gravity on the rotating system) and the temperature drift coefficient (describing the influence of temperature changes on the system). These parameters together determine the actual trajectory of the feature point during the rotation process.

[0156] Next, define the formula for calculating the ideal position of the feature point at any time t. This formula takes into account the direction and position of the rotation axis, the initial position of the feature point, and the angular velocity. Specifically, the ideal position of the feature point can be obtained by rotating the initial position around the rotation axis by a specific angle, which is equal to the product of the angular velocity and the time t.

[0157] In order to more accurately describe the actual situation, a gravity deformation correction term is introduced. This correction term is a function of the gravity deformation parameter, the gravity acceleration vector, and the angular velocity. Gravity deformation will cause a slight bend in the rotation axis, thus affecting the actual trajectory of the feature point. At the same time, a temperature drift correction term is introduced, which is a function of the temperature drift coefficient, the temperature difference between the current moment and the initial moment. Temperature changes will cause thermal expansion or contraction of the system, which will affect the position of the feature point.

[0158] Taking into account the ideal position calculation formula, gravity deformation correction term and temperature drift correction term, the actual position calculation formula of the feature point at time t is defined. This formula combines the ideal position with the two correction terms to more accurately describe the motion trajectory of the feature point in the actual environment.

[0159] Construct an objective function, which is the sum of the squares of the errors between the observed positions of all feature points at all times and the calculated actual positions. This objective function quantifies the difference between the model predictions and the actual observations, and is what needs to be minimized during the optimization process.

[0160] The robust Levenberg-Marquardt algorithm is used to solve the optimization model. First, the parameters and damping factors are initialized. Reasonable initial estimates can be given based on prior knowledge. Then the current residual (the difference between the observed value and the predicted value) and the Jacobian matrix (the partial derivative of the objective function with respect to each parameter) are calculated. In order to improve the robustness of the algorithm, the Huber loss function is used to define the weights, which can reduce the impact of outliers. The weighted matrix product is calculated to form the coefficient matrix and right-hand side vector of the incremental equation. The incremental equation is solved to obtain the parameter update amount and the parameter value is updated. The damping factor is dynamically adjusted according to the change in the objective function value: if the objective function decreases, the damping factor is reduced; otherwise, the damping factor is increased. This process is repeated until convergence (such as the parameter change is less than the threshold) or the maximum number of iterations is reached.

[0161] Based on the optimal rotation axis parameters, a dynamically compensated rotation axis coordinate system is established. First, the optimized rotation axis direction is defined as the Z axis. Then the vector perpendicular to both the Z axis and the gravity direction is calculated as the X axis, which ensures that the X axis is in the horizontal plane. The Y axis is obtained by the cross product of the X axis and the Z axis, thus forming an orthogonal coordinate system. Considering the influence of gravity deformation and temperature drift, the origin of the coordinate system is defined. Finally, the rotation matrix and translation vector are constructed based on these axis vectors and the origin.

[0162] Using the constructed rotation matrix and translation vector, the coordinates of any point in the world coordinate system can be transformed into the coordinate system of the dynamically compensated rotation axis. This conversion process involves first translating the point to move the point to the origin of the new coordinate system, and then applying the rotation transformation.

[0163] In practical applications, this method can be used for calibration and error compensation of high-precision rotating platforms. For example, in the test of a spacecraft attitude control system, a high-precision two-axis turntable was used. The nominal accuracy of the turntable is 0.001 degrees, but it is affected by gravity and temperature changes in actual work. 20 reflective markers were arranged on the surface of the turntable, and 8 high-speed cameras (1000fps, 2048x2048 resolution) were used to track the markers. The measurement time was 10 minutes, and the temperature changed by about 2°C during the measurement.

[0164] The above method is applied for data processing and model optimization. Initial parameter settings: the rotation axis direction vector is [0,0, 1], the coordinates of a point on the axis are [0, 0, 0], the gravity deformation parameter is 1e-6 m / rad^2, and the temperature drift coefficient is 1e-5m / °C. The maximum number of iterations of the Levenberg-Marquardt algorithm is set to 100, and the convergence threshold is 1e-8.

[0165] After optimization, the results show that the actual direction vector of the rotation axis is [0.000123, -0.000234, 0.999999], the coordinates of a point on the axis are [0.00567, -0.00321, 0.000789] (unit: meter), the gravity deformation parameter is 8.76e-7 m / rad^2, and the temperature drift coefficient is 1.23e-5 m / °C. The root mean square error of the final model is 0.005mm, which is 82% lower than the error (0.028mm) when the gravity and temperature effects are not considered.

[0166] The optimized parameters are used to establish a dynamically compensated rotation axis coordinate system, and the measured marker point coordinates are converted to this coordinate system. The converted coordinates show that the maximum deviation of the marker point relative to the ideal circular trajectory during rotation is 0.012 mm, while the maximum deviation without compensation is 0.078 mm.

[0167] In an optional embodiment,

[0168] Calculating the six-degree-of-freedom posture transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotation axis coordinate system at each angular position; extracting the actual rotation angle, axial play, radial runout and tilt error of the rotation axis based on the posture transformation matrix sequence includes:

[0169] Calculate the six-degree-of-freedom pose transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotation axis coordinate system at each angular position. The calculation process includes:

[0170] Obtaining the coordinates of the feature points on the polyhedron calibration target in the dynamically compensated rotation axis coordinate system and the target coordinate system;

[0171] Calculate the centroid of two sets of coordinate points and obtain the centroid coordinates;

[0172] Construct the covariance matrix and perform singular value decomposition;

[0173] Calculate the rotation matrix based on the singular value decomposition result, and determine whether correction is needed based on the determinant of the rotation matrix;

[0174] Calculate the translation vector and finally get the six-degree-of-freedom posture transformation matrix;

[0175] Based on the pose transformation matrix sequence, the error parameters of the rotation axis are extracted, including:

[0176] Calculate the actual rotation angle through the rotation matrix in two adjacent posture transformation matrices;

[0177] Compare the components of the translation vectors of two adjacent positions in the Z-axis direction to obtain the axial movement;

[0178] Calculate the displacement of the translation vector on the XY plane to obtain the radial runout;

[0179] Analyze the change of the third column of the rotation matrix with the angle, use the least squares method to fit the conical surface, and solve the direction of the average rotation axis;

[0180] The tilt error is obtained by calculating the angle between the average rotation axis direction and the ideal Z axis.

[0181] Exemplarily, the detailed technical contents of calculating the six-degree-of-freedom posture transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotation axis coordinate system at each angular position, and extracting the actual rotation angle, axial play, radial runout and tilt error of the rotation axis based on the posture transformation matrix sequence are as follows:

[0182] First, for each angular position, the six-DOF pose transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotation axis coordinate system is calculated. This process begins by obtaining the coordinates of the feature points on the polyhedron calibration target in two coordinate systems. Typically, multiple highly reflective spherical markers are arranged on the target, and the three-dimensional coordinates of these points can be obtained by a high-precision measurement system (such as an optical tracking system).

[0183] After obtaining the coordinates, calculate the centroid of the two sets of coordinate points. The centroid calculation method is to average the coordinates of all points. Then, subtract the corresponding centroid coordinates from the coordinates of each point to obtain the centroid coordinates. The purpose of this step is to transform the problem into a form that only involves rotation, simplifying subsequent calculations.

[0184] Next, a covariance matrix is ​​constructed. This matrix reflects the correlation between the two sets of de-centering coordinates. Specifically, the de-centering coordinate matrix in the dynamic compensation rotation axis coordinate system is multiplied with the de-centering coordinate matrix in the target coordinate system to obtain a 3x3 covariance matrix.

[0185] Perform singular value decomposition on the covariance matrix. This step decomposes the covariance matrix into the product of three matrices: the left singular vector matrix, the singular value diagonal matrix, and the transpose of the right singular vector matrix. Singular value decomposition is a key step in solving the rotation matrix.

[0186] Based on the results of the singular value decomposition, calculate the rotation matrix. The specific method is to multiply the right singular vector matrix by the transpose of the left singular vector matrix. Then, check the determinant of the rotation matrix. If the determinant is -1, it means that the transformation is a reflection rather than a rotation. You need to invert the third column of the right singular vector matrix and recalculate the rotation matrix.

[0187] Finally, calculate the translation vector. The translation vector is equal to the center of mass coordinates in the dynamic compensation rotation axis coordinate system minus the product of the rotation matrix and the center of mass coordinates in the target coordinate system. Combining the calculated rotation matrix and translation vector, the six-degree-of-freedom posture transformation matrix can be obtained.

[0188] After obtaining the sequence of posture transformation matrices, the error parameters of the rotation axis can be extracted. First, the actual rotation angle is calculated by the rotation matrices in the two adjacent posture transformation matrices. The specific method is to multiply the transpose of the previous rotation matrix with the next rotation matrix. The trace of the obtained matrix (the sum of the main diagonal elements) is related to the cosine of the rotation angle. The actual rotation angle can be obtained through the inverse cosine function.

[0189] The axial play is calculated by comparing the difference in the components of the translation vectors in the Z-axis direction between two adjacent positions. This reflects the amount of axial movement during the rotation process.

[0190] The radial runout is obtained by calculating the displacement of the translation vector in the XY plane. Specifically, the X and Y components of the translation vector are taken and the square root of the sum of their squares is calculated, which is the radial runout.

[0191] In order to analyze the tilt error of the rotation axis, it is necessary to observe the change of the third column of the rotation matrix with the angle. Ideally, if the rotation axis completely coincides with the Z axis, this column of vectors should remain unchanged. In practice, this column of vectors will describe a conical surface in space. Using the least squares method to fit this conical surface, the average rotation axis direction can be solved.

[0192] Finally, the angle between the average rotation axis and the ideal Z axis ([0,0,1]) is calculated, which is the tilt error. This angle can be calculated by the dot product of the two vectors and their modulus.

[0193] In practical applications, a high-precision large turntable calibration is taken as an example. The turntable has a diameter of 5 meters and a nominal accuracy of 0.0001 degrees. A cubic target (side length 50 cm) with 20 spherical markers is used for measurement. Data is collected every 10 degrees within a 360-degree range, for a total of 36 positions. 16 high-speed cameras (2000fps, 4096x3072 resolution) are used to track the target movement, and the measurement accuracy reaches 0.01mm.

[0194] Applying the above method to data processing, the following results are obtained:

[0195] 1. Actual rotation angle: average error is 0.00008 degrees, maximum error is 0.00023 degrees.

[0196] 2. Axial movement: mean value is 0.003mm, standard deviation is 0.001mm.

[0197] 3. Radial runout: Maximum value is 0.015mm, RMS value is 0.008mm.

[0198] 4. Tilt error: 0.00032 degrees.

[0199] These results show that the actual performance of the rotary table is basically consistent with its nominal accuracy, but there are slight deviations in some aspects. In particular, the radial runout is relatively large and may require further adjustment. This high-precision error analysis is crucial for the calibration and optimization of large precision equipment, and can guide subsequent adjustment and compensation strategies, thereby improving the accuracy and reliability of the overall system.

[0200] This method is not only applicable to turntable calibration, but can also be extended to the accuracy assessment of various rotating equipment, such as the rotary axis of CNC machine tools, the azimuth and pitch mechanisms of large antennas, etc. By accurately quantifying various error parameters, it can provide an important basis for the precise control and error compensation of the equipment, and ultimately improve the performance and reliability of the entire system.

[0201] In an optional embodiment,

[0202] Based on the periodic error components, a rotating axis geometric error model considering the mechanical structure characteristics is established; the random error components are modeled using a neural network algorithm to construct an adaptive error prediction model; the geometric error model and the adaptive error prediction model are combined to form a comprehensive error compensation model; in actual 3D level measurement, according to the current angle of the rotating axis and environmental parameters, the coordinate transformation is corrected in real time using the comprehensive error compensation model, including:

[0203] A rotating axis geometric error model considering the mechanical structure characteristics is established. The model includes:

[0204] Fourier series are used to describe the periodic error components of axial play, radial runout and tilt errors. The Bayesian information criterion is used to determine the optimal harmonic order to balance model complexity and fitting accuracy.

[0205] Construct an adaptive error prediction model, which includes:

[0206] Use long short-term memory network to model random error components; use current angle, environmental parameters and historical error data as network input; use multi-layer long short-term memory network structure, including dropout layer to improve generalization ability; use mean square error as loss function, and use Adam algorithm for optimization;

[0207] A comprehensive error compensation model is formed, which includes:

[0208] Combine the outputs of the geometric error model and the adaptive error prediction model; generate a comprehensive error estimate taking into account the parameters of current angle, temperature, and humidity;

[0209] Correction of coordinate transformation in real time during actual 3D level measurement. The process includes:

[0210] The ideal coordinate transformation matrix is ​​calculated according to the current angle of the rotation axis; the error compensation matrix is ​​calculated using the comprehensive error compensation model; the ideal coordinate transformation matrix is ​​multiplied by the error compensation matrix to obtain a corrected coordinate transformation matrix; and the coordinate transformation of the measurement point is performed using the corrected coordinate transformation matrix.

[0211] For example, based on the periodic error components, a rotating axis geometric error model considering the mechanical structure characteristics is established, and the random error components are modeled in combination with the neural network algorithm, and finally a comprehensive error compensation model is formed for real-time coordinate transformation correction in actual 3D level measurement. The detailed technical content of this process is as follows:

[0212] First, a rotating axis geometric error model considering the mechanical structure characteristics is established. This model is mainly used to describe the periodic error components of axial play, radial runout and tilt error. Fourier series is used to represent these periodic errors. Fourier series can decompose complex periodic functions into a series of simple sine and cosine functions. For each error type, a Fourier series model is established.

[0213] The Bayesian Information Criterion (BIC) is used to determine the optimal order of the Fourier series. BIC can strike a balance between model complexity and fitting accuracy to avoid overfitting problems. The specific approach is to gradually increase the order of the Fourier series from a low order and calculate the BIC value at each order. The order that minimizes the BIC value is selected as the optimal order.

[0214] For example, for a large precision turntable, analysis found that the optimal order of axial play is 3, the optimal order of radial runout is 4, and the optimal order of tilt error is 2. This means that axial play can be well described by 3 harmonic components, while radial runout requires 4 harmonic components and tilt error only requires 2 harmonic components.

[0215] Next, an adaptive error prediction model is built to handle random error components. Here, a long short-term memory (LSTM) network is used, which is a special recurrent neural network that can effectively capture long-term dependencies. The input of the LSTM network includes the current angle, environmental parameters (such as temperature and humidity), and historical error data.

[0216] The network structure adopts a multi-layer LSTM design, which specifically includes: an input layer that receives the above input data; two LSTM hidden layers, each containing 64 neurons; a dropout layer is added between the LSTM layers, and the dropout rate is set to 0.2 to prevent overfitting; and finally a fully connected output layer for predicting the error value.

[0217] The network training uses mean square error (MSE) as the loss function and uses the Adam optimization algorithm for parameter update. The learning rate is initially set to 0.001, and the learning rate decay strategy is adopted, decaying by 10% every 50 epochs. During the training process, the early stopping strategy is used, and the training is stopped when the loss on the validation set has not improved for 10 consecutive epochs.

[0218] In actual application, 10,000 sets of historical data were used for training, of which 80% were used for training sets and 20% for validation sets. After 200 epochs of training, the mean square error of the model on the validation set dropped to 0.0015 mm², indicating that the model has good predictive ability.

[0219] When forming the comprehensive error compensation model, the outputs of the geometric error model and the adaptive error prediction model are combined. The geometric error model provides an estimate of the periodic error, while the adaptive error prediction model provides an estimate of the random error. These two error estimates are added together, and the effects of parameters such as the current angle, temperature, and humidity are considered to finally generate a comprehensive error estimate.

[0220] In actual 3D level measurement, the process of real-time correction of coordinate transformation is as follows: First, the ideal coordinate transformation matrix is ​​calculated based on the current angle of the rotation axis. This matrix represents the ideal transformation of the object from one coordinate system to another without error.

[0221] Then, the error compensation matrix is ​​calculated using the comprehensive error compensation model. This matrix contains the correction values ​​for axial play, radial runout, and tilt errors. The ideal coordinate transformation matrix is ​​multiplied by the error compensation matrix to obtain the corrected coordinate transformation matrix. This corrected matrix comprehensively considers the ideal transformation and actual errors, and can more accurately describe the actual coordinate transformation relationship.

[0222] Finally, the coordinate transformation matrix is ​​used to transform the coordinates of the measurement points. In this way, the measurement results include compensation for various error sources, improving the measurement accuracy.

[0223] This method was applied to the 3D scanning measurement of a large wind turbine blade, and achieved remarkable results. The wind turbine blade is 80 meters long and uses a large turntable for rotational scanning. The original measurement data showed that the position error of the blade end reached ±5mm. After applying the comprehensive error compensation model, the position error was reduced to ±0.8mm, and the accuracy was improved by about 84%.

[0224] This comprehensive error compensation method is not only applicable to large turntables, but can also be extended to various precision rotating equipment, such as the azimuth system of large telescopes, the tracking system of satellite antennas, etc. By compensating various error sources in real time, the measurement accuracy and reliability of the system can be significantly improved, providing strong support for high-precision large-scale measurements.

[0225] In addition, an important advantage of this method is its adaptability. Over time and with the use of the equipment, the error characteristics may change. The adaptive error prediction model is able to continuously learn new data and update its prediction ability, thereby maintaining long-term high-precision performance. This is especially important for large equipment that operates for a long time, which can reduce the frequency of maintenance and extend the effective service life of the equipment.

[0226] Another embodiment of the present application is as follows:

[0227] Solving the equation of the rotation axis based on the plane equation

[0228] Assume a plane rotates around an axis, and measure The space plane equations are:

[0229] ;

[0230] Here the coefficients of the plane equation are is the unit normal vector of each plane, and the normal vectors have the same direction, all pointing to a certain point on the axis of rotation (such as the location of the measured object). is the signed distance from the origin to the plane.

[0231] The plane rotates around the axis, any point on the axis The perpendicular signed distances to each plane should be equal, set , then the directed distance satisfies the equation system:

[0232] ;

[0233] Writing the system of equations in matrix form:

[0234] ;

[0235] set up:

[0236] ,but:

[0237] ; Multiply both sides by ,have to:

[0238] ;

[0239] because is any point on the axis, so this equation is the axis of rotation with respect to the parameter The equation of the line, is the linear direction of the axis of rotation, , and the axis of rotation passes through the space point .Depend on:

[0240] ;

[0241] You can also calculate the unit vector in the direction of the axis of rotation ,here is the normalization factor, .

[0242] The above solution condition is that the matrix There exists an inverse matrix.

[0243] Furthermore, assuming that the rotation axis is the z-axis, the signed distance between the plane and the origin is +1, and the angle between the plane and the z-axis is 30 degrees (normal upward), the four plane equations of the plane rotating around the axis are:

[0244] , then let:

[0245] , then:

[0246] ;

[0247] Then the coordinates of any point on the axis line should meet the following requirements:

[0248] ;

[0249] Therefore, the unit vector of the axis of rotation is: , passing through the space point .

[0250] Now, give each coefficient of the plane equation a random perturbation, the coefficient The perturbation size is a uniformly distributed random number (-1,1) multiplied by a percentage, the directed distance The perturbation size is The uniformly distributed random number is multiplied by a percentage. The vector error refers to the spatial distance between the vector solution value and the ideal value.

[0251] From the above content, we can know that the axis of rotation can be solved by using the plane equation method. The solution error is proportional to the measurement error of the plane equation coefficient (proportional to when solving the inverse matrix). ). If the linear equation of the axis is solved by fitting a circle with the same name points, the solution error is proportional to the measurement error of the same name points. However, the coefficient of the plane equation is the result of fitting a multi-point point cloud, and the error is generally much smaller than the measurement error of a single spatial point. Therefore, the accuracy of calculating the axis using the plane equation method is generally higher. The condition for this method to be solved is that the matrix There exists an inverse matrix.

[0252] When the angle between the plane and the axis is 0 degrees and the plane is parallel to the axis, The rank of is 2, and the inverse matrix cannot be found. In this case, the rotation axis cannot be solved by this method. When the angle between the plane and the rotation axis is 90 degrees, and the plane is perpendicular to the rotation axis, The rank of is 1, the inverse matrix cannot be found, and the rotation axis cannot be solved by this method.

[0253] Figure 2 FIG. 1 is a schematic diagram of the structure of a spatial rotation axis calibration system for 3D level measurement according to an embodiment of the present invention. Figure 2 As shown, the system comprises:

[0254] The first unit is used to install a polyhedron calibration target on the rotation axis to be calibrated, wherein the polyhedron calibration target includes at least four feature planes that are not coplanarly distributed, and each feature plane is provided with at least three coded feature points that are not colinearly distributed; a plurality of high-resolution industrial cameras are arranged to ensure that at least two cameras simultaneously observe at least one feature plane of the polyhedron calibration target during the rotation of the rotation axis; the intrinsic and extrinsic parameters of the camera array are calibrated to establish a global coordinate system of the camera array; the rotation axis is driven to rotate one circle at a non-uniform angle interval, the rotation interval is determined according to the golden section ratio, and a multi-view image sequence of the polyhedron calibration target is collected at each angle position;

[0255] The second unit is used to perform adaptive threshold segmentation and sub-pixel edge extraction on the acquired image sequence, identify the encoded feature points, and extract the two-dimensional image coordinates of each feature point; reconstruct the two-dimensional image coordinates into high-precision three-dimensional space coordinates using the multi-view geometry principle and camera array parameters; based on the reconstructed three-dimensional space coordinates, construct a feature point trajectory fitting optimization model that considers gravity deformation and temperature drift, and the model takes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the gravity deformation parameters and the temperature drift coefficient as parameters to be optimized; use the robust Levenberg-Marquardt algorithm to solve the optimization model to obtain the optimal rotation axis parameters; use the optimal rotation axis parameters to establish a dynamically compensated rotation axis coordinate system;

[0256] The third unit is used to calculate the six-degree-of-freedom posture transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotating axis coordinate system at each angular position; based on the posture transformation matrix sequence, the actual rotation angle, axial play, radial runout and tilt error of the rotating axis are extracted; the wavelet analysis method is used to perform time-frequency analysis on the extracted actual rotation angle, axial play, radial runout and tilt error of the rotating axis to identify periodic error components and random error components; based on the periodic error components, a rotating axis geometric error model considering the mechanical structure characteristics is established; the random error components are modeled using a neural network algorithm to construct an adaptive error prediction model; the geometric error model and the adaptive error prediction model are combined to form a comprehensive error compensation model; in actual 3D position measurement, the comprehensive error compensation model is used to correct the coordinate transformation in real time according to the current angle of the rotating axis and environmental parameters.

[0257] According to a third aspect of the embodiments of the present invention,

[0258] An electronic device is provided, comprising:

[0259] processor;

[0260] a memory for storing processor-executable instructions;

[0261] The processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0262] A fourth aspect of the embodiments of the present invention is:

[0263] A computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the aforementioned method is implemented.

[0264] The present invention may be a method, an apparatus, a system and / or a computer program product. The computer program product may include a computer-readable storage medium carrying computer-readable program instructions for executing various aspects of the present invention.

[0265] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A spatial rotation axis calibration method for 3D level measurement, characterized in that: include: A polyhedron calibration target is installed on the rotation axis to be calibrated, wherein the polyhedron calibration target includes at least four feature planes that are not coplanarly distributed, and each feature plane is provided with at least three coded feature points that are not colinearly distributed; multiple high-resolution industrial cameras are arranged to ensure that at least two cameras simultaneously observe at least one feature plane of the polyhedron calibration target during the rotation of the rotation axis; the intrinsic and extrinsic parameters of the camera array are calibrated to establish a global coordinate system of the camera array; the rotation axis is driven to rotate one circle at a non-uniform angle interval, and the rotation interval is determined according to the golden section ratio, and a multi-view image sequence of the polyhedron calibration target is collected at each angle position; Adaptive threshold segmentation and sub-pixel edge extraction are performed on the collected image sequence to identify the coded feature points and extract the two-dimensional image coordinates of each feature point; the two-dimensional image coordinates are reconstructed into high-precision three-dimensional space coordinates using the multi-view geometry principle and camera array parameters; based on the reconstructed three-dimensional space coordinates, a feature point trajectory fitting optimization model considering gravity deformation and temperature drift is constructed, and the model takes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the gravity deformation parameters and the temperature drift coefficient as parameters to be optimized; the optimization model is solved using the robust Levenberg-Marquardt algorithm to obtain the optimal rotation axis parameters; using the optimal rotation axis parameters, a dynamically compensated rotation axis coordinate system is established; Calculate the six-degree-of-freedom posture transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotating axis coordinate system at each angular position; extract the actual rotation angle, axial play, radial runout and tilt error of the rotating axis based on the posture transformation matrix sequence; use the wavelet analysis method to perform time-frequency analysis on the extracted actual rotation angle, axial play, radial runout and tilt error of the rotating axis to identify periodic error components and random error components; based on the periodic error components, establish a rotating axis geometric error model that takes into account the mechanical structure characteristics; use a neural network algorithm to model the random error components and construct an adaptive error prediction model; combine the geometric error model and the adaptive error prediction model to form a comprehensive error compensation model; in actual 3D position measurement, use the comprehensive error compensation model to correct the coordinate transformation in real time according to the current angle of the rotating axis and environmental parameters; A rotating axis geometric error model considering the mechanical structure characteristics is established. The model includes: Fourier series are used to describe the periodic error components of axial play, radial runout and tilt errors. The Bayesian information criterion is used to determine the optimal harmonic order to balance model complexity and fitting accuracy. Construct an adaptive error prediction model, which includes: Use long short-term memory network to model random error components; use current angle, environmental parameters and historical error data as network input; use multi-layer long short-term memory network structure, including dropout layer to improve generalization ability; use mean square error as loss function, and use Adam algorithm for optimization; A comprehensive error compensation model is formed, which includes: Combine the outputs of the geometric error model and the adaptive error prediction model; generate a comprehensive error estimate taking into account the parameters of current angle, temperature, and humidity; Correction of coordinate transformation in real time during actual 3D level measurement. The process includes: The ideal coordinate transformation matrix is ​​calculated according to the current angle of the rotation axis; the error compensation matrix is ​​calculated using the comprehensive error compensation model; the ideal coordinate transformation matrix is ​​multiplied by the error compensation matrix to obtain a corrected coordinate transformation matrix; and the coordinate transformation of the measurement point is performed using the corrected coordinate transformation matrix.

2. The method according to claim 1, characterized in that: Adaptive threshold segmentation and sub-pixel edge extraction are performed on the acquired image sequence to identify the coded feature points and extract the two-dimensional image coordinates of each feature point; using the multi-view geometry principle and camera array parameters, the two-dimensional image coordinates are reconstructed into high-precision three-dimensional space coordinates, including: Performing adaptive threshold segmentation based on multi-scale analysis on the acquired multi-view image sequence, including performing wavelet decomposition on each image to obtain sub-images of different scales, applying the Otsu algorithm to each sub-image to obtain a local optimal threshold, obtaining a global optimal threshold by weighted fusion of thresholds of each scale, and using the global optimal threshold to segment the original image; Improved Zernike moment sub-pixel edge extraction is performed on the segmented image, including using the Canny algorithm to obtain the initial edge, constructing a circular window of variable size in the neighborhood of the edge point, calculating the Zernike moment up to the 7th order in the window, fitting the edge curve equation using the least squares method, and solving the equation to obtain the sub-pixel edge coordinates; Identify the coded feature points of the two-dimensional code mark designed based on the pseudo-random sequence in the extracted edges, the identification process includes detecting the positioning pattern using Hough transform, determining the two-dimensional code area based on the positioning pattern, extracting the code element using adaptive binarization, applying perspective transform to correct the distortion, decoding and verifying to obtain the unique identification of the feature point, and extracting the two-dimensional image coordinates of each coded feature point; Reconstructing the two-dimensional image coordinates into three-dimensional space coordinates by using an improved bundle adjustment algorithm, wherein the improved bundle adjustment algorithm comprises using a Schur complement method to decompose the optimization problem into two sub-problems of camera parameters and feature point coordinates, applying a preconditioned conjugate gradient method to solve the camera parameter sub-problem, and using a sparse Cholesky decomposition and parallel calculation to solve the feature point coordinate sub-problem; Establish a distortion correction model, including using the Brown model to describe radial distortion and tangential distortion, and using thermal imager data to establish a dynamic compensation model for the effect of temperature on camera intrinsic parameters; Substituting the reconstructed three-dimensional space coordinates into the distortion correction model to obtain corrected high-precision three-dimensional space coordinates; The design of the two-dimensional code mark includes generating an M sequence with a length of 63 as a basic code, rearranging the M sequence into a 9x7 matrix, adding positioning patterns at the four corners of the matrix, and using the BCH (63, 51) code to perform error correction encoding on the data.

3. The method according to claim 2, characterized in that Reconstructing the two-dimensional image coordinates into three-dimensional space coordinates using an improved bundle adjustment algorithm includes: The multi-view 3D reconstruction problem is modeled as a nonlinear least squares problem, with the objective function being the sum of squares of the reprojection error. The SIFT feature matching and PnP algorithm are used to obtain the initial camera pose and feature point 3D coordinate estimation. The Schur complement method is used to decompose the optimization problem into camera parameter sub-problems and feature point coordinate sub-problems, and the camera parameter increment is obtained using the Schur complement formula. The preconditioned conjugate gradient method is applied to solve the camera parameter subproblem, including: Construct a preconditioner matrix, using incomplete Cholesky decomposition as a preconditioner; Initialize residual and search direction; Iterate the solution until convergence or the maximum number of iterations is reached; The sparse Cholesky decomposition and parallel calculation are used for the feature point coordinate sub-problem, including: Construct the Hessian matrix of the feature point coordinates; Perform symbolic analysis to determine filling patterns and elimination order; Use multithreading to parallelize numerical decomposition; Use the decomposition results to solve the feature point coordinate increments; Combine the camera parameter increment and feature point coordinate increment to update the current estimate; Introducing an adaptive damping factor strategy, including: Initialize the damping factor; Calculate the actual reduction amount for the current iteration; Dynamically adjust the damping factor according to the actual reduction amount; Determine whether to accept the result of this iteration. If so, update the estimated value. Otherwise, increase the damping factor and solve again.

4. The method according to claim 1, characterized in that: Based on the reconstructed three-dimensional spatial coordinates, a feature point trajectory fitting optimization model considering gravity deformation and temperature drift is constructed. The model takes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the gravity deformation parameters and the temperature drift coefficient as parameters to be optimized; the robust Levenberg-Marquardt algorithm is used to solve the optimization model to obtain the optimal rotation axis parameters; using the optimal rotation axis parameters, the dynamic compensation rotation axis coordinate system is established, including: Based on the reconstructed three-dimensional space coordinates, a feature point trajectory fitting optimization model is constructed. The model takes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the gravity deformation parameter and the temperature drift coefficient as parameters to be optimized. Define the ideal position calculation formula of the feature point at time t, which includes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the initial position and angular velocity of the feature point; A gravity deformation correction term is introduced, which is a function of the gravity deformation parameter, the gravity acceleration vector and the angular velocity; at the same time, a temperature drift correction term is introduced, which is a function of the temperature drift coefficient, the current temperature and the initial time; Based on the ideal position calculation formula, gravity deformation correction term and temperature drift correction term, define the actual position calculation formula of the feature point at time t; Construct the objective function, which is the sum of squared errors between the observed positions of all feature points at all times and the calculated actual positions; The robust Levenberg-Marquardt algorithm is used to solve the optimization model, including: initializing parameters and damping factors, calculating the current residual and Jacobian matrix, defining weights using the Huber loss function and calculating weighted matrix products, solving incremental equations and updating parameters, dynamically adjusting the damping factor according to changes in the objective function value, and iterating the optimization until convergence or the maximum number of iterations is reached; Based on the optimal rotation axis parameters, a dynamically compensated rotation axis coordinate system is established, including: defining the rotation axis direction as the Z axis, calculating the X axis perpendicular to the Z axis and the gravity direction, obtaining the Y axis by cross-producting the X axis and the Z axis, defining the origin considering the influence of gravity deformation and temperature drift, and finally constructing the rotation matrix and translation vector; The constructed rotation matrix and translation vector are used to transform the coordinates of any point in the world coordinate system into the dynamic compensation rotation axis coordinate system.

5. The method according to claim 1, characterized in that Calculating a six-degree-of-freedom pose transformation matrix of the polyhedron calibration target relative to the rotation axis coordinate system of the dynamic compensation at each angular position; Based on the posture transformation matrix sequence, the actual rotation angle, axial play, radial runout and tilt error of the rotating axis are extracted, including: Calculate the six-degree-of-freedom pose transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotation axis coordinate system at each angular position. The calculation process includes: Obtaining the coordinates of the feature points on the polyhedron calibration target in the dynamically compensated rotation axis coordinate system and the target coordinate system; Calculate the centroid of two sets of coordinate points and obtain the centroid coordinates; Construct the covariance matrix and perform singular value decomposition; Calculate the rotation matrix based on the singular value decomposition result, and determine whether correction is needed based on the determinant of the rotation matrix; Calculate the translation vector and finally get the six-degree-of-freedom posture transformation matrix; Based on the pose transformation matrix sequence, the error parameters of the rotation axis are extracted, including: Calculate the actual rotation angle through the rotation matrix in two adjacent posture transformation matrices; Compare the components of the translation vectors of two adjacent positions in the Z-axis direction to obtain the axial movement; Calculate the displacement of the translation vector on the XY plane to obtain the radial runout; Analyze the change of the third column of the rotation matrix with the angle, use the least squares method to fit the conical surface, and solve the direction of the average rotation axis; The tilt error is obtained by calculating the angle between the average rotation axis direction and the ideal Z axis.

6. A spatial rotation axis calibration system for 3D level measurement, used to implement the method described in any one of claims 1 to 5, characterized in that: include: The first unit is used to install a polyhedron calibration target on the rotation axis to be calibrated, wherein the polyhedron calibration target includes at least four feature planes that are not coplanarly distributed, and each feature plane is provided with at least three coded feature points that are not colinearly distributed; multiple high-resolution industrial cameras are arranged to ensure that at least two cameras simultaneously observe at least one feature plane of the polyhedron calibration target during the rotation of the rotation axis; the intrinsic and extrinsic parameters of the camera array are calibrated to establish a global coordinate system of the camera array; the rotation axis is driven to rotate one circle at a non-uniform angle interval, the rotation interval is determined according to the golden section ratio, and a multi-view image sequence of the polyhedron calibration target is collected at each angle position; The second unit is used to perform adaptive threshold segmentation and sub-pixel edge extraction on the acquired image sequence, identify the encoded feature points, and extract the two-dimensional image coordinates of each feature point; reconstruct the two-dimensional image coordinates into high-precision three-dimensional space coordinates using the multi-view geometry principle and camera array parameters; based on the reconstructed three-dimensional space coordinates, construct a feature point trajectory fitting optimization model that considers gravity deformation and temperature drift, and the model takes the rotation axis direction vector, the three-dimensional coordinates of a point on the rotation axis, the gravity deformation parameters and the temperature drift coefficient as parameters to be optimized; use the robust Levenberg-Marquardt algorithm to solve the optimization model to obtain the optimal rotation axis parameters; use the optimal rotation axis parameters to establish a dynamically compensated rotation axis coordinate system; The third unit is used to calculate the six-degree-of-freedom posture transformation matrix of the polyhedron calibration target relative to the dynamically compensated rotating axis coordinate system at each angular position; based on the posture transformation matrix sequence, the actual rotation angle, axial play, radial runout and tilt error of the rotating axis are extracted; the wavelet analysis method is used to perform time-frequency analysis on the extracted actual rotation angle, axial play, radial runout and tilt error of the rotating axis to identify periodic error components and random error components; based on the periodic error components, a rotating axis geometric error model considering the mechanical structure characteristics is established; the random error components are modeled using a neural network algorithm to construct an adaptive error prediction model; the geometric error model and the adaptive error prediction model are combined to form a comprehensive error compensation model; in actual 3D position measurement, the comprehensive error compensation model is used to correct the coordinate transformation in real time according to the current angle of the rotating axis and environmental parameters.

7. An electronic device, characterized in that: include: processor; a memory for storing processor-executable instructions; The processor is configured to call the instructions stored in the memory to execute the method described in any one of claims 1 to 5.

8. A computer-readable storage medium having computer program instructions stored thereon, characterized in that: When the computer program instructions are executed by a processor, the method according to any one of claims 1 to 5 is implemented.

Images