Method and apparatus for automatic calibration of mobile LiDAR systems

Abstract

A method and apparatus for the automatic calibration method described in this application provides an integrated framework for performing a reliable and objective estimation of IMU-LiDAR latency and boresight angles. This method, based on the estimation of calibration parameters through the resolution of observation equations is able to deliver boresight and latency estimates as well as their precision. A new calibration method for the boresight method angles between a LiDAR and an IMU, based on an automatic data selection algorithm, followed by the adjustment of bore sight angles. This method, called LIBAC (LiDAR-IMU Boresight Automatic Calibration), takes in input overlapping survey strips following a sample line pattern over a regular slope. First, construct a boresight error observability criterion, used to select automatically the most sensitive soundings to boresight errors. From these soundings, adjust the boresight angle 3D, thus taking into account the coupling between angles. From a statistical analysis of the adjustment results, we derive the boresight precision.

Claims

1. A method of automatically calibrating a mobile LiDAR system comprising the steps of: providing a remotely controlled UAV equipped with a global network satellite (GNSS) system, an inertial measurement unit (IMU), and a light detection and ranging (LiDAR) device; calibrating LiDAR-IMU roll, pitch and yaw boresight angles and lever arms between the positioning reference point (PRP) and the optical center (OC) of the LiDAR to eliminate systematic errors resulting from misalignment between the LiDAR and IMU measurement frames; activating the LiDAR to acquire natural terrain data via scanning natural terrain through overlapping predetermined strips; subdividing, based on the natural terrain data, the strips into small surface elements whenever they can be fitted into a common planar surface; employing an automatic planar surface element selection process most relevant for boresight estimation, adjusting the selection points to planar surface elements to optimizing boresight angles, and storing resultant data in a system memory device.

2. The method of automatically calibrating a mobile LiDAR system of claim 1, further comprising the step of applying a boresight sensitivity criteria to select a surface element most relevant for boresight estimation.

3. The method of automatically calibrating a mobile LiDAR system of claim 1, wherein the step of scanning preselected points from natural terrain or man-made structures though overlapping strips comprises serially traversing at least two overlapping and opposite parallel survey lines wherein the outer beam of the second line overlaps points scanned at the nadir of the first line.

4. The method of automatically calibrating the mobile LiDAR system of claim 1, further comprising the step of repositioning said remotely controlled aerial UAV aircraft to alternate selected points of natural terrain and repeating said steps.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The present invention will now be described, by way of example, with reference to the accompanying drawings, in which:

(2) FIG. 1, is a perspective view of a quadrotor type UAV aircraft programmed and operated remotely by a ground based operator for LiDAR scanning of natural earth surfaces;

(3) FIG. 2, is a perspective view of the quadrotor type UAV aircraft of FIG. 1 carrying a lightweight LiDAR pod system illustrated on an enlarged scale;

(4) FIG. 3, is a schematic of a LiDAR pod system of FIG. 2 including a LiDAR network, an INS (or IMU) and a GPS receiver;

(5) FIG. 4, is a block schematic diagram of the basic systems of the UAV aircraft with LiDAR of FIG. 2;

(6) FIGS. 5A and 5B, represent the line pattern of fully automated LiDAR-IMU Boresight Automatic Calibration (LIBAC) system, formed by two overlapping and opposite survey lines, wherein the outer beam of the second line should overlap the points scanned at the Nadir of the first line;

(7) FIGS. 6A, 6B and 6C, represent the line pattern of a fully automated LiDAR-IMU Latency Automatic Calibration (LILAC) system, formed by one survey line over relatively flat and smooth terrain, wherein the UAV should perform this line with maximum allowable roll oscillations, generating high roll velocities;

(8) FIG. 7 represents a logic flow chart of the LiDAR IMU Boresight Automatic Calibration, or LIBAC, depicting the first phase of the method consisting in geo-referencing the point cloud from the calibration strips and to compute for each point a Combined Standard Measurement Uncertainty value taking into account errors from latency, LiDAR grazing angle, IMU, positioning, boresight, lever-arms, wherein each point from the point cloud is thus supplemented with its a priori uncertainty; and

(9) FIG. 8 represents a logic flow chart of the LiDAR IMU Latency Automatic Calibration, or LILAC, depicting a first calibration method devoted to IMU-LiDAR latency determination from neutral terrain data.

(10) Although the drawings represent embodiments of the present invention, the drawings are not necessarily to scale and certain features may be exaggerated in order to illustrate and explain the present invention. The exemplification set forth herein illustrates an embodiment of the invention, in one form, and such exemplifications are not to be construed as limiting the scope of the invention in any manner.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

(11) The present specification is organized as follows: First discussed is a general framework for system and mounting parameters of LiDAR systems. Then presented is the mathematical formulation of point geo-referencing from LiDAR systems data together with the problem setting of boresight calibration on natural surfaces in the second section. It will be shown that one of the key point of boresight calibration on natural surfaces lies on appropriate data selection methods.

(12) LiDAR System Parameter Calibration

(13) A typical LiDAR survey system consists of a positioning system giving the position of a Reference Point (RP), an IMU giving its orientation with respect to a Local Astronomic Frame (LAF) and the LiDAR itself delivering acoustic ranges along a swath oriented perpendicularly to the survey vessel trajectory. Let us mention that in practice, the IMU gives orientation in a Local Geodetic Frame (LGF), as most tactical grade IMU used in airborne of UAV surveying are unable to distinguish the (LAF) from the (LGF). Therefore, we shall denote the (LGF) frame by (n) (navigation frame) to avoid confusion between the geodetic and astronomical frames.

