Joint calibration method for double laser radars of mobile laser scanning system

Abstract

The invention provides a mobile laser scanning system dual-laser radar joint calibration method, a dual-laser radar mobile laser scanning system comprises a master laser radar and a slave laser radar, and systematic elimination of errors of the mobile laser scanning system is realized through a staged calibration strategy. In the first stage, a main laser radar is calibrated independently, a constraint equation is constructed through homonymy points of point cloud matching between air routes, and a placement angle parameter of the main laser radar is calibrated; and in the second stage, the calibration result of the master laser radar is taken as a reference, six-degree-of-freedom parameter collaborative calibration is carried out on three lever arm offsets and three placement angles of the slave laser radar, and finally the joint calibration of the dual laser radars is realized. According to the method, a specific calibration object, a professional calibration field and the like are not needed, the calibration process is fully automatically completed, small errors between different laser radar data in a single laser radar independent calibration method can be effectively eliminated, and the precision of the mobile laser scanning system is improved.

Description

Technical Field

[0001] This invention relates to the field of lidar installation parameter calibration and fully automatic calibration technology for mobile laser scanning systems, specifically to a method for joint calibration of dual lidars in a mobile laser scanning system. Background Technology

[0002] Mobile LiDAR Scanning (MLS) is a technology that uses laser scanners mounted on mobile platforms such as vehicles to acquire 3D point cloud data of the ground and objects. Due to its efficient 3D geospatial data acquisition capabilities, mobile LiDAR scanning systems have become an important technological means for geospatial information acquisition. Dual-scanner mobile LiDAR scanning systems, with their wider scanning angle, higher point cloud density, and stronger geometric stability, are widely used in intelligent transportation and road asset management, urban 3D modeling, autonomous driving, disaster monitoring, and environmental surveys. Mobile LiDAR scanning systems integrate sensors such as LiDAR, Global Navigation Satellite System (GNSS), and Inertial Measurement Unit (IMU) into the vehicle platform, which are rigidly mounted on the system. During installation, the origins of the LiDAR and IMU coordinate systems are almost never the same; this offset is called the lever arm offset. Simultaneously, the coordinate axes of the LiDAR and IMU are also difficult to perfectly align; the deviation angle between the coordinate axes is called the boresight angle. The presence of lever arm offset and placement angle significantly reduces the accuracy of data acquired by a mobile laser scanning system. Therefore, it is necessary to calibrate the lever arm offset and placement angle, i.e., to calibrate the placement parameters. Lever arm offset can usually be obtained through direct measurement or from design drawings. However, the placement angle cannot be directly measured and requires manual adjustment or calibration using algorithms. Manual adjustment is a time-consuming, labor-intensive method that highly depends on the operator's skills and experience, making it difficult to guarantee calibration accuracy. Using automated algorithms for placement angle calibration, on the other hand, achieves fully automatic parameter calibration without human intervention, significantly improving calibration efficiency and accuracy.

[0003] Researchers have proposed various methods for calibrating single-lidar radars in mobile laser scanning systems. These methods can be broadly categorized into non-rigorous and rigorous approaches. Non-rigorous methods estimate rotation and translation parameters between overlapping LiDAR strips by using point- or feature-based approaches that consider reference data information. They directly extract corresponding points from the overlapping areas of the strips and then estimate the rotation matrix R and translation vector t by iteratively minimizing the differences between corresponding features of the overlapping LiDAR strips. Latypov, in his 2002 paper "Estimating relative lidar accuracy information from overlapping flight lines," proposed a plane-based method to evaluate the differences between overlapping LiDAR strips, obtaining transformation parameters by comparing LiDAR data from the overlapping areas without relying on ground control points. Rigorous methods consider the causes of system errors by rigorously modeling point cloud errors by taking into account placement parameters and trajectories. In their 2006 paper, "Rigorous approach to bore-sight self-calibration in airborne laser scanning," Skaloud and Lichti proposed modeling system placement parameters using direct geolocation orientation equations. They then calibrated the placement parameters by jointly adjusting the planar parameters and placement parameters using plane constraints. However, this method is only suitable for urban areas and regions with abundant surfaces, such as man-made structures. In their 2019 paper, "Automatic data selection and boresight adjustment of LiDAR systems," Keyetieu and Seube proposed an automatic placement angle calibration method (LIBAC). This method uses overlapping measurement strips with simple linear distributions and regular slopes as input. By constructing a placement angle error observability criterion, it automatically filters points sensitive to placement angle errors and then completes the placement angle calibration based on these points. In his 2023 publication, "Research on GNSS / IMU / Laser / Vision Fusion SLAM Theory and Method", Yan Chao proposed a calibration method for the arm offset and placement angle between a GNSS / IMU and a multi-line lidar based on a specially designed calibration board. The method involves matching the positions of the four center points acquired by the lidar with the positions of the four center points measured by the total station, thereby solving for the arm offset and placement angle.

[0004] In a dual-LiDAR mobile laser scanning system, the system placement parameters can be calibrated by independently calibrating each LiDAR. However, this method has some drawbacks: (1) Independent calibration ignores the relative geometric constraints between the two LiDARs, which can easily lead to offsets or misalignments in the overlapping area during point cloud fusion, resulting in a lack of global consistency within the system. (2) Independent calibration requires repeated operation procedures, which is not only inefficient but also involves a lot of human intervention, easily introducing operational errors. In contrast, joint calibration of the two LiDARs makes full use of the constraints between the scanning data of different LiDARs. Furthermore, joint calibration of the placement parameters of the LiDARs can simultaneously optimize the arm offset and placement angle parameters of the two LiDARs within a unified framework, making full use of common-view information, suppressing error propagation, and improving calibration efficiency and system stability. Therefore, compared with independent calibration, joint calibration can better ensure the geometric consistency of the dual-LiDAR system and improve the quality of the acquired three-dimensional geospatial data. In their 2013 paper, "Pairwise LIDAR Calibration Using Multi-Type 3D Geometric Features in Natural Scene," Mengwen He et al. proposed a pairwise calibration algorithm based on multi-type geometric features. This algorithm first calibrates one LiDAR, then extracts features such as points, lines, surfaces, and secondary arrays from the point cloud of each LiDAR. By matching these multi-type features, the transformation parameters between the two sensors are estimated. This method requires sufficiently rich and reasonably distributed geometric features and is not suitable for scenes with simple structures and uniform textures (such as outdoor roads or open grasslands). In their 2020 paper, "Robust extrinsic calibration for arbitrarily configured dual 3D LiDARs using a single planar board," Kim et al. proposed a dual LiDAR relative placement parameter calibration algorithm. Its core is to use a single planar board with reflective strips and measurement data from three different orientations to achieve accurate estimation of the relative pose between the two LiDARs through geometric calculations. This method uses multi-line 3D LiDARs and is not applicable to single-line 2D LiDARs. Summary of the Invention

[0005] To address these issues, this invention proposes a joint calibration method for dual lidars in a mobile laser scanning system. This method requires no specific calibration objects or specialized calibration sites, and the calibration process is fully automated. It effectively eliminates the minor errors between different lidar data present in independent calibration methods for individual lidars, thereby improving the accuracy of the mobile laser scanning system.

[0006] This invention addresses the joint calibration problem of dual lidar systems in mobile laser scanning systems. It proposes a joint calibration method for dual lidar systems in mobile laser scanning systems that employs a master-slave progressive optimization algorithm. The dual lidar mobile laser scanning system includes a master lidar and a slave lidar. The method of this invention achieves systematic elimination of errors in the mobile laser scanning system through a staged calibration strategy.

[0007] The dual-lidar joint calibration method for a mobile laser scanning system proposed in this invention includes the following steps:

[0008] Step S1: In the first stage, the main lidar is independently calibrated. Constraint equations are constructed by matching corresponding points along different flight paths. The installation angle parameters of the main lidar are then solved by minimizing the distance between corresponding points along different flight paths. Specific steps include:

[0009] Step S11, Direct geolocation and orientation of the main lidar: Based on the initial placement parameters of the main lidar, the original point cloud data obtained after scanning by the main lidar, the POS data provided by GNSS and IMU, and the direct geolocation and orientation equation of the main lidar, the point cloud of the main lidar in the world mapping coordinate system is obtained.

[0010] Step S12, main lidar data preprocessing: preprocess the point cloud data in the world mapping coordinate system obtained after scanning by the main lidar; the data preprocessing steps include point cloud resampling and denoising, point cloud segmentation, and obtaining the route number of each point in the main lidar point cloud data based on POS data.

[0011] Step S13, Matching corresponding points of the main lidar: For the block-shaped point cloud obtained in step S12, the VGICP point cloud registration algorithm is used to match corresponding points of the point clouds of different routes within the same block of the main lidar.

