Extrinsic LiDAR/Ground Calibration Method Using 3D Geometrical Plane-Based Estimation

This paper details a new extrinsic calibration method for scanning laser rangefinder that is precisely focused on the geometrical ground plane-based estimation. This method is also efficient in the challenging experimental configuration of a high angle of inclination of the LiDAR. In this configuration, the calibration of the LiDAR sensor is a key problem that can be be found in various domains and in particular to guarantee the efficiency of ground surface object detection. The proposed extrinsic calibration method can be summarized by the following procedure steps: fitting ground plane, extrinsic parameters estimation (3D orientation angles and altitude), and extrinsic parameters optimization. Finally, the results are presented in terms of precision and robustness against the variation of LiDAR’s orientation and range accuracy, respectively, showing the stability and the accuracy of the proposed extrinsic calibration method, which was validated through numerical simulation and real data to prove the method performance.


Overview
With the evolution of technology, the 3D intelligent sensors have posed great challenges in signal processing, especially in their outstanding acquisition performance even in rough environment. Moving to road networks maintenance and transportation safety, the responsibility imposes itself in detecting and locating the road distortion (cracking, patching, potholes, rutting, shoving, etc.). The literature review in [1] presents different automated detection experiments and extensive research conducted on pavement adversity in recent years. The work shows the importance and the incredible progress of 3D sensors compared with the other sensors, especially the laser profiler that is characterized by its high precision measurement capability, high spatial resolution and acquisition flexibility.
In a related context, road defects present a big danger on the traffic circulation ending up with possible traffic accidents. Some traffic accidents result from the presence of disabilities or small obstacles on the roads. In fact, this is one of the major problems that people suffer from in their daily lives. The key problem in this research concerns the characterization of the road surface by detection, localization, and tracking the presence of potentially dangerous areas and road defects using a 3D LiDAR sensor. Although 2D LiDAR sensors can provide 3D data, they require the use of an additional instrument in the form of a tilt unit.
Basically, LiDAR sensor operation relies on the calibration process which improves the defect detection and other study procedure. Two different types of calibration exist: intrinsic and extrinsic calibration. The intrinsic calibration considers the modeling between the beam creation and measurement of the environment to estimate the sensor-environment relationship in terms of the sensor internal parameters. On the other hand, the extrinsic calibration considers the determination of the relationship between the sensor frame and the world reference frame by rotation and translation transformation.

Related Works
Numerous authors investigated intrinsic and extrinsic calibration methods in LiDAR sensors. An intrinsic calibration technique is presented in [18]. The calibration process is based on an optimization method, where the calibration pattern is a wide planar wall on a flat surface scanned using Velodyne HDL-64E. In addition, Glennie and Lichti [19] presented a static calibration technique to derive an optimal solution for the laser's intrinsic calibration parameters by a planar feature-based least squares in advantage of minimal constrained network. The study in [20] shows a correlation between the internal operating temperature of the LiDAR and the Laser scanner ranging error (intrinsic parameter). The calibration process considers a planar calibration approach to estimate the internal parameters for Velodyne VLP- 16. On the other hand, an extrinsic calibration technique is presented in [21]. In this technique, a flat plane is used for the calibration and an algorithm based on the inequality of two symmetric rays in azimuth with respect to the origin is proposed. This inequality is due to the shift angle of the center line. Another extrinsic calibration technique is presented in [22], where the authors worked on a 2D laser scanner and on the rotating platform to extract the rotation axis and radius using point-plane constraint. The Levenberg-Marquardt optimization method is applied in the two above extrinsic calibration methods to solve the non-linear least squares function problem.
In [23], a numerical algorithm is presented to compute both of the intrinsic and extrinsic parameters by minimizing the systematic errors due to the geometric calibration factors. Another approach introduced in [24] computes the intrinsic and extrinsic parameters of LiDAR sensor (Velodyne HDL-64E) by unsupervised calibration for each of multi-laser beams. An optimization function seeking to minimize the point-to-plane iterated closest point is then proposed.

