VALIDATION OF LIDAR SURVEY DATA BY COMPARISON OF SEVERAL UNCERTAINTY MODELS

This paper introduce a new method for validating the precision of an airborne or a mobile LiDAR data set. The proposed method is based on the knowledge of an a Combined Standard Measurement Uncertainty (CSMU) model which describes LiDAR point covariance matrix and thus uncertainty ellipsoid. The model we consider includes timing errors and most importantly the incidence of the LiDAR beam. After describing the relationship between the beam incidence and other variable uncertainty (especially attitude uncertainty), we show that we can construct a CSMU model giving the covariance of each oint as a function of the relative geometry between the LiDAR beam and the point normal. The validation method we propose consist in comparing the CSMU model (predictive a priori uncertainty) t the Standard Deviation Alog the Surface Normal (SDASN), for all set of quasi planr segments of the point cloud. Whenever the a posteriori (i.e; observed by the SDASN) level of uncertainty is greater than a priori (i.e; expected) level of uncertainty, the point fails the validation test. We illustrate this approach on a dataset acquired by a Microdrones mdLiDAR1000 system.


INTRODUCTION
The validation of LiDAR data is an area of interest for most organization that order surveys and apply quality assurance processes to check the consistency of the data against contractual clauses. The most classical quality assurance procedure consists in comparing LiDAR data on Ground Control Points to check the accuracy of the georeferenced LiDAR data in a geodetic system and with respect to a given projection. This procedure is based on the identification of a geometric target or an intensity sensitive target (photogrammetry like chess mark GCP) and checks the LiDAR georeferenced points on a particular location.
However, checking GCP only gives a local accuracy information which may not be representative of the whole survey area. To overcome this difficulty, GCP sampling can be done, but multiplying GCPs significantly affects the survey cost.
Another method can be applied in checking the overlaps between survey strips, like in a strip adjustment approach. Strip adjustment purpose is to correct the platform trajectory and attitude to minimize the strip inconsistencies within the overlap area. However, strip adjustment is an active method that corrects navigation parameter without identifying areas of conflicts within the data set.
The method we propose analyses the whole data set and do not rely on the presence of overlaps between LiDAR strip. In comparing two uncertainty models, we can detect areas where abnormal errors occurred without the need for LiDAR data overlap.

Survey data validation methods
The use of Ground Control Point is a classical method to detect inconsistencies between the LiDAR data point cloud and a geo-detic point monumented either by an intensity sensitive mark (like in photogrammetry) or a geometrical target. In [1], the author defines an hexagonal retro-reflexive target made of 6 reflectors. The target can be identified in the LiDAR point cloud and an accurate estimate of its centre can be done by matching the local point cloud data on a geometric model of the target. In comparing the target centre coordinates with the actual coordinates of the GCP, the local accuracy can be obtained.
More recently, in [3], a generalization of GCP analysis by planes intersection was proposed as well as a method to select appropriate planes to get a reliable virtual GCP.
The purpose of the LiDAR survey data validation we propose is to determine if the level of observed uncertainty due to systematic errors and random errors is compatible with a model of the a priori uncertainty of the system. We shall not consider here the presence of outliers and will focus on the contribution of systematic and random error to the uncertainty level of LiDAR data.

