AERIAL-TRIANGULATION AIDED BORESIGHT CALIBRATION FOR A LOW-COST UAV-LIDAR SYSTEM

The Laser-IMU boresight calibration is the precondition for an Unmanned Aerial Vehicle (UAV)-Light Detection and Ranging (LiDAR) system (ULS). The existing methods achieve good performance for calibrating ULSs with high-precision positioning and orientation systems (POS) (e.g., APX-15), in which, the systematic errors of the high-precision POS can be ignored, only the boresight parameters are estimated. However, these methods have difficulties in calibrating the low-cost ULSs with low-precision POS. To overcome the impact of the systematic errors of the low-precision POS on boresight calibration, an aerial-triangulation aided boresight calibration is proposed in this paper. It simultaneously estimates the laser-IMU boresight angles and system states (e.g. trajectory) by setting the point clouds derived from aerial-triangulation (AT point clouds) as the reference. Firstly, the planar voxels from the AT point clouds are extracted, due to the fact that they are more reliable in AT point clouds. Secondly, raw laser observations are matched with the extracted planar voxels to establish laser matching observations. Thirdly, a Dynamic Network (DN) is built using the GNSS observations, inertial observations, and laser matching observations to simultaneously optimize the initial laser-IMU boresight angles and the system states. All the sensor observations involved in the ULS are modeled with proper error models, which are essential for analyzing and refining the data quality of the low-cost ULS. The proposed method was tested to calibrate a low-cost ULS, KylinCloudII, in a calibration field. It showed that the average distance between the laser point clouds and the referenced AT point clouds was decreased from 2.560m (RMSE = 3.88m) to 0.08m (RMSE = 0.99m).