[0012] Step S14, Main LiDAR installation parameter calibration: Based on the pairs of corresponding points between different routes of the main LiDAR matched in step S13, the installation angle parameters are solved by minimizing the distance between the corresponding points using the LM least squares optimization algorithm.

[0013] Step S15: Repeat steps S11 to S14 until the average change in the point-to-surface distance between all corresponding points before and after calibration is less than the preset threshold, and obtain the optimal placement angle parameters of the main lidar.

[0014] Step S2, the second stage, uses the point cloud data of the calibrated main lidar as a benchmark to perform 6-DOF parameter collaborative calibration of the offsets of the three arms and the three mounting angles of the slave lidar, ultimately achieving joint calibration of the two lidars; the specific steps include:

[0015] Step S21, direct geolocation and orientation of the main lidar and the slave lidar: Based on the original point cloud data obtained after scanning by the main lidar, the original point cloud data obtained after scanning by the slave lidar, the POS data provided by GNSS and IMU, the initial placement parameters of the slave lidar, and the calibration results of the placement parameters of the main lidar, the point cloud of the main lidar and the slave lidar in the world mapping coordinate system is obtained through direct geolocation and orientation.

[0016] Step S22, Data preprocessing of main lidar and slave lidar: The point cloud data in the world mapping coordinate system obtained after scanning by the main lidar and slave lidar is preprocessed. The data preprocessing steps include point cloud resampling and denoising, point cloud segmentation, and obtaining the route number of each point in the slave lidar point cloud data based on POS data.

[0017] Step S23, Matching corresponding points between the main lidar and the slave lidar: For the block-based point cloud obtained in step S22, match corresponding points of each route in the slave lidar point cloud with the main lidar point cloud in the same block to obtain corresponding point pairs between the main lidar point cloud and the slave lidar point cloud.

[0018] Step S24, calibration of lidar placement parameters: by minimizing the distance between corresponding points obtained in step S23, the LM least squares optimization algorithm is used to solve for the lever offset parameters and placement angle parameters of the lidar.

[0019] Step S25: Repeat steps S21 to S24 until the average change in the point-to-surface distance between all corresponding points before and after calibration is less than a preset threshold, and obtain the optimal arm offset parameters and placement angle parameters from the lidar.

[0020] Furthermore, the specific steps of step S11 include:

[0021] Step S111: Calculate the coordinates of the scanning point in the main lidar coordinate system using the distance and direction angle recorded during lidar scanning. The calculation formula is shown in equation (1) below: (1) in, For scan points Spatial rectangular coordinates in the main lidar coordinate system For scanning distance, It is the direction angle. This is the transpose of the coordinate vector;

[0022] Step S112, the dual-LiDAR mobile laser scanning system includes a LiDAR coordinate system, an IMU coordinate system, a local horizontal coordinate system, and a world mapping coordinate system, based on the scanning points. The relationship between the main lidar and various coordinate systems is established, and the direct geolocation orientation equation is shown in equation (2) below: (2) in, Main lidar scanning point Coordinates in the world cartographic coordinate system This is the rotation matrix between the local horizontal coordinate system and the world cartographic coordinate system. Latitude, longitude, and geodetic height are provided for GNSS / IMU, where W is the world cartographic coordinate system; This is the rotation matrix between the IMU coordinate system and the local horizontal coordinate system. Roll angle, pitch angle, and yaw angle provided for GNSS / IMU; The rotation matrix from the master lidar coordinate system to the IMU coordinate system. The placement angle between the main lidar and the IMU system. For scan points Spatial rectangular coordinates in the main lidar coordinate system; This is the offset vector of the main lidar in the IMU coordinate system. , , The lever arm offset between the main lidar and the IMU system; This represents the position of the IMU coordinate system origin in the world cartographic coordinate system.

[0023] Step S113: Solve the installation angle of the main lidar by minimizing the distance between corresponding points on different routes in the main lidar point cloud data. Obtain the initial installation angle of the main lidar from the system design drawings, and calculate the rotation matrix from the main lidar coordinate system to the IMU coordinate system using the initial installation angle. The calculation formula is shown in the following formula (3): (3) Due to the existence of placement angle error, it is necessary to... Perform calibration, and record the calibration amount as follows: , , The error caused by the placement angle is corrected by the calculated corrected rotation matrix; the corrected rotation matrix is ​​shown in equation (4) below: (4) in, This represents the corrected rotation matrix. This represents the rotation matrix correction amount. , , These represent the installation angles of the main lidar. The correction amount, that is, the placement angle parameter that the main lidar needs to be calibrated;

[0024] Substituting the corrected rotation matrix shown in equation (4) into the direct geolocation orientation equation of the main lidar shown in equation (2), we obtain the corrected direct geolocation orientation equation of the main lidar, as shown in equation (5) below: (5) in, This represents the corrected direct geolocation orientation equation of the main lidar.

[0025] Furthermore, the specific steps of step S13 include:

[0026] Step S131, using and These represent the point cloud set to be registered and the reference point cloud set, respectively. Given the number of points in the point cloud to be registered and the reference point cloud; perform a nearest neighbor search to make... ,in, Indicates the registration point. Indicates the registration point The nearest reference point , ; This is the transformation matrix between the point cloud to be registered and the reference point cloud;

[0027] Assuming the surface containing the laser point cloud follows a Gaussian distribution, i.e. , It follows a Gaussian distribution. , Points and points The mean of the Gaussian distribution, , Points and points The variance of the Gaussian distribution is given; then the registration error is defined as: (9) in, Mean with the mean The error between;

[0028] Based on the properties of the Gaussian distribution, we know that... It follows a Gaussian distribution as follows: (10) in, Point With point The error between;

[0029] Step S132, calculate the registration point. Rather than in the reference point cloud set The distance in is less than neighborhood points The sum of the distances between: (12) in, Representing neighborhood points, Representing neighborhood points The mean of the Gaussian distribution, Point The mean of the Gaussian distribution Rather than in the reference point cloud set The distance in is less than The mean of the Gaussian distribution of the neighborhood points The sum of the distances between them Indicates the neighborhood radius;

[0030] Indicates the registration point Rather than in the reference point cloud set The sum of the distances between points in the neighborhood of r is less than the sum of the distances between those points. The distribution is as follows: (13) (14) (15) in, Indicates error The mean of the Gaussian distribution, Indicates error The variance of the Gaussian distribution, This represents the variance of the Gaussian distribution of the j-th neighborhood point. Point The variance of the Gaussian distribution;

[0031] Step S133: Obtain the optimal transformation matrix by maximizing the logarithm of the likelihood probability. : (16)

[0032] To calculate equation (16) more efficiently, equation (16) is transformed into: (17) in, For point In reference point cloud set The number of neighboring points in the array;

[0033] Step S134: Use the optimal transformation matrix for registration parameters. Registration point cloud set The points in the cloud are transformed, and for each transformed point, a reference point cloud is created. Find its nearest neighbor and take the nearest neighbor pair with a distance less than the threshold as the same point.

[0034] Furthermore, the specific steps of step S14 include:

[0035] Step S141, define the same-name point pairs in the point clouds of different routes of the main lidar as ( , ),in and If the index is a point, the placement angle of the main lidar is estimated as follows: (18) in, This represents the parameter vector of the main lidar. , The ) indicates the calibrated installation angle of the main lidar. This indicates the angular rotation of the X-axis of the calibrated master lidar coordinate system relative to the X-axis of the IMU coordinate system. This indicates the angular rotation of the Y-axis of the calibrated master lidar coordinate system relative to the Y-axis of the IMU coordinate system. N represents the angular rotation of the Z-axis of the calibrated master lidar coordinate system relative to the Z-axis of the IMU coordinate system, and N1 is the total number of pairs of points with the same name in the master lidar point cloud data. Let be the residual function of the o-th pair of corresponding points, defined as: (19) in, Point Time The vector, For point The normal vector of the plane; Represents the residual function. This represents the vector dot product, therefore Actual representation point Time The distance from a point on the plane to the plane; (20) in, Point Time ;

[0036] Step S142, using the residual function of corresponding point pairs in the main lidar point cloud data Parameter vector of the main lidar Differentiate to construct the Jacobian matrix: (twenty one) in, Representing the residual function Jacobian matrix, For the derivative sign;

[0037] Step S143: According to formula (5), the corrected direct geolocation orientation equation of the main lidar can be derived. The derivative is: (twenty two) Due to the rotation matrix correction It is a small-angle rotation matrix, which can be approximated as: (twenty three) in, It is the identity matrix;

[0038] Therefore, the derivative term can be linearized as follows: (twenty four) in, This represents the X-coordinate of the point cloud in the IMU coordinate system, calculated using the initial placement angle. This represents the Y-coordinate of the point cloud in the IMU coordinate system, calculated using the initial placement angle. This represents the Z-coordinate of the point cloud in the IMU coordinate system, calculated using the initial placement angle.

