Abstract
The calibration of the light detection and ranging (LiDAR) system is critical to ensure the accuracy of point data. In this paper, the lever-arm measurement of airborne LiDAR system (ALS) was realized by photogrammetry. An automatic iterative boresight calibration method based on approximate corresponding points (CPs) matching was proposed to correct the boresight misalignment. It was based on iterative closest point (ICP) registration algorithm with a normal space sampling strategy, and approximate CPs were obtained by establishing filter rules. The experimental results showed that the absolute accuracy of the calibrated ALS reached 7.13 cm when the flight altitude was 100 m, meeting the accuracy requirements.
1. Introduction
ALS is widely used in topographic mapping, digital city, power line inspection, and other fields. A laser scanner, a global position system (GPS), and an inertial navigation system (INS) were integrated into ALS. It calculates the distance between the scanner and the target by the time interval from transmitting laser pulse to receiving echo signal. Then, the 3D spatial coordinates of the target are calculated by combining the position and attitude of the laser scanner observed by the GPS and the inertial measurement unit (IMU) [1]. As a multisensor integrated system, the accuracy of laser point cloud data acquired by ALS is affected by various factors such as GPS positioning error, IMU angle measurement error, system integration error, and laser ranging error [2]. Although the manufacturer calibrates ALS equipment before leaving the factory, the transportation, disassembly, and installation processes will inevitably impact the different components of the airborne LiDAR system. Boresight (angular misalignment between the mounting axes of the laser scanner and IMU) and lever-arm (physical offset from the laser scanner to GNSS antenna) can occur in the system [3]. Boresight misalignment and lever-arm offsets are the primary error sources that bias the LiDAR point cloud positioning [4]. Their presence causes noncoincidence and misalignment between LiDAR point cloud strips, which seriously affects the overall accuracy and precision of the data [5, 6]. Given the above, high-precision geospatial applications of UAVs especially require calibration of the boresight misalignment and lever-arm offsets [7].
In the local coordinate system of IMU, the lever-arm offset is constant, usually measured by a straightedge or obtained from the design drawings. The boresight angles cannot be measured directly, so they should be calibrated by manual calibration (MC) or least-squares adjustment with the collected point cloud [8]. Manual calibration is performed by laying out suitable routes to obtain laser point clouds of characteristic features such as flat roads and gable-roofed buildings and then gradually computing the three error angles according to the empirical formula. The boresight angles are computed iteratively until the point clouds obtained from different strips overlap well [9]. Due to the limited number of point clouds selected for each profile, this method is susceptible to factors such as flight quality and characteristic ground conditions. Therefore, MC is time-consuming and challenging to ensure calibration accuracy.
The least-squares adjustment-based calibration methods can be divided into data-driven methods and strict model-driven methods, depending on whether the boresight angles are directly represented by the geometric reference equation (8). The methods described in Mass (2000) [10] and Ayman et al. (2010) [11] are data-driven approaches that use only the position information of georeferenced LiDAR points to reduce the interstrip variance. This approach does not require raw data, such as GNSS/INS measurements, which are not always available to the user. However, these methods are considered less rigorous because arbitrary transformations cannot compensate for all the biases associated with point cloud geo-alignment. The rigorous model-driven approach expresses the apparent axis angle explicitly in the georeferenced equations. The boresight angles are estimated by constraining the set of feature points (such as connection points/control points and plane points) obtained from multiple overlapped strips or other reference data to fit the corresponding geometric model and minimizing the weighted sum of squares of residuals [12]. To quantify the discrepancies between strips, conjugate tie points, planar patches, and/or modelled surfaces is typically matched in overlapping point clouds from the strips. These differences are then used to determine unaligned criteria that are typically defined as distances [13]. Compared with the image data obtained from aerial photogrammetry, ALS point cloud data do not have real CPs due to the irregularity and dispersion of spatial distribution, making it challenging to adjust directly by the “connection point” constraint [9]. The parity model based on conjugate plane matching is the most used in boresight correction. Filin and Sagi (2003) [14] first recovered the boresight in an ALS system using natural and artificial surfaces. They used the least-squares adjustment (LAS) to estimate the boresight angles from surfaces with known parameters in different directions. Skloude and Lichti (2007) [4] improved the boresight calibration based on plane features by estimating the plane parameters and the boresight angles simultaneously instead of using fixed, known plane parameters. The improved method significantly enhanced the calibration performance. Identifying and selecting conjugate surfaces require rigorous preprocessing steps (e.g., region growing, principal component analysis, or RANSAC) and suitable sites, such as urban environments [15]. Point-to-patch matching methods [16, 17] advantageously present direct and automated communication. Once the correspondence is established, the boresight angles are represented as a parameter of the optimized target. LSA was then performed to estimate the correction parameters. Glira et al. (2016) [17] proposed a LiDAR strip adjustment method. Their method was able to correct for deviations in the LiDAR-GNSS/INS calibration parameters and other system errors such as laser ranging deviations or scale factors and biases in GNSS/INS observations. Their approach minimized the discrepancy between robustly selected point-to-plane correspondences from overlapping LiDAR strips by LSA. Recently, Keyetieu and Seube (2019) [18] discussed the optimal selection of strip-adjusted LiDAR observations. The authors relied on model measurement uncertainty of georeferenced LiDAR points to achieve the smallest LSA problem size.
Generally speaking, the current calibration method using line/plane features has specific requirements for the normal direction and distribution of line/plane features. In order to achieve higher calibration accuracy, the plane should be perpendicular to three coordinate axes as much as possible, and the distribution of corresponding features should be uniform as much as possible [19]. These conditions are only applicable to areas with rich geometric features in urban areas. It is tough to automatically extract precise lines and planes in complex environments, limiting these algorithms’ applicability.
Compared with the extraction and matching of line/plane features, point features are more readily available in natural scenes and easier to realize automation. The adjustment method based on point features has potential advantages. It will help improve the practicality and automation of the calibration method if the CPs can be utilized for correction parameters calculation. However, due to the nature of LiDAR data, automatic identification of CPs is unreliable, restricting the development of a strict adjustment model with CPs. The ICP algorithm [20, 21] provides a way to accurately align the source point cloud and target point cloud due to the small boresight angles in ALS.
In view of the above, a strict adjustment model for boresight based on ICP and approximate CPs was presented to solve the challenge of automatic selection and matching CPs between LiDAR strips. The proposed method eliminates the dependence on prelaid calibration markers or targets. It can operate automatically even in natural scenes without important geometric features, such as planes in multiple directions. Further, benefitting from the rigorous screening process for CPs, the proposal method can effectively reduce the probability of incorrect convergence.
2. Materials and Methods
2.1. Data Acquisition
We carried out three flight experiments at the China University of Mining and Technology (Figure 1). In the three experiments, the flight heights were 150 m, 100 m, and 80 m, respectively, and the flight speeds were 8 m/s, 6 m/s, and 6 m/s, respectively. Two of the four flight belts were parallel, and the other two were perpendicular to them. Ten ground control points were set up in the two experimental areas before the flight, and a GPS receiver was used to measure them. GNSS/IMU data was processed by PosPac (v8.5) software to obtain attitude and position of the aircraft and time matched with the ranging data of the laser scanner. After rotation and translation, the final multistrip initial LiDAR point cloud was obtained.
2.2. Sensor Payload
All sensors were integrated into a mature industrial-grade UAV, a WIND 4 (DJI Technology Co., Ltd., Shenzhen, China). The integrated sensors include a RoboSense RS-LiDAR-32, a coupled GNSS/IMU sensor with a multifrequency, multiconstellation GNSS receiver (Applanix APX-15), and an onboard computer. We fixed the GNSS/IMU sensor and the on-board computer in a protective shield and attached the laser scanner to the protection by screws. The antenna was fixed on the side of the UAV and was 10 cm above the propeller to ensure the uninterrupted reception of satellite signals (Figure 2). Our installation ensured that the three axes of the UAV coordinate system, the IMU coordinate system, and the laser scanner coordinate system were nearly parallel, reducing the difficulty and error of postprocessing.
2.3. Antenna to Laser Scanner Lever-Arm Offsets
Measuring the lever-arm offsets from GNSS antenna to laser scanner is a necessary step for ALS self-calibration. Kersting et al. (2012) [16] claimed that the vertical lever-arm component could not be estimated by looking at the strip-to-strip correspondence alone since errors in the vertical lever-arm parameters have the same effect regardless of the flight direction or flight altitude. The vertical lever-arm offset can only be computed when at least one vertical ground control point (GCP) is employed. This limitation was again emphasized by Ravi et al. (2018) [22]. An accurate measurement of the offset vector from the GNSS antenna to the laser scanner center is required to eliminate lever-arm errors. However, it is not easy to ensure the exact alignment of the measurement direction with the coordinate system using manual measurements. Therefore, we took 293 highly overlapping photos of the UAV from different angles and heights using a cell phone. The ALS was finely modelled with the help of the structure-from-motion (SfM) method. It can be rotated randomly in space to correct misalignment to achieve accurate measurements of the lever-arm offsets [3]. Due to the smooth and reflective nature of the UAV and sensor housing, the texture was not sufficiently clear at specific angles. To ensure 3D reconstruction accuracy, we set up a total of 23 high-contrast marker points on the UAV and ground. Ten of them were measured with steel ruler at different heights as control points to assist aerial triangulation. The overlapping images were then processed using Agisoft’s Photoscan (v1.2.6) (http://www.agisoft.com/) to obtain a dense ALS 3D point cloud (Figure 3). High-quality and gentle filtering patterns were set to create dense point clouds. The dense point cloud was scaled with the measured GCPs to produce an absolute size ALS model. This scaled point cloud was then imported into Cloudcompare (v2.11.3) (http://www.cloudcompare.org/) to measure the lever-arm offsets.
2.4. Boresight Self-Calibration
We proposed a strict boresight calibration method based on approximate CPs matching. This method consists of two parts: approximate CPs matching and estimation of boresight angles using LSA. Taking the point clouds scanned by two strips as an example, the calculation process of boresight angles is shown in Figure 4.
2.4.1. General Formulation
The calibration of the ALS parameters is based on the laser point positioning equation. First, the target coordinates in the laser scanner coordinate system are converted into IMU coordinate system by using the lever-arm offsets and boresight angles between the laser scanner and GNSS/IMU. Then, the IMU coordinate system coordinates are converted to a local horizontal coordinate system using the traverse roll, pitch, and heading angle provided by GNSS/IMU. Finally, the local coordinates are converted to WGS84-ECEF using the latitude, longitude, and ellipsoidal altitude supplied by GNSS/IMU. The georeferenced equation describing the conversion of the scan point from the laser scan frame to the WGS84 ECEF frame is given.
where is the positioning vector of the point in the WGS84 ECEF frame; is the positioning vector of the point in the laser scanner frame; is the rotation matrix from the laser scanner frame to the IMU frame; is the position of the laser scanner in the IMU frame; are the boresight angles between the laser scanner and the GNSS/MU system; is the lever-arm offsets vector between the laser scanner and the GNSS/IMU system; is the rotation matrix between the IMU frame and the navigation coordinate system; is the rotation matrix between the navigation coordinate system and the WGS84 ECEF frame; is the position of the UAV in the WGS84 ECEF frame; are the latitude, longitude, and ellipsoidal height provided by the GNSS/IMU; are the roll, pitch, and heading provided by the GNSS/IMU.
The lever-arm offsets can be measured directly from the scaled model with sufficient accuracy so that the error caused by the lever arm offset error is negligible. Under good GNSS conditions, the positioning error and IMU attitude error have little effect on the positioning accuracy, while the boresight angles mainly determine the total error of the point cloud. , , and are the boresight angles in the roll, pitch, and heading directions, respectively. Since the boresight angles are small, generally in the range of -3° to 3°, the geolocation equation can be described as follows.
where
If the point is scanned from the -th strip and its corresponding point in the -th strip is , the deviation between the points and can be calculated by the error equation (5).
Suppose corresponding points are found from different scan directions or different strips. In that case, the boresight correction parameters can be derived by minimizing the total deviation of the corresponding points as follows.
2.4.2. Approximate CP Strategy
Unlike grayscale images, there are no true CPs in the point clouds of different strips due to the randomness and discreteness of LiDAR points. This property of LiDAR points is the most critical reason limiting strip adjustment models based on CPs matching. Of course, some scholars have proposed the intersection of the feature line with the surface as the feature point at this stage. However, such methods require filtering and fitting the points, which is challenging to implement in complex LiDAR points.
ICP algorithm is the mainstream algorithm in the fine alignment stage of point clouds. For two sets of point clouds with good initial positions, the ICP algorithm can achieve the ideal alignment effect. Since the boresight angels of ALS is generally minor, the different strip point clouds collected in a short time can be obtained by rigid transformation when the UAV flight altitude is low. Therefore, in this paper, the following matching strategies for different strip approximation CPs are proposed based on the ICP algorithm. (1)Removal of points in the point cloud that are too far from the overall point cloud using distance filtering, which is a prerequisite for fine alignment using ICP(2)The point cloud on each strip is split into small pieces using timestamps at 5 s intervals. The time interval is chosen according to the flight speed and should not be too large so that the large-scale point clouds can be better processed and the corresponding point clouds can be better aligned accurately by ICP(3)After partitioning the point cloud of each strip into small sections, the 3D bounding boxes of all sections are calculated. Then, they are used to detect the overlapping regions to construct the overlapping relationship between the parts of the point cloud. The alignment is performed between the sections that overlap each other(4)A large amount of data in a dense point cloud makes effective alignment of overlapping regions difficult. The normal vector is an essential geometric feature of the surface, and the angle of the normal vector can be used to determine the surface change at that place. To improve the matching efficiency, we adopt the normal vector space sampling strategy [23] in the ICP variant algorithm for point cloud refinement (Figure 5) and then use the ICP algorithm to align the overlapping sections(5)After alignment, 5000 points () are sampled uniformly from the source strip, and the three points () with the closest Euclidean distance to each sampled point are selected from the other strips ()(6)If the maximum height difference of the three points in is greater than 30 cm, they are removed from and to obtain the and
(a)
(b)
Two north-south flight strips in each experiment were taken as an example. First, the initial point cloud of each strip was filtered by distance to eliminate the outliers. Then, the point cloud was divided into small pieces to calculate the overlapping area of the strip point cloud (Figure 6). The normal vector space of the point cloud was computed, and the streamlined point cloud was obtained by uniform sampling in the normal vector space. The ICP alignment was performed on the streamlined point cloud to get the aligned point clouds of the two strips. Five thousand points were randomly selected from the overlapping area of the source point cloud (Figures 5 (a), 5(c), and 5(e)). The sampling points selected in experiments 1-3 were filtered according to the strategy proposed in Section 2.4.2 to obtain , , and , where , , and (Figures 5(b), 5(d), and 5(e)). The three points closest to each point in were extracted from the target strip to form . Therefore, , , and in experiments 1-3 are 7035, 9882, and 10866, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
It is worth noting that the convergence of ICP alignment is crucial to obtain reliable alignment results. We adopted three steps to improve the accuracy of ICP alignment: first, after removing the outliers, the point cloud was divided into small blocks to ensure the consistency between the source point cloud and the target point cloud through a rigid transformation; second, the feature points were selected by uniform sampling in normal vector space. This sampling strategy can preserve the local features of the point cloud surface while extracting the thin point cloud, which is more favorable for accurate matching; third, all boresight angles were prealigned. The prealignment makes the initial value of the placement angle error close to zero, allowing the algorithm to converge efficiently.
In addition, several parameters in the approximate CPs filtering process need to be set reasonably according to the actual scenario. For example, we sampled 5000 points uniformly from the source point cloud, and each point took the three points with the closest Euclidean distance in the target point cloud to constitute a set of approximate CPs. To get more CPs, the readers can increase “5000” or “3” appropriately. However, we recommend selecting at least three nearest points because the positioning accuracy of the points is greatly affected by the random error generated by the scanning angle and measuring distance of the laser scanner. The selection of only one nearest neighbor point will introduce many mistakes to the settlement of the boresight error. As for the elevation difference threshold, it needs to be adjusted according to the point cloud density. We suggest that when the point cloud density is lower than 150/m2, the elevation difference threshold should be increased appropriately.
2.4.3. Boresight Angle Correction Parameter Estimation
Because of the boresight error, the CPs in the overlapping strips are not coincident. Therefore, the calibration procedure uses the CPs position consistency as a constraint to recover the Boresight angles correction parameters. When the systematic error is eliminated, all point clouds are restored to their true positions.
Assuming that the CPs of the two strips are and , respectively, the observation equation can be expressed as
According to equation (2), the laser point coordinates are expressed as a function of the system parameters to be solved. After linearizing the equation, the error equation (8) can be listed by setting the coordinate residuals as and solved according to the least-squares principle.
where is the coefficient matrix of the unknown, ; is the matrix of boresight angles, ; is a constant matrix, .
When the value of the expression is the most minor, the correction parameters of the boresight angles can be obtained according to equation (9).
In order to solve the boresight angles correction parameter, it is only necessary to determine the corresponding laser points on different strips to estimate the system parameters according to the least-squares principle. The number of unknowns to be solved in Eq. (9) is three, and three equations can be listed for each pair of CPs, so at least one pair of CPs is needed for the solution. Multiple sets of boresight angle correction parameters can be derived by applying the above LSA model using the approximate CPs filtered by Section 2.4.2. However, because some points are located at the edge of the scanned area, vegetation area, or are too much affected by random errors, the accuracy of the obtained correction parameters is not enough or even completely wrong.
Therefore, we counted the frequency of the calculated , , and according to the gradient of 0.1°, respectively, and recorded the most frequent group of each as , , and . All CPs, whose corresponding in , were selected from to form set . The boresight angles in were averaged to obtain the final correction parameters. The calibration parameter solution is not performed only once but repeated iteratively until the resulting boresight error is below the set accuracy threshold.
2.5. ALS Alignment Evaluation
We laid some checkpoints on the ground in the experimental area. Taking each checkpoint as the center, the corrected point cloud in the sphere with a radius of 20 cm was intercepted. The point cloud intercepted at each checkpoint is plane-fitted, and is defined as the orthogonal distance from the -th checkpoint to the corresponding fit plane. If the ALS error correction accuracy is high enough, then should be close to zero. Two metrics and are introduced to evaluate the ALS error correction accuracy [24]. is defined as the average of all calculated for each checkpoint. is the average of all , where is the standard deviation of the distribution of the set of points around the plane fitted to the point cloud intercepted by each checkpoint. evaluates the absolute accuracy of the alignment, but evaluates its relative accuracy.
3. Results and Discussion
3.1. Lever-Arm Correction
The lever-arm offsets between sensors were measured in a scaled LAS point cloud, as shown in Table 1. The accuracy of the lever-arm offsets estimation depends on the accuracy of the GCPs measurement, the inherent error of the SfM modelling, and the estimation accuracy of the sensor reference points within the point cloud. The GCPs were measured with a rigid ruler and a lead hammer. The uncertainty of the rigid ruler measurement was within ±0.1 cm, while the accuracy of the lead hammer measurement was somewhat worse, around ±0.2 cm. Therefore, the maximum uncertainty of the GCP measurement was considered to be ±0.3 cm. The uncertainty of GPCs combined with the uncertainty in the intrinsic modelling of the SfM was translated into uncertainty in the 3D model. Since the images were taken at a very short distance of less than 1.5 m, the SfM inherent uncertainty was considered within the measurement accuracy of the GCPs. Thus, the uncertainty in the ALS 3D model was considered to be within ±0.5 cm. When the uncertainty of the selected points on the ALS model point cloud was added, the final calculated total uncertainty associated with the lever-arm offsets estimate was ±0.6 cm. Manual measurement is difficult to ensure that the measurement direction of the straightedge is in the same direction as the axis of the IMU coordinate system, so the SfM-based lever arm measurement method is more accurate than the manual measurement method.
3.2. Boresight Angle Correction Parameter Estimation
The least-squares adjustment model was applied to the first filtering results of experiment 1-3 to obtain 7035, 9882, and 10866 boresight angle correction parameter matrices, respectively. Taking experiment 2 as an example, , , and were counted as frequencies according to a gradient of 0.1° (Figure 7). The sampling points corresponding to [, ,] in [, , ] were picked in to get (). Similarly, we got and . The boresight angle correction parameters for the first iteration of experiments 1-3 were derived by averaging , , and obtained from the approximately CPs in ,, and , respectively.
The above process was iterated until the boresight angles obtained are less than the set threshold value (Figure 8). All three groups of experiments reached a steady state after 6-7 iterations. The convergence speed of experiment 3 was the fastest among the three groups. This is due to the fact that experiment 3 has more approximate CPs involved in the calculation. It is also related to the low flight altitude and small deformation between strip point clouds in experiment 3. The final cumulative boresight angle errors were shown in Table 2.
(a)
(b)
(c)
We corrected the point cloud using the final boresight angles (Figures 9–11). The points scanned from different strips had different colors, and there were obvious deviations between the original point clouds of different strips before calibration. To more intuitively and qualitatively evaluate the correction effect, seven point cloud profiles before and after boresight calibration were visualized (Figures 9–11). As seen from the figure, the point clouds of different strips after calibration by our method had achieved complete overlap, and the accuracy of the point clouds had been greatly improved.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
Finally, we quantified and evaluated the corrected point clouds using the ground checkpoints arranged in advance and conducted comparison tests with the virtual tie point model-based method (VTPM) proposed by Jing et al. (2013) [25] and the MC. The results were listed in Tables 2 and 3. Table 2 showed that the final boresight angles of three experiments calculated by the three methods were close. Table 3 showed that the point cloud corrected by the two automatic methods was more accurate than MC. Considering the random errors of the laser rangefinder range, GPS positioning, and attitude, we therefore believed that the boresight angle calculation results met the expectation of error correction. Among the three sets of experiments, the worst point cloud accuracy in experiment 1 was mainly due to the high flight altitude. The higher flight amplified the attitude and ranging random errors, resulting in a more discrete point cloud. It was also illustrated by the fact that is larger than . In experiment 1, the point cloud accuracy obtained by the VTPM method was better than our method. This was due to the fact that a large percentage of the 5000 sampling points uniformly generated by our method in experiment 1 were located in the plant areas and the edge areas of buildings. Few CPs remained after sampling point filtering (<1000). For this scene, the number of sampling points can be increased appropriately to improve the accuracy of our method. In experiments 2 and 3, the accuracy of our proposed method was much less time consuming and better than that of VTPM, mainly because our CP filtering process removed most of the pseudocorresponding points that were incorrectly matched due to random errors, speeding up the iterative process and avoiding incorrect convergence.
4. Conclusions
The calibration of ALS is a critical step before performing data postprocessing and directly affects the quality of the point cloud. We used a cell phone to finely model the ALS using the SfM method to measure the lever-arm. In addition, we proposed an automated, fast iterative boresight correction algorithm based on approximate corresponding point matching. It overcomes the difficulty of accurately extracting corresponding points from LiDAR point cloud. It enables the potential point matching adjustment model to be successfully applied to the boresight angles correction. Our research results showed that the accuracy of the proposed algorithm was ideal, and we believed it could be better promoted in the future.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare no conflict of interest.
Acknowledgments
The research was supported by the projects funded by National Key R&D Program of China, grant number 2017YFE0119600; the Fundamental Research Funds for the Central Universities, grant number 2019BSCX26; and the Slovak Research and Development Agency, grant number APVV-SK-CN-RD-18-0015.