1. Introduction

Laser range sensors have been a crucial cornerstone in the development of mapping and navigation over the past two decades [1]. Nowadays, rotating laser scanners and solid-state LiDAR scanners can provide dense and reliable 3D perception at high frequencies [2]. The rapid iteration of related products has also propelled the advancement of the autonomous driving industry [3,4]. Applications of this nature demand high levels of perceptual robustness and accuracy, often achieved through the fusion of multiple sensor systems to capitalize on their complementary strengths.

As one of the most widely used type of LiDAR, mechanical LiDAR offers an extensive sensing range and commonly used in industrial sector. However, it has several drawbacks, including high cost, sparse scanning, a limited vertical field of view, and less effective obstacle avoidance for distant objects. In contrast, the Hybrid Solid-State LiDAR(SSL) is designed to concentrate on a focused scanning area, employing a non-repetitive scanning pattern, leading to the acquisition of denser data [5]. The temporal accumulation of point clouds typically facilitates the detection of over 90% of the FoV(Field of View) [6,7]. Unlike mechanical LiDAR, hybrid solid-state LiDAR has cost-effectiveness and feasibility in large-scale deployment. However, it is limited to perceiving information within a specific narrow angular range [5]. In practical applications, to address the limitations of data sparsity and sensitivity to occlusion, mechanical LiDARs are often employed in conjunction with hybrid solid-state LiDARs. This approach capitalizes on the broad horizontal field of view offered by mechanical LiDARs and the focused characteristics of hybrid solid-state LiDARs within regions of interest. As a result, it significantly enhances both perception range and local measurement density, facilitating the achievement of robust omni-directional perception. As shown in Fig. 1, illustrating the mounting configurations of the two LiDAR systems onto an autonomous driving vehicle, alongside visualizations elucidating the perceptual extents of their respective sensing domains.

 

Fig. 1. An autonomous vehicle equipped with both mechanical and hybrid solid-state LiDARs(a). (b) Top view and (c) side view of the FoV visualization for the two sensors, where yellow and blue rays simulate the laser ranging rays of the mechanical and hybrid solid-state LiDARs, respectively.

Download Full Size | PDF

For any form of multisensory fusion technology, precise calibration stands as a fundamental prerequisite for its efficacy. The amalgamation of these two sensors leverages the complementary attributes of their respective FoV and point cloud densities, rendering advantages in terms of perception. However, these very attributes pose significant challenges for calibration, as they are the primary reasons for intricacies in aligning the point clouds. In the context of point cloud registration, the prevailing majority of existing approaches primarily focus on aligning homogeneous and densely overlapping point clouds, thereby neglecting the alignment intricacies associated with point clouds characterized by low overlap and disparate density.

To address the aforementioned challenges, this study proposes a robust and efficient extrinsic calibration method for multi-LiDAR systems involving heterogeneous scanning modes. The main contributions of this study are outlined below:

2. Related works

In recent years, the calibration of multi-LiDAR systems has captured the interest of numerous researchers, leading to the emergence of diverse methodologies, primarily categorized into motion-based and feature-based approaches. This section will be divided into these two parts for discussion.

2.1 Motion-based calibration

Motion-based methods are referred to as hand-eye calibration in the field of Simultaneous Localization and Mapping (SLAM) [8]. Hand-eye calibration was initially introduced to map measurements from the sensor center to the robot’s workspace frame, enabling precise manipulation of the sensor by the robot [9]. The classic formulation of this method is ${AX = XB}$, where ${A}$ and ${B}$ represent the motions of two sensors, and ${X}$ represents the transformation between the two sensors. Accurately calculating the individual motions ${A}$ and ${B}$ ensures the precision of obtaining the calibration parameters ${X}$. Motion-based calibration methods have successfully overcome the limitations imposed by overlapping FoV and initial calibration parameters. Martinelli et al. [10] introduced an Augmented Kalman Filter (AKF) approach, which concurrently estimates robot configuration and parameters, embedding calibration within a Kalman-based SLAM framework. Building upon their work, Kümmerle et al. [11] incorporates calibration parameters into a graph-based optimization problem. However, in the 6DoF calibration parameter space, the observability of relative motion constraints is closely related to the motion trajectory. Therefore, the approximately planar motion of a vehicle cannot provide sufficient constraints for achieving 6-DoF calibration. Moreover, calibration accuracy is heavily reliant on accurate LiDAR pose estimation. Errors in extrinsic parameters gradually accumulate within the LiDAR map, resulting in insufficient global calibration consistency. To address this, Chang et al. [12] proposed a novel multi-LiDAR calibration method based on global pose graph optimization, where coarse matching is initially performed on independent scan graphs, followed by global pose graph optimization to refine extrinsic parameters. This approach ensures global consistency between the map and extrinsic parameters through a gradual refinement process. However, this method is subject to significant limitations regarding environmental conditions, motion speed and scale, and pose estimation accuracy, which ultimately affect its convenience and practicality.