Uncertainty models
The precision and accuracy of a point issued by a mobile LiDAR mapping system depend on the system components characteristics that can be propagated through a point georeferencing model. This model is the relationship used to generate LiDAR point clouds from raw sensors measurements, including the LiDAR itself, an Inertial Navigation System (INS) and a positioning system based on a Global Navigation Satellite System (GNSS). The derivation of points uncertainty by a standard deviation, as a function of all components and integration parameters standard deviation can be done by propagating variancecovariance matrices of each parameter submitted to uncertainty through a linearized approximation of the point geo-referencing model. This type of approach leads to a Combined Standard Measurement Uncertainty (CSMU) model, as defined in the International Vocabulary of Metrology (JCGM, 2008) by the Bureau International des Poids et Mesures (BIPM). In (JCGM, 2008), the Combined Standard Measurement Uncertainty (CSMU) of a measurement device is defined as "standard measurement uncertainty 1 that is obtained using the individual standard measurement uncertainties associated with the input quantities in a measurement model".
A simplified and decoupled model of CSMU can be found in (Baltsavias, 1999), and more comprehensive 3D models are developped in (Goulden, Hopkinson, 2010) and (Glennie, 2007). In (Shaer et al., 2007), a model including uncertainty due to the grazing angle and to the divergence of the LiDAR beam is taken into account. The energy level of the LiDAR beam is projected on its footprint, which defines an additional uncertainty in the map frame coordinates. This uncertainty level is added to the CSMU model defined by a variance propagation.
To be relevant, a CSMU model should take into account all sources of errors. Most models take into account the standard uncertainty of each sensors (position, attitude, LiDAR range and beam angle). Integration parameters (boresight angles, lever-arms) are also taken into account by using their standard deviations as long as they can be estimated and returned by calibration methods (see for instance (Skaloud, Litchi, 2006, Hebel, Uwe, 2012, Keyetieu, Seube, 2019).

POINT GEOREFERENCING MODEL
The geometry of a LiDAR system (see figure 1) can be defined by: • A Local Geodetic Frame (denoted by (n)) and a Terrestrial Reference Frame; • A Positioning Reference Point (PRP): this is the point which position is computed by the GNSS system.
• A Body frame associated to the IMU/INS (bI) and the LiDAR (bS); • The lever-arm vector a from the PRP to the optical center of the LiDAR.
• The three boresight angles that defines the mis-alignment between the IMU/INS body frame and the LiDAR body frame; In the LiDAR body frame (bS) = (X bS , Y bS , Z bS )), the LiDAR return vector is: where ρ is the returned distance and γ the beam angle.
We shall denote by C bI bS the boresight transformation matrix between the (bS) frame and the (bI) frame. This transformation can be estimated by specific techniques like the ones described in (Skaloud, Litchi, 2006), (Hebel, Uwe, 2012) or (Keyetieu, Seube, 2019). The vector C bI bS r bS is the LiDAR return coordinated in the IMU frame.
is the LiDAR body frame, and (bI = (X bI , Y bI , Z bI ))is the IMU/INS body frame.
We consider the location of the PRP, denoted by Pn. The position of the LiDAR return, denoted by Xn in the (n) frame is: We consider that time-stamping errors may produce a latency dt between the different sources of sensor information (Seube et al., 2012). And in the following model, we suppose that there is a latency between the LiDAR and the GNSS-IMU system (i.e; the INS system).

COMBINED STANDARD MEASUREMENT UNCERTAINTY MODEL
Let us denote the attitude by Ξ = ϕ, θ, ψ, the boresight by β = δϕ, δθ, δψ, the distance and beam angle by Υ = ρ, γ and lever arms by a = (ax, ay, az). In equation (2), we can separate the different terms, in order to analyze the effect of positioning, lever arms and LiDAR uncertainties. Let us denote by the term due to the ranging device and by the term due to lever arms.
From (2), we have: The covariance matrix of Xn can be written by: The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition) Figure 2. Range bias δρi due to an attitude error affecting the beam pointing direction. Here, we represented the effect of a roll error δϕ. This bias depends on the local terrain morphology.
The covariance terms of Pn, an and rn can be derived by applying the variance propagation law to linearized models. For the sake of simplicity we skip this classical step, but we mention that the Jacobian matrices involved in this computation may depend on the latency to capture timing errors.

Uncertainty due to the incidence angle of the LiDAR beam
The error equation (6) describes the ranging term covariance as a function of all variables Ξ, β, dt, Υ , but it does not relate the dependencies that may exist between the ranging system measurement parameters ρ, the terrain slope and all angular parameters. Indeed, any orientation error of the beam (due to pointing or attitude errors) may produce a range error due to the terrain geometry, as shown in figure (2).
Equation (6) describes the range error coordinated in the (n) frame due to uncertainties in all parameters. But it ignores the effect due to the LiDAR as the consequence of a beam pointing error. Indeed, whenever the attitude, the boresight or the launch angle are submitted to uncertainty, the LiDAR beam will be oriented along a different vector than the expected one. As a consequence, the actual LiDAR range measurement will be biased and this bias will depend on the incidence angle between the LiDAR beam and the terrain.
In order to take into account the LiDAR range bias due to the terrain morphology around the actual LiDAR point, we can describe the variation of the LiDAR ranging information due to other parameters variations. We shall distinguish two terms in the range error: where δρi is the term due to the beam incidence and the term δρm is the measurement error.
If we suppose that on a neighborhood of Xn, the terrain can be represented by a planar surface, then it can be shown that δρi is a linear function of δΞ, δβ, dt, δΥ, and the coefficients of the linear mapping depends on the normal vector to the LiDAR point. Note that this normal vector can be determined by PCA.

Ranging error uncertainty due to the beam incidence angle
Let us denote by un = rn ||rn|| , the unit vector of the LiDAR return, coordinated in the (bs) frame. The measured range ρm, the true rangeρ and the range uncertainty δρ are related by: For this analysis, let us suppose that δρm = 0. Indeed, this is already taken into account in the standard approach of uncertainty model. where δsn = δρi un. The covariance matrix of δsn is with un, the unit vector from the LiDAR optical center to the measured point.
As explained in section (4.1), the ranging error uncertainty depends on attitude, boresight or launch angle uncertainties. Hence, the variance of δρi can be expressed as follows: with σ 2 i , the variance of the parameter i.