INTRODUCTION
Unmanned Aerial Vehicle (UAV)-Light Detection and Ranging (LiDAR) Systems (ULS) are widely used in many applications, such as building information modeling (Roca et al., 2014), urban changes detection (Qin et al., 2016), power-line corridor inspection (Chen et al., 2018), forest inventory estimation (Jaakkola et al., 2010;Li et al., 2019b;Liu et al., 2018;Wallace et al., 2012). Nevertheless, the costs of high-end ULSs are still high and unacceptable for non-professional users, compared with those of UAV photogrammetric systems.
Developing low-cost ULSs is conducive to the promotion of the ULS point cloud based applications (Habib, 2017), so it has attracted extensive attention from the communities of academy and industry in recent years (Jaakkola et al., 2017(Jaakkola et al., , 2010Li et al., 2019a;Ravi et al., 2018a). Wallace et al., (2012) developed a low-cost ULS with structure from motion (SfM) aided orientation, which was used for forest inventory successfully. Ravi et al. (2018b) integrated a ULS with a low-cost multi-beam laser scanner and a high-end APX-15 POS, and calibrated the boresight parameters with plane fitting accuracy in terms of 0.02m. Li et al. (2019a) integrated a low-cost ULS system, named KylinCloud-I, which equipped with a low-end MEMSbased IMU, achieving plane fitting accuracy in terms of 0.05cm.
Boresight calibration, estimating the laser-IMU mounting angles, was a core step for mobile mapping systems (Glennie, 2013;Skaloud and Lichti, 2006). Filin (2003) calibrated the laser-IMU boresight angles utilizing natural and man-made surfaces for an airborne LiDAR system (ALS). Zhang et al. (2015) calibrated the laser-IMU boresight angles using aerial triangulated images * Corresponding author as control information by matching corresponding feature points extracted from images and ALS points. The laser-IMU boresight angles can also be self-calibrated without reference data. Skaloud and Lichti (2006) calibrated laser-IMU boresight angles of an ALS using plane features from different strips by estimating plane parameters and boresight angles simultaneously. Chan et al. (2013) calibrated laser-IMU boresight angles of a mobile laser scanning (MLS) system using multi-features (plane feature and catenary feature). Glennie (2013) simultaneously calibrated the laser-IMU boresight parameters and interior scanner parameters using planar features. The trajectory errors of ULSs equipped with high-end POSs could be ignored (Li et al., 2019;Ravi et al., 2018b), then the laser-IMU boresight angles were estimated using the corresponding geometry features selected from different strips. The calibration accuracy is analyzed according to the bias impact analysis (Ravi et al., 2018b). However, the trajectory errors of the low-cost ULSs equipped low-end POSs can not be ignored, posing difficulties for the existing calibration methods. This problem was partly solved by integrating a camera in a ULS to improve the trajectory accuracies with the aid of aerial-triangulation (Elbahnasawy et al., 2018;Li et al., 2019a). However, as far as a low-cost ULS without a camera, the laser-IMU boresight calibration problem was still an open question.
To handle the boresight calibration problem for low-cost ULSs, system states (e.g. trajectory) and boresight parameters should be estimated together (Qin and Shen, 2018;Yang and Shen, 2017). Dynamic Networks (DN) was used as the basic framework for the parameter estimation in this paper. DN was firstly proposed to simultaneous modeling and adjusting of raw inertial observations with optical and GNSS data streams (Rouzaud and Skaloud, 2011). It was a post-processing framework, which was superior to the smoother based method for the integration of optical sensors (Rouzaud and Skaloud, 2011). Based on DN, Cucci et al. (2017a) proposed a bundle adjustment with raw inertial observations to estimate the UAV trajectories and the camera-IMU boresight angles through a dedicated calibration flight.
In this paper, an aerial-triangulation aided boresight calibration method is proposed to calibrate the laser-IMU boresight angles for a low-cost ULS equipped with low-end POS. Setting the AT point clouds as the reference, the planar voxel in the AT points are extracted first and used for the following calibration process, for the reliability of the planar points in AT point clouds. Then the laser point clouds are matched with the extracted planar voxels according to the distance. At last, all the sensor observations (GNSS, IMU, and laser matching) are integrated into a DN to optimize the involved parameters. Different from the existing methods, the proposed method considers the systematic error of the low-end POS, and estimates the laser-IMU boresight and system states simultaneously. In this manner, the sensors' error models are properly modeled, which are essential for analyzing and refining the data quality of the lowcost ULS.
The remainder of this paper is organized as follows. In Section 2, the steps of the proposed method are elaborated in detail. In Section 3, experiments were undertaken to evaluate the performance of the proposed method, after which conclusion is drawn at the end of this paper.

THE PROPOSED CALIBRATION METHOD
The proposed aerial-triangulation aided boresight calibration mainly consists of three steps: (1) planar voxel extraction (Section 2.2); (2) establishing laser matching (Section 2.3); (3) building Dynamic Networks (Section 2.4), as illustrated in Fig.1. For the convenience of method description, the system description of the low-cost ULS is first detailed as follow:

System description of the low-cost ULS
The low-cost ULS, named KylinCLoud-II, is composed of a multibeam laser scanner, a GNSS receiver, and a low-cost MEMS-based IMU. The multibeam laser scanner, R-FANS16 1 , is with 250 m measurement range and 2 cm distance accuracy. The GNSS receiver is Novatel OEM-719 2 . The gyro in-run bias of the low-cost IMU is 6 deg/h. The whole system was mounted on a multi-rotor UAV, DJI M600 3 , produced by DJI.
Similar to the definition of the common ULS (Ravi et al., 2018b), three coordinate frames are involved in the proposed method, namely, the mapping frame , body frame , and 16 laser frames for each scan lines { ( ) , = 0, … ,15}, as illustrated in Fig.2. A laser point in mapping frame is calculated using the observations by Eq.1 where, ( ) is the point observed by the ℎ laser scan line; ( ) is boresight rotation matrix for the ℎ laser scan line; is the level-arm vector for the laser scanner.
( ) and ( ) are the exterior parameter at time for the system.

Planar voxel extraction
The AT point clouds are used as the reference data. Due to the fact that the AT point clouds usually contain a lot of noise, only the planar points in the raw AT point clouds are used. To save computation cost, an octree with a 5 m grid size is built for the AT point clouds. Then principal components analysis (PCA) is performed for points located in each voxel to obtain the three eigenvalues 1 > 2 > 3 . The planarity of each voxel is evaluated by Eq.2.
Only the voxels, whose planarity more than 0.8, are selected for the following steps. For each selected planar voxel, plane fitting is performed. A 4× 1 parameter vector 4×1 are involved in each plane.

Establishing laser matching
Once the planar voxels are extracted from the referenced AT point clouds, we establish the matching between the raw laser 3 www.dji.com/cn/matrice600  observation and the planar voxels. Using the direct referencing equation Eq.1, we project the raw laser observation into the mapping frame according to the initial laser-IMU boresight parameters and initial trajectory obtained by a loosely coupled GNSS/IMU integration via Inertial Explorer. For each projected laser point, we find a closest planar voxel for it. If the distance between the projected laser point and the closest planar voxel is smaller than 5 m, the laser point, and the closest planar voxel construct a laser matching. However, there may exist a lot of outliers in these selected matches, which will be handled in the solving step (Section 2.4.5). This laser matching means the laser point is located on the plane presented by the planar voxel.

Building Dynamic Networks
DN is first introduced in (Rouzaud and Skaloud, 2011), it could simultaneously model and adjust raw inertial observations with optical and GNSS data streams. Compared with the traditional Kalman filter basded integration (Zhu et al., 2019) or GNSS/IMU aided bundle adjustment (Yuan et al., 2009), DN is a rigorous approach for the complex time-and space-dependent measurements from different sensors. Two types of equations, namely, static-condition equation and dynamic-condition equation are involved in the DN. The static-condition equations are established using GNSS measurements and optical measurements (e.g. from camera, laser scanner, and etc.) by back-projecting the mesurements to the time-dependent system states. The dynamic-condition equations are established using gyroscope and accelerometer measurements by differential the adjacent system states. In this manner, all the mesurements are took into a tight integration framework. For more details of DN, we refer the reader to (Cucci et al., 2017). In this work, we use DN to fuse GNSS data, IMU data, and laser matching for the laser-IMU boresight estimation.

Unknowns
The unknows in this problem consist of (1) System state ( ) where, ( ) is a quaternion corresponding to ( ). ( ) and ( ) are the biases of the gyroscope and the accelerometer, respectively.

GNSS positioning and velocity
The raw GNSS observations are processed using Inertial 4 www.novatel.com/products/software/inertial-explorer Explorer software 4 via carrier-phase differential postprocessing. Then the position observation and velocity observation of the GNSS antenna are obtained. The GNSS positioning nominal equation is written by: where, is the level arm vector between the body frame and the GNSS antenna.
( ) is the position of the GNSS antenna. GNSS positioning error equation is written by: where, is the observation error for position observation . The GNSS velocity nominal equation is written by: where, ( ) is the velocity of the GNSS antenna. GNSS velocity error equation is written by: where, is the observation error for velocity observation .

Angular velocity
The angular velocity is observed by the gyroscope at 200 Hz. The nominal equation for the angular velocity is written by: where, Ω is the rotation rate of the body frame. is the Earth rotation rate.
is the rotation between Earth-centered-Earth-fixed (ECEF) frame to the mapping frame. Ω is the angular velocity in the body frame. Then the angular velocity error equation is written by: where, is the observation error for angular velocity .

Specific force
The specific force is observed by the accelerator at 200 Hz. The nominal equation for the specific force is written by: where, is the acceleration of the body frame presented in the body frame. ̇( ) is the differential value of the velocity ( ) . is the gravity in the mapping frame. Then the specific force error equation is written by: where, is the observation error for specific force .

Laser matching
The laser scanning points are matched with the closest planar voxels. For a laser matching between the laser point and a planar voxel, the nominal equation is written by: which means a laser point is located on the plane. Then the error equation is written by: where, is the plane fitting error obtained from PCA. Due to the fact that there are a lot of outliers in the selected matches merely according to the distance (5m), we added a robust kernel function, Huber loss (Welsch, 1977), to every laser matching factor when building the DN.

Differential equations
Besides the sensor observation equations listed above, differential equations of the system states are also important and are listed below. As for angular velocity presented in the mapping frame ( ) at time , it is obtained as follow: where ∆ is the time interval between the two adjacent states (0.005s in this study). As for acceleration presented in mapping frame ̇( ) at time , it is obtained as follow: As for velocity presented in mapping frame ( ) at time , it is obtained as follow:

Solving
To summarize, all the unknowns, observations, and differential equations involved in the proposed method are included in a DN as illustrated in Fig.3. Circles in Fig.3 represent the unknowns to be estimated. Rectangles in Fig.3 represents the observations detailed above.
Because the DN is the non-linear system, the initialization of the system is a key factor for the final result. As for the laser-IMU boresight angles, we use the values obtained by the design drawings as the initial values. As for the initial values of the system states, we first perform a loosely coupled GNSS/IMU integration via Inertial Explorer. Then the initial values of the system states are set according to the results of the integration results. Fig.3 Dynamic Networks of the unknowns and observations involved in the proposed method.

Data Collection
To calibrate the parameters for the low-cost ULS, AT point clouds and ULS data were collected in a calibration field in Wuhan, China as shown in Fig.4 (a). The calibration field mainly consisted of planar structures. We collected images using a UAV, DJI Inspire 1 5 , as shown in Fig.4 (d). The flight height was about 200 m. The forward overlap and the side overlap were all set to 80% to ensure the data quality. The camera used in this UAV system is DJI zenmuse-x3 6 with a 4096x2160 image resolution. The approximate ground resolution of the images was 3 cm. The images were then processed using Pix4dmapper 7 to generate the AT point clouds, which was illustrated in Fig.4 (b).
The ULS point cloud was shown in Fig.4(c). The flight of the ULS was planned using DJI Ground Station Pro (GSP) 8 . The flight height and flight speed of the low-cost ULS were 120 m and 5m/s, respectively. The calibration field was collected with four strips of ULS. The laser point cloud density was over 150 points/m 2 . The laser point clouds rendered by collecting time were illustrated in Fig.5.

Results of planar voxel extraction
The AT point clouds were used as the reference data. As the AT point clouds usually contain a lot of noise as shown in Fig.6 (a), only the planar points in the raw AT point clouds were exracted and used for the following processing. The voxelization of the AT point clouds with planarity was illustrated in Fig.6 (b). The voxels, whose planarity is over 0.8, are selected for the following processing. 5 www.dji.com/cn/inspire-1-pro-and-raw 6 www.dji.com/cn/zenmuse-x3 7 www.pix4d.com/product/pix4dmapper-photogrammetrysoftware 8 www.dji.com/cn/ground-station-pro

Results of calibration
The proposed method was performed to estimate the laser-IMU boresight values for the low-cost ULS system. The initial values and the correction values of the laser-IMU boresight were listed in Table.1. To evaluate the result of the calibration results, we computed the distance between the laser point and the corresponding nearest point in the AT point clouds as shown in Fig.7. As for the laser point clouds calculated using initial values [ Fig.7(a)], the average distance is 2.560 m (RMSE = 3.88 m). After applying the values obtained by the proposed method, the average distance is reduced to 0.08 m (RMSE = 0.99 m). The resulted laser point clouds and the AT point clouds were aligned accurately as shown in Fig.8, which demonstrated the effectiveness of the proposed method.

Comparison with other methods
The proposed aerial-triangulation aided boresight calibration is focusing on the calibration problem of the the low-cost ULS equipped with low-end POS. Compared with the exisiting calibration methods (Radhika and Tamer, 2018;Ravi et al., 2018a), which mainly rely on selecting corresponding features from different strips, the proposed method needs additional control information from the aerial-triangulated images. However, the existing methods mainly based on the assumption that the error of the trajectory from the high-end POS could be ignored. As illustrated in Fig.9, we compared LiDAR point clouds geo-referenced using different strategies. The first one, Fig.9(a), is direct geo-referencing using calibrated boresight parameters and trajectory from low-end GNSS/IMU integration. The second one, Fig.9(b), is direct geo-referencing using calibrated boresight parameters and trajectory from the DN. Obviusly, the data quality of the second one is higher than the first one, which is mainly resulted by poor trajectory from the low-end POS. Thus, the existing boresight calibration methods are not suitable for the low-cost ULS.

CONCLUSION
In this work, we proposed an aerial-triangulation aided boresight calibration method to calibrate the laser-IMU boresight parameters for a low-cost ULS. Based on the DN framework, the proposed method simultaneously estimates the boresight parameters and system states (e.g., trajectory) with the reference of a pre-processed AT point clouds. The experiment was carried out to evaluate the proposed method in a calibration field. The calibration results showed that the average distance between the laser point clouds and the referenced AT point clouds is reduced from 2.560m (RMSE = 3.88m) to 0.08m (RMSE = 0.99m), which demonstrated that the accuracy and effectiveness of the proposed method.  System calibration is the first and primary step of integrating a low-cost ULS. In the future, we will explore the data quality refinement of the KylinCloud-II based on the calibration results obtained by the proposed method.