2.2 Feature-based calibration

Feature-based methods utilize characteristics of point clouds to establish calibration constraints by associating geometric targets within the overlapping field of view (FoV), thereby creating calibration constraints. This approach leverages matched data from the overlapping FoV of multiple LiDARs to formulate registration constraints, and introduces matching algorithms ICP to estimate calibration parameters. The accuracy of data association is a critical step in this process, and therefore, different methods primarily differ in their variations of features and forms of association. Pusztai et al. [13] associated the corners extracted from calibration boxes in the overlapping field of view of multiple LiDARs. In [14], the authors employed segmentation and fitting techniques to estimate the center of a moving calibration sphere, which was used as registration features. Similar methods involving specialized setups can be found in literature, such as the use of retro-reflective tapes [15] and AprilTags [16]. However, these methods are constrained by various factors, such as the need for manual setup of calibration targets, specialized calibration environments or equipment, and the labor-intensive nature of the calibration process, which reduces calibration efficiency. In this context, leveraging features from natural scenes has emerged as an effective strategy to mitigate the reliance on manual targets.

To overcome the limitations posed by calibration targets or specific environments, [17–19] employ the detection of planar features in natural environments to establish accurate associations. [17] presented an approach for plane association within indoor structured scenes, followed by the imposition of coplanarity and perpendicularity constraints on the vertically extracted planes. Lu et al. [20] fitting the ground in a multi-LiDAR setup is conducted to calibrate the roll and pitch angles. This approach involves compressing the 3D point cloud into a 2D grid map and constructing a loss function based on the grid, thereby aligning the calibration yaw angle to obtain the pose.

However, for 16-channel Mechanical LiDARs, due to their sparse nature, the reliability and accuracy of extracting plane features may face challenges when dealing with distant or small planes. Furthermore, for heterogeneous LiDARs, the issue of low overlap will further reduce the available features, which is particularly disadvantageous for feature-based methods. To address the issues of significant data type discrepancies and low point cloud overlap, FRICP [21] introduces the application of the Welsch robust function during the iterative process. Additionally, it employs Anderson acceleration techniques in the iteration process to facilitate easier convergence of the optimization, and achieved state-of-the-art performance.

From the above discussions, it is evident that feature-based methods often impose certain environmental requirements and necessitate stable feature detection and association to accurately correlate data from different LiDARs. Among the existing methods, FRICP’s approach to addressing these issues is most similar to ours. We will compare the effectiveness of our method with it in subsequent experiments. Our approach differs from the FRICP or aforementioned conventional methods in the following two aspects: Firstly, we utilize the matching between point features and plane features, significantly increasing the number of feature matches. Sufficient constraints enable calibration to be completed after a single acquisition, eliminating the need for scene switching and multi-frame fusion or the presence of rich structured features in the scene. Secondly, during the feature matching process, we employ a dynamic residual weighting based on the GMM strategy and utilize an improved trust threshold iteration method, replacing the traditional point-to-plane Gaussian-Newton optimization. This enables a convenient, efficient, and robust calibration process.

3. Methdology

This chapter commences with an elucidation of the fundamental symbols employed in this study and their respective interpretations. Subsequently, the methodological framework is delineated, categorized into two distinct segments: feature extraction and optimization strategies.

