Methods for Estimating Relative Accuracy of ALS Data 1

The data acquired by airborne laser scanning is a set of measurement points with certain XYZ coordinates (the so-called cloud of points), located at places where laser pulses sent by the scanning system are refl ected. The coordinates of these points are determined based on the distance measured from the scanner system to the ground along the predetermined laser beam, and the scanning angle at which the beam is sent. The position and orientation of the scanner during the measurement is determined with the use of the GPS/INS system integrated with the scanner [7]. The coverage of the developed area with data from the ALS is done by executing a fl ight in multiple parallel strips that partly overlap, thus ensuring the continuity of the development. The accuracy of the designated coordinates of ALS points depends on the accuracy of observation (distance, scan angle, position and orientation) and on the mathematical model used to determine the coordinates of these points [2]. In connection with the occurrence of errors in the measurement devices and in the defi nition of the mathematical model, there exist some discrepancy in the data from the adjacent strips in overlapping area. Therefore, in order to improve the position of the measuring points, we use ALS data alignment methods by basing on the adjustment of data from adjacent strips and their best integration in the exterior coordinate system. The use of ALS data for the purposes of subsequent developments, such as the DTM (Digital Terrain Model) and DSM (Digital Surface Model) requires the determination of their accuracy. By measuring the diff erences between the adjacent strips, we can estimate the relative accuracy of the ALS, and by sett ing the position of the examined control objects in the exterior coordinate system independently, we can also estimate their absolute accuracy. However, the specifi c nature of data


Introduction
The data acquired by airborne laser scanning is a set of measurement points with certain XYZ coordinates (the so-called cloud of points), located at places where laser pulses sent by the scanning system are refl ected.The coordinates of these points are determined based on the distance measured from the scanner system to the ground along the predetermined laser beam, and the scanning angle at which the beam is sent.The position and orientation of the scanner during the measurement is determined with the use of the GPS/INS system integrated with the scanner [7].The coverage of the developed area with data from the ALS is done by executing a fl ight in multiple parallel strips that partly overlap, thus ensuring the continuity of the development.The accuracy of the designated coordinates of ALS points depends on the accuracy of observation (distance, scan angle, position and orientation) and on the mathematical model used to determine the coordinates of these points [2].In connection with the occurrence of errors in the measurement devices and in the defi nition of the mathematical model, there exist some discrepancy in the data from the adjacent strips in overlapping area.Therefore, in order to improve the position of the measuring points, we use ALS data alignment methods by basing on the adjustment of data from adjacent strips and their best integration in the exterior coordinate system.
The use of ALS data for the purposes of subsequent developments, such as the DTM (Digital Terrain Model) and DSM (Digital Surface Model) requires the determination of their accuracy.By measuring the diff erences between the adjacent strips, we can estimate the relative accuracy of the ALS, and by sett ing the position of the examined control objects in the exterior coordinate system independently, we can also estimate their absolute accuracy.However, the specifi c nature of data acquisition (collection of distributed measurement points) prevents a direct comparison of the same points on data from two adjacent strips.It is therefore necessary to use additional objects obtained from the cloud of points.This has resulted in the development of several independent methods of estimating the accuracy of ALS, based on a selection of diff erent types of control elements.
Chapters 2, 3 and 4 present the three main approaches, described in the literature, of estimating the relative accuracy of the LiDAR data while Chapter 5 presents a method developed and used in Poland during the ALS data control in the ISOK project -IT System of the Country's Protection against extreme hazards [9].