Proposed Method
In transport applications, many articles use LiDAR to detect and track objects of interest (vehicle, pedestrians, etc.) from 3D measurements. The LiDAR sensor is also used to detect the road, often in addition to camera sensors. In these applications, the idea is to have a thorough view of the driver's environment over the widest possible horizon. This therefore involves a LiDAR sensor with a low angle of inclination (horizontally oriented sensor).
The perspective of our paper is to propose a calibration method (and road plane estimation) that works under difficult experimental conditions (high angle of inclination). Indeed, we aim at developing a calibration method that allows to determine precisely the road plane in a very close vicinity of the vehicle. The idea in the long term is to detect road defects when driving on the road network. Although developed with this in mind (i.e., with a high degree of accuracy in determining the road plan), our method is general enough to be applicable in any wider operational context.
In the context of this study (road defects detection), the LiDAR sensor is rotated toward the ground to increase the points' density covering the defects by the multi elevation laser. This causes a complicated modification in the ground 3D view scene with respect to the LiDAR frame. Therefore, extrinsic calibration was adopted in order to transform the LiDAR frame into a global reference frame, thus modifying the ground impact points transformation into an understandable view scene.
To attain the above key objective, a first method was applied by Zaiter et al. [25] on simulation data, which was restricted to the Euler angles estimation. In this conference paper, we proposed a first approach to extrinsic calibration of LRF sensors. It has been developed for the Velodyne VLP-16 LiDAR and the theoretical approach have been evaluated on some simulation results.
This paper addresses a new flexible extrinsic calibration method as compared to the previous plane-based methods. The approach developed is generalized to all types of scanning laser rangefinders. It now presents an optimized estimation of all extrinsic calibration parameters (angles and height). This global method can be implemented on different LiDAR sensors (low-cost 3D and full 3D) with various range accuracy. In addition, the proposed technique out performs in high orientation scenarios, which is a very interesting and challenging task that aims to increase the points' density coverage. The proposed calibration method can be summarized by the following two-fold contributions: (1) ground plane model estimation; and (2) rotation transformation matrix estimation from world ground reference to LiDAR sensor frame. The 3D Euler's angles (sensor orientation) and the height (sensor altitude above the ground) are two essential extrinsic parameters required to calibrate the full 3D LiDAR sensors, in order to improve the capability of road defect detection as explained in Section 2.1. In addition, the problem is modeled by 4-DOF (degree of freedom) transformation, namely 3-DOF rotation and 1-DOF height, instead of 6-DOF transformation, namely 3-DOF rotation and 3-DOF translation. This modeling advantage provides the simplicity in the optimization process of the extrinsic parameters.
The structure of this paper is as follows. Section 2 presents the correlation between extrinsic parameters and the geometrical pattern reflection on the ground, the rotated multi-laser beams projection modeling on the ground and the associated measurement errors on the 3D points cloud position. Section 3 presents in detail the different steps of the proposed LiDAR/Ground Calibration Method. Then, the method is evaluated on LiDAR's synthetic and real data in Section 4.

LiDAR/Ground Geometrical Impact Modeling
The synthetic data are generated depending on the features of multi-laser rangefinder or 3D LiDAR sensor, where the environment impact points can be modeled as an intersection between the LiDAR laser beams and the environmental surrounding surfaces. In this work, the LiDAR sensor must be oriented toward the ground to study the road defects. Therefore, the LiDAR laser beams are represented as straight lines and the ground surface as a flat plane in a 3D frame.
Depending on the application situation, two concepts can represent the geometrical reflection model between the LiDAR sensor and the ground surface, as shown in Figure 1: • Practical orientation concept: The LiDAR laser beams (d) are supposed to be rotated and the ground's real plane (P re ) is a fixed horizontal plane, as shown in Figure 1b. • Scientific orientation concept: The LiDAR laser beams (d) are supposed to be fixed and the virtual horizontal ground surface (P H ) must be rotated by the LiDAR's inverse orientation in the practical concept, to get the real oblique ground plane (P re ) in LiDAR frame, as shown in Figure 1c.