[0039] Step S144, final residual function The Jacobian matrix is: (25) in, Point The matrix after linearization of the derivative terms, Point The matrix after linearizing the derivative terms;

[0040] Step S145: Then, establish the incremental equation for the main lidar and introduce the damping factor. : (26) in, yes The Jacobian matrix of dimension N1, where N1 is the total number of pairs of identical points in the main lidar point cloud data. Representing the Jacobian matrix transpose, Represents a diagonal matrix; It is a parameter increment. ; Represents the residual vector of the incremental equation of the main lidar;

[0041] Step S146, the parameter update iteration formula is: (27) in, This represents the updated parameter estimate. This represents the parameter estimate at the e-th iteration;

[0042] The LM least squares optimization algorithm converges to the optimal parameters through iteration. The iteration terminates when the change in residual is lower than a preset threshold or when the maximum number of iterations is reached.

[0043] Furthermore, the specific steps of step S15 include:

[0044] Repeat steps S111 to S146 until the average change in the point-to-surface distance between all corresponding points before and after calibration is less than a preset threshold, thus obtaining the optimal placement angle parameters of the main lidar. Main lidar arm offset Obtained through direct measurement or from design drawings.

[0045] Furthermore, the specific steps of step S21 include:

[0046] From the lidar calibration, the arm offset after the main lidar calibration is calculated. and placement corner The initial lever offset is fixed and cannot be changed; it is obtained from the lidar by direct measurement or from the design drawings. and initial placement angle The initial calibration coordinates of the lidar are obtained by directly geolocating and orienting the raw point cloud data from the lidar. ;

[0047] The coordinates are obtained from the lidar point cloud after direct geolocation and orientation, expressed as: (28) in, This indicates the coordinates of the scan point from the lidar within the lidar coordinate system. This represents the offset vector from the lidar in the IMU coordinate system; This represents the rotation matrix from the lidar coordinate system to the IMU coordinate system; The small-angle rotation matrix from the lidar coordinate system to the IMU coordinate system is approximately: (29) in, It is the identity matrix. , , These represent the angles from which the lidar is installed. The correction amount.

[0048] Furthermore, the specific steps of step S24 include:

[0049] Step S241, define the corresponding point pairs between the point cloud data of the master lidar and the slave lidar as... The calibration of the lidar installation parameters is estimated as follows: (30) in, This indicates the calibration parameters of the lidar. This indicates the offset of the lever arm from the X-axis of the lidar coordinate system relative to the X-axis of the IMU coordinate system after calibration. This indicates the offset of the lever arm from the Y-axis of the lidar coordinate system relative to the Y-axis of the IMU coordinate system after calibration. This indicates the offset of the lever arm from the Z-axis of the lidar coordinate system relative to the Z-axis of the IMU coordinate system after calibration. This represents the angular rotation from the X-axis of the lidar coordinate system relative to the X-axis of the IMU coordinate system after calibration. This represents the angular rotation from the Y-axis of the lidar coordinate system relative to the Y-axis of the IMU coordinate system after calibration. This represents the angular rotation from the Z-axis of the lidar coordinate system relative to the Z-axis of the IMU coordinate system after calibration. N represents the parameter vector from the lidar, and N2 is the number of corresponding point pairs between the main lidar and the lidar point cloud data. Let be the residual function of the q-th pair of corresponding points, defined as: (31) in, Point Time The vector, For point The normal vector of the plane at that point; Represents the residual function. , These represent the correction amounts for the arm offset and mounting angle of the laser radar, respectively; that is, the arm offset and mounting angle parameters that need to be calibrated for the laser radar. Actual representation point Time The distance from a point on the plane to the plane; (32) in, Point Time ;

[0050] Step S242: The LM least squares optimization algorithm is used to estimate the optimal placement parameters, based on the residual function of corresponding point pairs between the main and slave lidar point cloud data. For the parameter vector of the lidar Differentiate to construct the Jacobian matrix: (33) in, Representing the residual function The Jacobian matrix;

[0051] Step S243, since the lidar needs to simultaneously optimize 6 degrees of freedom parameters. Therefore, the residual function Partial derivatives need to be calculated for the lever arm offset parameter and the placement angle, respectively: (34) (35) in, These are the coordinates of the lidar point in the IMU coordinate system;

[0052] Combining the two types of derivatives, namely equations (34) and (35), yields the Jacobian matrix from the lidar: (36) in, Representing the residual function The Jacobian matrix;

[0053] Step S244, introduce the damping factor Construct incremental equations from lidar: (37) in, yes The Jacobian matrix is ​​N2, where N2 is the number of identical point pairs between the main and slave lidar point cloud data. Representing the Jacobian matrix transpose, Represents a diagonal matrix; It is a parameter increment. For 6 degrees of freedom parameters The parameter increment; This represents the residual vector from the incremental equation of the lidar.

[0054] Step S245, the parameter update iteration formula is: (38) in, This represents the updated parameter estimate. This represents the parameter estimate at the e-th iteration. For 6 degrees of freedom parameters The parameter increment.

[0055] Furthermore, the specific steps of step S25 include:

[0056] Repeat steps S211 to S245 until the average change in the point-to-surface distance between all corresponding points before and after calibration is less than a preset threshold, thus obtaining the optimal arm offset parameters from the lidar. and placement angle parameters .

[0057] The technical effects of this invention are as follows:

[0058] (1) This method can automatically complete the entire calibration process based on the original point cloud data, POS data and initial placement parameters of the lidar, without the need for specific calibration objects, professional calibration sites and control points; the calibration process is fully automatic, without the need for manual intervention and human-computer interaction, saving manpower and material costs.

[0059] (2) The joint calibration of the two lidars eliminates the error between the scanning data of the two lidars, and the system calibration accuracy is high.

[0060] (3) It is highly practical and flexible. Attached Figure Description

[0061] Figure 1 This is a flowchart illustrating the overall process of the dual-lidar joint calibration method for a mobile laser scanning system in this embodiment of the invention.

[0062] Figure 2 This is a visualization of eight sets of point cloud data (MLS1-MLS8) actually collected by multiple dual-LiDAR mobile laser scanning systems in an embodiment of the present invention.

[0063] Figure 3 The image shows a cross-sectional view of the master and slave lidar scan lines in the point cloud datasets MLS1-MLS4. x-1, x-2, and x-3 (x represents a, b, c, ..., h) represent the point clouds before calibration, after calibration using the joint calibration method of this invention, and after independent calibration of the two lidars, respectively. Different colors represent point clouds collected along different flight paths.

[0064] Figure 4This is a cross-sectional view of the master and slave lidar scan lines in the point cloud dataset MLS5-MLS8. x-1, x-2, and x-3 (x represents a, b, c...h) represent the point clouds before calibration, after calibration using the joint calibration method of this invention, and after independent calibration of the two lidars, respectively. Different colors represent point clouds collected along different flight paths.

[0065] Figure 5 This is a cross-sectional view of the point clouds scanned by the master and slave lidars in the point cloud datasets MLS1-MLS4. x-1, x-2, x-3 (x is a, b, c...h) represent the point clouds before calibration, after calibration by the joint calibration method of this invention, and after independent calibration of the two lidars, respectively. The point cloud scanned by the master lidar is displayed in blue after calibration, the point cloud calibrated by the slave lidar using the joint calibration method of this invention is green, and the point cloud calibrated independently by the slave lidar is red.

[0066] Figure 6 This is a cross-sectional view of the point cloud scanned by the master and slave lidars in the point cloud datasets MLS5-MLS8. x-1, x-2, x-3 (x is a, b, c...h) represent the point cloud before calibration, after calibration by the joint calibration method of this invention, and after independent calibration of the two lidars, respectively. The point cloud scanned by the master lidar is displayed in blue after calibration, the point cloud of the slave lidar calibrated by the joint calibration method of this invention is green, and the point cloud of the slave lidar calibrated independently is red. Detailed Implementation

[0067] To better understand the above-described objects, features, and advantages of the present invention, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. Many specific details are set forth in the following description to provide a thorough understanding of the invention; however, the invention may be practiced in other ways different from those described herein, and therefore, the invention is not limited to the specific embodiments disclosed below.

[0068] Example:

[0069] A dual-LiDAR mobile laser scanning system comprises two LiDARs: a master LiDAR and a slave LiDAR. The proposed dual-LiDAR joint calibration method for this mobile laser scanning system can be broadly divided into two stages: master LiDAR installation parameter calibration and slave LiDAR installation parameter calibration. The overall calibration process is as follows: Figure 1 As shown, it includes the following steps:

[0070] Step S1 involves independently calibrating the main lidar. Constraint equations are constructed using corresponding point pairs matched between point clouds along different flight paths. The installation angle parameters of the main lidar are then solved by minimizing the distance between corresponding points along different flight paths. Specific steps include:

[0071] Step S11, Direct Geolocation and Orientation of the Main LiDAR: Based on the initial installation parameters of the main LiDAR, the raw point cloud data obtained after scanning by the main LiDAR, the POS (Position and Orientation System) data provided by GNSS (Global Navigation Satellite System) and IMU (Inertial Measurement Unit), and the direct geolocation and orientation equation of the main LiDAR, the point cloud of the main LiDAR in the world mapping coordinate system is obtained; specifically:

[0072] The dual-laser radar mobile laser scanning system uses a single-line 2D laser radar, which scans in two dimensions. During the scanning process, the distance and direction angle between the laser pulse emitted by the laser radar and the scanning point are accurately recorded. The laser radar coordinate system consists of a scanning plane and a rotation axis. The coordinates of the scanning point in the main laser radar coordinate system can be calculated using the distance and direction angle recorded during the scanning process. The calculation formula is shown in equation (1) below: (1) in, For scan points Spatial rectangular coordinates in the main lidar coordinate system For scanning distance, It is the direction angle. This is the transpose of the coordinate vector.

[0073] The dual-LiDAR mobile laser scanning system is a multi-sensor integrated system involving multiple different coordinate systems, including the LiDAR coordinate system, the Inertial Measurement Unit (IMU) coordinate system, the local horizontal coordinate system, and the world cartographic coordinate system. The transformation and correlation between these coordinate systems are crucial for achieving direct geolocation and orientation. The direct geolocation and orientation equation of the main LiDAR is: (2) in, Main lidar scanning point Coordinates in the world cartographic coordinate system This is the rotation matrix between the local horizontal coordinate system and the world cartographic coordinate system. Latitude, longitude, and geodetic height are provided for GNSS / IMU, where W is the world cartographic coordinate system; This is the rotation matrix between the IMU coordinate system and the local horizontal coordinate system. Roll angle, pitch angle, and yaw angle provided for GNSS / IMU; The rotation matrix from the master lidar coordinate system to the IMU coordinate system. The placement angle between the main lidar and the IMU system. For scan points Spatial rectangular coordinates in the main lidar coordinate system; This is the offset vector of the main lidar in the IMU coordinate system. , , The lever arm offset between the main lidar and the IMU system; This represents the position of the origin of the IMU coordinate system in the world cartographic coordinate system.

[0074] In practice, dual-LiDAR mobile laser scanning systems suffer from various errors, including LiDAR ranging and angle measurement errors, GNSS positioning errors, IMU attitude errors, boom offset errors, and placement angle errors. These errors prevent point clouds obtained from scanning in different directions and strips from perfectly overlapping. The ranging and angle measurement errors of the high-precision LiDAR itself have a negligible impact on positioning accuracy at the millimeter level; boom offset can be directly measured or obtained from design drawings, and its error is also negligible; under good GNSS conditions, GNSS positioning errors and IMU attitude errors are small, therefore the total error of the point cloud is mainly dominated by the placement angle error.

[0075] Based on this, the installation angle of the main lidar can be solved by minimizing the distance between corresponding points on different routes in the main lidar point cloud data. The initial installation angle of the main lidar can be obtained from the system design drawings. The angular rotation of the main lidar coordinate system relative to the IMU coordinate system along the X, Y, and Z axes are denoted as follows: The rotation matrix from the main lidar coordinate system to the IMU coordinate system is calculated using the initial placement angle, and the calculation formula is shown in equation (3) below: (3) Due to the presence of placement angle errors, the rotation matrix is ​​not accurate enough and needs to be adjusted. Perform calibration. The calibration amount is denoted as... , , The error caused by the placement angle is corrected by the calculated corrected rotation matrix. The corrected rotation matrix is ​​shown in equation (4) below: (4) in, This represents the corrected rotation matrix. This represents the rotation matrix correction amount. , , These represent the installation angles of the main lidar. The correction amount, that is, the placement angle parameter that the main lidar needs to be calibrated;

[0076] Substituting the corrected rotation matrix shown in equation (4) into the direct geolocation orientation equation of the main lidar shown in equation (2), we obtain the corrected direct geolocation orientation equation of the main lidar, as shown in equation (5) below: (5) in, This represents the corrected direct geolocation orientation equation of the main lidar.

[0077] Step S12, Main LiDAR Data Preprocessing: The original point cloud data in the world mapping coordinate system obtained after scanning by the main LiDAR is preprocessed; the data preprocessing methods include point cloud resampling and denoising, point cloud segmentation, and obtaining the route number of each point in the main LiDAR point cloud data based on POS data; the specific steps are as follows:

[0078] Step S121, point cloud resampling:

[0079] Due to the high density and large amount of point cloud data, the efficiency of matching corresponding points between different flight routes is low. To improve the efficiency of matching corresponding points, the point cloud data is downsampled and the point cloud density is reduced by grid sampling.

[0080] Grid sampling can be performed in two steps:

[0081] (1) Mesh generation and spatial index construction:

[0082] Based on the preset grid cell size V g Each point in the original point cloud is mapped to the corresponding grid cell according to its spatial location, and a spatial index relationship between the point cloud and the grid is established. The calculation formula is shown in the following formula (6): (6) in,( , , ) represents the index of the grid cell. , , Represents the coordinates of the point cloud. This represents the minimum X-axis coordinate of a point in the point cloud. This represents the minimum Y-axis coordinate of a point in the point cloud. This represents the minimum Z-axis coordinate of a point in the point cloud.

[0083] (2) Filtering representative points within the grid:

[0084] For each non-empty grid cell, traverse all point clouds within that grid cell, calculate the Euclidean distance between each point and the geometric center of the grid cell, and select the point closest to the geometric center as the representative point of the current grid.

[0085] Step S122, Point Cloud Denoising:

[0086] The presence of point cloud noise can affect the accuracy of calibration. Therefore, in this embodiment, intensity filtering is used to denoise the point cloud. Intensity filtering is a key step in removing aerial noise. By analyzing the statistical difference between the reflection intensity of aerial noise (such as airborne debris) and the actual ground objects, noise removal is achieved. The specific method can be described as follows:

[0087] Let the point cloud dataset be , ,in, This represents the i-th point in the point cloud dataset. This represents the coordinates of the i-th point in the point cloud dataset. Let N represent the reflection intensity of the i-th point, and N represent the total number of points in the point cloud. Noise points typically exhibit two characteristics: (1) low-intensity outliers, corresponding to weak reflection signals; and (2) high-intensity outliers, originating from abnormal reflection interference. Establish an intensity filtering model: (7) in, Indicates the intensity filtering model. This represents the reflection intensity at the i-th point. This represents the minimum value of the point cloud reflection intensity. This represents the maximum value of the point cloud reflection intensity.

[0088] Step S123, point cloud segmentation:

[0089] To accelerate the matching speed of the VGICP algorithm, the point cloud is divided into blocks. For a given block size S... block Points in the global coordinate system will be stored in a location with row number r. num and column number c num In the corresponding block, line number r num and column number c num The calculation method is shown in equation (8). During the point cloud segmentation process, it is recorded whether the point was collected by the main lidar or by the slave lidar.

[0090] (8)

[0091] Step S124, Obtaining the route number:

[0092] To obtain corresponding points among point clouds scanned from different positions and viewpoints, this invention uses trajectory data (including time, X / Y / Z axis coordinates, roll angle, pitch angle, and yaw angle, etc.) to determine the flight path number for each point. The specific determination criteria are as follows:

[0093] Time interval criterion: If the time interval between the current trajectory record and the previous record exceeds 1 second, the system will generate a new route. This is because data outside the target area may be truncated during the calibration process. The time threshold can effectively identify discontinuous trajectory segments.

[0094] Yaw angle change criteria: If the difference between the yaw angle of the current trajectory record and the starting record of the current route exceeds 90 degrees, it is determined that a turn or U-turn has occurred, and a new route is generated.

[0095] After segmenting the trajectory data based on the above rules, the system assigns a corresponding route number to each point through spatiotemporal correlation and records it. This route number will serve as a key index for point cloud matching.

[0096] Step S13, matching corresponding points between main lidar flight lines based on the VGICP point cloud registration algorithm:

[0097] For the segmented point cloud obtained in step S123, the VGICP point cloud registration algorithm is used to match corresponding points between different routes within the same block of the main lidar. The VGICP point cloud registration algorithm is described in the existing technical literature "Koide, Kenji, et al. "Voxelized GICP for fast and accurate 3D point cloud registration." 2021 IEEE international conference on robotics and automation (ICRA). IEEE, 2021".