Method of Matching DEM Models
In order to estimate the relative accuracy of ALS data, it is necessary to determine the discrepancy between independent measurements of the same objects, located in overlapping areas between strips.In the described method, DEM (Digital Elevation Model) [8] is a model of surface with elevation set out in the regular rectangular grid, resulting from the interpolation of a cloud of ALS points.The shifts between the strips are determined by comparing independently generated DEM models from the two ALS point clouds originating from two adjacent strips and located in the overlap area between these strips.In this case, therefore, the single control object is the DEM data window of a fi xed size.Based on the applied matching algorithm we determine dX, dY and dZ shift between the corresponding portions of two terrain models.Initially we determine the approximate parameters of the value, then by using DEM windows displacement in X and Y we determine dZ altitude deviation between the models.As the point cloud derived from two adjacent strips of ALS data is acquired at diff erent angles, it is characterized by the occurrence of variously located occlusion areas.At the moment of the DEM model creation, these areas are burdened with an additional error resulting from the interpolation of data, and they should be excluded from the analysis.In view of the above assumptions, the altitude deviations between DEM models in the presented method are determined solely on the basis of the analysis of fl at areas.
When it comes to the method of locating DEM areas appropriate for altitude comparison, in [18] it has been proposed to extend the above method with the use of the so-called smoothness-mask.In this method, the analysis does not concern occlusion areas or those portions of DEM models for which there is a sudden change in altitude (as in the case of areas covered by vegetation, or on the edge of building outlines).Thus, we avoid falsifying the results of the analysis, since large altitude diff erences between terrain models in these areas are primarily the result of the interpolation of regular GRIDs from two independent sets of dispersed XYZ points.In the proposed method, the individual height of DEM grid points are determined by the moving plane interpolation method, consisting in seamlessly integrating the plane into a set of ALS points that are closest to a given grid point for which the altitude is determined.During this process, each DEM grid point is assigned additional parameters, such as σ d , s and ε to form a smoothness-mask in the following stage.The value of σ d (standard deviation) determines the accuracy of the height interpolation of a given DEM grid point, thus enabling the detection of areas with sudden changes in altitude.The parameter s determines the distance between the point of the GRID and the nearest ALS point, while ε represents the planar distance between the LiDAR points' centre of gravity used to determine the said plane, and a given point of DEM grid.Increased values of s and ε parameters symbolize occurrence in a given strip of an area for which there are no ALS measurement points (e.g.occurrence of occlusion areas).On the basis of the above parameter values, for each strip a smoothness-mask is formed as a GRID of points, corresponding to the previously created grid of the DEM model.This mask determines whether the DEM model grid point represents an area considered to be "smooth", i.e. which satisfi es the condition (1): where: DEM model matching is carried out with the use of the LSM method (Least Squares Matching).This method gives the best results, assuming that the DEM resolution is similar to the original ALS data density [1].During the matching process, the DEM model windows are moved in the directions X and Y, until the squared sum of dZ between these models is minimised.The analysis of altitude diff erences takes into account only those cells of the DEM model for which there was an occurrence of "smooth" areas on both of the adjacent strips.In this way, the shifts in dX, dY and dZ at specifi c control objects (DEM fragments) represent the accuracy of strip adjustment of the ALS data.
Given the main problem of estimating the relative accuracy of the ALS data, which consists in comparing the two clouds of irregularly spaced XYZ points, to determine the off sets between data from adjacent strips we use indirectly derived DEM models.The described method is therefore fraught with an additional error resulting from the interpolation of these models based on the source data that are ALS point clouds.