The coordinate system for point cloud of the hybrid-solid-state LiDAR is defined as $\left \{\mathcal {C}_\mathcal {H} \right \}$, and it is regarded as the reference point cloud. Correspondingly, $\left \{ \mathcal {C}_\mathcal {M} \right \}$ is defined as a 3D coordinate system of the Mechanical LiDAR, and the its point cloud is applied as the target frame. The key in registering low similarity point clouds is to match and align the homogeneous features in the source and target point clouds. Our method uses at least three linearly independent planes in the environment to calibrate the dual lidar system. Then the geometric structure of these planes is used to obtain the extrinsic parameters. A plane clustered from livox point cloud is denoted as $\mathbf {P}^{\mathcal {H}}_{i}$. Its coresponding coffeicient of the fitting plane is $\pi _{i}^{\mathcal {H}}=\left [n_i, \beta _i\right ] ^{\top }$, where $n_i$ is the unit normal vector of the plane. Distance of point $p$ and this plane can be obtained by:

The set of planes in the point cloud of the whole scene is represented as $\mathbf {\Pi }^\mathcal {H} =\{ \pi _{i}^{\mathcal {H}}|i=0,1,\ldots,N_\mathcal {H}-1 \}$. For point clouds acquired by Mechanical LiDAR in $\left \{ \mathcal {C}_\mathcal {M} \right \}$, use $p^\mathcal {M}$ denote the extracted plane feature points, its coresponding set is represented as $\boldsymbol{\mathcal{P}}^\mathcal {M} = \left \{ \begin {array}{c} p_{j}^{\mathcal {M}}|j=0,1,\ldots,N_\mathcal {M}-1 \end {array} \right \}$. The detailed annotation of parameters are listed in the Tab.1, and the visualization of each notation and its meaning are marked and shown in Fig. 2.

 

Fig. 2. (a) Our mobile platform. (b) Proposed method pipeline. (c) Extracted plane features from raw point clouds. (d) Calibration results. The Hybrid Solid-State LiDAR point cloud with uniform coloring represents $\mathbf {\Pi }^\mathcal {H}$. The point cloud from the mechanical LiDAR is color-coded based on the z-axis for better distinction. Among them, the bold points indicate the points that have been matched with their corresponding planes.

Download Full Size | PDF

Table 1. Annotation Table

View Table | View all tables in this article

Figure 2(b) illustrates pipeline of our calibration process. The platform acquires data from two scanning modes of LiDAR sensors, which undergo data preprocessing. Initially, feature extraction is performed, encompassing both planar features and point features, as elaborated in Section 3.1. The second calibration step involves point-to-plane distance-based feature matching, which is pivotal to the entire system. To this end, we introduce the Gaussian Mixture Model (GMM) outlier clustering algorithm and the Huber function, discussed in Section 3.2.1. Section 3.2.2 delves into the presentation of the Trust Threshold algorithm employed in the optimization iteration process.

3.1 Plane feature extraction

Our method adopted distinct feature extraction approaches for various types of LiDAR data. In the case of hybrid solid-state LiDAR point cloud data, characterized by its high point density and disordered nature, we directly extract fitted planes from the dense point cloud to serve as matching features. We employ a region growing method as [22] do to extract planes from the LiDAR point cloud. Specifically, we first estimate the normals of each point, and then segment the point cloud by clustering points with similar normals. For each cluster, the RANSAC [23] algorithm is applied to detect a plane. Then merge planes with similar parameters and perform plane fitting. More detailed explanation is not provided here as it is not the focus of this paper.

For the repetitive scanning mechanical LiDAR, due to the regular arrangement of its point cloud, The point clouds along the same scanning line possess identical horizontal angles; thus, if the points lie on a planar object, they will form a straight line, as Fig. 3 shows. Curvature between adjacent points is calcurated as Eq. (2):

 

Fig. 3. In the arrangement of scan points by mechanical LiDAR on planar and non-planar surfaces, it is evident that the sequence of scan points on a plane forms a linear arrangement, while scan points on non-planar surfaces exhibit a higher curvature.

Download Full Size | PDF

This spatial curvature describes the planarity of the sampling region surrounding $\mathbf {p}_i$. When the planarity is high, the 10 points within the front and rear windows of $\mathbf {p}_i$ form a straight line, and $c_i$ approaches 0. On the contrary, the magnitude is larger. Therefore, we employ a threshold $c_{th}$ to select planar feature points.

3.2 Nonlinear optimization