[0098] In this embodiment, the specific method for matching corresponding point pairs between point clouds of different flight paths using the voxelized generalized iterative nearest point algorithm (i.e., the VGICP point cloud registration algorithm) is as follows: and These represent the point cloud set to be registered and the reference point cloud set, respectively. Let the number of points in the point cloud to be registered and the reference point cloud be specified. In classic ICP registration, nearest neighbor search is used to achieve... ,in, Indicates the registration point. Indicates the registration point The nearest reference point , ; This represents the transformation matrix between the point cloud to be registered and the reference point cloud. GICP (Generalized-ICP) assumes that the surface containing the laser point cloud follows a Gaussian distribution, i.e. , It follows a Gaussian distribution. , Points and points The mean of the Gaussian distribution, , Points and points The variance of the Gaussian distribution. Under this assumption, the registration error can be defined as: (9) in, Mean with the mean The error between;

[0099] Based on the properties of the Gaussian distribution, we know that... It follows a Gaussian distribution as follows: (10) in, Point With point The error between;

[0100] GICP obtains the optimal transformation matrix by maximizing the logarithm of the likelihood probability. : (11) in, Represents the logarithmic function. Indicates error The probability, Indicates error transpose, Represents the transformation matrix transpose;

[0101] The VGICP point cloud registration algorithm extends formula (9) by not only calculating the points to be registered. Its nearest reference point Instead of calculating the distance between them, we calculate the registration point. Rather than in the reference point cloud set The distance in is less than neighborhood points The sum of the distances between: (12) in, Representing neighborhood points, Representing neighborhood points The mean of the Gaussian distribution, Point The mean of the Gaussian distribution Rather than in the reference point cloud set The distance in is less than The mean of the Gaussian distribution of the neighborhood points The sum of the distances between them Indicates the neighborhood radius;

[0102] Indicates the registration point Rather than in the reference point cloud set The sum of distances between neighboring points whose distance is less than r is given by formula (12). Formula (12) can be considered as a smoothing of the distribution of the target point, similar to formula (14). The distribution is as follows: (13) (14) (15) in, Indicates error The mean of the Gaussian distribution, Indicates error The variance of the Gaussian distribution, Indicates the first The variance of the Gaussian distribution of the neighborhood points, Point The variance of the Gaussian distribution;

[0103] Finally, the optimal transformation matrix is ​​obtained by maximizing the logarithm of the likelihood probability in equation (13). : (16)

[0104] To calculate equation (16) more efficiently, equation (16) is transformed into: (17) In equation (17), For point In reference point cloud set The number of neighboring points in the array.

[0105] The registration parameters (i.e., the optimal transformation matrix) between flight lines are obtained using the VGICP point cloud registration algorithm. After that, the registration parameters are used to register the set of point clouds to be registered. The points in the cloud are transformed, and for each transformed point, a reference point cloud is created. Find its nearest neighbor and take the nearest neighbor pair with a distance less than the threshold as the same point.

[0106] Step S14, Main lidar installation parameter calibration:

[0107] After finding corresponding points between different flight paths, the placement angle parameters are solved. The LM least squares optimization algorithm is used to obtain the optimal placement angle parameters by minimizing the distance between corresponding points, thereby achieving the calibration of the main lidar. The LM least squares optimization algorithm is described in the existing technical literature "Moré, Jorge J. "The Levenberg-Marquardtalgorithm: implementation and theory." Numerical analysis: proceedings of the biennial Conference held at Dundee, June 28–July 1, 1977. Berlin, Heidelberg: Springer Berlin Heidelberg, 2006".

[0108] In this embodiment, the pairs of points with the same name between different routes of the main lidar are defined as ( , ),in and If the index is a point, the placement angle of the main lidar can be estimated as follows: (18) in, This represents the parameter vector of the main lidar. , The ) indicates the calibrated installation angle of the main lidar. This indicates the angular rotation of the X-axis of the calibrated master lidar coordinate system relative to the X-axis of the IMU coordinate system. This indicates the angular rotation of the Y-axis of the calibrated master lidar coordinate system relative to the Y-axis of the IMU coordinate system. N represents the angular rotation of the Z-axis of the calibrated master lidar coordinate system relative to the Z-axis of the IMU coordinate system, and N1 is the total number of pairs of points with the same name in the master lidar point cloud data. Let be the residual function of the o-th pair of corresponding points, defined as: (19) in, Point Time The vector, For point The normal vector of the plane. Represents the residual function. This represents the vector dot product, therefore Actual representation point Time The distance from a point on the plane to the plane.

[0109] (20) in, Point Time The vector.

[0110] Since equation (19) has obvious nonlinear characteristics, this invention uses the LM least squares optimization algorithm to estimate the optimal placement angle parameters. Specifically, this requires using the residual function of corresponding point pairs in the main lidar point cloud data. Parameter vector of the main lidar Differentiate to construct the Jacobian matrix: (twenty one) in, Representing the residual function Jacobian matrix, For the derivative sign;

[0111] According to formula (5), the corrected direct geolocation orientation equation of the main lidar can be derived. The derivative is: (twenty two) Due to the rotation matrix correction It is a small-angle rotation matrix, which can be approximated as: (twenty three) in, It is the identity matrix;

[0112] Therefore, the derivative term of equation (22) can be linearized as follows: (twenty four) in, This represents the X-coordinate of the point cloud in the IMU coordinate system, calculated using the initial placement angle. This represents the Y-coordinate of the point cloud in the IMU coordinate system, calculated using the initial placement angle. This represents the Z-coordinate of the point cloud in the IMU coordinate system, calculated using the initial placement angle.

[0113] Final residual function The Jacobian matrix is: (25) in, Point The matrix after linearization of the derivative terms, Point The matrix after linearizing the derivative terms;

[0114] Then, the incremental equation for the main lidar is established, and a damping factor is introduced. : (26) in, yes The Jacobian matrix of dimension N1, where N1 is the total number of pairs of identical points in the main lidar point cloud data. Representing the Jacobian matrix transpose, Represents a diagonal matrix; It is a parameter increment. ; This represents the residual vector of the incremental equation of the main lidar.

[0115] The parameter update and iteration formula is: (27) in, This represents the updated parameter estimate. This represents the parameter estimate at the e-th iteration;

[0116] The LM least squares optimization algorithm converges to the optimal parameters through iteration. The iteration terminates when the change in residual is lower than a preset threshold or when the maximum number of iterations is reached.

[0117] Step S15: Calculate the average point-to-surface distance of corresponding points before and after the placement angle calibration. If the difference in the average point-to-surface distance of all corresponding points is less than a preset threshold... If the calibration parameters of the main lidar installation angle converge, the estimated installation angle is considered the final optimal installation angle parameter. If the difference in the average point-to-surface distance of all corresponding points is greater than or equal to a preset threshold... The estimated main lidar placement angle parameters are then used to re-geolocate and orient the point cloud within each block. After re-geolocating and orienting block by block, corresponding points between different flight paths are re-matched. The placement angle of the main lidar is then re-estimated using the new corresponding points. This process is iterated until the change in the average point-to-area distance of all corresponding points before and after calibration is less than a preset threshold.

[0118] The installation angle of the main lidar can be determined by following the above steps. And lever arm offset The parameters can be obtained through direct measurement or from the design drawings. At this point, the calibration of the main lidar installation parameters is complete.

[0119] Step S2 involves using the calibrated point cloud data of the master lidar as a reference to perform collaborative calibration of the slave lidar's 6-DOF parameters (3 lever offsets and 3 mounting angles), thereby calibrating the slave lidar's mounting parameters and ultimately achieving joint calibration of the two lidars. Specific steps include:

[0120] Step S21, Direct Geoorientation of Master and Slave LiDAR: Based on the raw point cloud data obtained after scanning by the master lidar, the raw point cloud data obtained after scanning by the slave lidar, the POS data provided by GNSS and IMU, the initial placement parameters of the slave lidar, and the calibration results of the master lidar placement parameters, the direct geoorientation equation of the slave lidar is determined and corrected to obtain the point cloud of the master and slave lidars in the world mapping coordinate system; specifically:

[0121] From the lidar calibration, the arm offset after the main lidar calibration is calculated. and placement corner Fixed. The initial lever offset from the lidar is obtained through direct measurement or from the design drawings. and initial placement angle The initial calibration coordinates of the lidar are obtained by directly geolocating and orienting the raw point cloud data from the lidar. .

[0122] The coordinates are obtained from the lidar point cloud after direct geolocation and orientation, expressed as: (28) in, This indicates the coordinates of the scan point from the lidar within the lidar coordinate system. This represents the offset vector from the lidar in the IMU coordinate system; This represents the rotation matrix from the lidar coordinate system to the IMU coordinate system; The small-angle rotation matrix from the lidar coordinate system to the IMU coordinate system is approximately: (29) in, It is the identity matrix. , , These represent the angles from which the lidar is installed. The correction amount.