Expresssion of the CSMU model
The final CSMU model is defined by the covariance matrix of the positioning, lever-arm, LiDAR return, and ranging due to the incidence angle. The final model can be defined by the following equation It takes into account the range, the orientation, the terrain incidence angle coupling and the possible presence of timing errors.

LiDAR error model
In order to supplement the CSMU model with an accurate LiDAR range error model, some static tests were made on a fixed target (wall) at different distance. For these tests, we used a SICK LiDAR, mounted on a tripod and located at various distance from a concrete wall. In analysis the time series of each beam return, we fitted a quadratic model of both range and beam angle. The resulting model is plotted in Figure (3).

The CSMU ellipsoid
We can associate to each raw data the Combined Standard Measurement Uncertainty (CSMU) ellipsoid associated to the point. To use the CSMU, we also need to define the incidence angle of the beam with the terrain (or planar cell). Thus, from the knowledge of the normal and raw data set, we can compute the CSMU ellipsoid. Definition 1 The ellipsoid of uncertainty of a point Mi is the set of points M such that: where p is a probability level and Cov is the covariance matrix representing the CSMU.

STANDARD DEVIATION ALONG THE SURFACE NORMAL
To analyze the a posteriori point cloud produced by a LiDAR system, we can evaluate the standard deviation along surface normals on a regular grid.
Let us suppose that the point cloud has been segmented in planar areas, and that each of those has been divided into regular cells (in the local segment frame) over a grid. We consider a grid of a segmented planar area and the point cloud on a particular cell of the grid.
We can estimate the standard deviation of a cell point cloud by a PCA. Indeed the smallest eigenvalue (to which corresponds an eigenvector defining the normal vector ν to the cell) is the variance of the set of orthogonal distances of the point cloud along the normal vector. Thus, this quantity that we shall denote by σν is the variance of a 1D quantity. This quantity is the orthogonal distance of a point to the plane passing through the centroid C = (x,ȳ,z) of the cell point cloud, and orthogonal to the normal vector ν.
For a point M = (x, y, z), the orthogonal distance to the cell The property of the PCA we use here is the following: Proposition 1 Consider a point cloud of points Mi. The smallest eigenvalue λmin as given by the Principal Component Analysis of the set C = {Mi}i is the variance of the orthogonal distance to the plane Π passing through the centroid C of the point cloud and orthogonal to the eigenvector − → ν associated to λmin: Var(d(Mi, Π)) = λmin (9) This result can be used to estimate the standard deviation of the point cloud along the PCA fitted plane normal, and we shall call it the Standard Deviation Along the Surface Normal (SDASN).

Definition 2
The SDASN of a cell Cj from a point cloud is where λmin(Cj) is the smallest eigenvalue of the PCA of the cell Cj.

Selection of points
The comparison of the CSMU (a priori uncertainty) and of the SDASN (a posteriori uncertainty) cannot be done in all regions of a point cloud. It should be irrelevant in vegetation areas or on objects having irregular shapes. Indeed, we need a minimum level or regularity in the point cloud (to estimate properly the local geometry a points) to conduct a fair comparison of the two types of uncertainties. Note that these two models (CSMU and SDASN) require the knowledge of the normal vector to each point. Thus we should at least be sure that a local PCA can achieve properly this task.
For this purpose, we segment the point cloud in quasi planar areas using the approach presented in (Welhan et al., 2015). This apporach is a region growing method that uses a similarity criterion between point based on the local normal and the orthogonal distance between points. In tuning these parameters we can segment the point cloud with a more or less strict planarity criterion, enabling the clustering of quasi-planar regions. Indeed, we seek for a point cloud with neither major irregularities nor vegetation. Segmenting the point cloud in relaxing the similarity parameters seems to be a good option.

Comparison of uncertainty models
For each segment found in the point cloud the estimate of local normal to point can be considered as reliable. Therefore, we can conduct the comparison of uncertainty model over each segment. To do so, we define a grid on each local segment (local plane) and we can conduct the analysis on each segment grid cell. Note that this approach enable us to achieve the uncertainty model comparison on all segment whatever the orientation of their normal is. In particular, this enable us to analyze vertical segments.
Whenever a segment of points is found, we can defined he CSMU ellipsoid for each point and compare it to the SDASN on each grid cell. The SDASN as computed on the grid cell defines the level of a posteriori uncertainty that we obtain in using the redundancy of points within the cell. The a priori uncertainty, as defined by the CSMU defines the uncertainty of isolated points.
In general, adding redundancy should decrease the level of uncertainty of a cell since all points of the cell contribute to the measurement of a piece of planar area.
Thus the CSMU uncertainty (uncertainty prediction for a single point) will be consistent with the observed SDASN (uncertainty observation of a local group of points with a planar segment) if the CSMU ellipsoid is "greater" than the SDASN level. From a geometric point of view, the extrema value of the CSMU should lie outside of the SDASN boundaries for each point.