Method of Matching Points and TIN Model
Another method of estimating the relative accuracy of the ALS data is based on a comparison of dispersed cloud points in one strip with a TIN surfaces (Triangulated Irregular Network) stretched on a cloud of points from the adjacent strip [12].In this approach, the control object shall be deemed the two pieces of clouds of points, selected from the adjacent strips in overlapping area.These two sets of data are moved relative to each other in three directions, so that the sum of squared diff erences in their altitudes is as small as possible.The diff erence in altitude is determined between the source point of ALS from one strip and the point obtained from bilinear interpolation of the TIN triangle(which corresponds to a given location) formed from the point of the adjacent strip.To estimate the parameters of shifts in X, Y and Z we use the least squares adjustment (LSA) method, and the observation equations are formed independently for source points from both strips in the selected portion of data.Similarly to the method based on the comparison of DEM models, the inclusion of occlusion areas and areas with a sudden change in altitude in the analysis can cause the occurrence of coarse errors during data matching, and consequently lead to the erroneous determination of shifts between them.In order to eliminate this eff ect, when determining the diff erences between the objects we do not consider the points occurring in this type of areas, using the method for matching LSM together with some strong estimation techniques described in detail in [11].
The determination of shifts between the adjacent ALS strips in three directions, using the described method of matching points with the TIN model is justifi ed only when analysing the area that meets the criterion of adequate diversifi cation in terms of altitude.The selected portion of data should have surfaces located in three different non-coplanar directions.Therefore, in the case of fl at areas, with insuffi cient drops of terrain, where there are no buildings that would allow for proper analysis, we use constraints for data matching.Then, it is reasonable to designate only the altitude shifts between corresponding sets of points.In these areas, it may be helpful to have information about the intensity of the signal refl ection, obtained during the measurement of points with ALS.Since this information corresponds exactly to the location of the point, it can be used to determine an additional horizontal shift between the adjacent strips.Of course, the proper span of acquired refl ectance intensity in the analysed area is then necessary, which boils down to the presence of well-identifi able objects in the intensity image.Figure 1 shows exemplary ALS data, illustrating an area of slight diff erences in altitude in comparison with its intensity image, which allows for the interpretation of situational information, such as road edges or road strips.
The method of matching data from two adjacent strips with the use of the LSM method and allowing for additional information about the intensity of signal refl ection for the ALS data is described in [13].
A similar method of evaluating the accuracy of relative georeference of ALS data is presented in [4] and [5].It is also based on matching two pieces of point clouds, where one strip is represented by XYZ source points, while the second strip is represented by the TIN surfaces created on the basis of its ALS points.This approach determines transformation parameters between point q i of strip S 2 and the corresponding point q i ′ on the triangular patch S p , which has been formed based on the points from strip S 1 , as shown in Figure 2.With enough data on the pairs of points and TIN patches, we can determine transformation parameters between successive strips.This approach takes into account three displacement parameters between the data of the adjacent strips: X T , Y T , Z T , and a rotation angle ϕ -around the fl ight direction.The patches of TIN in strip S 1 , corresponding to given points of strip S 2 , are determined by an iterative method of matching data similar to the ICP (Iterative Closest Point).Initially, the triangular patch S p with the shortest normal distance to point q i ′, after allowing for the approximate transformation parameters, is considered the patch corresponding to point q i .For points in strips specifi ed in the same coordinate system, one can assume initial values of the transformation parameters (X T , Y T , Z T and ϕ) equal to zero.On the basis of this a) b) pre-defi ned matching of ALS points and TIN patches, we determine the correct transformation parameters between points q i and q i ′.Then, using the above calculations, we refi ne the matching of corresponding points with the patches of the TIN model.This process is repeated until there is no change in the values of the estimated transformation parameters between successive iteration.However, one should remember to assume the permitt ed normal distance between an ALS point and a given triangular patch of the TIN, in order to exclude the areas with sudden changes in altitude, which could skew the results of the analysis.This method is primarily designed to assess the relative accuracy of ALS data, in order to determine occurring in this data random and systematic errors.The resulting transformation parameters can be further used to remove the detected shifts between the strips in the analyzed point cloud.
In [14] and [15] presented diff erent approach based on the acquisition of surface from ALS points.In the overlapping area, the points from one strip are segmented (using method specifi ed in [3]) in order to obtain planar planes.In the next step, the points of the adjacent strip belonging to the identifi ed segments are selected.The discrepancies are measured between the segments from one strip and the surfaces obtained by fi tt ing them to selected points from other strip.The authors [14] and [15] proposes the measurement of the diff erences between the two surfaces and not between the surfaces and the points as in previously described approaches.

