ANALYSIS OF SYSTEMATIC ERRORS OF MOBILE LiDAR SYSTEMS: A SIMULATION APPROACH

The systematic error analysis of the mobile LiDAR system (MLS) is always a challenging task in real-world situations. This challenge is mainly due to the mixture of systematic errors with non-systematic errors. To tackle this issue, in this paper, we introduce a conceptual model of an MLS simulator. The main advantage of the simulation-based approach is the full control over the erroneous systematic and non-systematic parameters that affect an MLS’s output. In the proposed simulation approach, we only take into account systematic errors that affect the simulated georeferenced point cloud. These systematic errors are as follows, POS-LiDAR boresight angles, POSLiDAR leverarms, range offset, and scan angle offset. To simplify our analysis, we concentrate only on modeling the effects of systematic errors on planar targets and we focus solely on the terrestrial platform. Based on an independent analysis performed on each of the eight systematic errors of an MLS, to obtain strong visibility over systematic errors of an MLS, we suggest two planar targets of 1m x 1m dimensions with vertical and inclined orientations and a five-line pattern for MLS, two parallel and three side-looking passages. The proposed configuration generates an ideal input point cloud for the detection of systematic errors (except for the ZLeverarm error) and ultimately it will lead to the proper input data for calibration of a terrestrial MLS. To validate our methodology, with an in-house assembled terrestrial MLS, we scanned a set of planar targets with three different orientations (vertical, inclined, and horizontal). This real-data validation test illustrated that with only two out of three planar targets (vertical and inclined) and with five out of six passages (two parallel to the planar targets and three side-looking passages), we will obtain expected visibility over the systematic errors of a terrestrial MLS, which approves the results with the simulation data.


INTRODUCTION
In the last four decades, mobile LiDAR systems (MLS) have evolved from a cutting-edge, expensive, and unreachable geomatics technology into a more user-friendly and accessible surveying technique for the acquisition of georeferenced point clouds. We categorize them based on their platforms as terrestrial, marine, and aerial MLS. MLS has enabled geomatics professionals to generate millions of georeferenced points rapidly and at a lower cost than other surveying techniques. The product of these systems can be used in various fields and applications, such as 3D city modeling, autonomous vehicle, and virtual reality (Vosselman and Maas, 2010;Shan and Toth, 2009).
An MLS consists of two main components: a position and orientation system (POS) and a LiDAR scanner (Ackermann, 1999;Wehr and Lohr, 1999;Shan and Toth, 2009;Vosselman and Maas, 2010). The POS combines a GNSS antenna with an inertial measurement unit (IMU). However, the interconnection between these components can be affected by systematic and non-systematic errors degrading the quality of the final georeferenced points.
The main systematic errors are POS-LiDAR boresight angles, POS-LiDAR leverarms, range offset, and scan angle offset. These errors have a systematic, repetitive, and constant effect on the data. The non-systematic errors consist of two categories; blunders or gross errors, such as GNSS multipath error, which affects the platform trajectory, and random errors such as noise originated from the sensors (LiDAR or POS) that affects the point cloud. It is relatively difficult to separate these two kinds * Corresponding author of errors to study them independently without the influence of one on the other. For example, to generate an accurate and precise point cloud, we have to acquire trajectory data to be able to georeference the points. Thus, the georeferenced points will inevitably influence the trajectory errors and as a result, the impact of other systematic errors will not be really clear on the georeferenced points, and therefore, the analysis of such errors will be complex.
In recent years, in the geomatics world, there are two major tendencies in MLS simulation-based approach. The first type is a radiometric simulation approach based on radiometric and spectral interaction of LiDAR photons with environmental complex scenes like forest canopy (Brown et al., 2005;Cifuentes et al. 2018;Kukko and Hyyppä, 2007). The second type is a geometric simulation approach in which a direct georeferencing mathematical model constitutes the base of the simulator (Friess, 2006;Lohani et al., 2006;Kim et al., 2009;Heinz et al., 2015). Lohani (2006) introduced simulation of aerial altimetry MLS to generate specific terrain surface with objects on top of that surface. In Friess (2006), simulated data are used as input data to test the proposed boresight angles estimation algorithms without any discussion on the details of the simulated data production (Skaloud and Lichti, 2006). Kim (2009) proposed a geometric model of the LiDAR sensor to simulate aerial mobile LiDAR data of a ray-tracing algorithm and to facilitate the assessment of data quality and the development of the data processing algorithms. Heinz (2015) generates simulated planar targets with a robotic simulation toolbox called (V-REP). Most of the geometric-based simulation approaches produce simulated data from airborne MLS with a nadir look on scenes (Friess, 2006;Kim et al., 2009). Moreover, as we can see, the most popular geometric feature in the geometric-based simulation approach is a planar target. Indeed, a planar target shows a simple mathematical model and can be used to demonstrate the effect of MLS systematic errors very well (Filin, 2003). In various MLS systematic error estimation approaches, for geometric features of their algorithms, they consider in-situ planar targets like building roofs and facades (Glennie et al., 2016;Skaloud and Lichti, 2006;Zeng et al., 2018;Heinz et al., 2015).
In this paper, the objective is to analyze the effects of the systematic errors of an MLS on simulated planar targets. For this purpose, we introduce in detail, the concept of an MLS simulator that generates simulated point clouds on a pre-defined hypothetical planar target. The simulator is designed in such a way that we can choose to have the impact of each of the systematic errors of an MLS of the simulated point cloud or not. This ability allows us to study the systematic errors separately and independently without the influence of other systematic and non-systematic errors. After the generation of various simulated point clouds, we present the systematic-error visibility criteria that enable us to find out the best configuration in terms of planar target orientation and also MLS passages for each systematic error separately.
The subsequent sections of the paper are organized as follows. In section 2, the conceptual model of the proposed simulator and the visibility criteria index are presented. In section 3, the results of the analysis of the systematic errors for the simulated data are discussed. In section 4, the validation procedure with the realdata point clouds is presented. Finally, conclusions and future works are recommended in section 5.