[0123] Step S22: Data preprocessing of master and slave lidar: The scan point cloud data of the master and slave lidar in the world mapping coordinate system are preprocessed. The data preprocessing steps include point cloud resampling and denoising, point cloud segmentation, and point cloud route acquisition. The principle of the preprocessing steps is the same as that of step S12.

[0124] Step S23: Matching corresponding points between the main and secondary lidar: Matching corresponding points of each route in the secondary lidar point cloud with the main lidar point cloud at the same location on the same route. The matching principle is the same as in step S13, which is to use the voxelized generalized iterative nearest point algorithm (i.e. VGICP point cloud registration algorithm) to perform block-by-block matching to obtain a set of corresponding point pairs between the main lidar and the secondary lidar.

[0125] Step S24, Calibrate the lidar installation parameters:

[0126] The VGICP point cloud registration algorithm is used to obtain corresponding point pairs after direct geolocation and orientation from the raw point cloud data of the master and slave lidars. Subsequently, the calibration of the lidar installation parameters can be estimated as follows: (30) in, This indicates the calibration parameters of the lidar. This indicates the offset of the lever arm from the X-axis of the lidar coordinate system relative to the X-axis of the IMU coordinate system after calibration. This indicates the offset of the lever arm from the Y-axis of the lidar coordinate system relative to the Y-axis of the IMU coordinate system after calibration. This indicates the offset of the lever arm from the Z-axis of the lidar coordinate system relative to the Z-axis of the IMU coordinate system after calibration. This represents the angular rotation from the X-axis of the lidar coordinate system relative to the X-axis of the IMU coordinate system after calibration. This represents the angular rotation from the Y-axis of the lidar coordinate system relative to the Y-axis of the IMU coordinate system after calibration. This represents the angular rotation from the Z-axis of the lidar coordinate system relative to the Z-axis of the IMU coordinate system after calibration. N2 represents the parameter vector from the lidar, and N2 is the number of corresponding point pairs between the main lidar and the lidar point cloud data. Let be the residual function between the corresponding point pairs in the q-th group of master and slave lidar point cloud data, defined as: (31) in, Point Time The vector, For point The normal vector of the plane. Represents the residual function. , These represent the correction amounts for the arm offset and mounting angle of the laser radar, respectively; that is, the arm offset and mounting angle parameters that need to be calibrated for the laser radar. Actual representation point Time The distance from a point on the plane to the plane.

[0127] (32) in, Point Time The vector.

[0128] The LM least squares optimization algorithm is used to estimate the optimal placement parameters. Specifically, this requires using the residual function of corresponding point pairs between the main and slave lidar point cloud data. For the parameter vector of the lidar Differentiate to construct the Jacobian matrix: (33) in, Representing the residual function The Jacobian matrix;

[0129] Because the lidar requires simultaneous optimization of 6 degrees of freedom parameters. Therefore, the residual function Partial derivatives need to be calculated for the lever arm offset parameter and the placement angle, respectively: (34) (35) in, These are the coordinates of the lidar point in the IMU coordinate system.

[0130] Combining the two types of derivatives, namely equations (34) and (35), yields the Jacobian matrix from the lidar: (36) in, Representing the residual function The Jacobian matrix.

[0131] Introducing damping factor Construct incremental equations from lidar: (37) in, yes The Jacobian matrix is ​​N2, where N2 is the number of identical point pairs between the main and slave lidar point cloud data. Representing the Jacobian matrix transpose, Represents a diagonal matrix; It is a parameter increment. For 6 degrees of freedom parameters The parameter increment; This represents the residual vector from the incremental equation of the lidar.

[0132] The parameter update and iteration formula is: (38) in, This represents the updated parameter estimate. This represents the parameter estimate at the e-th iteration. For 6 degrees of freedom parameters The parameter increment.

[0133] Step S25: Calculate the average point-to-surface distance between the main lidar and the corresponding points in the point cloud at the same location along the same route of the secondary lidar before and after joint calibration. If the difference in the average point-to-surface distance of all corresponding points is less than a preset threshold... If the calibration of the slave lidar placement parameters is successful, the estimated placement parameters are considered the final optimal placement parameters. Otherwise, the estimated slave lidar placement parameters are used to re-geolocate and orient the point cloud within each block, and after re-geolocating and orienting block by block, the corresponding points between the master and slave lidar point clouds are re-matched. The slave lidar placement parameters are then re-estimated using the new corresponding points. This process is iterated until the change in the average point-to-surface distance of all corresponding points before and after calibration is less than a preset threshold.

[0134] Compared to primary lidar calibration, secondary lidar calibration is unique in three aspects: First, the optimization parameters are expanded from three placement angles to six degrees of freedom. Second, the corresponding points originate from point clouds of both the primary and secondary lidars along the same flight path, rather than overlapping areas across flight paths, thus requiring higher point cloud matching accuracy. Third, the placement parameters after primary lidar calibration serve as ground truth inputs, making the secondary lidar calibration results directly dependent on the calibration accuracy of the primary lidar. Through this joint calibration framework, the geometric deviations of the dual lidar system are systematically eliminated, providing reliable data support for 3D reconstruction and measurement in complex scenarios.

[0135] To verify the accuracy and effectiveness of the system calibration, this embodiment selected eight sets of point cloud data (MLS1-MLS8) after independent and joint calibration of the dual-LiDAR mobile laser scanning system. Figure 2 The eight sets of urban point cloud data collected are displayed. Figures 3 to 4 The comparison of point cloud effects before calibration, joint calibration, and independent calibration is shown. Different colors represent point clouds collected from different flight paths. Among them, (a-1)~(h-1) are point clouds before calibration, (a-2)~(h-2) are point clouds after calibration by the joint calibration method of this invention, and (a-3)~(h-3) are point clouds calibrated independently. Figure 3(a-2), (b-2), (c-2), and (g-2) in the figure indicate that the joint calibration method of the present invention can effectively reduce the wall thickness; Figure 4 The values ​​(d-2), (f-2), and (h-2) in the figure show that the calibration accuracy at sharp feature locations is basically equivalent to that of independent calibration. Figure 4 (a-2) and (c-2) in the figure demonstrate that the method has stability when dealing with sharp features of heterogeneous datasets. Figure 4 The values ​​(b-2), (e-2), (f-2), and (g-2) in the figure show that there was no significant difference between the joint calibration and independent calibration results in the target point cloud calibration task. Figures 3-4 It can be seen that the joint calibration method of the present invention can effectively complete the calibration and reduce the error between the flight paths on eight different datasets.

[0136] Figures 5-6 The accuracy improvement effects of joint calibration and independent calibration on the point clouds of master and slave lidar were further compared. (a-1)~(h-1) represent the point clouds before calibration, (a-2)~(h-2) represent the point clouds after calibration using the joint calibration method of this invention, and (a-3)~(h-3) represent the point clouds after independent calibration. The point cloud scanned by the master lidar is shown in blue, the point cloud of the slave lidar after joint calibration is shown in green, and the point cloud of independent calibration is shown in red. The results show that in most cases, the difference between the two methods is not significant. Figure 5 (b-2) and (b-3) show that under specific conditions, this method can better maintain calibration stability and avoid point cloud discrepancies. Figure 6 (b-2) and (b-3) in the figure indicate that in some scenarios, joint calibration can make the point cloud spacing after master-slave lidar calibration smaller than that after independent calibration. Figures 3-6 The joint calibration method of the present invention was verified to have the following effect: while maintaining the accuracy between independent calibration lines, the joint calibration can achieve higher calibration accuracy between master and slave lidars under specific conditions.

[0137] To quantitatively verify and evaluate the effectiveness and performance of the proposed joint calibration method for dual lidar placement parameters, this embodiment performs quantitative analysis by calculating the mean (mean) and root mean square error (RMSE) of the point-to-surface distances of matched corresponding points before and after calibration. Specifically, given N3 sets of matched corresponding points, where the point-to-surface distance of the nth corresponding point is defined as... The formulas for calculating Mean and RMSE are as follows: (39) in, Indicates the average error; (40) in, This represents the root mean square error.

[0138] If the calibration is successful, the average point-to-surface distance (Mean) and root mean square error (RMSE) of matched corresponding points should be small, and the smaller the better. To better evaluate the accuracy of the point cloud before and after calibration, this embodiment calculates the Mean and RMSE for matched corresponding points between different scan lines of the same lidar and between master and slave lidar point clouds. The calculation results are shown in Tables 1 and 2. The Mean and RMSE of matched corresponding points between master and slave lidar point clouds are denoted as Mean1 and RMSE1, respectively, while the Mean and RMSE of matched corresponding points between different scan lines of the same lidar are denoted as Mean2 and RMSE2.