To establish a point-to-plane feature matching pair, first, a plane feature point from the source set $\mathbf {p}_i\in \boldsymbol{\mathcal{P}}^\mathcal {M}$ is transformed to reference viewpoint $\left \{ \mathcal {C}_\mathcal {H} \right \}$ using the transformation matrix $\mathbf {T}$ . Then the closest plane $\mathbf {\hat {\pi }}_{i}^{\top }\in \mathbf {\Pi }^\mathcal {H}$ is assigned as the corresponding matching pair for $\mathbf {p}_i$. The minimal distance to the corresponding matching plane is calculated as the residual of $\mathbf {p}_i$:

In contrast to the majority of existing point-to-point [24–27], point-to-plane [21,28], and plane-to-plane methods [29], our approach diverges by focusing on identifying the nearest plane instead of seeking the closest point within the target set. In this way, we obviate the need to construct a kd-tree for pair searching and calculate residuals for all target point clouds. In each iteration, every point $\mathbf {p}\in \boldsymbol{\mathcal{P}}^\mathcal {M}$ is associated with a unique matching plane, while each plane $\mathbf {\pi }\in \mathbf {\Pi }^\mathcal {H}$ may be associated with multiple points from the source set $\boldsymbol{\mathcal{P}}^\mathcal {M}$. The global optimization procedure can be formulated as follows:

This typical non-linear optimization problem can be addressed through the utilization of iterative gradient descent methodology. Nonetheless, a notable challenge arises from the acquisition of dependable initial poses essential for iterative optimization. Inaccurate matches during the iteration process can yield detrimental consequences, resulting in suboptimal result. A simple gradient descent scheme cannot achieve the global optimal results. Our method employs the Majorization-Minimization (MM) algorithm framework. This approach alternates between matching plane queries and updating alignment poses. We re-match features after each optimization iteration and subsequently utilize the updated matching associations to further optimize the minimum distance. The primary alignment process involves the following two steps:

The above two steps are performed iteratively. Moreover, we employed the alignment step as a sub-iteration optimization. Thus, the entire optimization process can be regarded as a Majorization-Minimization (MM) algorithm: a successive upper bound minimization.

Let the objective function in the optimization described by Eq. (4) be $f\left ( \mathbf {T} \right )$. After the correspondence step, we obtained the pose of the current iteration $\mathbf {T}^{\left ( k \right )}$ and the updated surrogate function $g\left ( \left. \mathbf {T} \right |\mathbf {T}^{\left ( k \right )} \right )$, i.e. the objective function in Eq. (6). The surrogate functions satisfy two properties for $\forall \mathbf {T}$: the dominance condition: $g\left ( \left. \mathbf {T} \right |\mathbf {T}^{\left ( k \right )} \right ) \geqslant f\left ( \mathbf {T} \right )$, and the tangent condition: $g\left ( \left. \mathbf {T}^{\left ( k \right )} \right |\mathbf {T}^{\left ( k \right )} \right ) =f\left ( \mathbf {T}^{\left ( k \right )} \right )$. The gradient direction of the agent function at $\mathbf {T}^{\left ( k \right )}$ is the same as that of the objective function. Thus we have:

We repeat the correspondence step to update the surrogate function and repeat the alignment step to turn the non-convex optimization problem into a smooth sub-problem. Such a iteration scheme makes $f\left ( \mathbf {T} \right )$ converge to a global optimal solution. A reasonable construction of $\mathcal {F}$ in Eq. (6) can effectively improve the robustness and efficiency of the optimization, enabling the optimization to be less sensitive to the initial value and accelerating the convergence. The iteration quit when residual difference of two adjacent iterations is less than the threshold $\epsilon$, or the iteration exceeds maximum number of times.

Traditional Point-to-Point based ICP algorithms are usually apply $\mathcal {F}\left ( \cdot \right )$ as a $\ell _2$-norm metric, such that the problem have a closed solution via SVD, while this method is considered noise-sensitive. We use the nonlinear optimization method instead, and apply a robust cost function with adaptive weights, this part is introduced in the following Section 3.2.1.

3.2.1 Robust weighted cost via GMM

Due to the low global similarity of their point clouds in this study, the align optimization may encounter a large number of mismatched pairs and outliers. Thus, the greatest challenge is to improve the align robust by reducing the interference from these pairs. Some ICP variants improve robustness by excluding point pairs with large deviations between positions or normal in the alignment step [30]. Research in [31] indicated that such methods can be difficult to tune and inefficient for iteration. Our method designed a Normal-distribution-based soft weight to penalize point pairs with large differences in position, and the weight function is adaptively adjusted according to the distribution of the residual in each iteration. As the distance of associated point gets farther, the corresponding cost gradually decreases along with the weight, so as to reduce the impact of outliers and non-overlapping regions.