METHODOLOGY
This section introduces the conceptual model of the simulator that enables us to study the effect of systematic errors of an MLS in a simulated environment. Figure 1 illustrates this conceptual model and different parts of the methodology that generates an MLS simulated point cloud. As we can see in Figure 1, two kinds of input data enter the direct georeferencing module function: the positioning and orientation system (POS) as observations and LiDAR range with systematic error effects as the result of LiDAR range estimation function. In section 2.1, we describe the mathematical model of direct georeferencing and its components.

Direct Georeferencing Mathematical Model
Equation (1) presents the mathematical model of the direct georeferencing of an MLS (Ackermann, 1999;Wehr and Lohr, 1999;Shan and Toth, 2009;Vosselman and Maas, 2010). The output of the Equation (1) model, ( ) ′ , is the 3D position vector of the georeferenced point in a relative coordinate system such as the Local Geodetic Frame (LGF).
In Equation (1), we have two transformations between three frames, from the LiDAR frame to the POS frame and from the POS frame to the LGF frame. Here also, for simplifying our simulation design, we consider that both the GNSS antenna phase center and the IMU gravity center are located at the POS reference center. The reason why we choose to georeference the data in a relative coordinate system such as LGF instead of a projection system like UTM is that the survey and the area that we consider for the analysis of systematic errors such as boresight angles normally covers a small area. Thus, for the sake of simplifying calculations, we avoid the use of a projection coordinate system and remain in a relative local coordinate frame.
The following sections explain each part of the direct georeferencing mathematical model by taking into account their respective contribution to the proposed simulation design and concept.

Trajectory Position vector ( )
In the simulator design, the positions of the trajectory, ( ) ′ which is a 3D position vector for each point of the trajectory in the Local Geodetic Frame (LGF) system are represented by a sequence of points. Figure 2 shows a sample of trajectory points (in black) and simulated LiDAR points (in red) on a hypothetical planar target. We consider that the trajectory follows a straight line at a fixed height from the ground.

Attitude Angles Transformation Matrix ( )
The transformation matrix consists of three attitudes angles roll (r), pitch (p), and heading(h) that form the rotation matrix between the gravity center in the POS frame and the reference point of the LGF system. In the proposed simulator design, we consider that roll and pitch angles are zero, and solely, the heading angle determines the orientation of the platform. For example, for having round-trip trajectories parallel to the planar target, we consider that in the first passage, the heading angle is 0°, and in the second passage, the heading angle is 180°. In our simulation, we can change the heading angle to produce nonparallel points w.r.t the planar target, for example, 45°, 225°, 315°, and 135° heading. Figure 3 illustrates the simulated trajectory of the six passages w.r.t the planar target. In this paper, we consider only these trajectories to simplify the analysis of systematic errors. Our objective is to consider a simple platform (vehicle) motion scenario and analyze the effect of the systematic error on the planar targets based on the proposed configurations.

