Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Using Spatial Structure Analysis of Hyperspectral Imaging Data and Fourier Transformed Infrared Analysis to Determine Bioactivity of Surface Pesticide Treatment
Next Article in Special Issue
Monitoring Automotive Particulate Matter Emissions with LiDAR: A Review
Previous Article in Journal
Digital Northern Great Plains: A Web-Based System Delivering Near Real Time Remote Sensing Data for Precision Agriculture
Previous Article in Special Issue
Ground Filtering Algorithms for Airborne LiDAR Data: A Review of Critical Issues
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Alternative Methodologies for LiDAR System Calibration

Department of Geomatics Engineering, The University of Calgary, 2500 University Drive NW, T2N 1N4, Calgary, AB, Canada
*
Author to whom correspondence should be addressed.
Remote Sens. 2010, 2(3), 874-907; https://doi.org/10.3390/rs2030874
Submission received: 26 January 2010 / Revised: 8 March 2010 / Accepted: 16 March 2010 / Published: 23 March 2010
(This article belongs to the Special Issue LiDAR)

Abstract

:
Over the last few years, LiDAR has become a popular technology for the direct acquisition of topographic information. In spite of the increasing utilization of this technology in several applications, its accuracy potential has not been fully explored. Most of current LiDAR calibration techniques are based on empirical and proprietary procedures that demand the system’s raw measurements, which may not be always available to the end-user. As a result, we can still observe systematic discrepancies between conjugate surface elements in overlapping LiDAR strips. In this paper, two alternative calibration procedures that overcome the existing limitations are introduced. The first procedure, denoted as “Simplified method”, makes use of the LiDAR point cloud from parallel LiDAR strips acquired by a steady platform (e.g., fixed wing aircraft) over an area with moderately varying elevation. The second procedure, denoted as “Quasi-rigorous method”, can deal with non-parallel strips, but requires time-tagged LiDAR point cloud and navigation data (trajectory position only) acquired by a steady platform. With the widespread adoption of LAS format and easy access to trajectory information, this data requirement is not a problem. The proposed methods can be applied in any type of terrain coverage without the need for control surfaces and are relatively easy to implement. Therefore, they can be used in every flight mission if needed. Besides, the proposed procedures require minimal interaction from the user, which can be completely eliminated after minor extension of the suggested procedure.

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, X G , can be derived through the summation of three vectors ( X o , P G and ρ ) after applying the appropriate rotations: Ryaw, pitch, roll, RΔωφΔκ and Rα, β. In this equation, X o is the vector from the origin of the ground reference frame to the origin of the IMU coordinate system, P G 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 X o 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.
Remotesensing 02 00874 g001
X G = X o + R y a w ,   p i t c h ,   r o l l   P G + R y a w ,   p i t c h ,   r o l l R Δ ω , Δ φ , Δ κ R S α α , S β β [ 0 0 ρ ]
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.

2. Proposed Calibration Methods

In this section, two new calibration procedures for the determination of the biases in the system parameters are presented. The first introduced method, denoted as “Simplified method”, utilizes only the LiDAR point cloud and deals with parallel strips. The second proposed method, denoted as “Quasi-rigorous method” can deal with both parallel and non-parallel strips. However, this method requires time-tagged LiDAR point cloud and the trajectory position data. In the experimental results section, a comparative analysis between the two proposed calibration procedures will be presented.

2.1. Simplified Method

In the Simplified method, the system biases are estimated using the detected discrepancies between overlapping LiDAR strips. More specifically, this calibration method consists of a two-step procedure: first, the discrepancies between parallel overlapping strips are determined; then, the system biases are estimated using the detected discrepancies between the strips.
The proposed method will be explained in the following subsections. First, the mathematical model relating the discrepancies between parallel overlapping strips and the system biases will be derived. Such analysis will also lead to the optimum flight configuration for reliable estimation of biases in the system parameters. The mathematical derivation will utilize a simplified LiDAR equation, which is based on some reasonable assumptions. Finally, a methodology for detecting the discrepancies and evaluating the transformation parameters relating overlapping strips will be presented.

2.1.1. Mathematical model

In the Simplified method, the following assumptions are considered in the mathematical derivation: (a) we are dealing with a linear scanner with the mirror scanning in the across-flight direction, (b) flying directions are parallel to the positive and negative directions of the Y-axis of the ground coordinate system, (c) the flight lines follow a straight-line trajectory with constant attitude, (d) the LiDAR system is almost vertical (i.e., Ryaw, pitch, roll ≈ Identity matrix for a system flying along the positive direction of the Y-axis), (e) the LiDAR system has relatively small boresight angles, (f) the mapped area is comprised of a relatively flat terrain, where the height variations are much smaller than the flying height above ground, and (g) the Y-axis of the ground coordinate system is defined half-way between the overlapping strips at the ground level. The convention used for the laser scanner and IMU body frame coordinate systems is right-forward-up (right-handed). Such assumptions simplify the LiDAR geometric model as represented by Equation 1 to the form in Equations 2. Note that the scan angle (β) and the lateral distance (x) for a given point are defined relative to the flight trajectory. Therefore, one should use the appropriate signs when dealing with two flight lines, which are flown in opposite directions. Based on the above assumptions, one can conclude that the Simplified method requires the availability of parallel flight lines captured by a steady platform over an area with moderately varying elevation. If the strips are not flown along the Y-axis direction, then the dataset just need to be rotated in order to use the equations that will be derived in this section. An alternative way, which is implemented in this paper, would be to rotate/recalculate the detected discrepancies between overlapping strips to correspond to a coordinate system where the flight direction is parallel to the Y-axis. The implementation details are provided in the experimental results section as well as the impact of deviations from the above listed assumptions.
X G X o + [ Δ X Δ Y Δ Z ] + [ 1 Δ κ Δ φ Δ κ 1 Δ ω Δ φ Δ ω 1 ] [ ( ρ + Δ ρ )   sin ( S   β ) 0 ( ρ + Δ ρ )   cos ( S   β ) ] = X o + [ Δ X Δ Y Δ Z ] + [ 1 Δ κ Δ φ Δ κ 1 Δ ω Δ φ Δ ω 1 ]   [ x 0 H ]
X G X o + [ Δ X Δ Y Δ Z ] + [ 1 Δ κ Δ φ Δ κ 1 Δ ω Δ φ Δ ω 1 ]   [ ( ρ + Δ ρ )   sin ( S   β ) 0 ( ρ + Δ ρ )   cos ( S   β ) ] = X o + [ Δ X Δ Y Δ Z ] + [ 1 Δ κ Δ φ Δ κ 1 Δ ω Δ φ Δ ω 1 ]   [ x 0 H ]
where,
  • Equation 2a is valid for a system flying along the positive direction of the Y-axis (this flight line will be denoted as the forward strip),
  • Equation 2b is valid for a system flying along the negative direction of the Y-axis (this flight line will be denoted as the backward strip),
  • X, ∆Y, ∆Z are the components of the lever-arm offset vector P G ,
  • S is the scale factor for the mirror angle β (this scale factor should be unity for a bias-free system),
  • H is the flying height above ground, and
  • x is the lateral distance between the LiDAR point in question and the projection of the flight trajectory onto the ground.