[0139] Table 1 shows the RMSE1 and RMSE2 values ​​before and after calibration using the independent calibration method and the joint calibration method. It can be seen that the average RMSE1 after using the joint calibration method reaches 1.686 cm, which is better than the 1.720 cm of the independent calibration method. The average RMSE2 of the two calibration methods are very close (2.490 cm for the independent method and 2.495 cm for the joint method). In the eight datasets, the joint calibration method performs better in the RMSE1 metric for master-slave LiDAR scan point clouds in six datasets (MLS2, MLS4, MLS5, MLS6, MLS7, and MLS8). In the RMSE2 metric, both methods perform better in four datasets each. The results in Table 1 indicate that the joint calibration method has higher accuracy between master-slave LiDAR point clouds than the independent calibration method, while the accuracy of the two methods is comparable between point clouds of different scan lines of the same LiDAR.

[0140] Table 2 shows the Mean1 and Mean2 values ​​before and after calibration using the independent and joint calibration methods. As can be seen from Table 2, the average Mean1 for the eight datasets using the joint calibration method is 1.015 cm, which is better than the 1.126 cm of the independent calibration method. Unlike RMSE2, the average Mean2 of the joint calibration method (1.610 cm) is slightly better than the independent method (1.702 cm). When using Mean1 as the evaluation metric, the joint calibration method performs better in seven of the eight datasets (MLS1, MLS2, MLS4, MLS5, MLS6, MLS7, MLS8), and performs similarly in the remaining one (MLS3). However, when using Mean2 as the evaluation metric, the independent calibration method performs better in six datasets (MLS1, MLS2, MLS5, MLS6, MLS7, MLS8). The results in Table 2 indicate that the joint calibration method can achieve higher accuracy between master-slave LiDAR scanned point clouds.

[0141]

[0142]

[0143] The verification and evaluation of this embodiment show that the present invention can significantly improve the calibration performance of a dual-LiDAR mobile laser scanning system.

[0144] (1) For the mean distance from point to surface, this embodiment achieves a relative accuracy of about 1.1 cm between master and slave lidar point clouds and a relative accuracy of about 1.6 cm between point clouds of different scan lines;

[0145] (2) Regarding the root mean square error (RMSE) of the distance from point to surface, this embodiment can achieve a relative accuracy of about 1.7 cm between the master and slave lidar point clouds and a relative accuracy of about 2.5 cm between the point clouds of different scan lines;

[0146] (3) Compared with the independent calibration method, the joint calibration method has better relative accuracy between master and slave lidar point clouds, while maintaining comparable calibration accuracy between point clouds of different scan lines.

[0147] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for joint calibration of dual lidar systems in a mobile laser scanning system, characterized in that, Includes the following steps: Step S1: In the first stage, the main lidar is independently calibrated. Constraint equations are constructed by matching corresponding points between flight routes. The installation angle parameters of the main lidar are solved by minimizing the distance between corresponding points between different flight routes. The specific steps include: Step S11, Direct geolocation and orientation of the main lidar: Based on the initial placement parameters of the main lidar, the original point cloud data obtained after scanning by the main lidar, the POS data provided by GNSS and IMU, and the direct geolocation and orientation equation of the main lidar, the point cloud of the main lidar in the world mapping coordinate system is obtained. Step S12, main lidar data preprocessing: preprocess the point cloud data in the world mapping coordinate system obtained after scanning by the main lidar; the data preprocessing steps include point cloud resampling and denoising, point cloud segmentation, and obtaining the route number of each point in the main lidar point cloud data based on POS data. Step S13, Matching corresponding points of the main lidar: For the block-shaped point cloud obtained in step S12, the VGICP point cloud registration algorithm is used to match corresponding points of the point clouds of different routes within the same block of the main lidar. Step S14, Main LiDAR installation parameter calibration: Based on the pairs of corresponding points between different routes of the main LiDAR matched in step S13, the installation angle parameters are solved by minimizing the distance between the corresponding points using the LM least squares optimization algorithm. Step S15: Repeat steps S11 to S14 until the average change in the point-to-surface distance between all corresponding points before and after calibration is less than the preset threshold, and obtain the optimal placement angle parameters of the main lidar. Step S2, the second stage, uses the point cloud data of the calibrated main lidar as a benchmark to perform 6-DOF parameter collaborative calibration of the offsets of the three arms and the three mounting angles of the slave lidar, ultimately achieving joint calibration of the two lidars; the specific steps include: Step S21, direct geolocation and orientation of the main lidar and the slave lidar: Based on the original point cloud data obtained after scanning by the main lidar, the original point cloud data obtained after scanning by the slave lidar, the POS data provided by GNSS and IMU, the initial placement parameters of the slave lidar, and the calibration results of the placement parameters of the main lidar, the point cloud of the main lidar and the slave lidar in the world mapping coordinate system is obtained through direct geolocation and orientation. Step S22, Data preprocessing of main lidar and slave lidar: The point cloud data in the world mapping coordinate system obtained after scanning by the main lidar and slave lidar is preprocessed. The data preprocessing steps include point cloud resampling and denoising, point cloud segmentation, and obtaining the route number of each point in the slave lidar point cloud data based on POS data. Step S23, Matching corresponding points between the main lidar and the slave lidar: For the block-based point cloud obtained in step S22, match corresponding points of each route in the slave lidar point cloud with the main lidar point cloud in the same block to obtain corresponding point pairs between the main lidar point cloud and the slave lidar point cloud. Step S24, calibration of lidar placement parameters: by minimizing the distance between corresponding points obtained in step S23, the LM least squares optimization algorithm is used to solve for the lever offset parameters and placement angle parameters of the lidar. Step S25: Repeat steps S21 to S24 until the average change in the point-to-surface distance between all corresponding points before and after calibration is less than a preset threshold, and obtain the optimal arm offset parameters and placement angle parameters from the lidar.

2. The dual-lidar joint calibration method for a mobile laser scanning system according to claim 1, characterized in that, The specific steps of step S11 include: Step S111: Calculate the coordinates of the scanning point in the main lidar coordinate system using the distance and direction angle recorded during lidar scanning. The calculation formula is shown in equation (1) below: （1） in, For scan points Spatial rectangular coordinates in the main lidar coordinate system For scanning distance, It is the direction angle. This is the transpose of the coordinate vector; Step S112, the dual-LiDAR mobile laser scanning system includes a LiDAR coordinate system, an IMU coordinate system, a local horizontal coordinate system, and a world mapping coordinate system, based on the scanning points. The relationship between the main lidar and various coordinate systems is established, and the direct geolocation orientation equation is shown in equation (2) below: （2） in, Main lidar scanning point Coordinates in the world cartographic coordinate system This is the rotation matrix between the local horizontal coordinate system and the world cartographic coordinate system. Latitude, longitude, and geodetic height are provided for GNSS / IMU, where W is the world cartographic coordinate system; This is the rotation matrix between the IMU coordinate system and the local horizontal coordinate system. Roll angle, pitch angle, and yaw angle provided for GNSS / IMU; The rotation matrix from the master lidar coordinate system to the IMU coordinate system. The placement angle between the main lidar and the IMU system. For scan points Spatial rectangular coordinates in the main lidar coordinate system; This is the offset vector of the main lidar in the IMU coordinate system. , , The lever arm offset between the main lidar and the IMU system; This represents the position of the IMU coordinate system origin in the world cartographic coordinate system. Step S113: Solve the installation angle of the main lidar by minimizing the distance between corresponding points on different routes in the main lidar point cloud data. Obtain the initial installation angle of the main lidar from the system design drawings, and calculate the rotation matrix from the main lidar coordinate system to the IMU coordinate system using the initial installation angle. The calculation formula is shown in the following formula (3): （3） Due to the existence of placement angle error, it is necessary to... Perform calibration, and record the calibration amount as follows: , , The error caused by the placement angle is corrected by the calculated corrected rotation matrix; the corrected rotation matrix is ​​shown in equation (4) below: （4） in, This represents the corrected rotation matrix. This represents the rotation matrix correction amount. , , These represent the installation angles of the main lidar. The correction amount, that is, the placement angle parameter that the main lidar needs to be calibrated; Substituting the corrected rotation matrix shown in equation (4) into the direct geolocation orientation equation of the main lidar shown in equation (2), we obtain the corrected direct geolocation orientation equation of the main lidar, as shown in equation (5) below: （5） in, This represents the corrected direct geolocation orientation equation of the main lidar.