(14) The following frames will be used in the framework of this paper: The Local Geodetic Frame (LGF), that will be denoted by (n) and called the navigation frame; The IMU body frame, denoted by (bI); The LiDAR body frame, denoted by (bS).

(15) The main objective of this invention is to design a calibration method for estimating the frame transformation from the (bS) frame to the (bI) frame, denoted by C_(bS){circumflex over ( )}bI which depends on three boresight angles, denoted by &p the boresight roll angle, δθ, the boresight pitch angle and δψ, the boresight yaw angle. Refer Paragraph [00199]. We shall denote hereafter by on C_(F1){circumflex over ( )}F2 the direction cosine matrix corresponding to the transformation from frame F1 to F2. Refer Paragraph [00200]. The LiDAR system geometrical parameters that need to be known by the user are the boresight angles (or mounting angles), the lever arms between the position Reference Point (RP) and the Optical Center of the LiDAR, as illustrated in FIG. 3, which constitutes a typical LiDAR system composed of a LiDAR, an INS (or IMU), a GPS receiver. The LiDAR frame is (bS), the IMU frame is (bI). The lever arm between the LiDAR and positioning Reference Point (RP) is a.sub.bI. The LiDAR ranging returns, coordinated in the (bS) frame, are denoted by r.sub.bS.

(16) Lever arms and boresight angles could be jointly estimated using the same adjustment method. However, data sets to be used for boresight calibration and for lever arms calibration should be distinct, since the effect of these two systematic errors cannot generally be observed simultaneously. Therefore, the methodology we propose is based on two distinct boresight and lever arm calibration processes, based on data selection procedures dedicated to each parameter. However, this paper will only focus on the estimation of boresight angles.

(17) LiDAR System Geo-Referencing Model

(18) Geo-referencing is a combination of raw data from GNSS, IMU and LiDAR to provide points, coordinated in the navigation frame. Geo-referencing of a given LiDAR return can be done as follows: From the beam angles, compute the launch vector expressed in the LiDAR (bS) frame; From IMU attitude data, coordinate the launch vector in the (n) frame; Using the beam travel time construct the path from the OC of the LiDAR to the terrain; Finally, using the LiDAR system PRP position, coordinate the point in the LGF.

(19) Denoting by X.sub.n(t)=(x, y, z).sup.τ.sub.n a geo-referenced point in the (n) frame at time t, we have:
X.sub.n(t)=P.sub.n(t)+C_(bI){circumflex over ( )}n(t−dt)(C_(bS){circumflex over ( )}BIr_bS(t)+a_bI) Refer Paragraph [00201].  (1)

(20) where, P.sub.n(t) is the position delivered by the GNSS receiver in frame (n), C_(bI){circumflex over ( )}n is the coordinate transformation from the IMU body frame to the navigation frame (which can be parameterized using Euler angles (φ, θ, ψ), denoting pitch, roll and yaw, respectively), the LiDAR return r.sub.bS, coordinated in the LiDAR frame (bS), the lever-arm vector coordinated in the IMU frame a.sub.bI and the boresight coordinate transformation C_(bS){circumflex over ( )}bI. Refer Paragraphs [00201 and 00199].

(21) In equation (1), t denotes the reference time from the GPS, which is supposed to be propagated to the IMU through a distributed time and message synchronization system. Thus, the OPS system and the IMU are supposed to share the same time base.