The LiDAR point coordinates ( X G ) , as presented in Equations 2a,b, are function of the system parameters ( x ) and measurements ( l ) (Equation 3) and represent the true point coordinates ( X T r u e ) . In the presence of biases in the system parameters, the LiDAR point coordinates will become biased ( X B i a s e d ) and will be function of the system parameters and measurements as well as the biases in the system parameters ( δ x ) , as expressed by Equation 4. In this work, as shown in Equation 4, we will investigate the impact of biases in the lever-arm offset components ( δ Δ X , δ Δ Y , δ Δ Z ), biases in the boresight angles ( δ Δ ω , δ Δ φ , δ Δ κ ), constant bias in the measured ranges ( δ Δ ρ ) , and constant scale bias in the scan mirror angles ( δ S ) . Equation 4 can be linearized with respect to the system parameters using Taylor series expansion, yielding the form in Equations 5 and 6, after ignoring second and higher order terms. The term f x represents the partial derivatives with respect to the system parameters, while the term f x   δ x represents the impact of the system biases onto the derived point cloud ( δ X G ) .
X G = X T r u e = f ( x , l )
where,
  • x = ( Δ X ,   Δ Y ,   Δ Z ,   Δ ω ,   Δ φ ,   Δ κ ,   Δ ρ ,   S ) , and
  • l = ( X o ,   y a w ,   p i t c h ,   r o l l ,   β ,   ρ )