3. The method for joint calibration of dual lidar in a mobile laser scanning system according to claim 2, characterized in that, The specific steps of step S13 include: Step S131, using and These represent the point cloud set to be registered and the reference point cloud set, respectively. Given the number of points in the point cloud to be registered and the reference point cloud; perform a nearest neighbor search to make... ,in, Indicates the registration point. Indicates the registration point The nearest reference point , ; This is the transformation matrix between the point cloud to be registered and the reference point cloud; Assuming the surface containing the laser point cloud follows a Gaussian distribution, i.e. , It follows a Gaussian distribution. , Points and points The mean of the Gaussian distribution, , Points and points The variance of the Gaussian distribution is given; then the registration error is defined as: (9) in, Mean with the mean The error between; Based on the properties of the Gaussian distribution, we know that... It follows a Gaussian distribution as follows: (10) in, Point With point The error between; Step S132, calculate the registration point. Rather than in the reference point cloud set The distance in is less than neighborhood points The sum of the distances between: (12) in, Representing neighborhood points, Representing neighborhood points The mean of the Gaussian distribution, Point The mean of the Gaussian distribution Rather than in the reference point cloud set The distance in is less than The mean of the Gaussian distribution of the neighborhood points The sum of distances between them, where r represents the neighborhood radius; Indicates the registration point Rather than in the reference point cloud set The sum of the distances between points in the neighborhood of r is less than the sum of the distances between those points. The distribution is as follows: (13) (14) (15) in, Indicates error The mean of the Gaussian distribution, Indicates error The variance of the Gaussian distribution, Indicates the first The variance of the Gaussian distribution of the neighborhood points, Point The variance of the Gaussian distribution; Step S133: Obtain the optimal transformation matrix by maximizing the logarithm of the likelihood probability. : (16) To calculate equation (16) more efficiently, equation (16) is transformed into: (17) in, For point In reference point cloud set The number of neighboring points in the array; Step S134: Use the optimal transformation matrix for registration parameters. Registration point cloud set The points in the cloud are transformed, and for each transformed point, a reference point cloud is created. Find its nearest neighbor and take the nearest neighbor pair with a distance less than the threshold as the same point.

4. The dual-lidar joint calibration method for a mobile laser scanning system according to claim 3, characterized in that, The specific steps of step S14 include: Step S141, define the same-name point pairs in the point clouds of different routes of the main lidar as ( , ),in and If the index is a point, the placement angle of the main lidar is estimated as follows: (18) in, This represents the parameter vector of the main lidar. , The ) indicates the calibrated installation angle of the main lidar. This indicates the angular rotation of the X-axis of the calibrated master lidar coordinate system relative to the X-axis of the IMU coordinate system. This indicates the angular rotation of the Y-axis of the calibrated master lidar coordinate system relative to the Y-axis of the IMU coordinate system. N represents the angular rotation of the Z-axis of the calibrated master lidar coordinate system relative to the Z-axis of the IMU coordinate system, and N1 is the total number of pairs of points with the same name in the master lidar point cloud data. Let be the residual function of the o-th pair of corresponding points, defined as: (19) in, Point Time The vector, For point The normal vector of the plane; Represents the residual function. This represents the vector dot product, therefore Actual representation point Time The distance from a point on the plane to the plane itself; (20) in, Point Time ; Step S142, using the residual function of corresponding point pairs in the main lidar point cloud data Parameter vector of the main lidar Differentiate to construct the Jacobian matrix: （21） in, Representing the residual function Jacobian matrix, For the derivative sign; Step S143: According to formula (5), the corrected direct geolocation orientation equation of the main lidar can be derived. The derivative is: （22） Due to the rotation matrix correction It is a small-angle rotation matrix, which can be approximated as: （23） in, It is the identity matrix; Therefore, the derivative term can be linearized as follows: （24） in, This represents the X-coordinate of the point cloud in the IMU coordinate system, calculated using the initial placement angle. This represents the Y-coordinate of the point cloud in the IMU coordinate system, calculated using the initial placement angle. This represents the Z-coordinate of the point cloud in the IMU coordinate system, calculated using the initial placement angle; Step S144, final residual function The Jacobian matrix is: （25） in, Point The matrix after linearization of the derivative terms, Point The matrix after linearization of the derivative terms; Step S145: Then, establish the incremental equation for the main lidar and introduce the damping factor. : （26） in, yes The Jacobian matrix of dimension N1, where N1 is the total number of pairs of identical points in the main lidar point cloud data. Representing the Jacobian matrix transpose, Represents a diagonal matrix; It is a parameter increment. ; Represents the residual vector of the incremental equation of the main lidar; Step S146, the parameter update iteration formula is: （27） in, This represents the updated parameter estimate. This represents the parameter estimate at the e-th iteration; The LM least squares optimization algorithm converges to the optimal parameters through iteration. The iteration terminates when the change in residual is lower than a preset threshold or when the maximum number of iterations is reached.

5. The dual-lidar joint calibration method for a mobile laser scanning system according to claim 4, characterized in that, The specific steps of step S15 include: Repeat steps S111 to S146 until the average change in the point-to-surface distance between all corresponding points before and after calibration is less than a preset threshold, thus obtaining the optimal placement angle parameters of the main lidar. The offset of the main lidar arm Obtained through direct measurement or from design drawings.

6. The method for joint calibration of dual lidar in a mobile laser scanning system according to claim 5, characterized in that, The specific steps of step S21 include: From the lidar calibration, the arm offset after the main lidar calibration is calculated. and placement corner The initial lever offset is fixed and cannot be changed; it is obtained from the lidar by direct measurement or from the design drawings. and initial placement angle The initial calibration coordinates of the lidar are obtained by directly geolocating and orienting the raw point cloud data from the lidar. ; The coordinates are obtained from the lidar point cloud after direct geolocation and orientation, expressed as: （28） in, This indicates the coordinates of the scan point from the lidar within the lidar coordinate system. This represents the offset vector from the lidar in the IMU coordinate system; This represents the rotation matrix from the lidar coordinate system to the IMU coordinate system; The small-angle rotation matrix from the lidar coordinate system to the IMU coordinate system is approximately: （29） in, It is the identity matrix. , , These represent the angles from which the lidar is installed. The correction amount.

7. The method for joint calibration of dual lidar systems in a mobile laser scanning system according to claim 6, characterized in that, The specific steps of step S24 include: Step S241, define the corresponding point pairs between the point cloud data of the master lidar and the slave lidar as... The calibration of the lidar installation parameters is estimated as follows: (30) in, This indicates the calibration parameters of the lidar. This indicates the offset of the lever arm from the X-axis of the lidar coordinate system relative to the X-axis of the IMU coordinate system after calibration. This indicates the offset of the lever arm from the Y-axis of the lidar coordinate system relative to the Y-axis of the IMU coordinate system after calibration. This indicates the offset of the lever arm from the Z-axis of the lidar coordinate system relative to the Z-axis of the IMU coordinate system after calibration. This represents the angular rotation from the X-axis of the lidar coordinate system relative to the X-axis of the IMU coordinate system after calibration. This represents the angular rotation from the Y-axis of the lidar coordinate system relative to the Y-axis of the IMU coordinate system after calibration. This represents the angular rotation from the Z-axis of the lidar coordinate system relative to the Z-axis of the IMU coordinate system after calibration. N represents the parameter vector from the lidar, and N2 is the number of corresponding point pairs between the main lidar and the lidar point cloud data. Let be the residual function of the q-th pair of corresponding points, defined as: (31) in, Point Time The vector, For point The normal vector of the plane at that point; Represents the residual function. , These represent the correction amounts for the arm offset and mounting angle of the laser radar, respectively; that is, the arm offset and mounting angle parameters that need to be calibrated for the laser radar. Actual representation point Time The distance from a point on the plane to the plane; (32) in, Point Time ; Step S242: The LM least squares optimization algorithm is used to estimate the optimal placement parameters, based on the residual function of corresponding point pairs between the main and slave lidar point cloud data. For the parameter vector of the lidar Differentiate to construct the Jacobian matrix: (33) in, Representing the residual function Jacobian matrix; Step S243, since the lidar needs to simultaneously optimize 6 degrees of freedom parameters. Therefore, the residual function Partial derivatives need to be calculated for the lever arm offset parameter and the placement angle, respectively: （34） （35） in, These are the coordinates of the lidar point in the IMU coordinate system; Combining the two types of derivatives, namely equations (34) and (35), yields the Jacobian matrix from the lidar: （36） in, Representing the residual function Jacobian matrix; Step S244, introduce the damping factor Construct incremental equations from lidar: （37） in, yes The Jacobian matrix is ​​N2, where N2 is the number of identical point pairs between the main and slave lidar point cloud data. Representing the Jacobian matrix transpose, Represents a diagonal matrix; It is a parameter increment. For 6 degrees of freedom parameters The parameter increment; This represents the residual vector from the incremental equation of the lidar. Step S245, the parameter update iteration formula is: （38） in, This represents the updated parameter estimate. This represents the parameter estimate at the e-th iteration. For 6 degrees of freedom parameters The parameter increment.

8. The method for joint calibration of dual lidar in a mobile laser scanning system according to claim 7, characterized in that, The specific steps of step S25 include: Repeat steps S211 to S244 until the average change in the point-to-surface distance between all corresponding points before and after calibration is less than a preset threshold, thus obtaining the optimal arm offset parameters from the lidar. and placement angle parameters .

Images