Extrinsic Parameters vs. Practical Concept
The four extrinsic parameters (altitude and orientation angles) affect the research goals, where the orientation parameter is the strongest influence factor on the ground points distribution process, as shown in Figure 2a,b. The proposed calibration method must satisfy two contradictory conditions in relation to the final research objectives: • Goal: The plane-based extrinsic calibration needs large sparsity area to improve the plane estimation, which requires high altitude and low orientation angles. • Constraint: The stated finality of road surface object detection needs high density points to improve the capability of defect coverage points, which requires low altitude and high orientation angles.
Therefore, a trade-off is needed to optimize the extrinsic parameters (altitude and orientation angles), providing the suitable coverage points distribution over the ground. Four geometric view patterns are summarized in three cases depending on the variation of pitch angle φ x with respect to the LiDAR's vertical field of view "vFOV", as shown in Figure 2: The study case in this problem focuses on the hyperbolas case, as shown in Figure 2e, in order to increase the points density on the ground.

LiDAR Laser Beams and Oblique Ground Surface Intersection
In the following part, the scientific orientation concept is chosen to model the reflections of the laser beams (d) on the real oblique ground surface (P re ). Therefore, the LiDAR is supposed to be fixed and the parametric equations of the fixed straight lines (d) in LiDAR frame are given by: where α and β describe the azimuth and the elevation angles of each laser beam and t is the parameter of the parametric representation. The virtual horizontal ground plane (P H ) must be rotated by the 3D Euler's angles ψ z , θ y , φ x , so that the equation of the rotated real ground plane (P re ) is expressed as a function of the horizontal ground plane (P H ) with height h and rotational matrix R z,y,x (ψ z , θ y , φ x ). This transformation is expressed as: The parametric and Cartesian equations of the horizontal ground plane (P H ) are expressed as follows: The rotational matrix R z,y,x (ψ z , θ y , φ x ) is expressed as: Therefore, the Cartesian coordinates of the real points cloud c re obtained from the intersection between the fixed straight lines (d) and the rotated real plane (P re ) are expressed as:

Error Modeling in Polar and Cartesian Coordinates
In this section, the systematic and random errors w ρ , w α , w β are taken in consideration as a source of error to represent the modeling of the additive white Gaussian noise for each polar coordinates ρ, α, β of the real points cloud c re in Equation (7), where ρ w , α w , β w are the real measurements of the range ρ, azimuth α and elevation β, respectively, for each reflecting point.
The 3D transformation from polar coordinates ρ w , α w , β w to Cartesian coordinates x w , y w , z w of the ground noisy points cloud c w is given by: Then, the standard deviation of the error can be derived from polar to Cartesian parameters in Equation (9), assuming that σ x w , σ y w , σ z w , σ ρ w , σ α w , σ β w are, respectively, the standard deviations of the added noise on x w , y w , z w , ρ w , α w , β w . The terms of the standard deviations σ ρ w , σ α w , σ β w with a power higher than two can be neglected in this derivation to obtain this approximation: Glennie et al. [20] showed in particular that the error of a scanning LiDAR sensor is mainly manifested over range. In this type of sensor, angles are not directly measured, but the error is mainly related to the reproducibility of the measurement for a given angle. The hypothesis of neglecting the scanning angle error is a very common assumption in the field of LiDAR detection: it is part of the manufacturers' specifications and is commonly used in the literature. This is particularly related to the very small influence of the angle reproducibility errors on the range measurement of the object of interest.
In this study, we then focus on the range error σ ρ w and neglect the azimuth and elevation errors σ α w , σ β w , respectively, in the simulation data as given by the constructor. The transformation in Equation (9) is then simplified as:

LiDAR/Ground Extrinsic Calibration Method
In multi-sensor applications, data acquired from the different sensors must be fused in one common reference frame. In this application, the calibration of LiDAR frame scans is necessary to merge them in one world reference frame, in order to increase the points density coverage on the ground, which facilitates the road defect detection. Therefore, the extrinsic calibration aims to model the relationship between the LiDAR frame and the world reference frame.
We thus propose the LiDAR/Ground Calibration Method (LGCM) presented in Figure 3. The method includes the following procedures: fitting ground plane by Least Squares estimator, rotation about axis by Rodrigues formula, Least Squares Conic Algorithm, and height estimation. The proposed method is supplemented by Levenberg-Marquardt optimization algorithm as opt-LGCM. The main role of LGCM procedure is to estimate the extrinsic parameters: the Euler's rotational angles ψ z , θ y , φ x and the height h. Then, the opt-LGCM is initialized by the estimated extrinsic parameters to optimize them. Finally, the distributed ground noisy points c w along the real plane (P re ) are rotated along the horizontal plane (P H ) by the optimized extrinsic parameters in the frame of fixed LiDAR. The proposed method consists mainly in two steps. The first, totally unsupervised step, consists of estimating a first value of the steering angles. This first estimate is then used as a basis for the optimization step, which seeks the best orientation parameters. The developed method is therefore totally unsupervised and does not require a priori knowledge of the orientation of the sensor by a tilt unit for example.

Fitting Plane
The first step aims to fit an estimated plane (P est ) with the rotated ground noisy points c w . The Least Squares estimator is used to obtain the normal vector of the plane (P est ).
The equation of the estimated plane (P est ) in the LiDAR frame is expressed by: where A, B, and D are the plane parameters and w is an additive white Gaussian noise with standard deviation σ w . Therefore, Equation (11) of the estimated plane (P est ) can be written in linear form as: where where N is the number of reflected points. The solution of Least Squares estimator for this linear model is expressed as:

Rotation about Axis
Rodrigues formula is an efficient rotation transformation that computes the rotation matrix R rod , which rotates a vector into another vector in 3D frame around a fixed axis vector −−→ Axis by rotational angle η [26]. Therefore, after having estimated the parameter vector of the oblique estimated plane (P est ) according to Section 3.1, the next step is to compute the rotational matrix R rod from the normal vector − → n 1 of the oblique estimated plane (P est ) to the normal vector − → n 2 of the horizontal plane (P H ) that is parallel to X L Y L -plane with height −h (cf. Figure 1c). The objective of this step is to use Rodrigues formula in order to estimate the first two Euler's angles pitchφ x , rollθ y , and the first partial yaw angleΨ z1 -due to the incomplete calibration in yaw rotation, which is solved by the next step-from Rodrigues Matrix R rod .
Assuming that − → n 1 , − → n 2 , and −−→ Axis are expressed as: , the Rodrigues rotation formula R rod can be then written as: where Now, by using Equation (15), Then, the Rodrigues matrix R rod provides the computation ofΨ z1 ,θ y , andφ x as expressed in the equations below:Ψ where ij represents the matrix element index of (R rod ) ij . As a graphical result, the ground noisy points c w are rotated by Rodrigues matrix R rod to the distributed points c H1 along the horizontal plane (P H ) by Equation (19), as shown in Figure 4.

Yaw Angle Estimation
After rotating the noisy points c w to the horizontal points c H1 , the second partial yaw angleΨ z2 is estimated by the efficient Algorithm 1 that we propose in Figure 5 to rotate the points c H1 to c H2 about z-axis. This algorithm is called Least Squares Conic Algorithm (LSCA), which takes advantage of the center s characteristic of the geometrical impact patterns (hyperbolas, parabolas, and circles) formed by the points c H1 , as shown in Figure 5. The aim of this part is to compute yaw angleψ z from the partial anglesΨ z1 andΨ z2 , as shown in Equation (20).