X B i a s e d = f ( x + δ x ,   l )
where,
δ x = ( δ Δ X ,   δ Δ Y ,   δ Δ Z ,   δ Δ ω ,   δ Δ φ ,   δ Δ κ ,   δ Δ ρ ,   δ S )
X B i a s e d f ( x ,   l ) + f x   δ x = X T r u e + f x   δ x = X T r u e + [ δ X G δ Y G δ Z G ]
X B i a s e d X T r u e + [ δ X G δ Y G δ Z G ] δ Δ X ,   δ Δ Y ,   δ Δ Z + [ δ X G δ Y G δ Z G ] δ Δ ω ,   δ Δ φ ,   δ Δ κ + [ δ X G δ Y G δ Z G ] δ Δ ρ + [ δ X G δ Y G δ Z G ]   δ S = X T r u e + [ δ X G δ Y G δ Z G ] T o t a l
By differentiating Equations 2 with respect to the system parameters, Equation 6 would lead to the form in Equation 7. The multiple signs in this equation indicate the impact for the forward and backward strips (with the top sign referring to the forward strip). So far, we discussed the impact of various biases in the LiDAR system parameters on the derived point cloud. Using the derived expression, we would like to derive the mathematical relationship between conjugate points in overlapping strips, which are flown in the same or opposite directions.
[ δ X G δ Y G δ Z G ] T o t a l = [ X B i a s e d X T r u e Y B i a s e d Y T r u e Z B i a s e d Z T r u e ] [ ± δ Δ X H   δ Δ φ sin ( β )   δ Δ ρ H   β   δ S ± δ Δ Y   ± H   δ Δ ω ±   x   δ Δ κ δ Δ Z   x   δ Δ φ   cos ( β )   δ Δ ρ x   β   δ S ]
Due to the presence of various systematic errors, the bias-contaminated coordinates of conjugate points in overlapping strips will show systematic discrepancies. The mathematical relationship between these points can be derived by rewriting Equation 7 for the two overlapping strips and subtracting the resulting equations from each other. An example of such a relationship for two strips, which are flown in opposite directions, is shown in Equations 8, 9, and 10. It should be noted that Equation 8 refers to the forward strip (denoted by the subscript A) while Equation 9 refers to the backward strip (denoted by the subscript B).
[ X A X True Y A Y True z A z True ] [ δ Δ X H δ Δ φ sin ( β A )   δ Δ ρ H β A   δ S δ Δ Y + H δ Δ ω + x A   δ Δ κ δ Δ Z x A   δ Δ φ   cos ( β A )   δ Δ ρ x A   β A   δ S ]
[ X B X True Y B Y True z B z True ] [ δ Δ X + H δ Δ φ + sin ( β B )   δ Δ ρ + H β B   δ S δ Δ Y H δ Δ ω x B   δ Δ κ δ Δ Z x B   δ Δ φ   cos ( β B )   δ Δ ρ x B   β B   δ S ]
[ X A X B Y A Y B z A z B ] [ 2 δ Δ X   2 H δ Δ φ [ sin ( β A ) + sin ( β B ) ]   δ Δ ρ   H ( β A + β B )   δ S 2 δ Δ Y +   2 H δ Δ ω + ( x A + x B )   δ Δ κ ( x A x B )   δ Δ φ   [ cos ( β A ) cos ( β B ) ]   δ Δ ρ ( x A   β A x B   β B )   δ S ]
Equation 10 can be simplified by assuming that in the scan angle range of [±25°], which is the typical scan angle range for operational LiDAR systems, the sine and tangent functions are almost equivalent. Therefore, the term [ sin ( β A ) + sin ( β B ) ] δ Δ ρ in Equation 10 can be simplified to [ tan ( β A ) + tan ( β B ) ] δ Δ ρ = [ x A / H + ( x B ) / H ] δ Δ ρ = δ Δ ρ [ D / H ] , since ( x A + x B ) = ( D ) = D , where D refers to the lateral distance between the two flight lines in question (see Figure 2a). Also, one can assume that the cosine function in this range is almost constant. As a result, the term [ cos ( β A ) cos ( β B ) ]   δ Δ ρ can be reduced to zero. Finally, assuming that β A x A / H and β B x B / H , the term ( x A   β A x B   β B )   δ S can be simplified to ( x A 2 x B 2 ) H   δ S = ( x A x B ) β T δ S , where ( β T D / H ) refers to the total scan angle between the two flight lines (i.e., the scan angle from one flight line to an object point, which is vertically below the second flight line, refer to Figure 2a). Such simplifications would lead to Equation 11.
[ X A X B Y A Y B Z A Z B ] [ 2   δ Δ X 2   H   δ Δ φ D / H   δ Δ ρ H   β T   δ S 2   δ Δ Y + 2   H   δ Δ ω D   δ Δ κ ( x A x B )   δ Δ φ ( x A x B )   β T   δ S ] = [ 2   δ Δ X 2   H   δ Δ φ D / H   δ Δ ρ H   δ θ 2   δ Δ Y + 2   H   δ Δ ω D   δ Δ κ ( x A x B )   δ Δ φ ( x A x B )   δ θ ]
where, δθ = βTδS.
Equation 11 can be rewritten after mathematical manipulation to produce the form in Equation 12, where the coordinates of a given point in the forward strip is expressed in terms of the coordinates of the corresponding point in the backward strip. As it can be seen in Figure 2, the coordinates of the involved points are referred to a local coordinate system whose Y-axis is half-way between the two strips. The mathematical manipulation leading to the transition from Equation 11 to Equation 12 is graphically illustrated in Figure 2b. As it can be seen in this figure, the height discrepancy between conjugate points in overlapping strips (last row in Equation 11) is the result of a rotation angle around the flight direction (i.e., roll angle equivalent to 2δφ + 2δθ).
[ X A Y A Z A ] = [ 2   δ Δ X 2   H   δ Δ φ D / H   δ Δ ρ H   δ θ 2   δ Δ Y + 2   H   δ Δ ω D   δ Δ κ 0 ] + R ( 2 δ Δ φ + 2 δ θ ) [ X B Y B Z B ]
Figure 2. (a) Observed object point in overlapping strips flown in opposite directions. (b) The introduced tilt across the flight direction by the boresight roll bias and the scan angle scale bias.
Figure 2. (a) Observed object point in overlapping strips flown in opposite directions. (b) The introduced tilt across the flight direction by the boresight roll bias and the scan angle scale bias.
Remotesensing 02 00874 g002
For two flight lines, which are flown in opposite directions with 100% overlap, the expression in Equation 12 would reduce to the form in Equation 13 (for such flight lines, D and βT are equal to zero).
[ X A Y A Z A ] = [ 2   δ Δ X 2   H   δ Δ φ 2   δ Δ Y + 2   H   δ Δ ω 0 ] + R 2 δ Δ φ   [ X B Y B Z B ]
In a similar fashion, one can derive Equation 14, which expresses the mathematical relationship between the bias-contaminated point coordinates of the same object point in two overlapping strips that are flown in the same direction.
[ X A Y A Z A ] = [ D / H   δ Δ ρ H   δ θ   D   δ Δ κ D   δ Δ φ ] + R 2   δ θ   [ X B Y B Z B ]
Equations 12, 13, and 14 reveal the possibility of identifying the presence of systematic errors in the system parameters by evaluating the discrepancies between conjugate points in overlapping strips, which are flown in the same or opposite directions. Moreover, using these equations, it is possible to determine the flight configuration that maximizes the impact of systematic errors. For example, as it can be seen in Equation 14, parallel strips with the least amount of necessary overlap for identifying conjugate surface elements are useful for magnifying the boresight yaw and roll biases as well as biases in the range and mirror angle scale (this is caused by the large lateral distance D and total scan angle βT). In addition, higher flying heights are optimal for magnifying the boresight pitch and roll biases (Equations 12 and 13). Also, closer inspection of these equations would allow for the determination of an optimal flight configuration design, which decouples various systematic errors. For example, working with four strips which are captured from two flying heights in opposite directions with 100% overlap are optimal for the recovery of the planimetric lever-arm offsets as well as the boresight pitch and roll biases (Equation 13). In addition, two flight lines, which are flown in the same direction with the least overlap possible, are optimal for the recovery of the boresight yaw and roll biases, range bias, and mirror angle scale bias (Equation 14). Only a vertical bias in the lever-arm offset components cannot be detected by observing discrepancies between conjugate surface elements in overlapping strips. Such inability is caused by the fact that a vertical bias in the lever-arm offset components produces the same effect regardless of the flying direction, flying height, or scan angle. The inability of estimating the vertical bias in the lever arm components is not critical given the fact that the vertical lever-arm offset component can be determined with a very high accuracy by field survey. This has been confirmed by the fact that the quality of LiDAR data in the vertical direction is known to be much better when compared to the quality in the XY directions.
In summary, the discrepancies between parallel strips in the presence of the studied biases can be modeled by 3 shifts (XT, YT, ZT) and a rotation angle around the flight line (ϕ). The relationship between the system biases and the discrepancies (XT, YT, ZT, ϕ) between conjugate bias-contaminated points in two flight lines, which are flown in opposite directions is given by Equation 15. It should be noted that this equation would reduce to the form in Equation 16 for flight lines with 100% overlap. For two flight lines flown in the same direction, the relationship between the system biases and the discrepancies between the strips is expressed by Equation 17. The multiple signs ( , ± ) in these equations depend on the relationship between forward and backward strips, with the top sign used when the forward strip is to the right of the backward strip.
[ X T Y T ϕ ] = [ 2   δ Δ X 2   H   δ Δ φ D / H   δ Δ ρ H   δ θ 2   δ Δ Y + 2   H   δ Δ ω D   δ Δ κ 2 δ Δ φ ± 2 δ θ ]
[ X T Y T ϕ ] = [ 2   δ Δ X 2   H   δ Δ φ 2   δ Δ Y + 2   H   δ Δ ω 2 δ Δ φ ]
[ X T Y T Z T ϕ ] = [ D / H   δ Δ ρ   H   δ θ D   δ Δ κ ± D   δ Δ φ ± 2 δ θ ]
The derived mathematical relationship between the biases in the system parameters and the discrepancies among overlapping strips assumes that we can identify conjugate points in overlapping strips, which is not the case for LiDAR surfaces. A discrepancy detection procedure, which can deal with LiDAR surfaces, will be presented in the next section.

2.1.2. Discrepancy detection methodology

In order to come up with a methodology for the determination of the discrepancies between overlapping strips, one must address the following questions:
  • What are the appropriate primitives, which can be used to identify conjugate surface elements in overlapping strips comprised of irregular sets of non-conjugate points?
  • What is the possibility of automatic derivation of these primitives?
  • What is the possibility of automated identification of conjugate primitives in overlapping strips?
  • What is the appropriate similarity measure, which utilizes the involved primitives and the defined transformation function to describe the correspondence of conjugate primitives in overlapping strips?
The following subsection presents the answers to the above questions as they are related to the selected primitives.

2.1.2.1. Primitives extraction