(22) The dependency of the calibration parameters on LIDAR points geo-referencing is described by equation (1), among them are: dt, then latency between the IMU and the LiDAR (it is to be noticed that in most LiDAR systems, latency between GNSS and the LiDAR impact can be considered as negligible, but latency between the LiDAR and IMU is not [Seube et al., 2012]; C_(bS){circumflex over ( )}bI, the boresight coordinate transformation. [Refer Paragraph 00199]; a.sub.BI, the lever-arms; The LiDAR range and beam angle, affecting the term r.sub.bS.

(23) Described here is an estimation procedure of the boresight coordinate transformation C_(bS){circumflex over ( )}bI and the latency between the IMU and the LiDAR. [Refer Paragraph 00199].

(24) General Principle of the Boresight and Latency UAV Calibration

(25) The calibration method is devoted to aerial UAVs and is in part automated, without supervision of the user. The LiDAR system is integrated on the UAV, including mission planning. The integrated hardware system is composed of a LiDAR sensor, a GNSS-INS unit, an embedded computer and a camera.

(26) “Boresight angles” being defined as misalignment between the INS frame and the LiDAR frame, and “Latency” being defined as time delay between the time-tag of INS and LiDAR data, the procedure for calibrating the boresight and latency parameters is as follows: 1. The UAV takes-off and the users flies it to a location where the terrain or man-made structures has a minimum slope of 10%; 2. Whenever the calibration procedure is launched by the user from the control panel, the UAV is making a 360° turn at a fixed point while scanning the terrain with the LiDAR system; 3. From this data set, the terrain steepest slope is determined; 4. The UAV computes a LIBAC line pattern in starting by a calibration line along the steepest slope. It flies the opposite line to make the calibration line 2; 5. The UAV computes the outer beam location of the right side of the first calibration line, and the calibration line 3 is done in making the Nadir to match the outer beams locations; 6. The UAV is flying the opposite line 4; 7. The UAV transmits the data acquired during the 4 lines to the ground; 8. The UAV is flying again line 3 with the highest possible roll oscillations; 9. The UAV transmits line 3 with roll oscillation to the ground;

(27) The calibration is done using post-processed GNSS reference point positions and therefore in post-processing. Whenever the post-processing is performed: 10. The IMU-LIDAR latency calibration is achieved by applying LILAC; 11. The data set acquired for the boresight calibration is corrected from the estimated IMU-LIDAR latency; 12. Boresight calibration is performed using LIBAC using the dataset corrected for latency.

(28) Once the calibration is done, the ground computer georeferenced the LiDAR dataset corresponding to the 5 survey lines; and computes the uncertainty along surface elements normal on a grid which step size is defined by the user. The result is the uncertainty map of the calibrated LiDAR system.

(29) The calibration is based on two methods: LILAC for latency calibration and LIBAC for boresight estimation.

(30) The calibration procedure consists in: 1. Acquiring data along a survey line with roll motion over a relatively flat surface. This line is fed to LILAC, which output a latency estimate. This latency estimate is then set in a geo-referencing software to correct the point cloud for IMU-LiDAR latency. 2. Acquiring data along 4 survey lines, following the line pattern described in FIGS. 5A and 5B. Then, boresight calibration can occur, in performing a line pattern consisting in two quasi parallel lines over a slope. The outcomes of LIBAC are three boresight angles between the IMU frame and the LiDAR frame and the associated standard deviations.

(31) LILAC (LiDAR IMU Latency Automatic Calibration)

(32) LILAC is a calibration method devoted to IMU-LiDAR latency determination from natural terrain data. The estimation is performed by adjusting a set of points on a quadratic surface, meanwhile determining the latency and the surface parameters. Selected data used for the calibration are done using a dynamic criterion (high angular velocities and low angular acceleration).

(33) Referring to FIGS. 6A-6C. LILAC needs only one line on a relatively flat and smooth terrain. Since latency is strongly related to angular velocities, it is advised to have the highest variation of attitude while running the calibration lines.

(34) Line Pattern for LILAC

(35) LILAC requires one calibration line without any overlap on a smooth terrain, but with relatively high roll dynamics, as shown in FIGS. 6A-6C. The line length range is 20-30 meters at a velocity of 1-2 m/s.

(36) To adjust the IMU to LiDAR latency, we derive a parametric model depending explicitly on dt. Considering a first order approximation of the frame transformation matrix C_(bI){circumflex over ( )}n (t), equation (1) [Refer Paragraph 00202] can be re-written by:
X_n(t)=P_n(t)+C_(bI){circumflex over ( )}n(t)(Id−Ω_(n/bI){circumflex over ( )}bI)(C_bS){circumflex over ( )}bIr_bS(t−dt)+a_bI). Refer Paragraph [00203].

(37) Where Ω_(n/I){circumflex over ( )}bI is the angular velocity skew-symmetric matrix between the IMU frame and the navigation frame, coordinated in the IMU frame. [Refer Paragraph 00204].

(38) The next step is to determine the IMU/INS latency by the following rule: dt should be adjusted such that all soundings from a supposed smooth and regular seafloor should belong to the same quadratic surface. The plane or quadratic equation model that will be used is defined as follows:
Z=s(X,Y,R)

(39) where s is surface cartesian equation depending on 3 (for a plane) or 6 (for a quadratic surface) parameters (denoted by R), and X, Y, Z are the coordinates of georeferenced soundings.

(40) Points selected for calibration should have the two following characteristics: High angular rate to select data that are the most sensitive to latency; Small angular acceleration to comply with the assumptions of the first order Taylor expansion of C_(bI){circumflex over ( )}n(t−dt) [Refer Paragraph 00205]; and To lie on a supposed smooth and regular seafloor.

(41) To satisfy these three conditions, a data selection module was implemented. This module checks these assumptions and feeds an adjustment module with appropriate sounding and raw data from the MBES systems.

(42) The estimation problem will therefore determine both the surface parameters and the latency dt of interest.

(43) Any selected sounding should verify the above equation and the geo-referencing equation. In other words, if
X_n=(.square-solid.(x_n@y_n@z_n))=P_n(t)+C_(bI){circumflex over ( )}n(t)(Id−Ω_(n/bI){circumflex over ( )}bI)(C_(bS){circumflex over ( )}bIr_bS(t−dt)+a_bI) Refer Paragraph [00206]:

(44) then we should have the following observation equation.
Z.sub.n=ƒ(x.sub.n,y.sub.n,R)

(45) One can readily see that this last equation can be re-written as a function of the unknown parameters
g(R,dt)=0

(46) and that we have one equation per point. We can therefore solve this by applying a least square method to the linearized problem, following an iterative least square procedure.

(47) In the present approach, the least square weight matrix requires the knowledge of the Combined Standard Measurement Uncertainty of each sounding. As the weight matrix is determined according to the model, it can be used to have an estimate of the precision of the estimated parameters. This procedure is conducted following Least Squares statistical tools.

(48) Data selection is a requirement to feed LILAC with data satisfying the assumptions. Selected points correspond to high angular rate and low angular acceleration. This computation is done by using the direct output of IMU angular rates. To be more explicit, the norms of angular rate and angular acceleration associated to roll and pitch are computed. Then data corresponding simultaneously to high norm of angular rate and low norm of angular acceleration are selected for adjustment.

(49) LIBAC (LiDAR IMU Boresight Automatic Calibration)

(50) Line Pattern for LIBAC