Algorithm 1: Least Squares Conic Algorithm.
Input: x,y,z,α,β of the distributed points c H1 Output:Ψ z2 1 Fit the lines (l) and (l ) that pass through the points at each ζ = 10 • consecutive azimuth by Least Squares estimator. The solution of Least Squares estimator for linear model: Compute the coordinates of the intersection points s of each two symmetric lines of (l) and (l ). Assume that: (l) : y =m 1 x +b 1 (l ) : y =m 2 x +b 2 Therefore, the intersection points s of the straight lines (l) and (l ) are computed as follows: Finally, compute the angleΨ z2 formed by the fitting line (v) and y-axis: Therefore, the third Euler's angle of rotation (yaw angle)ψ z is computed as follows: Finally, the points c H1 are rotated to points c H2 by an angle −Ψ z2 around z-axis, as shown in Figure 6.

Height Estimation
At the end of LGCM approach, a suitable way to estimate the height is to compute the altitude mean of the points c H2 , due to the ground geometrical model used in this paper. Therefore, the estimated height is then expressed as:

Extrinsic Parameters Optimization
Levenberg-Marquardt algorithm is an optimization algorithm that combines gradient descent and Gauss-Newton methods [27]. In addition, it is a very efficient technique to find the minima and it performs well on most non-linear functions. Therefore, the role of opt-LGCM is to optimize the extrinsic parameters ψ z , θ y , φ x , h by Levenberg-Marquardt algorithm, which is initialized by the estimated extrinsic parametersψ z ,θ y ,φ x ,ĥ to obtain the optimized extrinsic parametersψ z ,θ y ,φ x ,ĥ opt , in order to minimize the mean square error mse that represents the square difference between the position of noisy points c w and the optimized position of points c opt in Equation (24). The optimized points c opt represent the intersection between all the LiDAR beams (d) and the optimized plane (P opt ), which is the rotation of the horizontal plane (P H ) by the new optimized Euler's anglesψ z ,θ y ,φ x and the optimized heightĥ opt in Equation (25). In other words, the importance of the above procedure is to get the optimized heightĥ opt and the optimized Euler's anglesψ z ,θ y ,φ x that rotate in the inverse ordering orientation the horizontal plane (P H ), to fit the noisy points c w that are distributed along the oblique real plane (P re ) with minimum mse on the position.
The non-linear function mse is expressed by: The optimized points c opt represent the intersection between the straight lines (d) and the rotated optimized plane (P opt ), where the plane (P opt ) is the rotation of the fixed ground horizontal plane (P H ) of heightĥ opt by −ψ z , −θ y , −φ x based on R x,y,z rotation matrix, as shown in Equation (25): where the rotation matrix R x,y,z is the reverse of R z,y,x .
Finally, the Cartesian coordinates of rotated optimized points cloud c opt in the reverse orientation sense are estimated byψ z ,θ y ,φ x ,ĥ opt , as expressed below: (c opt ) :

Experimental Results
The proposed calibration method LGCM was applied on two types of data: simulation data were obtained by the modeling mentioned in Section 2 and real data acquisition from the Velodyne VLP-16 LiDAR. The most important features of Velodyne VLP-16 LiDAR for modeling are shown in Table 1. The extrinsic calibration results are presented in terms of precision and robustness. According to our application, the precision shows the stability of the method with respect to the variation of pitch angle φ x toward the ground, while the robustness shows the method strength with respect to the variation of range accuracy σ ρ of the measurements. Therefore, the evaluation parameters of the results focus on the point cloud features of the real points c re on the real plane (P re ), noisy points c w distributed along the real plane (P re ), estimated points c est on the estimated plane (P est ) obtained by LGCM, and the optimized points c opt on the optimized plane (P opt ) obtained by opt-LGCM, as described below: • The real height h, estimated heightĥ, and the optimized heightĥ opt .