Target Point Position Vector ( )
The vector describes the position of the target point w.r.t the optical center of the LiDAR scanner. This vector consists of two parameters, the (range) distance from the optical center of the LiDAR scanner to the target point and (scan angle), which can be between 0° and 360°. The direct georeferencing mathematical model uses the polar coordinates representation for the vector. In the proposed simulation design, we directly insert systematic errors such as the boresight angle, leverarm errors, range offset, and scan angle offset, inside the range estimation function. In section 2.2, we will describe in detail how we insert systematic errors into the range estimation function. In Figure 4, a schematic LiDAR scanner in the simulation design is illustrated with its corresponding spherical representation. Figure 4 -Local LiDAR frame with spherical parameters representation

POS-LiDAR Leverarm Vector ( )
In the proposed simulator, to simplify the design, we consider the POS-LiDAR leverarm vector, ( ) ′ , as a translation vector between the LiDAR frame and the POS frame. Thus, the value for the leverarm vector is zero, and we only consider the leverarm components to have errors in X, Y, and Z directions.

Boresight Angles Transformation Matrix ( )
The transformation matrix consists of three boresight angles (dr, dp, dh) between LiDAR and POS frames. These imposed boresight angles affect the calculated range of LiDAR scanner w.r.t a specific target. Due to the effect of these parameters on the range of georeferenced data, we expect that generated point clouds with imposed POS-LiDAR boresight angles will be placed at a wrong location or may even be placed outside the hypothetical planar surface.

LiDAR Range Estimation Function
Among the presented direct georeferencing parameters of Equation (1), as described in section 2.1.3, the range ( ) which is one of the parameters of vector has a significant role in modeling POS-LiDAR systematic errors of an MLS. By taking into account the misaligned MLS due to systematic errors such as boresight angles and leverarms, we calculate the range of the LiDAR point (Morin and El-Sheimy, 2002). To generate a simulated point cloud with inserted systematic errors to the system, we have to be able to implicitly introduce systematic errors such as POS-LiDAR boresight angles, POS-LiDAR leverarms, LiDAR range offset, and LiDAR scan angle offset to the simulator algorithm and produce data affected by these misalignments created by the MLS systematic errors.
In the first step, as shown in Equation (2), we calculate the georeferenced point without considering any of the systematic errors, In a second step, we introduce a function that calculates explicitly the range ( ) of the emitted point from an MLS with imposed systematic errors. The concept of the MLS systematic errors is exactly like in a real-world situation where we do not have any a priori knowledge about these systematic errors. Hence, we can generate misaligned points by adding systematic errors to boresight angles, leverarm vectors, range, and scan angle. After georeferencing the point cloud, we can observe a data misalignment coming from these POS-LiDAR systematic errors. Figure 5 illustrates the procedure to calculate the LiDAR range in the simulation model. Figure 5 -Flowchart of the range estimation function As mentioned before, the geometric feature that we use in our analysis is a planar target. Thus, we study how systematic errors of an MLS impact points constrained to lie over a planar target. Equation (3) represents the mathematical model for the planar target.
In Equation (3), the ( , , ) are the coordinates of the georeferenced point obtained with Equation (1). The plane parameters ( , , , ) define the hypothetical planar target where parameters , , and are the normalized (unit length) normal vector of the planar target and the parameter is the negative orthogonal distance of the planar target w.r.t the origin of the LGF coordinate system (Skaloud and Lichti, 2006).
Thus, each georeferenced point generated by the simulator is constrained to lie on the hypothetical planar target. Based on this assumption, we calculate the range of each point defined as the distance between the mobile LiDAR system and the planar target. This range variable is part of the vector of the direct georeferencing model in Equation (1).
To define the explicit function for range based on the georeferencing mathematical model, we use the symbolic language. In this paper, to generate the explicit function of the range, we have used the MATLAB symbolic toolbox.
Thus, the explicit range function will be a function of all input parameters described in Figure 5 (blue cases), as follows:

Point to Plane Belonging Function
One of the crucial parts of the simulator design is to verify if a georeferenced point belongs to a planar surface or not. Using a mathematical model that defines an infinite plane, we first have to define four corner points of our planar target w.r.t. the system frame. In our proposed method as illustrated in Figure 6, if a point is inside a rectangle, the summation of the areas of the four triangles that the point makes with the four vertices of the rectangle must be equal to the area of the rectangle (Toma, 2012). Thus, the point "M" belongs to the rectangle " 1 2 3 4" because the summation of the areas of the following four triangles: ∆ 1 2, ∆ 2 3, ∆ 3 4 and ∆ 4 1 is equal to the area of the rectangle " 1 2 3 4". If the point is outside the rectangle, the above summation will produce a larger value than the area of the rectangle, and the algorithm will reject it.

Systematic-Error Visibility Criteria Function
Until now, we explained how to generate a point cloud, which lies on a planar target using an MLS with inserted systematic errors. With this capability, we can produce various simulated point clouds with their specific configurations. We define the term "Configuration" as the setup that leads to the final point cloud which is as follows, 1. The orientation of the planar target 2. The passage of the MLS 3. The dimension of the planar target 4. The height difference between MLS and planar target In this research, we consider parameters (1) and (2) of the configuration as variables and parameters (3) and (4) as constants.
Based on various configurations, we will have simulated point clouds on the planar target. To quantify the impact of systematic errors on the plane's point cloud, we introduce two types of "systematic-error visibility criteria". Type I, which makes a raw selection from the configurations, is a quantitative value based on two parameters, planarity and centroidal distance, and Type II, which makes a fine selection from the retained configurations, consists of two parameters; angular deviation and signed orthogonal distance.

Type I -Planarity and Centroidal Distance
As mentioned beforehand, Type (I) systematic-error visibility criterion consists of two parameters, planarity and centroidal distance which makes a raw selection from all the available configurations. We consider the eigenvalues of the covariance matrix of the simulated point cloud to calculate the planarity parameters (West et al., 2004;Pauly et al., 2003;Blomley et al., 2016;Gross et al., 2007). On the other hand, the centroidal distance calculates the distance between the centroid of the hypothetical planar target and the centroid of the simulated planar target and is considered as a dispersion indicator. The combination of these two parameters into one quantitative value gives us the Type (I) systematic-error visibility criteria, as we can see in Equation (4) In Equation 4, the Centroidal_Dist is normalized between [0, 1]. The more the point cloud geometric form is near a plane shape (planarity) and at the same time it is near the hypothetical plane ( near zero), the more the value of Type (I) criteria is near 1. In our raw analysis, these selected configurations (Type (I) criteria near 1) will be kept for a finer selection. In Figure 7, two samples for strong (near 1) and weak (near 0) Type (I) systematic-error visibility criterion are introduced.
(a) (b) Figure 7 -(a) Strong configuration Type (I) systematic-error visibility criterion (b) Weak configuration Type (I) systematic-error visibility criterion As we can see schematically in Figure 7(a), the Type (I) systematic-error visibility criterion is strong. Thus, we keep this configuration for further analysis with Type (II) criteria. On the other hand, in Figure 7(b), the Type (I) systematic-error visibility criterion is weak and we eliminate this configuration.

Type II -Angular Deviation and Signed Orthogonal Distance
In this section, we consider Type (II) systematic-error visibility criteria, which make a finer selection from the retained configurations. The normal vector between the planar surfaces is a very strong and simple indicator that shows the angular deviation between the hypothetical plane and the plane generated by the simulator. The signed orthogonal distance also shows the shortest distance between the centroid of the simulated planar surface point cloud and the hypothetical planar target. As shown in Figure 8, a systematic error of the roll boresight angle in parallel round-trip passages w.r.t the planar target has an impact on the generated point clouds when compared to the hypothetical planar surface in terms of (a) angular deviation and (b) orthogonal distance drift. For Type (II) criteria, we have to make a comparison between various passages and then choose the configuration which has a bigger angular deviation difference and also higher orthogonal distance.
Thus, Type (I) and Type (II) systematic-error visibility criteria make a raw and a fine selection respectively from configurations that ultimately generate the simulated point clouds.
Configurations that there Type (I) and Type (II) criteria value attain the expected threshold will be considered as an ideal input for the detection of systematic errors and consequently chosen for calibration of an MLS. In the next section, these two types of systematic-error visibility criteria will be used as an indicator to select the best planar target configurations.
(a) (b) Figure 8 -Type (II) systematic-error visibility criteria for roll systematic error analysis (a) angular deviation (b) signed orthogonal distance