As it has been mentioned earlier, one cannot assume point-to-point correspondence in overlapping strips due to the irregular and sparse nature of the LiDAR points. Instead of points, one can use linear and areal features, which can be identified in overlapping strips. Such primitives, however, would require preprocessing of the LiDAR point cloud to extract areal and linear features (e.g., segmentation, plane fitting, and neighboring plane intersection). Moreover, linear and areal features can be reliably extracted only in urban areas. Such a restriction would limit the practicality of the calibration procedure. In this research, we aim at selecting primitives, which can be derived with minimal preprocessing of the original LiDAR points. Moreover, the selected primitives should be reliably derived in any type of environment (e.g., urban and rural areas). To satisfy these objectives, we chose to represent one strip using the original points, while the other strip is represented by triangular patches, which can be derived from a Triangular Irregular Network (TIN) generation procedure. TIN generation from a given set of points is an automated procedure, which is available in several software packages. As an example, Figure 3 illustrates the case where the strip denoted by (S2) is represented by a set of points while the other strip denoted by (S1) is represented by a set of triangular patches. Due to the high density of the LiDAR data as well as the relatively smooth characteristics of terrain and man-made structures, using TIN patches to describe the physical surface is quite acceptable. It is quite obvious that there are some exceptions where the TIN patches would not represent the physical surface. These exceptions are caused by the irregular and sparse distribution of the LiDAR points. In vegetation areas, for example, planar patches through the LiDAR points do not represent the vegetation physical surface. Another instance of such a deviation will take place at building boundaries where the TIN patches are connecting points along the building rooftop and neighboring ground. In such a case, the TIN patches do not represent the building walls since the LiDAR point cloud would not necessarily coincide with the building boundaries and its footprint along the ground. Figure 4 illustrates these exceptions. In summary, although one cannot assume that there is point-to-point correspondence in overlapping strips, one can argue that we have point-to-patch correspondence as long as the patch represents the physical surface.
Figure 3. Conceptual basis of the proposed point-to-patch correspondence procedure.
Figure 3. Conceptual basis of the proposed point-to-patch correspondence procedure.
Remotesensing 02 00874 g003
The discrepancy detection procedure can be carried out using the entire overlap area between neighboring strips or it can be conducted using an aggregated set of local areas to speed up the process. In this research, an aggregated set of local areas will be utilized to speed up the process. Regardless of utilizing the whole overlap region or the aggregated set of local areas, the correspondence between conjugate point-patch pairs should be established. The identification of the point-patch pairs as well as the estimation of the transformation parameters are explained in the following subsection.
Figure 4. Exceptions where the TIN patches will not represent the physical surface (highlighted in yellow).
Figure 4. Exceptions where the TIN patches will not represent the physical surface (highlighted in yellow).
Remotesensing 02 00874 g004

2.1.2.2. Similarity measure and primitives matching

This section explains the proposed procedure for establishing the correspondence between conjugate point-patch pairs and their utilization for the estimation of the parameters of the transformation function. As shown in Figure 3, if the point qi in S2 belongs to the triangular patch Sp represented by the vertices S p a , S p b , and S p c in S1, then this point should coincide with that patch after applying the appropriate transformation function, which is assumed to be defined by three shifts and three rotation angles (although the considered systematic errors in this research would lead to three shifts and one rotation across the flight direction, we can consider a general transformation involving three rotation angles – the reason for such a generalization will be discussed in the experimental results section). In other words, the volume of the pyramid whose vertices are qi, S p a , S p b , and S p c should be zero. Such a volume can be mathematically described by Equation 18. Using multiple point-patch pairs, we can estimate the transformation parameters which satisfy the volume constraints. The solution to the defined constraints after the linearization shown in Equation 19 can be derived through Equation 20. For reliable estimation of the transformation parameters, the utilized patches should have a balanced distribution of slopes and aspects.
| X q i Y q i Z q i 1 X p a Y p a Z p a 1 X p b Y p b Z p b 1 X p c Y p c Z p c 1 | = 0
where,
[ X q i X q i X q i ] = [ X T Y T Z T   ] + R   [ X q i Y q i Z q i ]
w = A   x + B   e     e   ~   ( 0 ,   Σ Y )
where,
-
x are the corrections to the approximate values of the unknown parameters (XT, YT, ZT, ω, ϕ, κ),
-
A is the partial derivative matrix w.r.t. the unknown parameters,
-
B is the partial derivative matrix w.r.t. the points’ coordinates,
-
w is the estimated determinant using the approximate values for the unknown parameters and points’ coordinates, and
-
ΣY is the a-priori variance-covariance matrix of the points’ coordinates.
x ^ = [ A T ( B   Σ Y   B T ) 1   A ] 1   A T ( B   Σ Y   B T ) 1   w B   e ˜ = w A   x ^ σ ^ o 2 = e ˜ T   B T   ( B   Σ Y   B T ) 1   B   e ˜ / r e d u n d a n c y   a posteriori variance factor Σ ˜ Y = σ ^ o 2   Σ Y    a-posteriori variance-covariance matrix of the point cloud coordinates
After introducing the necessary similarity measure, we need to propose a matching methodology for establishing the correspondences between points in S2 and patches in S1. The approach proposed by [19] deals with the point-to-patch matching problem in the presence of significant rotations and shifts between two overlapping surfaces. For such a case, the Modified Iterated Hough Transform (MIHT) is used to sequentially estimate the transformation parameters through a voting scheme in an accumulator array. For the current work, the correspondence is performed in an iterated manner, using the Iterative Closest Patch (ICPatch) procedure, which corresponds to the Iterative Closest Point (ICP) approach [20] after some modification. The initial correspondence between qi and Sp can be established by determining the patch Sp with the shortest normal distance to the point qi, which corresponds to qi after applying an initial approximation of the transformation parameters. Since we are dealing with small biases in the system parameters, which will lead to discrepancies between overlapping strips in the order of few decimeters to few meters in the worst cases, we can consider zero shifts and rotations as the initial parameters. To be considered a correct match, the shortest normal distance should be less than a given threshold to avoid situations where the triangular patch does not represent the physical surface (e.g., in vegetation and building boundaries as shown in Figure 4). Moreover, the projection of qi onto the patch Sp should be located inside the patch. The initial correspondences can then be used to derive an updated estimate of the transformation parameters, which can be used to derive a new set of correspondences. The new correspondences can be utilized to derive a better estimate of the transformation parameters. Such a procedure is repeated until convergence, where there is no significant change in the estimated transformation parameters from one iteration to the next.
Once the transformation parameters relating parallel overlapping strips are determined, the biases in the mounting parameters can be estimated using Equations 15 and 17 for strips flown in opposite directions and in the same direction, respectively. In these equations, the transformation parameters are expressed as a linear combination of the biases in the LiDAR system parameters. The resulting linear system can be then solved using a least-squares adjustment procedure to derive an estimate of the systematic biases in the data acquisition system.

2.2. Quasi-Rigorous Method

In contrast to the Simplified calibration procedure, the Quasi-rigorous calibration procedure does not have restrictions in terms of the flight configuration and terrain elevations. In other words, this calibration procedure can deal with non-parallel strips over terrain with significant elevation variation regardless of the flying height. This method only assumes that we are dealing with a linear scanner and that the LiDAR system unit is almost vertical. The biases in the system mounting parameters are estimated in a single-step procedure while reducing discrepancies between conjugate surface elements in overlapping strips.
This section starts by deriving the mathematical model relating conjugate points in overlapping strips in the presence of systematic biases in the system parameters. Then, the similarity measure, which incorporates the primitives together with the established mathematical relationship to describe their correspondence, is introduced.

2.2.1. Mathematical model