•
The standard deviation σ d w/i of the noisy points c w orthogonal Euclidean distance with respect to the real plane (P re ), the estimated plane (P est ) and the optimized plane (P opt ).
where x w , y w , z w are the Cartesian coordinates of the noisy points c w , A i , B i , C i , D i are the coefficient parameters of the planes, i = {re, est, opt}, and N is the number of impact points.

•
The standard deviation σ ρ re /ρ i of the real points c re range difference with respect to the noisy points c w , the estimated points c est , and the optimized points c opt .
where i = {w, est, opt} and N is the number of impact points.

•
The standard deviation σ ρ w /ρ i of the noisy points c w range difference with respect to the real points c re , the estimated points c est , and the optimized points c opt .
where i = {re, est, opt} and N is the number of impact points.

•
The gain in performance PF i that describes the range accuracy enhancement obtained from the Levenberg-Marquardt optimization algorithm which is defined as: where σ ρ w /ρ re is the LiDAR range accuracy and i = {est, opt}.

Simulation Data Results
Using the simulation data, the setups used to validate the proposed calibration method are separated in two categories:

Standard Deviation σ d w/i in Terms of Precision and Robustness
After analyzing Figure 7a, the increasing of standard deviation σ d w/i along the planes is due to the orientation effect of the LiDAR by the pitch angle φ x on σ d w/i . Thus, as pitch angle φ x tends to 90 • , the standard deviation σ d w/i tends to the LiDAR range accuracy σ ρ w /ρ re . In Figure 7b, the increasing of the standard deviation σ d w/i is due to increasing of LiDAR range accuracy σ ρ w /ρ re . Moreover, Equation (34) describes the relation of σ d w/i with φ x and σ ρ w /ρ re ,which proves the increasing of σ d w/i .
where ϕ k is the elevation angle of each VLP-16 LiDAR laser, and k = {1, 2, ..., 16} represents the laser index. In Figure 7c,d, we can see that the standard deviation σ d w/opt is closer to the the standard deviation σ d w/re than the standard deviation σ d w/est . This shows that the optimized plane (P opt ) is better fit to the real plane (P re ) than the estimated plane (P est ). In terms of precision and robustness, Figure 8 shows the increasing behavior of the range standard deviations σ ρ re /ρ est and σ ρ w /ρ est after the LGCM calibration, due to: • The increase of pitch angle φ x on positive and negative sides decreases the sparsity of impact points on the ground. This leads to decrease the precision of plane fitting estimation, as shown in Figure 8a,c.
• The increase of LiDAR range accuracy σ ρ w /ρ re decreases the precision of plane fitting estimation, as shown in Figure 8b,d.
The standard deviation σ ρ re /ρ opt is lower than the standard deviation σ ρ re /ρ est , as shown in Figure 8a,b, which indicate how the optimized points c opt are closer to the real points c re than the estimated points c est . On the other hand, the standard deviation σ ρ w /ρ opt is closer to LiDAR range accuracy σ ρ w /ρ re than the standard deviation σ ρ w /ρ est , as shown in Figure 8c,d, which indicate the equality of the noisy points c w range distribution along the real plane (P re ) and the optimized plane (P opt ).
The negligible standard deviation σ ρ re /ρ opt in Figure 8a,b and the coincident standard deviations σ ρ w /ρ opt and σ ρ w /ρ re in Figure 8c,d prove the similarity of the real plane (P re ) and the optimized plane (P opt ) compare to the estimated plane (P est ). (d) Variation of σ ρ w /ρ i with respect to σ ρ w /ρ re Figure 8. The variation of σ ρ re /ρ i and σ ρ w /ρ i in terms of precision and robustness.