After the corresponding step, each point $\mathbf {p}_i\in \boldsymbol{\mathcal{P}}^\mathcal {M}$ forms a matching pair with its closest plane. For each plane $\mathbf {\pi }_s\in \mathbf {\Pi }^\mathcal {H}$, the inliers, which is actually belong to the plane, are spatially close to each other, thus the histogram of the residuals of these points shows less discrete. We assume that it follows a Gaussian distribution $\left \{ r_{in} \right \} _s\sim \mathcal {N}\left ( \mu _{in},\varSigma _{in} \right )$. The outliers may originally belong to another spatial plane or outside of the common view area, or simply noise points. The degree of dispersion and spatial position of the outliers are also completely different compared to the inliers. Accordingly, we assume that the outlier follows another Gaussian distribution $\left \{ r_{out} \right \} _s\sim \mathcal {N}\left ( \mu _{out},\varSigma _{out} \right )$. Therefore, the align residual $\left \{ r_1,r_2,\ldots,r_{M_s} \right \} _s=\left \{ r_{in} \right \} _s\cup \left \{ r_{out} \right \} _s$ by the points matching with $\mathbf {\pi }_s$ conforms to a two-component Gaussian Mixture distribution Model (GMM), and its probability density function is defined as Eq. (8), expressed as the weighted sum of the Gaussian densities of the two components:

The parameter estimation is conducted using Expectation-Maximization(EM) method in two steps: Expectation step: Initialize the current distribution parameters $\mu _k,\varSigma _k$ by its value in the latest align optimization, and update it by the results of the previous EM iteration. Calculate the probability of generating the current data $\left \{ r_1,r_2,\ldots,r_{M_s} \right \}$ by each distribution component:

Maximization step: Update the parameters $\alpha _k,\mu _k,\varSigma _k$ using the following formula:

Repeat the first two steps of iteration until the value of the likelihood function converges, the residual distribution parameters $\mu _{in}, \mu _{out}, \varSigma _{in}, \varSigma _{out}$ are obtained.

After GMM fitting, the cutoff thresholds for both inliers and outliers are obtained. The weight function for each point residual is defined as follows:

The weight function can be viewed as a zero-mean Gaussian distribution with a peak value of 1 and a variance of $\nu _s$. According to the 3 sigma principle, the inlier threshold range is $r_{in}<\mu _{in}+3\sigma _{in}$. The weight function is truncated when the residual greater than this value, that is $3\nu _s=\mu _{in}+3\sigma _{in}$. Thus we have $\nu _s=\left ( \mu _{in}+3\sigma _{in} \right ) /3$. For ease of understanding, an example of the residual histogram of $\mathbf {\pi }_s$ and the fitted GMM probability density curve as shown in Fig. 4(a), and the weight function is shown as Fig. 4(b). We can see that the inlier residual distribution in the initial optimization are more discrete. After one iteration, the distribution demonstrated a strong trend toward decreasing variance and mean value. The inliers are more clearly distinguished from the outliers through the distribution characteristics. The inliers are more clearly distinguished from the outliers through the distribution characteristics. After that, the new GMM fitting is applied to update residual distribution and the cutoff threshold $\nu _s$, such that the weight function can be adaptively adjusted to make the optimization converge faster and more robust. Accordingly, in our method, Eq. (6) in alignment step can be written as:

 

Fig. 4. The superscript of the formula in the figure represents the number of current iterations, (a) the red and blue curves and histograms represent the residual distribution of the 0th and 1st iteration optimization and the fitted Gaussian Mixture Model curve (b) The curve red and blue is the weight function used in the 0th and 1st iterations

Download Full Size | PDF

3.2.2 Trust region based iterative approximation

In classical nonlinear-optimization-based ICP methods, the Gauss-Newton method is usually applied in the Alignment step. However, due to its linear delta in iteration, the convergence rate may be slow. Moreover, this scheme is relatively sensitive to incorrect matching and outliers, especially when encounters a large number of mismatch pairs induced by the bad initial pose or small overlapping areas, the Gauss-Newton method may not guarantee a global convergence. The trust region based optimization method are applied in our alignment step. It can avoid erroneous iterations caused by singular values, and adaptively adjust the convergence rate and direction to conduct the optimization toward the global minimum. We have proved through experiments 3 that this method can effectively accelerate the convergence and reduce iterations times.