(51) LIBAC needs four calibration strips lines over a smooth slope, as shown in FIGS. 5A and 5B. The line pattern of LiBAC formed by two overlapping and opposite survey lines, wherein the outer beam of the second line should overlap the points scanned at the Nadir of the first line. The length of each survey line is typically 40-50 meters. The height is 10-20 meters. The slope is at least 10%. The end points of the slope should not be exceeded.

(52) Principle of LIBAC

(53) The first phase of the method consists in geo-referencing the point cloud from the calibration strips and to compute for each point a Combined Standard Measurement Uncertainty value taking into accounts errors from latency, LiDAR grazing angle, IMU, positioning, boresight, lever-arms. Each point from the point cloud is thus supplemented with its a priori uncertainty.

(54) The second phase consists in automatically finding planar surface elements from each calibration strip. This is done by a quad-tree algorithm that defines recursively a grid on which the LiDAR points can be modeled by a planar surface. For each cell of the quad-tree, if there is enough density of point originating for at least 2 different strips, if the point cloud is isotropic, and if those points can be modelled by a planar surface thanks to a Deming Least Square algorithm, then the cell is tagged as eligible for calibration.

(55) The third phase consists in computing a boresight sensitivity criterion, which measures the variation of the cell's centroid elevation with respect to boresight. This sensitivity criterion is used to select the number of cells on which the boresight adjustment could be performed.

(56) The last phase adjusts all points from each selected cell to belong to the same parametric (planar or quadratic) surface elements by optimizing boresight angles. The model used to perform this adjustment is not a priori linearized with respect to boresight angles. In other words, we do not consider the boresight transformation matrix to be a micro-rotation represented by a skew-symmetric matrix. We keep the non-linear representation of boresight and we make the adjustment by an iterative least square adjustment using this non-linear model.

(57) Data Selection for Boresight Adjustment

(58) The purpose of data selection is twofold:

(59) First, to extract planar areas from the calibration strips. Indeed, we shall see in the next section that the observation equation we use for boresight adjustment express the fact that any point X.sub.n (as defined by equation (1)) from an overlapping survey line and belonging to a planar surface element, should verify s(X.sub.n)=0 where s(x, y, z)=z−ax−by−c is a plane equation, or a quadratic equation s(x, y, z)=z−ax.sup.2−by.sup.2−cxy−dx−ey−ƒ.

(60) Secondly, to select the planar surface elements having the highest sensitivity to boresight errors. Indeed, the adjustment should be performed from data exhibiting the highest boresight error, in order to maximize the input data systematic error to noise ratio.

(61) Detection of Planar Surface Elements

(62) The present approach follows [Hebel et al. 2012] and [Skaloud 2007], but in addition, we propagate the point uncertainty to the parametric surface parameters uncertainty, and we test the normalized residual to verify the correct geometric modelling of each surface element.

(63) It consists in finding planar or quadratic surface elements from the terrain, for which we can observe the maximum effect of systematic error due to boresight. As we are constrained to find planar areas in natural surfaces we adopted a variable resolution approach based on a quad-tree decomposition. To decide if a quad-tree surface element is a plane, we use a Deming Least-Square (DLS) plane fitting method [Deming, 1943], [Moniot, 2009]. We use this method as it enables to take into accounts the propagated uncertainty of the LiDAR system on all points. Indeed, LiDAR returns are subjected to a certain level of uncertainty that should be considered by the plane fitting method. While we compute soundings thanks to the point geo-referencing model (1), we also compute a Combined Standard Measurement Uncertainty (CSMU) estimator. Using the CSMU of each point in input of the DLS method gives more reliable plan coefficient estimates than using a classical Total Least Square or a Principal Component Analysis estimator without weighted uncertainties.

(64) The quad-tree subdivision process termination test takes into accounts the number of points within the surface elements and tests the presence of different survey lines (at least two survey strips are required).

(65) When the subdivision of overlapping strips in planar areas is done, we look for the best surface elements to be used for boresight angles estimation. To do so, we construct a sensitivity criterion which computes the relevance of each surface element to boresight estimation. Points from the selected surface elements will be used for boresight adjustment. We remove surface elements produced by irregular point clouds (like trees, canopy) by letting the user to check the presence of irregularities and associated doubtful surface elements.

(66) Boresight Sensitivity Criterion for Surface Elements

(67) One of the key elements of the calibration procedure is the computation of boresight sensitivity indexes for each planar element, as detected by the quad-tree analysis.

(68) This section describes the selection process of the most relevant surface elements for boresight estimation. This phase is essential to minimize the size of the underlying boresight adjustment (It shall become apparent that the size of the boresight estimation problem is 3 P+3, where P is the number of selected surface elements). In addition to this, we should guarantee that over the selected surface elements, the errors between overlapping data due to boresight angles are maximum.

(69) On each surface element selected by the quad-tree process, the point cloud includes data from several overlapping strips. Let us consider a given planar surface element. We define C), the elevation difference between the centroid of the plane fitted with points from all survey strips and the centroid of the plane fitted with points from a strip j, namely;
C.sub.j=Z.sub.j−aX−bY−c  (2)
where
Z.sub.j=a.sub.jX+b.sub.jY+c.sub.j  (3)

(70) is the elevation computed with the plane coefficients a.sub.j, b.sub.j, c.sub.j adjusted with point data from strip j, and X, Y are the coordinates of the center of the surface element.