Height Recovering in Terms of Precision and Robustness
In terms of precision and robustness, Figure 9 highlights the recovering of the height parameter and how the optimized heightĥ opt is closer to the real height h than the estimated heightĥ, which presents the height optimization importance and the strength of the Levenberg-Marquardt optimization algorithm.  Figure 10 shows the gain in performance of the optimized plane points c opt against the estimated plane points c est distributed by the noisy points c w compared to the LiDAR range accuracy σ ρ w /ρ re as expressed in Equation (33), with respect to the variation of pitch angle φ x and LiDAR range accuracy σ ρ w /ρ re . Moreover, the negligibility of the method performance PF opt appears after the optimization, which means that the standard deviation σ ρ w /ρ opt after optimization is closer to the standard deviation σ ρ w /ρ est before optimization with respect to the LiDAR range accuracy σ ρ w /ρ re . In addition, it presents the recovering of noisy points c w range distribution along the real plane (P re ) after the optimization algorithm, taking advantage of maintaining the standard deviation σ ρ w /ρ opt value as negligible. The gain feature PF i proves again the better fit between the optimized plane (P opt ) and the real plane (P re ) rather than the estimated plane (P est ).

Real Data Results
The 3D point cloud acquisitions were obtained using a multi-lasers rangefinder VLP-16 LiDAR mounted on a vehicle. To obtain a telemetric information about the ground surface and to achieve the application goal, the VLP-16 LiDAR was rotated toward the ground direction with a pitch angle φ x 70 • , and it was at a height of h 1.05 m above the ground surface. The real setup is shown in Figure 11. The proposed method was applied to two different acquisitions: • Acquisition 1: The vehicle was at rest on the road. • Acquisition 2: The vehicle was moving at a slow speed on the road.
Standard Deviation σ ρ w /ρ i per LiDAR Frames In the absence of real plane (P re ) when using real data, the results focus on the range distribution of the noisy points c w along the estimated plane (P est ) and the optimized plane (P opt ). It is clear that the standard deviations σ ρ w /ρ opt curve is lower than the σ ρ w /ρ est in the two acquisitions, as shown in Figure 12. The optimization algorithm is thus proved to be more efficient for real data as well in decreasing the range distribution of the noisy points c w along the fitting planes.

Results Discussion
In general, the results prove the efficiency of the optimization algorithm, which is represented by the optimized plane (P opt ), versus the estimated plane (P est ), compared with the real plane (P re ), in terms of precision and robustness. On the other hand, the convergence of the optimization algorithm is granted automatically by the suitable initialization parameters, namely the estimated Euler's angleŝ ψ z ,θ y ,φ x and the estimated heightĥ, which are computed in Stage 1 (LGCM) to obtain the estimated plane (P est ), and then optimized by Levenberg-Marquardt optimization algorithm (opt-LGCM) in Stage 2 to get the optimized Euler's anglesψ z ,θ y ,φ x and the optimized heightĥ opt to obtain the optimized plane (P opt ). Finally, the results show the strength and the method performance in terms of precision and robustness against the variation of pitch angle φ x and LiDAR range accuracy σ ρ w /ρ re , respectively, in order to achieve the application's aim, as shown in Figure 13.

Conclusions
A new extrinsic LiDAR/Ground calibration method for 3D LiDARs is presented in this paper. The solution relies on plane-based modeling of the ground, which allows the estimation of the LiDAR's orientation and altitude using Rodrigues formula, Least Squares Conic Algorithm for yaw angle estimation and height estimation. The proposed method (LGCM) is extended to an optimized derivation (opt-LGCM) using the Levenberg-Marquardt algorithm and is shown to be a suitable solution to LiDAR/Ground calibration problem. It is implemented on synthetic and real LiDAR telemetric data. The results show the performance in terms of precision and robustness against the variation of LiDAR's orientation and range accuracy, respectively, proving the stability and the accuracy of the proposed calibration method.

Conflicts of Interest:
The authors declare no conflict of interest.