The derivative of the matching residual with respect to the pose is:

Simultaneously, the pose Jacobian matrix under the current iteration can be obtained, and the Levenberg-Marquardt method is used to solve the iteration increment.

Algorithm 1. Extrinsic Optimization Algorithm.

View Table | View all tables in this article

4. Experiment

This chapter employs two approaches, simulation and real-world experiments, to validate the accuracy and effectiveness of our algorithm. In the simulation experiments, we constructed a Gazebo simulation environment, comprising common scene models such as buildings and street vehicles, to realistically simulate the sampling data of the LiDARs in a real-world scenario as closely as possible. In the practical experiments, we conducted data collection on an Automated Guided Vehicle (AGV) equipped with the target types of LiDAR sensors. Both simulation and real-world experiments utilize two types of LiDAR sensors with identical models and parameters: the mechanical LiDAR Velodyne VLP16, featuring 16 scan channels covering a 360$^{\circ }$ horizontal field of view and a 30$^{\circ }$ vertical field of view, with a maximum measurement range of 100 meters; and the hybrid Solid-state LiDAR DJI Livox Avia, which offers a non-repetitive scan with a 70.4$^{\circ }$ horizontal field of view and a 77.2$^{\circ }$ vertical field of view, capable of measuring distances up to 450 meters.

To assess the efficacy of our proposed method, we undertook a comparative analysis with a variant of the Iterative Closest Point (ICP) algorithm, known as FRICP [21]. FRICP stands out as a remarkable variant of the ICP algorithm, showcasing state-of-the-art performance. Opting for this method provides us with a solid benchmark, allowing us to perform comparative tests against the traditional point-to-plane matching ICP, thereby offering a clearer and more comprehensive evaluation of our approach.

4.1 Simulation experiment

4.1.1 Establishment of simulation environment

Due to the challenges in obtaining ground truth for the extrinsic parameters between LiDARs in real-world scenarios, evaluating the calibration accuracy becomes difficult. To address this, our experiments employed Gazebo as the simulator. As one of the premier 3D simulators in the robotics community, Gazebo [32] has the capability to synthesize data from various sensors, including the two different types of LiDAR data used in this study. Leveraging its powerful physics engine, Gazebo allows for the rapid construction of diverse scenarios and the simulation of various sources of errors during data acquisition, thus providing a close approximation to real-world data collection and enabling the retrieval of nearly perfect ground truth information. As illustrated in Fig. 5(a), we have established a scenario for LiDAR calibration comprising common objects such as buildings and street vehicles.

 

Fig. 5. The autonomous vehicle equipped with LiDARs from two perspectives in the street scene, wherein the blue rays represent the laser scanning range of the hybrid solid-state LiDAR.

Download Full Size | PDF

Since our method necessitates the acquisition of plane features, we designed three different scenarios based on the richness of planar structures in the scenes. The visualized point cloud results are illustrated in the Fig. 6.

 

Fig. 6. (a) color-encoded point clouds of scene with abundant plane features, (b) scene with a moderate amount of plane features, and (c) scene with a limited number of plane features.

Download Full Size | PDF

Considering the limitations imposed by vehicle width in vehicular applications, the testing extrinsic parameter positions for both laser radars are restricted to within 0.5 meters. Taking into account that the hybrid-solid-state LiDAR has a Field of View (FoV) of 70.4$^{\circ }$ horizontally and 77.2$^{\circ }$ vertically, the ranges for roll angle and pitch angle are constrained to within 5$^{\circ }$. This configuration ensures that approximately 40% to 50% of the FoV within AVIA is perceived by the mechanical LiDAR VLP16, resulting in a point cloud overlapping ratio. Under this ratio, both the proposed method and the FRICP used for comparison maintain accuracy and robustness [21,33,34]. The six group of ground truth extrinsic parameters, transformed into translation and Euler angles, are presented as Tab.2.

Table 2. Extrinsic parameter configurations

View Table | View all tables in this article