(71) To define the sensitivity criterion, we need to express C.sub.j, as a function of boresight angles. To do so, we introduce the virtual vector from the averaged position of OC positions P_n{circumflex over ( )}j [Refer Paragraph 00207] for strip j (which corresponds to LiDAR point within the given surface element) and the centroid of the surface element X_n{circumflex over ( )}j. Refer Paragraph [00208]. Note that this last point depends on j as its elevation is computed using the coefficients of the plane fitted using strip j point data, as given by equation (3).

(72) Therefore, X_n{circumflex over ( )}j=(X, Y, Z) and we can also write this point using an approximation of the geo-referencing equation
X_n{circumflex over ( )}j=P_n{circumflex over ( )}j+C.sup.−_bI{circumflex over ( )}n_j(C_(bS){circumflex over ( )}bIr_rS{circumflex over ( )}j(t−dt)) Refer Paragraph 00209].  (4)

(73) where C.sup.−_bI{circumflex over ( )}n_j is the average coordinate transformation matrix defined by the LiDAR system attitudes while the surface element is scanned by the LiDAR. [Refer Paragraph 00210]. To express the virtual LiDAR return we approximate it without boresight, as follows:
r_bS{circumflex over ( )}j=C.sup.−_bI_j(X_n{circumflex over ( )}j−P_n{circumflex over ( )}j) [Refer Paragraph 00211].   (5)

(74) Then, using equation (4) and (5), we can construct the difference of elevation in (2) as the function of boresight. The sensitivity criterion we use is the min-max error of the gradient of C.sub.j with respect to the boresight angles (i.e.); the variation of C.sub.j due to boresight angles), over all survey strips j. Indeed, a surface element is sensitive to boresight whenever the planar surfaces fitted from points of strip j have a significant elevation difference due to boresight variations.

(75) The global sensitivity criterion is an average of the roll, pitch and yaw sensitivity criteria. The method can be used for calibrating one angle: in this case the sensitivity criterion is limited to the one relative to this angle.

(76) Adjustment of Boresight Angles

(77) In this section, we explain how the data selected by planar surface elements are used for the estimation of the boresight angles. The principle of the method is to adjust the boresight angles in such a way that all the points belonging to selected surface elements satisfy the same plane or quadratic surface equation. Let X.sub.n=(x.sub.n, y.sub.n, z.sub.n).sup.T be a point from a planar or quadratic surface element (p), satisfying the following cartesian equation:
Z.sub.n=s(X.sub.n,Y.sub.p,Q)  (6)

(78) Where Q is the set of parameter defining the parametric surface (plane of quadratic surface).

(79) We shall denote by P the set of planar or quadratic surface elements selected for calibration by the data selection method. For a given point, all elements of equation (1), except the boresight angles δφ, δθ, δψ and the lever arms a.sub.bI are given by the LiDAR survey system sensors (positioning, IMU, LiDAR). In the following, we will suppose that we already know the lever arms.

(80) We write equation (1) as a function of the unknown variables δφ, δθ, δγ and (Q.sub.p).sub.p∈P, the rest of the parameters P.sub.n, C_bI{circumflex over ( )}n, r.sub.bS being measured or computed from the knowledge of LiDAR survey raw data, or assumed to be already known, like the lever arm a.sub.bI. Refer Paragraph [00202].

(81) From this, one can readily show that the left-hand side of equation (6) is a non-linear function denoted by ƒ:
Z.sub.n−s(X.sub.n,Y.sub.p,Q.sub.p)=ƒ(δφ,δθ,δψ,Q.sub.p)

(82) and which depend on the following variables: 3 boresight angles, and planar surface element coefficients (3 for a plane, 6 for a quadratic surface).

(83) The problem of boresight calibration is to find the boresight angles δφ, δθ, δψ and the plane equation parameters (Q.sub.p).sub.p∈P such that for all p∈P:
ƒ(δφ,δθ,δψ,Q.sub.P)=0  (7)

(84) A significant difference between our method and the methods presented in [Hebel et al. 2012] and [Skaloud 2007] lies in the fact that we do not consider that the boresight matrix is a linear operator (a skew symmetric matrix representing a micro-rotation). Indeed, we do not linearize the boresight direction transformation matrix as C_bS{circumflex over ( )}bI=I+Ω, where Ω is the skew-symmetric approximation of the micro. Refer Paragraph [00199]. Notice that this approximation is done in existing LiDAR calibration methods and may lead to significant boresight calibration errors.

(85) The data selection process is a very important component of the boresight calibration process. Its goal is indeed to select a relatively small number of planar surface elements for which the sensitivity to boresight errors is maximum.

(86) The boresight adjustment problem is a collection of equations of the type of equation (7), that can be solved by applying an iterative weighted least square approach. From the statistical analysis of the least square solution, we can get the boresight precision as the variance of the adjusted parameter of the least square problem.

Conclusion

(87) The automatic calibration method described in this application provides an integrated framework for performing a reliable and objective estimation of IMU-LiDAR latency and boresight angles. This method, based on the estimation of calibration parameters through the resolution of observation equations is able to deliver boresight and latency estimates as well as their precision.

(88) Referring to the drawing figures, and particularly to FIG. 1, a preferred embodiment of a LiDAR boresight automatic calibration apparatus and method in accordance with the present invention is illustrated. A quadrotor type UAV 10 can be operated remotely by a ground based operator 12 manipulating a base station 14 for programming and operating the UAV 10 as well as uploading and downloading data. The base station 14 is interconnected with the UAV 10 by a multi-channel radio frequency (RF) link 16. Mission planning information including flight control and operating instructions for the UAV 10 are either stored in controller memory of the base station 14 or a controller (not illustrated) within the drone aircraft 10. The UAV 10 carries an easily reconfigurable mounted instrument or operating package such as a LiDAR pod 18 which is cable connected and powered by the UAV 10.