RESULTS AND DISCUSSIONS
In the analysis of systematic errors of an MLS, the objective is to find out the best configuration (which in this paper is considered the combination of Plane orientation and MLS passage) to obtain Type (I) and (II) systematic-error visibility criteria for each of the eight MLS systematic errors in the expected threshold. In the terrestrial MLS simulator, we consider that the height of the system is fixed and is 'h = 2.5 m' from the ground and the planar target is on the ground at a lower height (less than h). We consider 19 orientations of the planar target w.r.t the terrestrial MLS simulator which varies between -90 o and 90 o with 10 o intervals as illustrated in Figure 9.

Systematic-Error Visibility Criteria Analysis
In this section, we analyze each of the eight systematic errors. With the Type (I) criteria analysis, we keep only the configurations within the expected threshold of [0.8, 1]. This expected threshold includes the configurations that produce point clouds with acceptable planar shape (planarity value near 1) and with proximity to the hypothetical plane (the centroidal distance around 0). Once we have selected the configurations with Type (I) criterion, we then use the Type (II) criteria, which will make a finer selection from retained configurations and will result in a single plane's orientation with a combination of all the necessary MLS passages. The result of Type (II) criteria of each terrestrial MLS systematic error represents the best possible configuration to generate point cloud with the best visibility for that specific error. In the following, a complete analysis is performed on the pitch systematic error. Based on the same procedure, we analyze the other seven MLS systematic errors and we summarize the best configurations for each of them.
To demonstrate the procedure, we take into account the impact of the pitch boresight angle of 2° inserted to the MLS on the hypothetical planar target. Figure 10 presents the values for the Type (I) criterion and Figure 11 the values for the Type (II) criteria for all the six MLS passages and the 19 plane orientations. The blue box in these two figures is considered the plane's orientation that the Type (I) and Type (II) criteria are within the expected threshold.
Passage 0° 180° 45° 225° 315° 135° Color Based on the Type (I) criterion, Based on the possible configurations presented in Table 1, we perform the fine selection with the Type (II) criteria. As shown in Figure 11, the best visibility configuration for the pitch boresight angle systematic error is one of the eight passage possibilities introduced in Table 1, like (180° -Orange color, 135° -Green color, 315° -Light blue color) passages and 0° orientation (vertical plane). This configuration produces the maximum value of both the signed orthogonal distance and the angular deviation when using a vertical plane (0° orientation) which accordingly results in better visibility of pitch boresight angle systematic error. While analysing all MLS systematic errors, we recognized that there is a correlation between the scan angle offset error and the roll boresight angle systematic error. This is due to the same rotational effect of two systematic errors on the planar target. Also, we did not suggest any possible configuration for the Z-Leverarm systematic error. Due to the nature of the Z-Leverarm error, which is a translation and shift in the Z direction in a similar fashion for all the points of all passages. Thus, it is impossible to observe the Z-Leverarm error with our actual planar target configurations (Leslar et al. 2014;Skaloud and Lichti, 2006). Therefore, we suggest the measurement of the Z-Leverarm with precise and accurate technique, like the total station method, then enter the measured value in the georeferencing model like a constant.
(a) Based on our analysis results, the majority of systematic errors can be assessed with a planar target of 0° orientations, which represents a vertical plane and the remaining systematic errors with an inclined planar target between -60° and -40°. For the line pattern passages, we have to choose the combination of necessary passages that covers all the systematic errors.