For each scenario, data collection was carried out using six different ground truth extrinsic poses listed in Tab.2. We performed a comparative analysis to compare our proposed method with the conventional ICP approach. To provide a clear demonstration, we conducted testing in moderately challenging scenes under the aforementioned six poses. We assessed the error of $\mathrm {T}(\mathrm {t}^{*},\mathrm {R}^{*})$ for each iteration, including rotational and translational errors, with the calculation as follows:

4.1.2 Evaluation of calibration results on synthetic Data

In the preceding section, we elaborated on the test data and methodology. For the sake of both representativeness and readability, we initially present six test results obtained from the most challenging scenario depicted in Fig. 6(c). The error curve plots depicting the iterative optimization process for our method and the traditional approach are shown in Figure Fig. 7.

 

Fig. 7. The error curves during the registration iterations of our method and traditional point-to-plane matching ICP on simulated data.

Download Full Size | PDF

The results from the graphs demonstrate that our method significantly outperforms the traditional optimization approach in terms of iterative speed and efficiency, primarily attributed to the following two aspects.

Firstly, in terms of iterative robustness, the nature of ICP introduces numerous outlier points in each iteration, which can significantly undermine the effectiveness of the iterative process. At times, this can even lead the iterations to optimize in the opposite direction. As the current point-to-plane ICP optimization methods have not incorporated extensive adaptability improvements with respect to matching errors, the compared method in our study, namely FRICP, applies nonlinear Welsch weights to different-sized residuals. However, this weight adjustment does not adequately accommodate and distinguish between inlier and outlier matching points, resulting in a notable inefficiency during the early iterations, particularly when the initial estimate is poor. Consequently, the rate of error reduction remains suboptimal. Our method benefits from the weight adjustment based on Gaussian Mixture Model (GMM), as discussed in Section 3.2.1. This approach enables dynamic adjustment of residual weights according to their distribution during the iterations. Consequently, the presence of outlier points is effectively attenuated in the matching process, reducing their adverse impact on the iteration process.

Secondly, during the iterative process, the compared method often exhibit distinct phases of low efficiency, characterized by a gradual plateau in the error curve. In contrast, our approach maintains a consistently uniform slope in the error curve throughout the iterative process, a benefit attributed to the optimization of our iterative strategy. This is because the iterative process of ICP involves a cycle between optimization and re-matching. In cases where FRICP and other traditional ICP optimization encounters a low rate of parameter variation, it triggers a re-matching process, discarding the previous erroneous optimization direction and continuing with optimization based on the new match with lower outlier rate. In contrast, fundamentally, our approach employs a trust-region-based method, dynamically adjusting the iterative direction and step size based on the current error result. This significantly reduces occurrences of incorrect iterative directions or convergence stalls.

To validate the aforementioned perspectives and demonstrate the efficacy of the approach based on Gaussian Mixture Model (GMM) weights for outlier elimination, the residuals from the test results were clustered. The truncation threshold, as described in Section 3.2.1, was set to $\left ( \mu _{in}+3\sigma _{in} \right )$. For each iteration, the residuals of each iteration’s outcome were labeled as inliers and outliers, and depicted using a box plot format to clearly illustrate their distribution characteristics, as shown in Fig. 8. During the iterative optimization process, all matched points adhere to Gaussian weights. Inliers, which are points with residuals smaller than the truncation distance $\nu _s=\left ( \mu _{in}+3\sigma _{in} \right ) /3$, carry a higher optimization weight. Their residual values converge to smaller values after only a few iterations and are more tightly clustered. On the other hand, outliers stabilize over multiple iterations, yet due to the influence of initial matching, their distribution appears scattered. Consequently, their mean and distribution variance both increase. The distribution in Fig. 8 aligns with our algorithm’s expectations.

 

Fig. 8. Inlier and outlier residual distribution during the registration iteration of our method.

Download Full Size | PDF

We present the final optimized error results of all data in the form of a scatter plot, as depicted in Fig. 9, which visually demonstrates the final error outcomes of different methods in an intuitive manner. Our method exhibits a mean translation error of 9.3003e-05m, outperforming FRICP’s 1.7939e-04m. In terms of rotational error, our method achieves a value of 1.6e-3 degree, surpassing FRICP’s 4e-3 degree. These improvements indicate a minimum enhancement of 48% in accuracy compared to FRICP.

 