(89) Varied configurations of UAVs 10 such as depicted herein can be employed in modified form to practice the present invention.

(90) The communications network can also include continuous access to a Global Navigation Satellite system (GNSS) (not illustrated) either directly to the UAV 10 or relayed through one or more local GNSS base stations 22 for recording data during flight. The base stations 22 are also interconnected directly with the UAV 10 and the base station 14 by additional multi-channel radio frequency (RF) links 24 and 26, respectively.

(91) The LiDAR system has an emitter which produces a sequence of outgoing pulses 20 of coherent collimated light that is transmitted in a given direction towards a predetermined ground based natural or man-made surfaces.

(92) As best seen in FIG. 2, the quadrotor type UAV 10 consists of a center pod 32 containing a battery, motor and control electronics. Resilient landing struts 34, as well as the LiDAR operating package 18 extend below the lower surface 36 of the center pod 32. The gimbal or fixed mount 46 is illustrated in an extremely simplified, schematic form, it being understood that it enables at least three degrees of rotational freedom (i.e., yaw, pitch and roll) of the operating package 30 with respect to the center pod 32 of the drone aircraft 10. Furthermore, it can enable independent bidirectional translation along three (e.g., X, Y and Z) axes of the operating package 18 with respect to the center pod 32 of the UAV 10. Four circumferentially arranged arms with motor pods 38 extend outwardly from the center pod 32, each supporting a rotor 40.

(93) Referring to FIG. 3, a typical LiDAR system 28 includes a Lidar 18, an INS (or IMU) and a GPS receiver 42. The LiDAR frame is (bS) and the IMU frame is (bI) based. The lever arm between the LiDAR and positioning Reference Point (RP) is a.sub.bI. The LiDAR ranging returns, coordinated in the (bS) frame, are denoted by r.sub.bS.

(94) Referring to FIG. 4, a block schematic diagram of the basic systems of the UAV aircraft with LiDAR of FIG. 2 include a laser 48, a scanner and associated optics 50, a photodetector and receiver electronics 52, and a position and navigation system 54 are each interfaced with a controller 56, including active and passive memory and a power supply.

(95) Referring to FIGS. 5A and 5B, boresight calibration is defined as capturing the angles between a LiDAR and the IMU body frames of reference. The goal of the calibration is to correct the misalignment between the two sensors. The applicant has developed a boresight calibration tool that can be used in the natural environment to facilitate flight operations.

(96) Flight planning and area specifications include: defining alignment waypoints, employing four (4) scan lines, providing a slope of 20-30% with a length of 40 to 50 m in length with no vegetation on the ground.