The following assumptions are considered in the mathematical derivation of the Quasi-rigorous method: (a) we are dealing with a linear scanner, (b) the LiDAR system is almost vertical (i.e., pitch and roll angles are almost zero), and (c) the LiDAR system has relatively small boresight angles. Such assumptions simplify the LiDAR geometric model as represented by Equation 1 to the form in Equation 21. The convention for the laser scanner and IMU body frame coordinate systems is the same as the one used in Simplified method, right-forward-up (right-handed).
X G X o + [ cos κ sin κ 0 sin κ cos κ 0 0 0 1 ]   [ Δ X Δ Y Δ Z ] + [ cos κ sin κ 0 sin κ cos κ 0 0 0 1 ]   [ 1 Δ κ Δ φ Δ κ 1 Δ ω Δ φ Δ ω 1 ]   [ ( ρ + Δ ρ )   sin ( S   β ) 0 ( ρ + Δ ρ )   cos ( S   β ) ] = X o + [ cos κ sin κ 0 sin κ cos κ 0 0 0 1 ]   [ Δ X Δ Y Δ Z ] + [ cos κ sin κ 0 sin κ cos κ 0 0 0 1 ]   [ 1 Δ κ Δ φ Δ κ 1 Δ ω Δ φ Δ ω 1 ]   [ x 0 z ]
where,
Δ X , Δ Y , Δ Z are the components of the lever-arm offset vector P G ,
S is the scale factor for the mirror angle β (this scale factor should be unity for a bias-free system),
z is the vertical coordinate of the laser point with respect to the laser unit coordinate system, and
x is the lateral coordinate of the laser point with respect to the laser unit coordinate system, which is the lateral distance (with the appropriate sign) between the LiDAR point in question and the projection of the flight trajectory onto the ground.
The same mathematical derivation presented in section 2.1.1, through Equations 3 – 6, will be utilized here. The only difference is that the LiDAR mathematical model expressed by Equation 21 is now considered. By differentiating Equation 21 with respect to the system parameters, the form in Equation 22 can be derived.
[ δ X G δ Y G δ Z G ] T o t a l = [ X B i a s e d X T r u e Y B i a s e d Y T r u e Z B i a s e d Z T r u e ] = [ cos κ δ Δ X sin κ δ Δ Y + sin κ z δ Δ ω + cos κ z δ Δ φ sin κ δ Δ X + cos κ δ Δ Y c os κ z δ Δ ω + sin κ z δ Δ φ δ Δ Z x δ Δ φ ] + + [ sin κ x δ Δ κ cos κ sin ( S β ) δ Δ ρ + cos κ z β δ S cos κ x δ Δ κ sin κ sin ( S β ) δ Δ ρ + sin κ z β δ S cos ( S β ) δ Δ ρ x β δ S ]
By rewriting Equation 22 for two overlapping strips (A and B) and subtracting the resulting equations from each other, we can derive the mathematical relationship between conjugate points in overlapping strips. The result is shown in Equation 23. In a similar fashion, the mathematical relationship between control points, if available, and the LiDAR points can be derived by rewriting Equation 22 for the control surface (A) and the LiDAR strip (B) and subtracting the resulting equations from each other. Since the term f / x δ x will be equal to zero for the control surface, the discrepancy between the control and LiDAR points will reduce to the form in Equation 24. It is important to mention that Equation 24 can be written for all LiDAR strips which overlap with the control surface.
Equations 23 and 24 are the final linear observation equations, when dealing with overlapping strips and control surfaces, respectively. These equations allow us to recover the biases in the system parameters ( δ Δ X , δ Δ Y , δ Δ z , δ Δ ω , δ Δ φ , δ Δ κ , δ Δ ρ , δ S ) . As already mentioned in the Simplified method, the vertical bias in the lever-arm offset parameters ( δ Δ Z ) cannot be detected when using only overlapping strips. Such inability is caused by the fact that a vertical bias in the lever-arm offset parameters produces the same effect regardless of the flying direction, flying height, or scan angle. To estimate such bias, control information would be needed. In the absence of control data, the inability of estimating ( δ Δ Z ) is not critical since this quantity can be obtained with high accuracy from field survey. It should be noted that vertical control will contribute towards the estimation of systematic errors that produce vertical discrepancies between the LiDAR points and the control surface. Therefore, vertical control contributes towards the estimation of the vertical lever-am offset component ( δ Δ Z ) , the range ( δ Δ ρ ) , the roll ( δ Δ φ ) and the mirror angle scale ( δ S ) biases (refer to the third row in Equation 24). However, we cannot recover ( δ Δ Z ) and ( δ Δ ρ ) simultaneously due to the high correlation between these parameters. As can be observed in the third line in Equation 24, the vertical lever-arm offset bias ( δ Δ Z ) and the impact of the range bias ( cos ( S β ) δ Δ ρ ) are highly correlated in the range of the utilized scan angle (i.e., cos ( S β ) 1 for β [ 25 ° , + 25 ° ] ). To avoid correlation between ( δ Δ φ ) and ( δ S ) , well distributed data in the across flight direction should be used. The use of full control data over sloped surfaces will contribute towards the estimation of all parameters.
Once the biases are recovered, we can reconstruct the corrected point cloud using Equation 25. In this equation, the terms ( δ Δ ^ X , δ Δ ^ Y , δ Δ ^ z , δ Δ ^ ω , δ Δ ^ φ , δ Δ ^ κ , δ Δ ρ ^ , δ S ^ ) are the estimated biases in the system parameters.
[ X A Y A z A ] Biased [ X B Y B z B ] Biased = [ ( cos κ A cos κ B ) δ Δ X   ( sin κ A sin κ B ) δ Δ Y ( sin κ A sin κ B ) δ Δ X + ( cos κ A cos κ B ) δ Δ Y 0 ] + + [ ( sin κ A   z A sin κ B   z B ) δ Δ ω + ( cos κ A   z A cos κ B   z B ) δ Δ φ ( sin κ A   x A sin κ B   x B ) δ Δ κ ( cos κ A   z A cos κ B   z B ) δ Δ ω + ( sin κ A   z A sin κ B   z B ) δ Δ φ + ( cos κ A   x A cos κ B   x B ) δ Δ κ ( x A x B ) δ Δ φ ] + + [ [ cos κ A   sin ( S β A ) cos κ B   sin ( S β B ) ]   δ Δ ρ + ( cos κ A   z A   β A cos κ B   z B   β B ) δ S [ sin κ A   sin ( S β A ) sin κ B   sin ( S β B ) ]   δ Δ ρ + ( sin κ A   z A   β A sin κ B   z B   β B ) δ S [ cos ( S β A ) cos ( S β B ) ]   δ Δ ρ ( x A   β A x B   β B ) δ S ]
[ X A Y A z A ] ControlSurface [ X B Y B z B ] Biased = [   cos κ B δ Δ X + sin κ B δ Δ Y sin κ B δ Δ X   cos κ B δ Δ Y δ Δ z ] + + [   sin κ B   z B δ Δ ω cos κ B   z B δ Δ φ + sin κ B   x B δ Δ κ cos κ B   z B δ Δ ω sin κ B   z B   δ Δ φ   cos κ B   x B δ Δ κ x B δ Δ φ ] + [ cos κ B   sin ( S β B )   δ Δ ρ   cos κ B   z B   β B δ S sin κ B   sin ( S β B ) δ Δ ρ   sin κ B   z B   β B δ S cos ( S β B ) δ Δ ρ + x B   β B δ S ]
[ X Y z ] Corrected = [ X Y z ] Biased [ cos κ   δ Δ ^ X sin κ   δ Δ ^ Y + sin κ   z   δ Δ ^ ω + cos κ   z   δ Δ ^ φ sin κ   δ Δ ^ X + cos κ   δ Δ ^ Y cos κ   z   δ Δ ^ ω + sin κ   z   δ Δ ^ φ δ Δ ^ Z x   δ Δ ^ φ ] [ sin κ   x   δ Δ ^ κ cos κ   sin ( S β ) δ Δ ρ ^ + cos κ   z   β δ S ^ cos κ   x   δ Δ ^ κ sin κ   sin ( S β ) δ Δ ρ ^ + sin κ   z   β δ S ^ cos ( S β ) δ Δ ρ ^ x   β δ S ^ ]
The procedure for estimating the necessary quantities ( x , z ,  and  k ) presented in Equations 23 and 24 using the available data (time-tagged point cloud and trajectory positions) is as follows:
(I) For a LiDAR point mapped at time (t), we identify trajectory positions within a certain time interval (t – ∆t, t + ∆t);
(II) Then, a straight line is fitted through the selected trajectory positions to come up with a local estimate of the trajectory. After defining the local trajectory, the necessary quantities can be estimated as follows:
  • x, which is the lateral coordinate of the laser point with respect to the laser unit coordinate system, can be determined by computing the normal distance (with the appropriate sign) between the LiDAR point and the interpolated trajectory data. The intersection of the normal from the LiDAR point to the interpolated trajectory will define the position of the trajectory at time t;
  • z, which is the vertical coordinate of the laser point with respect to the laser unit coordinate system, can be determine by subtracting the elevation of the laser firing point (H) at time t, given by the interpolated flight trajectory, from the LiDAR point elevation (Z), i.e., z = Z − H; and
  • κ, which is the trajectory heading, can be computed once we have the local estimate of the trajectory and its direction (defined by the neighbouring trajectory positions);
  • β, which is the scan angle, can be computed by simple trigonometric operation using the estimated lateral distance and the trajectory height above ground.