Fig. 9. The scatter plot of registration errors between our method and the traditional point-to-plane matching ICP on the simulated dataset, with translation error on the horizontal axis and rotation error on the vertical axis. Data points that are closer to the origin indicate better results.

Download Full Size | PDF

4.2 Real-world experiment

In the real-world experiments, two LiDAR sensors of the same as those used in the simulation experiments were employed. Additionally, three different scenes with varying levels of environmental plane complexity were tested, as illustrated in Fig. 10. To enhance point cloud density, data from the hybrid Solid-state LiDAR, featuring non-repetitive scans, were time-integrated over 1 second in static conditions. Meanwhile, the repetitive scanning LiDAR acquired one frame of point cloud data for matching with the integrated data from the VLP-16 point cloud.

 

Fig. 10. Calibration scenes.

Download Full Size | PDF

Due to the unavailability of accurate ground truth extrinsics in real-world settings, an assessment of data stability was conducted through multiple measurement comparisons. Similarly, a comparative analysis between the proposed algorithm and FRICP was performed. The distribution of converged extrinsic parameters, including translation (x,y,z) and Euler angles (roll, pitch, yaw), for all scene configurations is depicted in the Fig. 11.

 

Fig. 11. Distribution of converged extrinsic values for 3 scene settings

Download Full Size | PDF

The proposed method demonstrates calibration outcomes with mean translational standard deviations of 0.0082, 0.0133, and 0.0140 meters for low, medium, and high difficulty calibrations, respectively. The standard deviations of Euler angles are 0.31, 0.37, and 0.38 degree, respectively. The heat map in Fig. 12 effectively illustrates the calibration result standard deviations across various scenarios. The results reveals that, in each of the aforementioned scenarios, our method exhibits noticeably lower variance and fewer outliers compared to the results obtained with FRICP. Particularly, when the scenes encompass non-planar features like trees and glass, our method maintains a relatively stable outcome distribution, demonstrating greater resilience against environmental noise and interference.

 

Fig. 12. Distribution of converged extrinsic values for 3 scene settings

Download Full Size | PDF

In real-world street scenarios, the alignment of point clouds before and after algorithm calibration is demonstrated in Fig. 13. The larger-sized point cloud, which appears more sparse, represents data acquired by the mechanical LiDAR, color-coded by laser return intensity. The smaller-sized point cloud, denser in nature, corresponds to data collected using the hybrid Solid-state LiDAR, and different colors denote points extracted from distinct planes. From the presentation in Fig. 13(a) and (b), it is evident that even in noisy and limited co-visibility scenarios, the calibration alignment exhibits excellent performance. Further details in sub figures (c), (d) and (e) specifically illustrate the alignment efficacy of the frontal, lateral, and top-down views of building exteriors, respectively.

 

Fig. 13. The raw (a) and aligned-fused (b) point cloud data captured from the streets, along with details of partially aligned point clouds(c),(d),(e).

Download Full Size | PDF

5. Conclusion and discussion

This study presents an automatic extrinsic calibration method for hybrid solid-state LiDAR and mechanical LiDAR systems. The proposed method eliminates the need for auxiliary calibration tools by utilizing structural information within the scene as the calibration medium, thereby achieving high-precision extrinsic calibration. The proposed method is applicable for aligning point cloud data from non-repetitive scanning LiDARs with repetitive scanning LiDARs. It effectively addresses challenges arising from differences in point cloud features and field of view (FoV) between these two types of LiDARs. This method proves to be particularly effective in applications requiring advanced precision calibration, such as autonomous driving perception and large-scale high-precision map reconstruction. Given the LiDAR ranging accuracy of only 2 to 3 cm, the proposed approach achieves millimeter-level positional accuracy and sub-degree rotational accuracy. The primary limitation of this method lies in its reliance on environmental structural features. Additionally, due to the necessity of plane extraction, it imposes point cloud density requirements on hybrid solid-state LiDAR data, often necessitating accumulation of point clouds over a certain period of time. These limitations render our approach suitable solely for calibration and offline mapping purposes, rather than real-time point cloud registration and SLAM applications. In the future, we intend to enhance this method by incorporating a Gaussian Mixture Model (GMM)-based outlier removal technique for point-to-point matching, thereby reducing the reliance on environmental features. Furthermore, the application of multi-frame joint optimization could also elevate the precision of our algorithm.