(97) The required flight path requires the following steps; LiDAR 18 before take-off. 2. Initiation. 3. Scan one scan line above the slope in one direction. 4. Turn around 180°. 5. Scan one line of 40 m above the slope in an opposite direction of the first line. 6. Fly 50 m perpendicular from the line (to obtain 50% overlap with the first two scan lines. 7. Scan one line of 40 m above the slope in one direction. 8. Turn around 180°. 9. Scan one line above the slope in opposite direction of the first scan line. 10. Perform waypoint alignment. 11. Follow the usual workflow when you land with md4-1000 and mdLiDAR payload.

(98) FIGS. 6A, 6B and 6C represent the line pattern of a fully automated LiDAR-IMU Latency Automatic Calibration (LILAC) system, formed by one survey line over relatively flat and smooth terrain, wherein the UAV should perform this line with maximum allowable roll oscillations, generating high roll velocities.

(99) FIG. 7 is a flowchart of an exemplary process of a logic flow of the LiDAR IMU Latency Automatic Calibration, or LILAC, depicting a first calibration method devoted to IMU-LiDAR latency determination from neutral terrain data. The depicted flowchart is initiated at a start step 100 feeding a priori values of boresight angles logic block 102, which, in turn feeds an input data (GNSS, IMU/INS, LiDAR) logic block 104. The input data block 104 feeds a process merging of data block 106. A LILAC processor 108 feeds an estimated; latency logic block 110, which, also feeds the process merging of logic block 106, with the combined resultant feeding a selection of the most relevant points (QUADTREE, SENSITIVITY CRITERION) to a logic block 112. Logic block 112 feeds logic block 114 which effects computation of matrices (MODEL, OBSERVATION, WEIGHTED) for the associated least squares problem. Logic block 114 feeds a logic block 116 which provides estimation of boresight correction with iterative weighted least squares. Logic block 118 which tests whether the hypothesis test is validated. If not, Logic block 118 feeds a new weight matrix in logic block 120, which, in turn, feeds back to logic block 114, and the process steps are repeated. If the hypothesis test is validated in logic block 118. The results are fed to logic block 122 for the correction of initial boresight and to provide quality information. Thereafter logic block 122 is fed to end step 124, and the algorithm is completed.

(100) FIG. 8 is a flowchart of an exemplary process of a logic flow of the LiDAR IMU Boresight Automatic Calibration, or LIBAC, depicting the first phase of the method consisting in geo-referencing the point cloud from the calibration strips and to compute for each point a Combined Standard Measurement Uncertainty value taking into account errors from latency, LiDAR grazing angle, IMU, positioning, boresight, lever-arms, wherein each point from the point cloud is thus supplemented with its a priori uncertainty. The depicted flowchart is initiated at a start step 200 feeding a priori values of latency logic block 202, which, in turn feeds an input data (GNSS, IMU/INS, LASER) logic block 204. The input data block 204 feeds a data block 208 for selection of the most relevant points. An a priori boresight angles data block 206 also feeds the logic block 208. Logic block 208 feeds logic block 210 providing Computation of Matrices (MODEL, OBSERVATION, WEIGHTED) for the associated Least Squares Program. Logic Block 210 feeds a logic block 212 which provides an Estimation of Latency Correction with Iterative Weighted Least Squares. Logic Block 212 feeds a Logic Block 214 which determines if the Hypothesis Test is Validated. If “No”, Logic Block 214 feeds a Logic Block 216 which establishes a New Weight Matrix and feeds it to logic block 210. If “Yes”, logic block 214 feeds a logic block 218 which provided for Correction of Initial Latency and Provides Quality Information. Thereafter logic block 218 is fed to end step 220, and the algorithm is completed.

(101) The following documents are deemed to provide a fuller background disclosure of the inventions described herein and the manner of making and using same. Accordingly, each the below-listed documents are hereby incorporated into the specification hereof by reference. [Rabine Keyetieu, Nicolas Seube and Stéfanie Van-Wierts], unpublished, “Boresight and Latency Automatic Calibration of LiDAR Systems on Natural Surfaces”. [Rabine Keyetieu, Nicolas Seube and Stéfanie Van-Wierts], unpublished, “Boresight Automatic Calibration of LiDAR Systems on Natural Surfaces”. U.S. Patent Application No. 2016/0291134 A1 to Droz et al. entitled “Long Range Steerable LiDAR System”. U.S. Pat. No. 9,285,477 B1 to Smith et al. entitled “3D Depth Point Cloud from Timing Flight of 2D Scanned Light Beam Pulses”. U.S. Pat. No. 8,543,265 B2 to Ekhaguere et al. entitled “Systems and Methods for Unmanned Aerial Vehicle Navigation”. U.S. Pat. No. 7,969,558 to Hall entitled “High Definition LiDAR System”. [Barber et al., 2008] Barber, D., Mills, J., Smith-Voysey, S., 2008. Geometric validation of ground-based mobile laser scanning system. ISPRS Journal of Photogrammetry and Remote Sensing 63˜(1), 128-141. [Burman, 2000] Burman, H., 2000. Calibration and orientation of airborne image and laser scanner data using OPS and INS PhD Dissertation, Royal Institute of Technology Department of Geodesy and Photogrammetry, Stockholm, Sweden, April 2000. [Deming, 1943] Deming, W. E., Mills, Statistical Adjustment of Data. Wiley, New-York, 1943. [Filin, 2003] Filin, S., 2003. Recovery of systematic biases in laser altimetry data using natural surfaces. Photogrammetric Engineering and Remote Sensing 69, 1235-1242. [Filin et al., 2004] Filin, S., Vosselman, G., 2004. Adjustment of airborne laser altimetry strips. In: ISPRS Congress Istanbul, Proceedings of Commission Ill. [Friess 2006] Friess, P., 2006. Toward a Rigorous Methodology for Airborne Laser Mapping. Proceedings of the International Calibration and Validation Workshop EURO COW, Castelldefels, Spain. [Habib et al., 2010] Habib, A., Bang, K., Kersting, A., Chow, J., 2010. Alternative methodologies for lidar system calibration. Remote Sensing 2˜(3), 874-907. [Hebel et al. 2012] Hebel, M. and Uwe, S., Simultaneous Calibration of ALS Systems and Alignment of Multiview LiDAR Scans of Urban Areas, IEEE Transactions on Geoscience and Remote Sensing, Vol. 50, No. 6, June 2012, pp. 2364-2379. [Kumari et al., 2011] Kumari, P., Carter, W. E., Shrestha, R. L., 2011. Adjustment of systematic errors in als data through surface matching. Advances in Space Research 47, 1851-1864. [Moniot (2009)] Moniot, R. K., Deming least-squares fits to multiple hyperplanes. In Applied Numerical Mathematics, 59(2009), pp: 135-150. [Morin et al., 2002] Morin, K., Naser El-Sheimy, 2002. Post-mission adjustment methods of airborne laser scanning data. In: FIG. XXII Int. Congress, Washington D.C., Apr. 19-26, 2002. [Schenk, 2001] Schenk, T., 2001. Modeling and analyzing systematic errors of airborne laser scanners. Tech. rep., Department of Civil and Environmental Engineering and Geodetic Science. The Ohio State University, Columbus, Ohio [Seube et al., 2012] Seube, N, Picard, A. and Rondeau, M., 2012. A simple method to recover the latency time of tactical grade IMU systems, ISPRS Journal of Photogrammetry and Remote Sensing 74 (2012) 85-89. [Skaloud 2006] Skaloud, J., 2006. Reliability of direct geo-referencing: Phase 0. Tech. rep., Euro SDR commission 1: sensors, primary data, acquisition and geo-referencing. [Skaloud and Litchi, 2006] Skaloud, J., Litchi, D., 2006. Rigorous approach to boresight self-calibration in airborne laser scanning. ISPRS Journal of Photogrammetry & remote Sensing 61, 47-59. [Skaloud 2007] Skaloud, J., Shaer, P., 2006. Towards automated LiDAR boresight self-calibration. Proc. 5th Int. Symp. Mobile Mapping Technol., May 29-31, 2007.

(102) It is to be understood that the invention has been described with reference to specific embodiments and variations to provide the features and advantages previously described and that the embodiments are susceptible of modification as will be apparent to those skilled in the art.

(103) Furthermore, it is contemplated that many alternative, common inexpensive materials can be employed to construct the basis constituent components. Accordingly, the forgoing is not to be construed in a limiting sense.

(104) The invention has been described in an illustrative manner, and it is to be understood that the terminology, which has been used is intended to be in the nature of words of description rather than of limitation.

(105) Obviously, many modifications and variations of the present invention are possible in light of the above teachings. For example, various types of UAV can be employed. It is, therefore, to be understood that within the scope of the appended claims, wherein reference numerals are merely for illustrative purposes and convenience and are not in any way limiting, the invention, which is defined by the following claims as interpreted according to the principles of patent law, including the Doctrine of Equivalents, may be practiced otherwise than is specifically described.

Notes Re. Equation Equivalents

(106) $C\left(\mathrm{bS}\right)^\mathrm{bI}=C_{\mathrm{bS}}^{\mathrm{bI}}$$C\left(F1\right)^F2=C_{F1}^{F2}$$\mathrm{Xn}\left(t\right)=\mathrm{Pn}\left(t\right)+\mathrm{C\_}\left(\mathrm{bI}\right)^n\left(t-\mathrm{dt}\right)\left(\mathrm{C\_}\left(\mathrm{bS}\right)^\mathrm{bI}\mathrm{r\_bS}\left(t\right)+\mathrm{a\_bI}\right)=X_n\left(t\right)=P_n\left(t\right)+C_{\mathrm{bI}}^n\left(t-\mathrm{dt}\right)\left(C_{\mathrm{bS}}^{\mathrm{bI}}{r}_{\mathrm{bs}}\left(t\right)+a_{\mathrm{bI}}\right)$$\mathrm{C\_}\left(\mathrm{bI}\right)^n=C_{\mathrm{bI}}^n$$\mathrm{X\_n}\left(t\right)=\mathrm{P\_n}\left(t\right)+\mathrm{C\_}\left(\mathrm{bI}\right)^n\left(t\right)\left(\mathrm{Id}-\mathrm{\Omega \_}\left(n/\mathrm{bI}\right)^\mathrm{bI}\right)\left(\mathrm{C\_}\left(\mathrm{bS}\right)^\mathrm{bI}\mathrm{r\_bS}\left(t-\mathrm{dt}\right)+\mathrm{a\_bI}\right)=X_n\left(t\right)=P_n\left(t\right)+C_{\mathrm{bI}}^n\left(t\right)\left(\mathrm{Id}-\mathrm{dt}{\Omega }_{n/\mathrm{bI}}^{\mathrm{bI}}\right)\left(C_{\mathrm{bS}}^{\mathrm{bI}}{r}_{\mathrm{bs}}\left(t\right)+a_{\mathrm{bI}}\right)$$\mathrm{\Omega \_}\left(n/\mathrm{bI}\right)^\mathrm{bI}={\Omega }_{n/\mathrm{bI}}^{\mathrm{bI}}$$\mathrm{C\_}\left(\mathrm{bI}\right)^n\left(t-\mathrm{dt}\right)=C_{\mathrm{bI}}^n\left(t-\mathrm{dt}\right)$$\mathrm{X\_n}=\left(■\left(\mathrm{x\_n}@\mathrm{y\_n}@\mathrm{z\_n}\right)\right)=\mathrm{P\_n}\left(t\right)+\mathrm{C\_}\left(\mathrm{bI}\right)^n\left(t\right)\left(\mathrm{Id}-\mathrm{\Omega \_}\left(n/\mathrm{bI}\right)^\mathrm{bI}\right)\left(\mathrm{C\_}\left(\mathrm{bS}\right)^\mathrm{bI}\mathrm{r\_bS}\left(t-\mathrm{dt}\right)+\mathrm{a\_bI}\right)=X_n=\left(\begin{array}{c}{x}_n\\ y_n\\ z_n\end{array}\right)=P_n\left(t\right)+C_{\mathrm{bI}}^n\left(t\right)\left(\mathrm{Id}-\mathrm{dt}{\Omega }_{n/\mathrm{bI}}^{\mathrm{bI}}\right)\left(C_{\mathrm{bS}}^{\mathrm{bI}}{r}_{\mathrm{bS}}\left(t\right)+a_{\mathrm{bI}}\right)$$\mathrm{P\_n}^j=P_n^j$$\mathrm{X\_n}^j=X_n^j$$\mathrm{X\_n}^j=\mathrm{P\_n}^j+C^{-}\mathrm{\_bI}^nj\left(\mathrm{C\_}\left(\mathrm{bS}\right)^\mathrm{bI}\mathrm{r\_bS}^j\left(t-\mathrm{dt}\right)\right)=X_n^j=P_n^j+{\stackrel{\_}{C}}_{{\mathrm{bI}}_j}^n\left(C_{\mathrm{bS}}^{\mathrm{bI}}{r}_{\mathrm{bS}}^j\left(t-\mathrm{dt}\right)\right)$$C^{-}\mathrm{\_bI}^nj={\stackrel{\_}{C}}_{{\mathrm{bI}}_j}^n$$\mathrm{r\_bS}^j=C^{-}\mathrm{\_bI}^nj\left(\mathrm{X\_n}^j-\mathrm{P\_n}^j\right)=r_{\mathrm{bS}}^j={\stackrel{\_}{C}}_{{\mathrm{bI}}_j}^n\left(X_n^j-P_n^j\right)$