As already mentioned, there is no direct correspondence between points in overlapping LiDAR surfaces. Similar to the discrepancy detection procedure presented in the Simplified method, we will represent one strip using the original points, while the second strip will be represented by triangular patches, which can be derived from a Triangular Irregular Network (TIN) generation procedure. In addition to the previously listed advantages of utilizing points/TIN patches, such primitives provide a way of describing the surface while having an explicit link between the surface representation scheme and the LiDAR equation. In the next section, the point-to-patch matching procedure along with the similarity measure, which utilizes the matched primitives and the derived mathematical model to come up with an estimate of the biases in the system parameters, will be introduced.

2.2.2. Similarity measure and primitives matching

As abovementioned, we will represent one strip using the original points, denoted by S2, while the second strip, denoted by S1, is represented by triangular patches, which can be derived from a Triangular Irregular Network (TIN) generation procedure. When dealing with control information, the control data is represented by points and the overlapping LiDAR strips are represented by TIN patches. The correspondence between points in S2 and patches in S1 is established in an iterated manner, using the Iterative Closest Patch (ICPatch) procedure, which was already described in Section 2.1.2.2.
Differently from the discrepancy detection used in the Simplified method that utilizes the volume constraint as the similarity measure, a point-based similarity measure will be introduced for the Quasi-rigorous method. The triangular patch will be represented by one of its vertices together with the orientation of its surface normal. Since the patch vertex in S1 and the LiDAR point in S2 might not be conjugate, this fact should be considered by the implemented similarity measure. Therefore, a point-based similarity measure, which can deal with non-conjugate points, is proposed. More specifically, non-conjugate points will be used in the constraints in Equations 23 and 24 to derive an estimate of the biases in the system parameters. In order to compensate for fact that the points (e.g., patch vertex point X B in S 1 and point X A in S 2 ) are not conjugate, one can manipulate the weight matrix of the patch vertex point. For the vertex point X B of a triangular patch, zero weights will be assigned to that point along the patch plane. This weight restriction will ensure the minimization of the random component between conjugate point-patch pairs along the patch normal in overlapping strips, after applying the estimated system biases (as will be explained later in this section). First, a local coordinate system (UVW) with the U and V axes aligned along the triangular plane is defined. The relationship between the strip coordinate system (XYZ) and the local coordinate system (UVW) can be represented by Equation 26. The rotation matrix in that equation is defined using the orientation of the normal to the triangular planar patch including the vertex in question. The original weight matrix P X Y Z , as shown in Equation 27, is defined as the inverse of the variance-covariance matrix Σ X Y Z , which depends on the accuracy specification of the data acquisition system. Using the law of error propagation, the weight of that vertex point in the local coordinate system P U V W can be derived according to Equation 28. As it has been noted earlier, the weight along the triangle plane normal is the only useful information when working with non-conjugate points along corresponding conjugate point-patch pairs. Therefore, the weight matrix can be modified according to Equation 29. Finally, the modified weight matrix P X Y Z ' in the strip coordinate system can be derived according to Equation 30.
[ U V W ] = R [ X Y Z ]
P X Y Z = X Y Z 1
P U V W = R   P X Y Z R T = [ P U P U V P U W P V U P V P V W P W U P W V P W ]
P U V W ' = [ 0 0 0 0 0 0 0 0 P W ]
P X Y Z ' = R T P U V W ' R
The evaluation of the redundancy in the calibration procedure will be based on the rank of the weight matrix for the individual points. For a point defining a planar patch, whose weight matrix has been manipulated according to Equation 29, a contribution of one will be considered towards the redundancy computation (i.e., one effective equation per constraint).
As mentioned earlier, the proposed weight restriction in Equation 29 will ensure the minimization of the random component between conjugate point-patch pairs, after applying the estimated system biases. In the Least Squares Adjustment (LSA) procedure, the unknown parameters are estimated to minimize the weighted-squared-sum of the mis-closure vector e , as represented by Equation 31. It should be noted that the mis-closure vector e in Equation 31 is composed of a random ( e A B ) and a non-random ( d X ) component. The random component is introduced by the random noise in the system measurements while the non-random component is introduced by the fact that we are dealing with non-conjugate points along the point-patch pair. Equations 32–35 demonstrate that utilizing the modified weight matrix, Equation 30, the non-random component is eliminated. Since the vector d U , shown in Equations 33 and 34, represents the difference vector relative to the plane coordinate system, d W will equal to zero (the non-conjugate points lie on the same plane). Therefore, the term d U T P u v w d U equals to zero and the LSA target function gets the form in Equation 35, where we have only the random component. Finally, Equations 36 and 37 show that the LSA target function, which contains only the random component, minimizes the weighted sum of the squared random components between conjugate point-patch pairs after the manipulation of the weight matrix according to Equation 29.
Σ e T   P   e = min
where,
e = [ ( X A Y A Z A ) T r u e ( X B Y B Z B ) T r u e ] + [ ( e X   A e Y   A e Z   A ) r a n d o m ( e X   B e Y   A e Z   B ) r a n d o m ] = [ d X d Y d Z ] + [ e X   A B e Y   A B e Z   A B ] = d X + e A B
[ d X + e A B ] T P X Y Z [ d X + e A B ]   d X T P X Y Z d X + ( e A B ) T P X Y Z e A B = min
d X T R T P u v w R d X + ( e A B ) T P X Y Z e A B   d U T P u v w d U + ( e A B ) T P X Y Z e A B = min
d U T P u v w d U = [ d U d V d W ]   [ 0 0 0 0 0 0 0 0 P w ]   [ d U d V d W ] = P W   d W 2 = 0
LSA Target Function = ( e A B ) T P X Y Z e A B
Σ [ e X   A B e Y   A B e Z   A B ] P X Y Z ' [ e X   A B e Y   A B e Z   A B ] Σ [ e X   A B e Y   A B e Z   A B ] R T P U V W ' R   [ e X   A B e Y   A B e Z   A B ] = min
Σ [ e U   A B e V   A B e W   A B ] P U V W ' [ e U   A B e V   A B e W   A B ] Σ [ e U   A B e V   A B e W   A B ] [ 0 0 0 0 0 0 0 0 P W ] [ e U   A B e V   A B e W   A B ] = P W e W   A B 2 = min
Where [eU AB eV AB eW AB] represent the random components of the mis-closure vector in the UVW coordinate system (i.e., the planar patch local coordinate system with eW AB being along the plane normal). For reliable estimation of the system biases, the utilized patches should have a balanced distribution of slopes and aspects.
In summary, the Quasi-rigorous calibration procedure proceeds as follows:
  • The correspondence between points in S2 and patches in S1 is established using the Iterative Closest Patch (ICPatch) procedure;
  • For each conjugate point-patch pair, (e.g., patch vertex X B in S1 and point X A in S2), one can write the observation equations similar to those in Equations 23 and 24, when dealing with overlapping strips and control surfaces, respectively;
  • The initial correspondences are used to derive an updated estimate of the system biases. Then, the estimated system biases are used to reconstruct both surfaces using Equation 25. Since after reconstructing the LiDAR point cloud the correspondence between point-patch pairs might have changed, a new set of correspondences has to be established. The new correspondences are utilized to derive a better estimate of the system biases.
  • Such a procedure is repeated until the corrections to the estimated calibration parameters are almost zero.