Method of Matching Linear Elements
A separate approach to estimating the discrepancies between data in adjacent ALS strips is the use of the linear objects that are obtained on their basis.This method is mainly based on the edges designated from the intersection of two adjacent planes.The positions of these planes are represented by multiple ALS points so that their equations can be determined with high accuracy [6].In the fi rst phase, we use the ALS data segmentation methods to make it possible to automatically designate this type of linear objects.For this purpose, from the data of one strip is created a TIN surface, based on the analysis which sets of points are identifi ed, that have a degree of spatial coherence.Thus, as a result of the segmentation we obtain sets of points that represent a given plane, and by determining its parameters and identifying the adjacent planes, we can easily determine the equation of the line of intersection of two such planes.The most frequently used surfaces in this approach are roofs of buildings, and the lines of the intersection of these planes constitute the ridges of the roofs.
With the given equation of the line drawn from the intersection of two planes, we must still specify its end points, i.e. choose one fragment of the line which actually corresponds to the physical ridge of the roof.To this end, the points representing the analysed planes in a predetermined buff er are projected on the identifi ed line of intersection and the extreme points that are successively selected from it are used to designate the end point of the line [10], as shown in Figure 3.
The segmentation and determination of the intersections of the lines is performed independently for data from each of the ALS strips, followed by the identifi cation the lines corresponding to each other.Because of the way of acquiring data with the use of the airborne laser scanning method(diff erent scanning angle, positioning of planes, occurrence of occlusion areas etc.), one should take into account the relative position of the detected lines to determine whether they actually correspond to each other.For this purpose, we analyze their mutual distance along the normal vector, their mutual parallelism and mutual percent coverage [6].
Having obtained the corresponding linear objects on the basis of data acquired from two adjacent ALS strips, we can apply a number of approaches to measuring diff erences between them.In one method, the matching points are the end points of such lines [4,10,20].In this type of approach one should bear in mind that these points do not necessarily correspond to each other, and they can be shifted relative to each other along the line from which they originate.Because of that, when determining transformation parameters between such observations, one should extend the variance for such a point from the direction of the orientation of the line which it represents [4].Finally, the elements of transformation (shift and rotation) between successive strips are determined by using the least squares adjustment (LSA) method.
Another solution is based on a comparison of points from the intersection of linear objects.In [17] it was suggested to compare 2D and 3D points resulting from the intersection of the lines representing the ridges of roofs.2D points, on the basis of which we determine only horizontal shifts, arise from the intersection of the ridges of buildings located in short distances from each other.3D points, however, used to determine shifts in the three directions: X, Y and Z, are determined by the intersection of the ridge of the roof of one of the buildings and the surface of a roof of a neighbouring building.An example of a 2D point and two 3D points obtained from the two adjacent buildings is shown in Figure 4.

Fig. 3. Extraction of line end point
Source: [20] When designating this type of points independently based on the data from two adjacent strips, we calculate the diff erences in coordinates between them.The horizontal and vertical shifts obtained in this way, resulting from the shift and rotation between the corresponding strips, are used to estimate the relative accuracy of the ALS data.