EXPERIMENTS AND VALIDATION
To validate the proposed methodology, we implemented a realworld test with a mobile LiDAR system composed of a LiDAR scanner (Z+F 9012 profiler), an IMU (iXBlue ATLANS-C) and GNSS antenna (Septentrio) and a structure made of three planar targets as seen in Figure 12. The validation test was carried out on a parking lot where we installed the three planar (horizontal, vertical, and inclined) target assembly and collected data by following our proposed six-line pattern. The dimension of each planar target is 1m x 1m.
To be able to compare the results of this real-data validation test with those of the simulated data from section 3, we consider that the vertical planar target has an orientation between [-10°, 10°], the inclined planar target has an orientation between [-60°, -40°], and orientation between [70°, 90°] is considered a horizontal target. Also, we tried to follow the same line pattern presented in Figure 3 for the simulation. Based on the Type (I) criterion analysis, five out of six passages were within the accepted threshold of 0.8. Only passage P4, one of the side-looking passages, was removed because its Type (I) value was below the 0.8 thresholds. On the other hand, as the horizontal plane does not satisfy the Type (I) criterion for most of the passages, it was not used. Thus, we finally consider only the vertical and inclined planar targets and the five passages (P1, P2, P3, P5, and P6) in our experiment.
In the next step, we analyze the Type (II) criteria based on these selected configurations. For the orthogonal distance criteria with the vertical plane, the difference between the maximum value, obtained in passage P6 and the minimum value, obtained in passage P3 is 9.4 mm. With the inclined plane, the difference between the maximum value, obtained in passage P3 and the minimum value obtained in passage P2 is 49.6 mm. Thus, the inclined plane generates point clouds with the highest orthogonal distances, which demonstrates strong visibility for some of the systematic errors of the selected MLS. Based on the angular distance criteria, for the vertical plane, the difference between the maximum (P3) and minimum (P2) is 1.0892° and for the inclined plane, the difference between the maximum (P5) and the minimum (P3) is 0.3863°. Thus, the vertical plane generates point clouds with the highest angular deviations, which demonstrates strong visibility for some of the systematic errors of the selected MLS. We conclude that both planar targets (vertical and inclined) are essential to satisfy the Type (II) criteria and as a result to obtain strong visibility of all systematic errors of a terrestrial MLS.
Even though we tried as much as we could to reproduce in the real-data validation test the same configurations recommended in the simulation analysis, there are still significant differences between the real and simulation configurations which leads to different values for Type (II) criteria, especially for the orthogonal distance. However, these configurations produce point cloud data with the necessary deviation and drift to generate strong visibility of all systematic errors of a mobile LiDAR system.

CONCLUSIONS AND FUTURE WORKS
In this paper, we analyze eight systematic errors of a terrestrial MLS with a simulation approach. We introduced the concept, design, and methodology of the MLS simulator. The advantage of adopting this simulation approach is to be able to study each systematic error independent from other errors (e.g., trajectorybased errors and sensor noises). On the other hand, the simulator can generate point clouds on hypothetical planar targets with various configurations without any limitation. The term configuration in this paper mainly refers to the orientation of the planar target and the set of passages, known as a line pattern followed by an MLS. To detect the best configuration of the systematic errors of an MLS, we introduced two types of systematic-error visibility criteria, which implement a raw and a fine selection on the existing configurations. The results show that to have strong visibility on the systematic errors of a terrestrial MLS, we need to have two planar targets of 1m x 1m size with vertical and inclined orientations combine with a five passages line pattern. This result was confirmed with a real data acquisition test. The proposed configuration in this research will lead to relevant information for the calibration of a terrestrial MLS. Furthermore, the approach adopted in this paper shows that we can study the behavior (here, systematic errors) of a mobile LiDAR system (here, a terrestrial MLS), merely with a conceptual mathematical model, without the need to invest in expensive and time-consuming field tests.
In future works, we will consider more variables for the configuration parameters like the dimension of planar targets and also various platforms (aerial or marine) with multiple motion scenarios (position and orientation). Consequently, we will generate a lot of more possible configurations that will not be easy to analyze by using the approach proposed in this paper. Thus, we will adopt machine learning techniques to optimize the selection of the best configuration with strong visibility of the systematic errors of a mobile LiDAR system.

ACKNOWLEDGMENTS
This study was financed by PSR-SIIRI-953 and MITACS Acceleration in a research program that focuses on "Automatic calibration and performance analysis of a mobile LiDAR system" through a partnership between Université Laval and CIDCO (Centre Interdisciplinaire de Développement en Cartographie des Océans), an NPO company.