The main differences between the Quasi-rigorous and the Simplified methods can be summarized as follows:
  • The Simplified method utilizes the LiDAR point cloud only, while the Quasi-rigorous method utilizes time-tagged LiDAR point cloud and trajectory position data;
  • The Simplified method requires parallel flight lines with small pitch and roll angles and minor terrain elevation variations compared to the flying height above ground, while the Quasi-rigorous method only requires small pitch and roll angles;
  • The Simplified method is a two-step procedure, while the Quasi-rigorous method is a single-step procedure;
  • The Simplified method cannot handle the use of control information since the derived equations are based on the relationship between conjugate surface elements in overlapping strips. More specifically, the derived mathematical model relating conjugate surface elements in overlapping strips, which is a rigid body transformation (3 shifts and 1 rotation angle around the flight direction), is not valid for relating conjugate LiDAR and control surfaces. The Quasi-rigorous method, on the other hand, can handle control data through Equation 24.

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.
Remotesensing 02 00874 g005
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 OverlapDirection
(i) Strips 1&2100%Opposite directions
(ii) Strips 3&4100%Opposite directions
(iii) Strips 3&550%Same direction
(vi) Strips 1&670%Opposite directions
(v) Strips 5&750%Opposite directions
(vi) Strips 7&840%Opposite directions
(vii) Strips 2&670%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 ScenarioOverlapping 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 ( X T ' and R ω ϕ κ ' ) in the local coordinate system given the transformation parameters ( X T and R ω ϕ κ ) in the ground coordinate system and the heading angle Κ . Table 3 presents the recalculated transformation parameters (w.r.t. the local coordinate system).
X T ' = R Κ X T
R ω ϕ κ ' = R Κ R ω ϕ κ R Κ T
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.
X T ' (m) Y T ' (m) Z T ' (m) ω ' (sec) ϕ ' (sec) κ ' (sec)
Strips 1&2−0.251.27−0.0110.541.3467.68
Strips 3&4−0.010.520.02−2.59−9.7220.52
Strips 3&5−0.32−0.19−0.067.2096.486.48
Strips 2&6−0.420.150.01−4.6574.61−136.10
Strips 1&6−0.671.51−0.064.6891.0871.64
Strips 5&7−0.130.550.040.6812.96−5.40
Strips 7&80.480.820.061.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.
[ X T Y T ϕ ] S t r i p s   1 & 2 = [ 2   δ Δ X     2   H 2   δ Δ φ 2   δ Δ Y + 2   H 2   δ Δ ω 2 δ Δ φ ] = [ 0.25 1.27 1.34 ]
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 δ Δ X (m) δ Δ Y (m) δ Δ ω (") δ Δ φ (") δ Δ κ (") δ Δ ρ (m) δ S
I−0.07−0.1175−1800.260.000565
II−0.08−0.1175−1−39.50.220.000567
III−0.11−0.1086−16410.280.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 1Strip 3
I vs. III vs. IIIII vs. IIII vs. III vs. IIIII 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.010.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.020.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 δ Δ X (m) δ Δ Y (m) δ Δ ω (") δ Δ φ (") δ Δ κ (") δ Δ ρ (m) δ S
I−0.07−0.1179.6−2.652.50.590.00038
II−0.08−0.1280.8−4.817.30.000.00097
III−0.05−0.0984.3 −1.232.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 1Strip 3
I vs. III vs. IIIII vs. IIII vs. III vs. IIIII 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.130.03 ± 0.12−0.02 ± 0.000.00 ± 0.01−0.00 ± 0.000.00 ± 0.01
(0.13)(0.13)(0.02)(0.01)(0.00)(0.01)
Z (m) 0.62 ± 0.03 0.71 ± 0.04−0.09 ± 0.01 0.60 ± 0.01 0.69 ± 0.01−0.09 ± 0.00
( 0.62 )( 0.71 )(0.09)( 0.60 )( 0.69 )(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.
Remotesensing 02 00874 g006aRemotesensing 02 00874 g006b
Figure 7. Profiles location along overlapping strips 1&2.
Figure 7. Profiles location along overlapping strips 1&2.
Remotesensing 02 00874 g007
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.
Remotesensing 02 00874 g008
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 CalibrationAfter Calibration Simplified MethodAfter Calibration Quasi-rigorous Method
IIIIIIIIIIII
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.251.27-0.010.090.030.010.17-0.050.010.40 0.160.030.200.040.010.210.060.020.210.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.541.3467.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.671.51-0.060.25-0.310.030.51-0.150.030.52 -0.100.05 0.38-0.250.040.28-0.170.020.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.6891.0871.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.060.380.110.05-0.080.180.06-0.07 0.04 0.02 -0.080.020.05-0.080.090.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.2096.486.482.84-45.00-9.36-5.6317.56-5.76-5.04 40.32 -6.12 -5.7615.12-7.56-5.1568.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.

Acknowledgements

This work was supported by the Canadian GEOIDE NCE Network (SII-72) and the National Science and Engineering Council of Canada (Discovery Grant). The authors would like to thank LACTEC, Institute of Technology for the Development, Brazil for providing the LiDAR data and the valuable feedback.

References and Notes

  1. Vaughn, C.R.; Bufton, J.L.; Krabill, W.B.; Rabine, D.L. Georeferencing of airborne laser altimeter measurements. Int. J. Remote Sens. 1996, 17, 2185–2200. [Google Scholar] [CrossRef]
  2. Schenk, T. Modeling and Analyzing Systematic Errors in Airborne Laser Scanners; Technical Report in Photogrammetry No. 19; Ohio State University: Columbus, OH, USA, 2001; pp. 1–42. [Google Scholar]
  3. El-Sheimy, N.; Valeo, C.; Habib, A. Digital Terrain Modeling: Acquisition, Manipulation and Applications, 1st ed.; Artech House Remote Sensing Library: Boston, MA, USA, 2005; pp. 1–256. [Google Scholar]
  4. Kilian, J.; Haala, N.; Englich, M. Capture and evaluation of airborne laser scanner data. Int. Arch. Photogramm. Remote Sens. 1996, 31, 383–388. [Google Scholar]
  5. Crombaghs, M.; De Min, E.; Bruegelmann, R. On the adjustment of overlapping strips of laser altimeter height data. Int. Arch. Photogramm. Remote Sens. 2000, 33, 230–237. [Google Scholar]
  6. Maas, H.G. Method for measuring height and planimetry discrepancies in airborne laserscanner data. Photogramm. Eng. Remote Sens. 2002, 68, 933–940. [Google Scholar]
  7. Bretar, F.; Pierrot-Deseilligny, M.; Roux, M. Solving the strip adjustment problem of 3D airborne LiDAR data. In Proceedings of the IEEE IGARSS’04, Anchorage, AK, USA, 2004; pp. 4734–4737.
  8. Vosselman, G. Strip offset estimation using linear features. In Proceedings of 3rd International Workshop Mapping Geo-Surficial Processes Using Laser Altimetry, Columbus, OH, USA, 2002; pp. 1–9.
  9. Filin, S.; Vosselman, G. Adjustment of airborne laser altimetry strips. Int. Arch. Photogramm. Remote Sens. 2004, 35, 285–289. [Google Scholar]
  10. Pfeifer, N.; Elberink, S.O.; Filin, S. Automatic tie elements detection for laser scanner strip adjustment. Int. Arch. Photogramm. Remote Sens. 2005, 36, 1682–1750. [Google Scholar]
  11. Kager, H. Discrepancies between overlapping laser scanning strips- simultaneous fitting of aerial laser scanner strips. In Proceedings of the International Society for Photogrammetry and Remote Sensing XXth Congress, Istanbul, Turkey, 2004; Vol. 34, pp. 555–560.
  12. Filin, S. Calibration of spaceborne and airborne laser altimeters using natural surfaces. Ph.D. Dissertation, Department of Civil and Environmental Engineering and Geodetic Science, the Ohio-State University, Columbus, OH, USA, 2001; pp. 1–129. [Google Scholar]
  13. Skaloud, J.; Lichti, D. Rigorous approach to bore-sight self-calibration in airborne laser scanning. ISPRS J. Photogramm. Remote Sens. 2006, 61, 47–59. [Google Scholar] [CrossRef]
  14. Friess, P. Toward a rigorous methodology for airborne laser mapping. In Proceedings of EuroCOW, Castelldefels, Spain, 2006. [CD-ROM].
  15. Morin, K.W. Calibration of Airborne Laser Scanners. M.S. Thesis, University of Calgary, Alberta, Canada, 2002; pp. 1–125. [Google Scholar]
  16. Burman, H. Calibration and Orientation of Airborne Image and Laser Scanner Data Using GPS and INS. Ph.D. Dissertation, Royal Institute of Technology, Stockholm, Sweden, 2000; pp. 1–111. [Google Scholar]
  17. Toth, C.K. Calibrating airborne LiDAR systems. In Proceedings of ISPRS Commission II Symposium, Xi’an, China, 2002; pp. 475–480.
  18. Habib, A.; Bang, K.; Kersting, A.P.; Lee, D.C. Error budget of LiDAR systems and quality control of the derived data. Photogramm. Eng. Remote Sens. 2009, 75, 1093–1108. [Google Scholar] [CrossRef]
  19. Habib, A.; Cheng, R.; Kim, E.; Frayne, M.; Ronsky, J. Automatic surface matching for the registration of LiDAR data and MR imagery. ETRI J. 2006, 28, 162–173. [Google Scholar] [CrossRef]
  20. Zhang, Z. Iterative point matching for registration of free-form curves and surfaces. Int. J. Comput. Vis. 1994, 13, 119–152. [Google Scholar] [CrossRef]
  21. Habib, A.; Kersting, A.P.; Bang, K. Detecting systematic biases and GNSS/INS Drifts in LiDAR data. In Proceedings of ION 2009 International Technical Meeting, Anaheim, CA, USA, 2009.

Share and Cite

MDPI and ACS Style

Habib, A.; Bang, K.I.; Kersting, A.P.; Chow, J. Alternative Methodologies for LiDAR System Calibration. Remote Sens. 2010, 2, 874-907. https://doi.org/10.3390/rs2030874

AMA Style

Habib A, Bang KI, Kersting AP, Chow J. Alternative Methodologies for LiDAR System Calibration. Remote Sensing. 2010; 2(3):874-907. https://doi.org/10.3390/rs2030874

Chicago/Turabian Style

Habib, Ayman, Ki In Bang, Ana Paula Kersting, and Jacky Chow. 2010. "Alternative Methodologies for LiDAR System Calibration" Remote Sensing 2, no. 3: 874-907. https://doi.org/10.3390/rs2030874

Article Metrics

Back to TopTop