Method of Matching Roof Ridge Lines and Altitude Grids
This method consists in obtaining information from the comparison of both linear elements and planar, horizontal surfaces on the data of the adjacent LiDAR strips.It has been developed in order to estimate the relative accuracy of the data derived from ALS, obtained within the framework of the ISOK project, covering almost 70% of Poland.
In the suggested method, the control objects are both linear elements(mainly ridges of simple pitched roofs), used for the estimation of horizontal shifts, as well as height points measured at regular intervals, forming 6x6 points with 1m grid size (located on fl at, paved surfaces with a slight decrease in the area), used to estimate vertical shifts.In this way we can obtain complete information about both the horizontal and vertical discrepancies between the adjacent ALS strips.The control object in this case is a pair of possibly mutually perpendicular buildings, as the corresponding ridges only give us the possibility to obtain information about the shift in the direction perpendicular to them [19].In the local coordinate system of these objects, we examine the horizontal shifts between pairs of the corresponding lines, along the direction of their normal vector.Taking into account the possibility of deviation of the ridges relative to each other, the distance measurements for each pair of lines are carried out independently at three locations, as shown in Figure 5.
With the obtained average divergence Δx and Δy along two perpendicular directions, we determine the resultant horizontal displacement for a given control object, according to equation (2):  Source: [16] In less urbanized areas, the control objects may be linear elements derived from the intensity images.These types of objects constitute details that are well identifiable in the intensity image, e.g.road stripes, roadway edges or farmland borders.However, one should remember to choose two of these elements, which, like in the case of ridges, are positioned perpendicular to each other.
In the proposed method, the measurement of vertical diff erences between the adjacent strips is determined on fl at, smooth surfaces, such as paved roads or parking lots.For a set of grid points (5 m × 5 m) of a known XY coordinates, we determine the elevation by interpolation of ALS points, independently for the data obtained from the two adjacent strips.Based on the values of ΔZ displacement between the successive grid points, we calculate the average vertical error on the examined control object.
In the proposed method to determine horizontal and vertical shifts between strips we used separate control objects.It is possible, of course, to obtain this information using only linear elements, i.e. roof ridges (as in the method described in Chapter 4, matching linear elements).However, these objects are often not available in less urbanized areas, where the independent use of information from the image of refl ectance intensity and the measurement of diff erences in altitude gives an additional opportunity to assess the relative accuracy of the ALS data.Also, due to the nature of the subsequent use of the data acquired within the ISOK project (mainly fl ood protection), estimating the accuracy of altitude constitutes a key element.Therefore, it appears justifi able to locate objects for altitude control on fl at, horizontal surfaces, where the measurement of vertical discrepancies is not aff ected by a horizontal shift between the strips.
The average error in the control objects within the ISOK project in urban areas (for which the acquired ALS point clouds had a density of 12 points/m 2 ) should not exceed 0.60 m horizontally and 0.15 m vertically.These areas represent approximately 7% (13.7 thousand km 2 ) of the examined land.In turn, for other data with the acquired point cloud of minimum 4 points/m 2 the permissible average horizontal error is 0.75 m and the vertical error is 0.22 m.In addition, it is required that diff erences in 68% of the measured control objects within a single, consistently examined area, does not exceed the specifi ed limit values, in 95% of the control points double their value, while for none of the control objects can they be more than triple the average value of the maximum permissible errors.All data developed within the ISOK project and made available by the Centre of Geodetic and Cartographic Documentation (CODGiK) meets the above requirements.

Summary
Due to the rapid development of ALS, it is necessary to develop eff ective methods of controlling the quality of the data obtained with the use of this technology.It follows from the presented review of the literature that so far diff erent methods have been used to estimate the accuracy of the ALS data, and each approach has its advantages and disadvantages.The described methods have been developed primarily to determine the relative accuracy of the ALS data, but by determining the independent location of selected control objects in exterior coordinate system, one can also estimate their absolute accuracy.The assessment of the usefulness of a particular solution is highly infl uenced by the type of coverage of the analysed area.In the case of heavily urbanized areas, where there are a lot of occlusion areas and areas with sudden changes in altitude, roof ridge lines and roof patches constitute good control objects.However, in less urbanized areas, a bett er approach is to choose the method based on the comparison of points and TIN surfaces.As for a rough and rapid determination of the ALS data accuracy, it might be enough to compare DEM models generated on their basis, even with the application of a limitation regarding the determination of just vertical shifts between strips.
Within the ISOK project, where the data is located throughout Poland and includes both heavily urbanized areas and areas diffi cult to access (mountains), we have developed a method that uses linear objects as well as fl at surfaces to obtain information about horizontal and vertical accuracy of the received data.This method allows independent comparison of ridges of buildings, information on the intensity of the refl ected signal or fl at surfaces.It has proved to be the optimal solution for such a large project because it gives the possibility to adjust the selected type of the control object to the coverage of land in the analysed data.
the standard deviation, ε -the planar distance between point of DEM and the ALS points' centre of gravity, s -the distance between the point of DEM and the nearest ALS point, value of, and s.

Fig. 5 .
Fig. 5. Schematic location of measured distances d between pairs of corresponding ridges