1. Introduction
A LiDAR system is composed of a laser ranging and scanning unit and a position and orientation system (POS), which consists of an integrated differential global positioning system (DGPS) and an inertial measurement unit (IMU). The laser scanner measures distances from the sensor to the ground, while the integrated GPS/IMU observations provide the position and attitude information of the platform. The coordinates of the LiDAR points are the result of combining the derived measurements from each of its system components, as well as the mounting parameters relating such components. The relationship between the system measurements and parameters is embodied in the LiDAR equation [
1,
2,
3], Equation 1. As it can be seen in
Figure 1, the position of the laser point,
, can be derived through the summation of three vectors (
,
and
) after applying the appropriate rotations:
Ryaw, pitch, roll,
RΔω,ΔφΔκ and
Rα, β. In this equation,
is the vector from the origin of the ground reference frame to the origin of the IMU coordinate system,
is the offset between the laser unit and IMU coordinate systems (lever-arm offset vector), and
is the laser range vector whose magnitude is equivalent to the distance from the laser firing point to its footprint. It should be noted that
is derived through the GPS/INS integration process while considering the lever-arm offset vector between the IMU body frame and the phase centre of the GPS antenna. The term
Ryaw, pitch, roll stands for the rotation matrix relating the ground and IMU coordinate systems, which is derived through the GPS/INS integration process;
RΔω,Δφ,Δκ represents the rotation matrix relating the IMU and laser unit coordinate systems (which is defined by the boresight angles), and
Rα, β refers to the rotation matrix relating the laser unit and laser beam coordinate systems with
α and
β being the mirror scan angles. For a linear scanner, which is the focus of this research work, the mirror is rotated in one direction only leading to a zero
α angle. The involved quantities in the LiDAR equation are all measured during the acquisition process except for the boresight angles and lever-arm offset vector (mounting parameters), which are usually determined through a calibration procedure.
Figure 1.
Coordinate systems and involved quantities in the LiDAR equation.
Figure 1.
Coordinate systems and involved quantities in the LiDAR equation.
The calibration process is usually accomplished in several steps: (i) Laboratory calibration, (ii) Platform calibration, and (iii) In-flight calibration. In the laboratory calibration, which is conducted by the system manufacturer, the individual system components are calibrated. In addition, the eccentricity and misalignment between the laser mirror and the IMU as well as the eccentricity between the IMU and the sensor reference point are determined. In the platform calibration, the eccentricity between the sensor reference point and the GPS antenna is determined. The in-flight calibration utilizes a calibration test field composed of control surfaces for the estimation of the LiDAR system parameters. The observed discrepancies between the LiDAR-derived and control surfaces are used to refine the mounting parameters and biases in the system measurements (mirror angles and ranges). Current in-flight calibration methods have the following drawbacks: (i) They are time consuming and expensive; (ii) They are generally based on complicated and sequential calibration procedures; (iii) They require some effort for surveying the control surfaces; (iv) Some of the calibration methods involve manual and empirical procedures; (v) Some of the calibration methods require the availability of the LiDAR raw measurements such as ranges, mirror angles, as well as position and orientation information for each pulse; and (vi) There is no commonly accepted methodology since the calibration techniques are usually based on a manufacturer-provided software package and the expertise of the LiDAR data provider.
As a result of the non-transparent and sometimes empirical calibration procedures, collected LiDAR data might exhibit systematic discrepancies between conjugate surface elements in overlapping strips. The elimination and/or reduction of the effect of systematic errors in the system parameters have been the focus of the LiDAR research community in the past few years. The existing approaches can be classified into two main categories: system driven (calibration) and data driven (strip adjustment) methods. This categorization is mainly related to the nature of the utilized data and mathematical model. System driven (or calibration) methods are based on the physical sensor model relating the system measurements/parameters to the ground coordinates of the LiDAR points. These methods incorporate the system’s raw data or at least the trajectory and time-tagged point cloud for the estimation of biases in the system parameters with the help of the LiDAR equation. In this paper, the term “raw data” is used to denote all the quantities present in the LiDAR equation (
i.e., position and orientation information as well as the measured range and scan angle for each pulse). Moreover, the utilized sequence of rotation angles defining the system attitude and boresight angles has to be defined. The access to the system raw measurements is usually restricted to LiDAR system manufacturers. Such a restriction has triggered the development of methods that utilize the XYZ coordinates of the LiDAR point cloud only. These procedures are classified as data-driven (or strip adjustment) methods [
4,
5,
6,
7,
8,
9,
10,
11]. The major drawback of such methods is the mathematical model employed. The effects of systematic errors in the system parameters are usually modeled by an arbitrary transformation function between the laser strip coordinate system and the reference coordinate system. The utilized transformation function might not be appropriate depending on the nature of the inherent biases in the LiDAR system parameters.
Although existing calibration procedures are based on the appropriate mathematical model (
i.e., LiDAR equation), they require the access to the raw data [
12,
13,
14], which may not be available for all end-users, or at least the trajectory and time-tagged point cloud [
15,
16,
17]. In [
12], natural and man-made control surfaces, represented by a set of planar surfaces, are utilized to estimate the system calibration parameters. More specifically, the calibration parameters are estimated by enforcing the coincidence of the LiDAR points and the control surfaces they belong to. Without the availability of control surfaces, the proposed procedure in [
12] cannot be applied. The control requirement is circumvented in [
13,
14,
15,
16,
17], where the system biases are estimated using overlapping strips. The method proposed in [
13] estimates the calibration parameters by conditioning a group of points to lie on a common plane. The utilized planes are selected manually. The calibration procedure estimates the system parameters as well as the parameters describing the involved planes. The proposed procedure in [
13] can only be applied whenever planar surfaces with varying slopes and aspects are available. Thus, the procedure is restricted to LiDAR data covering urban areas. Moreover, the number of unknowns changes with the number of planes used in the calibration procedure. Similar to [
13], the proposed procedure in [
14] is based on planar surfaces, which are automatically derived through a segmentation procedure. Pre-processing of the LiDAR point cloud to extract the planar patches might negatively affect the quality of the calibration procedure if it is not properly implemented. The calibration methods proposed in [
15,
16,
17] utilize trajectory and time-tagged point cloud only. In spite of the relaxed data requirements in these methods, they have some limitations. For instance, in [
15] point primitives are utilized to establish the correspondence between overlapping strips. Due to the irregular nature of the LiDAR points, the identification of distinct and conjugate points in overlapping strips is quite difficult and not reliable. Moreover, this method assumes that the true ground coordinates can be derived by averaging the coordinates of tie points in overlapping strips. The averaging process will not lead to the correct estimate of the LiDAR point in the presence of some systematic errors. In [
16,
17], only biases in the boresight angles are estimated in the calibration procedure. Moreover, the parameters describing the utilized surface in [
16] are estimated within the calibration procedure. Thus, the number of estimated parameters depends on the extent of the area being utilized in the calibration procedure. In [
17], the boresight angles are estimated using observed discrepancies between conjugate surface elements in overlapping LiDAR strips. The discrepancies are obtained through a matching procedure that works on interpolated data. The implemented matching procedure does not lead to accurate estimate of the horizontal discrepancies when compared to the estimated vertical discrepancy. Therefore, the estimated biases in the boresight pitch and heading angles are of lower accuracy than the boresight roll angle since the former lead to horizontal discrepancies while the latter leads to horizontal as well as vertical discrepancies among overlapping strips.
In this paper, two new calibration procedures for the estimation of biases in the system parameters that overcome the limitations of existing calibration procedures are introduced. The two presented methods differ in terms of data requirement to satisfy the needs of several users. The first procedure, denoted as the “Simplified method”, is designed to work with the LiDAR point cloud only. In spite of the increasing adoption of the LAS format, some of the mapping agencies still require only the delivery of the XYZ coordinates of the LiDAR point cloud. Therefore, several users only have access to the point cloud coordinates. Although the Simplified method has the same data requirement as strip adjustment procedures, it properly models the impact of the present systematic errors on the LiDAR point cloud. The Simplified method provides a new concept of utilizing estimated discrepancies between conjugate surface elements in overlapping strips for the estimation of biases in the system parameters. In this regard, the proposed procedure can also be used as a quality control tool to evaluate the performance of the data acquisition system, as in [
18]. The simple nature of the utilized data in the Simplified method requires a strict flight configuration and terrain characteristics. More specifically, the Simplified method assumes the availability of parallel flight lines with small pitch and roll angles as well as minor terrain elevation variations compared to the flying height above ground.
The second procedure, denoted as the “Quasi-rigorous method”, is designed to work with time-tagged LiDAR point cloud and navigation data (trajectory position). In contrast to the position and orientation information for each pulse, which is needed for rigorous calibration procedures, the Quasi-rigorous method only requires a sample of the trajectory positions at a much lower rate. The trajectory attitude is not used in this approach since the utilization of such information would require the knowledge of the sequence of rotation angles defining the system attitude, which might not be always available. Another advantage of the Quasi-rigorous method is that it computes approximations of some of the system raw measurements (
i.e., position, heading, and scan mirror angle of the platform at the moment of pulse emission). Such characteristic makes the method ready for implementation with datasets where raw measurements are available. Compared to the Simplified method, the Quasi-rigorous method doesn’t require strict flight configuration since it can deal with non-parallel strips over terrain with significant terrain elevation variation regardless of the flying height. The only limitation of this method is having a system with small pitch and roll angles, which is quite realistic for steady platforms (e.g., fixed-wing aircraft). Although the Quasi-rigorous procedure can make use of control surfaces, the availability of such control is not mandatory for the calibration process. In the Quasi-rigorous method, biases in the system parameters are estimated while reducing discrepancies between conjugate surface elements in overlapping strips or control surfaces, if available. In terms of the nature of the needed data, the Quasi-rigorous method has similar requirement as some of the existing techniques [
15,
16,
17]. However, the proposed method overcomes the existing limitations of such procedures in terms of the modeled biases in the system parameters, the utilized primitives as well as the mechanism of utilizing them. In contrast to [
15], suitable primitives (point-patch pairs) are implemented in the Quasi-rigorous method. When compared to [
16], the unknown parameters do not depend on the extent of the area or the number of primitives being utilized in the Quasi-rigorous calibration procedure. Moreover, a robust matching procedure capable of determining accurate planimetric and vertical discrepancies between conjugate surface elements in overlapping strips is proposed, circumventing the above mentioned limitation in [
17].
In summary, the proposed methods have several advantages when compared with existing techniques such as (i) they demand data that can be easily available to the end-user; (ii) the utilized primitives do not impose any restrictions on the areas where we can conduct the calibration procedure (e.g., the methods do not require the availability of planar patches with varying slopes and aspects, which can only be present in LiDAR data over urban areas); (iii) they don’t require control information; (iv) they are based on suitable primitives (point-patch pairs) that do not require pre-processing of the LiDAR data; (v) they are relatively easy to implement; (vi) they can be implemented in every flight mission if needed; (vii) they can also be utilized as a "quality control" tool (especially the first method) to evaluate the quality of the LiDAR point cloud; and (viii) the proposed procedures are characterized by having a high level of automation (manual interaction is restricted to the selection of local areas, which will be utilized for the estimation of discrepancies between overlapping strips). The procedures can be fully automated through the development of a tool that isolates local areas characterized by facets with different slopes and aspects (e.g., using eigen value analysis).
The paper starts by investigating the appropriate mathematical model relating conjugate surface elements in overlapping strips in the presence of systematic errors in the system parameters. The mathematical derivation is based on a simplified LiDAR equation, which is obtained based on the considered assumptions in each method. Due to the irregular nature of the LiDAR points, appropriate primitives that can be extracted with a satisfactory level of automation (i.e., requiring minimum user interaction) are introduced for both methods. In addition, experimental results using simulated and real data are presented. Finally, the paper presents some conclusions and recommendations for future work.
3. Experimental Results
In this paper, we presented two new calibration procedures for the estimation of systematic errors in the LiDAR acquisition system. Experimental results using simulated data have proven the validity of the presented calibration procedures. More specifically, it has been verified that the proposed calibration procedures lead to an accurate estimate of the system parameters when using simulated LiDAR data that comply with the listed assumptions. Moreover, the simulation experimental results have shown that reasonable deviations from the listed assumptions would not significantly affect the quality of the estimated parameters. More specifically, the Simplified method produced good estimates when dealing with non parallel strips with 20° deviation from parallelism, pitch and roll angles of up to 5°, and variations in the terrain elevation of up to 10% relative to the flight height. The Quasi-rigorous method, on the other hand, can tolerate pitch and roll angles of up to 5°. Due to the space limitation, only the results from real data will be presented in this paper.
To perform the experiments using real data, a LiDAR dataset, which was captured by an Optech ALTM 2050 using the optimum flight configuration, was utilized. The optimum flight configuration, as mathematically demonstrated in
Section 2.1.1, consists of four strips which are captured from two flying heights in opposite directions with 100% overlap, and two flight lines, which are flown in the same direction with the least overlap possible. In addition to strips in the optimum configuration, some extra strips were acquired as well.
Figure 5 shows the characteristics of the acquired dataset and
Table 1 presents the involved overlapping strip pairs. In addition to testing the validity and comparing the performance of the proposed calibration procedures, we would like to investigate whether the calibration results are significantly different when adding more overlapping pairs to the minimum recommended optimum configuration.
Table 2 shows the tested scenarios. Test scenario “I” corresponds to the minimum optimal configuration, consisting of three overlapping pairs. Test scenario “II” adds one more overlapping pair to the minimum configuration, for a total of four overlapping pairs. Finally test scenario “III” adds three overlapping pairs to the minimum configuration, for a total of six overlapping pairs.
Figure 5.
LiDAR dataset configuration.
Figure 5.
LiDAR dataset configuration.
Table 1.
Characteristics of the involved overlapping pairs in the proposed calibration procedures.
Table 1.
Characteristics of the involved overlapping pairs in the proposed calibration procedures.
Overlapping Pairs | % of Overlap | Direction |
(i) Strips 1&2 | 100% | Opposite directions |
(ii) Strips 3&4 | 100% | Opposite directions |
(iii) Strips 3&5 | 50% | Same direction |
(vi) Strips 1&6 | 70% | Opposite directions |
(v) Strips 5&7 | 50% | Opposite directions |
(vi) Strips 7&8 | 40% | Opposite directions |
(vii) Strips 2&6 | 70% | Same direction |
Table 2.
List of the utilized overlapping pairs for the tested scenarios.
Table 2.
List of the utilized overlapping pairs for the tested scenarios.
Test Scenario | Overlapping Pairs |
I | (i), (ii), and (iii) |
II | (i), (ii), (iii), and (vii) |
III | (i), (ii), (iii), (iv), (v), and (vi) |
3.1. Simplified Method Results
As described in
Section 2.1, the Simplified method consists of a two-step procedure. First, the discrepancies between parallel overlapping strips are determined; then, the system biases are estimated based on the determined discrepancies between the strips using the derived equations in
Section 2.1.1.
Note that the flight lines are either in the NE-SW or SW-NE direction. The provided formulas in
Section 2.1.1 assume that the flight lines are parallel to the Y-Axis. Therefore, the estimated transformation parameters using the discrepancy detection procedure need to be recalculated to correspond to a coordinate system where the flight direction is parallel to the Y-Axis. Equations 38 and 39 are utilized to compute the transformation parameters (
and
) in the local coordinate system given the transformation parameters (
and
) in the ground coordinate system and the heading angle
.
Table 3 presents the recalculated transformation parameters (w.r.t. the local coordinate system).
Although the conducted analysis of the impact of the studied biases has proven that the effect on the derived point cloud will consist of three shifts (
XT,
YT,
ZT) and one rotation around the flight direction (
ϕ), we also computed the other two rotation angles in order to check for the presence of other biases, such as unmodeled biases in the system parameters and GPS/INS navigation errors. We can notice in
Table 3 a significant discrepancy in the along flight direction (Y direction) for the strips flown in opposite directions. This discrepancy is believed to be caused by biases in the boresight pitch angle. Also, significant values for the
angle for strips 1&2, 2&6 and 1&6 are evident. This can be explained by deviation from parallelism between overlapping strips.
Table 3.
Estimated transformation parameters w.r.t. the local coordinate system using the discrepancy detection procedure.
Table 3.
Estimated transformation parameters w.r.t. the local coordinate system using the discrepancy detection procedure.
| (m) | (m) | (m) | (sec) | (sec) | (sec) |
---|
Strips 1&2 | −0.25 | 1.27 | −0.01 | 10.54 | 1.34 | 67.68 |
Strips 3&4 | −0.01 | 0.52 | 0.02 | −2.59 | −9.72 | 20.52 |
Strips 3&5 | −0.32 | −0.19 | −0.06 | 7.20 | 96.48 | 6.48 |
Strips 2&6 | −0.42 | 0.15 | 0.01 | −4.65 | 74.61 | −136.10 |
Strips 1&6 | −0.67 | 1.51 | −0.06 | 4.68 | 91.08 | 71.64 |
Strips 5&7 | −0.13 | 0.55 | 0.04 | 0.68 | 12.96 | −5.40 |
Strips 7&8 | 0.48 | 0.82 | 0.06 | 1.62 | −166.32 | −2.84 |
In order to determine the magnitude of the biases in the system parameters, the estimated transformation parameters in
Table 3 are then expressed as a linear combination of the biases in the LiDAR system using Equation 15 for strips flown in opposite directions. It should be noted that for strips which are flown in opposite directions with 100% overlap (
i.e., strips 1&2 and 3&4), Equation 15 would reduce to the form in Equation 16. For strips flown in the same direction, Equation 17 should be used. One example of these linear equations is given in Equation 40 for strips 1&2. Finally, the resulting equations are used in a least-squares adjustment to derive an estimate of the systematic biases in the data acquisition system. The system biases were estimated using three different test scenarios (refer to
Table 2). The results obtained for the different test scenarios are reported in
Table 4. We can observe in this table that, the bias in the boresight pitch angle is the most significant one, followed by the bias in the boresight yaw angle.
Table 4.
Estimated system biases using the Simplified method for the 3 test scenarios in
Table 2.
Table 4.
Estimated system biases using the Simplified method for the 3 test scenarios in Table 2.
Test Scenario | (m) | (m) | (") | (") | (") | (m) | |
I | −0.07 | −0.11 | 75 | −1 | 80 | 0.26 | 0.000565 |
II | −0.08 | −0.11 | 75 | −1 | −39.5 | 0.22 | 0.000567 |
III | −0.11 | −0.10 | 86 | −16 | 41 | 0.28 | 0.000670 |
To test the equivalency of the estimated calibration parameters obtained using the three test scenarios, one can compare the numerical values of the individual parameters. However, such a procedure would not lead to a meaningful quantitative measure of their equivalency or deviation. Similar to [
21], a pair-wise comparison of the estimated calibration parameters is devised in this research. More specifically, the original LiDAR strip is reconstructed using two sets of calibration parameters. The coordinates of the LiDAR points after reconstruction using the two sets of calibration parameters are compared through an RMSE analysis. If the estimated RMSE values are within the predicted noise level in the point cloud coordinates, which is based on the law of error propagation using the expected noise in the system measurements and navigation data (from the system manufacturer specification), the estimated calibration parameters are deemed equivalent. The average noise level in the LiDAR point cloud in the XYZ directions is (0.27m, 0.29m, and 0.12m) from a flying height of 1000m and (0.45m, 0.48m, and 0.18m) from a flying height of 2000m.
Table 5 presents the RMSE analysis for the three test scenarios, using strip 1 (flight height of 2000 m) and strip 3 (flying height of 1,000 m). The values presented are equal or smaller than the noise level, indicating the equivalency of the calibration parameters.
Table 5.
RMSE equivalency analysis for the three test scenarios (Simplified method).
Table 5.
RMSE equivalency analysis for the three test scenarios (Simplified method).
| Strip 1 | Strip 3 |
| I vs. II | I vs. III | II vs. III | I vs. II | I vs. III | II vs. III |
Mean ± Std (RMSE) | X (m) | −0.08 ± 0.15 | 0.13 ± 0.04 | 0.21 ± 0.11 | −0.03 ± 0.10 | 0.06 ± 0.02 | 0.10 ± 0.08 |
(0.17) | (0.13) | (0.24) | (0.10) | (0.07) | (0.12) |
Y (m) | −0.02 ± 0.05 | −0.08 ± 0.05 | −0.07 ± 0.01 | 0.00 ± 0.03 | −0.02 ± 0.03 | −0.02 ± 0.00 |
(0.05) | (0.10) | (0.06) | (0.03) | (0.04) | (0.02) |
Z (m) | 0.04 ± 0.00 | −0.00 ± 0.02 | −0.04 ± 0.02 | 0.04 ± 0.00 | −0.01 ± 0.01 | −0.05 ± 0.01 |
(0.04) | (0.02) | (0.05) | (0.04) | (0.02) | (0.05) |
3.2. Quasi-Rigorous Method Results
As in the Simplified method, the system biases in the acquisition system were estimated using the three test scenarios. The estimated system biases are reported in
Table 6, which are compatible with the estimated system biases using the Simplified method (compare the reported results in
Table 4 and
Table 6) except for the range and the mirror angle scale biases. The reason for such incompatibility is explained later on in this paragraph. As presented in
section 3.1 for the Simplified method results, the equivalency of the estimated calibration parameters obtained using the three test scenarios were evaluated for the Quasi-rigorous method as well. Here again, a pair-wise RMSE comparison of the estimated calibration parameters was performed.
Table 7 presents the RMSE analysis for the three test scenarios, using strip 1 (flight height of 2,000 m) and strip 3 (flying height of 1,000 m). Comparing the estimated RMSE and predicted noise magnitude in the point cloud coordinates, one can conclude the equivalency of the calibration parameters except for the pair-wise analysis involving scenario I, where we can observe RMSE values in the Z direction larger than 60cm. For this test scenario, we can notice the same values for the mean differences, which signify that these RMSE values are due to constant shifts in the Z direction. Correlations between the boresight angles and lever-arm offset components, and between the range and the mirror angle scale biases could be observed in all tested scenarios. The correlation between the range and the mirror angle scale biases were found to be −0.98, −0.87, and −0.91 for scenarios I, II and III, respectively. The correlation between these parameters explains the difference between the estimated values for such biases using the Quasi-rigorous and the Simplified methods (
i.e., observed discrepancies between
Table 4 and
Table 6). The test scenario I presented the highest correlation, explaining the high RMSE values in the equivalency analysis involving scenario I. It will be shown in the results evaluation section that even though we have this constant shift introduced when using test scenario I, the obtained compatibility of the overlapping strips after the calibration procedure will be quite similar for all three test scenarios.
Table 6.
Estimated system biases using the Quasi-rigorous method for the 3 test scenarios.
Table 6.
Estimated system biases using the Quasi-rigorous method for the 3 test scenarios.
Test Scenarios | (m) | (m) | (") | (") | (") | (m) | |
I | −0.07 | −0.11 | 79.6 | −2.6 | 52.5 | 0.59 | 0.00038 |
II | −0.08 | −0.12 | 80.8 | −4.8 | 17.3 | 0.00 | 0.00097 |
III | −0.05 | −0.09 | 84.3 | −1.2 | 32.5 | −0.09 | 0.00103 |
Table 7.
RMSE equivalency analysis for the three configuration cases (Quasi-rigorous method).
Table 7.
RMSE equivalency analysis for the three configuration cases (Quasi-rigorous method).
| Strip 1 | Strip 3 |
| I vs. II | I vs. III | II vs. III | I vs. II | I vs. III | II vs. III |
Mean ± Std (RMSE) | X (m) | 0.01 ± 0.02 | 0.06 ± 0.01 | −0.05 ± 0.03 | 0.00 ± 0.03 | 0.05 ± 0.02 | −0.04 ± 0.01 |
(0.03) | (0.06) | (0.06) | (0.03) | (0.05) | (0.05) |
Y (m) | 0.02 ± 0.13 | 0.03 ± 0.12 | −0.02 ± 0.00 | 0.00 ± 0.01 | −0.00 ± 0.00 | 0.00 ± 0.01 |
(0.13) | (0.13) | (0.02) | (0.01) | (0.00) | (0.01) |
Z (m) | ± 0.03 | ± 0.04 | −0.09 ± 0.01 | ± 0.01 | ± 0.01 | −0.09 ± 0.00 |
() | () | (0.09) | () | () | (0.09) |
3.3. Results Evaluation
The calibration results will be evaluated qualitatively and quantitatively. The qualitative evaluation will be performed by visually comparing the quality of generated intensity images from the original and adjusted (reconstructed using the estimated system biases) point cloud to check any improvements in the clarity and definition of various objects. In addition, profiles will be generated using the original and adjusted point cloud to check for any improvements in the quality of fit between overlapping strips. The quantitative assessment, on the other hand, will be performed by applying the discrepancy detection procedure (described in
Section 2.1.2) to the overlapping strips, after reconstructing them using the estimated system biases (Equation 24). The evaluation is reported in the next subsections.
Qualitative Evaluation: The improvement in the quality of the generated intensity images is illustrated in
Figure 6.
Figure 6a shows the intensity image generated before the calibration procedure while
Figure 6b and
Figure 6c show the intensity images generated after the calibration procedure using the Simplified and the Quasi-rigorous methods, respectively. In both methods, configuration III was utilized. As it can be seen in this figure (refer to the circled areas in
Figure 6), enhancement in the feature definition, such as at airplanes and surface markings (the word “PIPER”), are visually clearer in the generated intensity images after the calibration procedure. This means that, besides the positional accuracy of the LiDAR point cloud, the calibration procedure also improves the data interpretation. We can observe from the visual analysis of the intensity images quite compatible results for the two proposed methods. To illustrate the improvement in the quality of fit between overlapping strips, profiles whose location are shown in
Figure 7 were drawn.
Figure 8 shows profiles along strips 1&2 before and after the calibration procedure. We can observe a significant improvement in the compatibility of these profiles after the calibration, especially in the profiles which are in (or almost in) the flight direction—Y direction, caused by improved estimate of the boresight pitch angle. Profile 4, on the other hand, does not present improvements since it is perpendicular to the flight direction, where no significant discrepancies were detected (refer to
Table 3).
Figure 6.
(a) Intensity image generated before the calibration procedure. (b) Intensity image generated after the calibration procedure using the Simplified method and test scenario III. (c) Intensity image generated after the calibration procedure using the Quasi-rigorous method and test scenario III.
Figure 6.
(a) Intensity image generated before the calibration procedure. (b) Intensity image generated after the calibration procedure using the Simplified method and test scenario III. (c) Intensity image generated after the calibration procedure using the Quasi-rigorous method and test scenario III.
Figure 7.
Profiles location along overlapping strips 1&2.
Figure 7.
Profiles location along overlapping strips 1&2.
Figure 8.
Profiles in overlapping strips 1&2 (refer to
Figure 7 for the location of these profiles) showing the degree of compatibility between the point cloud before and after the calibration procedure.
Figure 8.
Profiles in overlapping strips 1&2 (refer to
Figure 7 for the location of these profiles) showing the degree of compatibility between the point cloud before and after the calibration procedure.
Quantitative Evaluation: Table 8 reports the discrepancies before and after applying the calibration parameters determined using the Simplified and Quasi-rigorous methods, respectively. Also, these tables present the results for the three tested scenarios for some of the overlapping pairs. We can observe in the Quasi-rigorous method results that even though we got a high RMSE value in the Z direction (constant shift) when comparing the test scenario I with the other two tested scenarios (refer to the highlighted values in
Table 7), the obtained compatibility of the overlapping strips after the calibration procedure is quite similar for all three tested scenarios. Overall, we can observe quite compatible results between the two proposed methods, and between the three tested scenarios. By comparing the discrepancies after applying the calibration parameters and the discrepancies before the calibration, we can see an improvement in the quality of fit between the strips, especially in Y direction (flight direction) for the strips flown in opposite directions (e.g., strips 1&2 and strips 1&6). The expected outcome would be the absence of discrepancies between the strips after applying the calibration parameters. Remaining discrepancies after point cloud reconstruction might be attributed to errors in the GPS/INS navigation data. In this regard, it is important to emphasize that this work is focusing on the implementation of system calibration procedures. More specifically, this work is focusing on detecting biases in the system parameters and reducing their impact on the derived point cloud. This does not mean that the adjusted point cloud will be error-free. Other errors, such navigation errors, will still impact the quality of the derived point cloud as well as the estimated transformation parameters between adjusted overlapping strips. These errors cannot be included in a system calibration procedure since they are mission or strip-dependent (
i.e., depend on the mission time, GPS constellation, distance to GPS base station,
etc.).
Table 8.
Discrepancies between some of the overlapping strip pairs before and after applying the calibration parameters obtained using real data, the Simplified and Quasi-rigorous methods, and the 3 test scenarios.
Table 8.
Discrepancies between some of the overlapping strip pairs before and after applying the calibration parameters obtained using real data, the Simplified and Quasi-rigorous methods, and the 3 test scenarios.
Before Calibration | After Calibration Simplified Method | After Calibration Quasi-rigorous Method |
---|
| I | II | III | I | II | III |
---|
Strips 1&2 |
XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) |
-0.25 | 1.27 | -0.01 | 0.09 | 0.03 | 0.01 | 0.17 | -0.05 | 0.01 | 0.40 | 0.16 | 0.03 | 0.20 | 0.04 | 0.01 | 0.21 | 0.06 | 0.02 | 0.21 | 0.18 | 0.01 |
ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) |
10.54 | 1.34 | 67.68 | -5.76 | -2.38 | -69.48 | -5.96 | -3.50 | -68.40 | -6.48 | -31.68 | -69.48 | -6.48 | -6.48 | -57.24 | -5.82 | -8.18 | -71.28 | -5.76 | -0.90 | -70.56 |
Strips 1&6 |
XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) |
-0.67 | 1.51 | -0.06 | 0.25 | -0.31 | 0.03 | 0.51 | -0.15 | 0.03 | 0.52 | -0.10 | 0.05 | 0.38 | -0.25 | 0.04 | 0.28 | -0.17 | 0.02 | 0.25 | -0.07 | 0.01 |
ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) |
4.68 | 91.08 | 71.64 | -1.84 | -37.08 | -73.8 | -39.01 | -1.79 | -73.08 | -1.69 | -57.96 | -72.72 | -1.76 | -48.96 | -73.80 | -1.62 | -16.34 | -73.44 | -1.73 | -5.40 | -76.68 |
Strips 3&5 |
XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) | XT (m) | YT (m) | ZT (m) |
-0.32 | -0.19 | -0.06 | 0.38 | 0.11 | 0.05 | -0.08 | 0.18 | 0.06 | -0.07 | 0.04 | 0.02 | -0.08 | 0.02 | 0.05 | -0.08 | 0.09 | 0.05 | -0.07 | 0.06 | 0.05 |
ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) | ω (sec) | ϕ (sec) | κ (sec) |
7.20 | 96.48 | 6.48 | 2.84 | -45.00 | -9.36 | -5.63 | 17.56 | -5.76 | -5.04 | 40.32 | -6.12 | -5.76 | 15.12 | -7.56 | -5.15 | 68.25 | -6.84 | -5.04 | 71.64 | -6.48 |
4. Conclusions and Recommendations for Future Work
In this paper, two new calibration procedures for the estimation of biases in the LiDAR system parameters are introduced. The key the advantages of the proposed methods are the nature of the required data, the utilized primitives as well as the mechanism of using such primitives for the estimation of residual biases in the system parameters.
The first proposed calibration procedure, denoted as the Simplified method, was designed to satisfy the needs of the users that have only access to the LiDAR point cloud. The loose data requirement mandates a strict flight configuration and terrain characteristics for implementing such a method. More specifically, it demands parallel overlapping LiDAR strips acquired by a linear scanner using a steady platform over an area with moderately varying terrain elevation. This method consists of a two-step procedure: first, the discrepancies between parallel overlapping strips are determined; then, the system biases are estimated using the estimated discrepancies between the strips. In this regard, the proposed procedure can be used as a quality control tool to evaluate the performance of the data acquisition system (i.e., inaccurate system calibration or navigation errors would be manifested in the form of significant discrepancies between overlapping strips). Different from existing strip adjustment procedures, which rely on the LiDAR point cloud only, the Simplified method properly models the impact of systematic errors in the system parameters on the point cloud coordinates.
The second presented calibration procedure is denoted as the Quasi-rigorous method due the fact that few reasonable assumptions are undertaken. This method only assumes that we are dealing with a linear scanner and that the LiDAR system unit is almost vertical, which is quite realistic for flight missions using a steady platform. This method requires time-tagged point cloud and trajectory position data. In contrast to the position and orientation information requirement for each pulse in the rigorous calibration, the Quasi-rigorous procedure only requires a sample of the trajectory positions at a much lower rate. Access to this type of data is not a concern. Moreover, since this calibration procedure derives approximations of some of the system raw measurements (i.e., position, heading, and scan mirror angle of the platform at the moment of pulse emission), the proposed procedure is ready for implementation with datasets where raw measurements are available. In contrast to some existing procedures, the unknown parameters do not change with the extent of the area or the number of primitives being utilized in the calibration procedure. Moreover, control information is not essential for the system calibration. Different from the Simplified method, the Quasi-rigorous method can deal with non-parallel strips and can handle heading variations and varying terrain elevation. Similar to the Simplified method, the system parameters are estimated while reducing discrepancies between conjugate surface elements in overlapping strips.
The Simplified and Quasi-rigorous procedures are both based on appropriate primitives (point-patch pairs) that do not require pre-processing of the LiDAR data. In contrast to some of the existing procedures, which rely on the availability of large planar patches with varying slopes and aspects and restrict the calibration to datasets covering urban areas, the utilized primitives can be present in any type of terrain coverage. Moreover, the presented methods are characterized by having a high level of automation. Manual interaction is restricted to the selection of local areas, which are utilized for the estimation of discrepancies between overlapping strips. The procedures can be fully automated through the development of a tool that isolates local areas characterized by facets with different slopes and aspects (e.g., using eigen value analysis). The proposed methods are easy to implement and can be used in every flight mission if needed. Finally, the developed procedures provide a detailed analysis of the optimum flight configuration for reliable estimation of residual biases in the system parameters.
The performance of the developed calibration procedures has been verified using experimental results from real datasets. The results have shown the feasibility of the proposed calibration procedures in operational environments. It was shown that collected LiDAR data might exhibit significant incompatibilities due to insufficient calibration procedures. The calibration procedures were evaluated qualitatively and quantitatively. The qualitative evaluation of the calibration procedures was performed by visually comparing the quality of intensity images generated from the original and adjusted (reconstructed using the estimated system biases) point cloud to check for any improvements in the clarity and definition of various objects. In addition, profiles were generated using the point cloud before and after the calibration procedure to check any improvements in the quality of fit between overlapping strips. The quantitative evaluation was performed by computing the discrepancies between overlapping strips after reconstructing the LiDAR point cloud using the estimated biases in the system parameters. The qualitative assessment has demonstrated a significant improvement in visual appearance of the generated intensity images from multiple strips, as well as the quality of fit between overlapping strips (visually checked in drawn profiles) before and after the calibration procedure. Although the quantitative evaluation has shown considerable improvement in the quality of fit between the strips (e.g., the discrepancy along the flight direction between strips 1&2 and strips 1&6 was reduced more than 1.20m and 1.40m, respectively, after the calibration procedure), it could be observed that there are still remaining discrepancies between the strips after the calibration procedure. It is important to clarify that this work has focused on the implementation of system calibration procedures. More specifically, this work focused on detecting biases in the system parameters and reducing their impact of on the derived point cloud. This does not mean that the adjusted point cloud will be error-free. Other errors, such navigation errors, will still impact the quality of the derived point cloud as well as the estimated transformation parameters between adjusted overlapping strips. These errors cannot be included in a system calibration procedure since they are mission-dependent. Finally, the experimental results have shown that the estimated biases in the system parameters do not change by using additional flight lines beyond what has been established as the optimal flight configuration.
In this paper, only the relative accuracy of the adjusted LiDAR point cloud has been demonstrated. Current work is dealing with verifying the impact of the calibration procedures on the absolute accuracy of the LiDAR data. Future work will focus on diagnosing the source of the remaining discrepancies and extending the error model if needed. In addition, effort will be given towards the operability of the calibration procedures by automating the selection of the areas in the overlapping region as well as by improving the computational speed of the proposed methods (e.g., parallel processing). Moreover, the incorporation of control information in the calibration method will be investigated.