Maximum distance of an ellipsoid from a plane
We consider here the CSMU ellipsoid E, represented by its cartesian equation, and whose center belongs to the cell plane.
To compare the CSMU ellipsoid with the SDASN, we have to consider the maximum distance from the plane to the ellipsoid. This is to say that we need to solve the following optimization problem: max where E is the CSMU ellipsoid, Π is the local planar segment to which belong the point M , ν is the local normal vector to the plane Π and Cov is the CSMU covariance matrix.

The SDASN-CSMU comparison criterion
For a given cell C | of a planar segment, we know the standard deviation of the orthogonal distances to the plan Π, which is the square root of the smallest eigenvalue of the covariance matrix of points, as given by equation (10). This quantity is the SDASN of the cell j, denoted by SDASN(C | ).
On the other hand, we know the maximum distance of the CSMU ellipsoid of all points lying on the plane Π of the cell C | .
Definition 3 We shall say that a cell Cj is validated whenever the following inequality holds: where E is the CSMU ellipsoid.
The interpretation of this definition is the following: If the CSMU ellipsoid is "larger" than the SDASN means that individual points have a greater uncertainty than the one formed by a subset of point forming the cell plan. In other words, the consistency of the underlying plan over the cell is smaller than the uncertainty of each individual points. This definition matches with the notion of cell validation we were seeking for.

NUMERICAL RESULTS
We made a test on a dataset acquired at the HIEG-VD with a Microdrones mdLiDAR1000 UAV system. This system comprises a md4-1000 UAV on which is integrated a SICK LiDAR and an Applanix APX15 INS. The CSMU model includes the model of the range uncertainty of the SiCK LiDAR. The uncertainty level of the IMU (attitudes) and of the GNSS (position of the PRP) is determined from the Applanix data sheet. The CSMU model was shown to be consistent with the Microdrones mdLIDAR1000 performance level. Indeed, the (1-σ) precision of individual points measured by this system is about We can observe a straight line of non validated point, tagged in red. All segments are represented in various levels of green 5.5cm which was shown to be consistent with the prediction of our model. Once the CSMU model is validated from test data, we can conduct the analysis on any point cloud acquired by the mdLiDAR1000.
We first segmented the data set from the HEIG-VD building which had the effect of removing all vegetation and irregular shapes. Only quasi planar segment were kept to conduct the survey data validation analysis. Figure (4) shows the point cloud segmented in planar areas. The different level of green indicate the several segments that were found by the region growing algorithm.
After griding each segment (inducing the vertical segments corresponding to walls), we applied the validation test as defined by inequality 12 In each circle point, we colorized in green validated points and in red non validated ones. On a different view (see Figure (5)), we can see that the non validated points are not lying on the planar surface. Actually, we checked that these points (forming a straight line in the horizontal plane) are from a gutter located below the ground and partially hidden by a metallic grid. It appears that some LiDAR points reached the bottom of the gutter due to the relatively low beam divergence and multiple returns capabilities of the mdLiDAR1000.
These points were not validated by our algorithm as the observed SDASN is not consistent with the a priori CSMU level. This illustrates the relevance of our approach since the points tagged as non validated can be presented to the user to further analysis or validation.
Another possible application of this validation system is the detection of systematic (lever-arm or boresight errors): Whenever a mismatch between some part of a point cloud (not horizontal to enable the effect of boresight or lever-arms errors to be observable), the SDASN level becomes not consistent with the a priori CSMU model.

CONCLUSION
We presented an new CSMU model which includes timing effect and LiDAR beams incidence angle. We studied a particular UAV LiDAR system (the mdLiDAR1000 from Microdrones) in doing static tests of the LiDAR to provide a Standard Deviation model depending on range and beam angle.
A new methodology to validate airborne LiDAR data was proposed. It is based on the comparison of the uncertainty level that we observe from the point cloud and the level of uncertainty we can expect from the system we use (i.e; the predictive CSMU model). Whenever the a priori CSMU model is lower than the observed local uncertainty model, we conclude that the system gathered point with an abnormal level of uncertainty and we tag it as non validated.
This approach enabled us to detect an inconsistency in the data set that should not have been detected otherwise. Moreover, further analysis of non validated point can reveal non regular behavior of the system or, in case of the survey realized with the mdLiDAR1000 non regular features (the underground gutter), as all point were validated.
Another possible application of this method can be the detection of the presence of a systematic error within the LiDAR system. Indeed, any calibration error can create local discrepancies between survey lines and thus make the SDASN to be greater than the CSMU. In this case, the validation tool will return unvalidated points in areas where systematic errors are observable. This could facilitate the data analysis for clients of LiDAR surveys in order to validate the relative consistency of point cloud delivered by contractors.