REGISTRATION OF UAV DATA AND ALS DATA USING POINT TO DEM DISTANCES FOR BATHYMETRIC CHANGE DETECTION

This paper shows a method to register point clouds from images of UAV-mounted airborne cameras as well as airborne laser scanner data. The focus is a general technique which does rely neither on linear or planar structures nor on the point cloud density. Therefore, the proposed approach is also suitable for rural areas and water bodies captured via different sensor configurations. This approach is based on a regular 2.5D grid generated from the segmented ground points of the 3D point cloud. It is assumed that initial values for the registration are already estimated, e.g. by measured exterior orientation parameters with the UAV mounted GNSS and IMU. These initial parameters are finely tuned by minimizing the distances between the 3D points of a target point cloud to the generated grid of the source point cloud in an iteration process. To eliminate outliers (e.g., vegetation points) a threshold for the distances is defined dynamically at each iteration step, which filters ground points during the registration. The achieved accuracy of the registration is up to 0.4 m in translation and up to 0.3 degrees in rotation, by using a raster size of the DEM of 2 m. Considering the ground sampling distance of the airborne data which is up to 0.4 m between the scan lines, this result is comparable to the result achieved by an ICP algorithm, but the proposed approach does not rely on point densities and is therefore able to solve registrations where the ICP have difficulties.


INTRODUCTION
Most registration approaches are based either on human made linear structures, planar structures or point correspondences.But these have the disadvantages that human made structures mainly exist in urban areas and direct point registrations, i.e.ICP (Iterative Closest Point), are limited to point clouds with almost the same point density and distribution.Therefore, this paper focusses on a more general technique which does rely neither on linear or planar structures nor the point cloud density.The main application for this registration will be the change detection of river ground points.By looking on river areas it is obvious that the river mostly crosses rural areas.So the registration of bathymetric measurements in river areas must mostly deal with rural areas.Furthermore, for the change detection information is needed on the accuracy of the transformed coordinates so that it is possible to distinguish real geometric changes from data noise or registration errors.Therefore, an accurate registration is needed which also considers measurement errors and take them into account for determining the variance of the resulting transformation parameters.
The proposed method is based on a regular 2.5D grid generated from the segmented ground points of the 3D point cloud.It is assumed that initial values for the registration are already estimated, e.g. by measurements of the UAV mounted GNSS and IMU sensors.Figure 1 shows an overview of the steps for registration from the two input point clouds to the aligned point clouds.
The input source data is an airborne laser scanner point cloud where the ground points are segmented based on a voxel grid and the pulse class of the laser return.Like descripted in a previous Figure 1.Sheme of registration: The source point cloud is generalised with a regular raster DEM and the target Point cloud is aligned using a minimisation of point to raster distances publication of the authors (Boerner et al., 2017), it is assumed that ground points correspond to the lowest voxels inside the grid including mainly last pulses.The 2.5D grid is generated from the ground points with a least squares fit of a regular raster with smoothing criteria.This grid has the advantage that each raster cell holds a reconstruction accuracy and that the point to grid distance can be calculated very fast.The registration is solved with a least squares approach based on weighted point to raster distances by minimizing this point to raster distance for the whole point cloud.The reconstruction accuracy of each raster cell is used to filter false observations by not considering observations which belong to raster cells with low accuracy.This is checked by setting a threshold on the variance of the heights of the points inside a raster cell.Furthermore, this variance of the raster cell height is also used to adjust the observation weights.Therefore, the statistical model considers not only the accuracy of the 3D coordinates of each point but also the reconstruction accuracy of the raster.This reconstruction accuracy is low (high variance) in areas with less 3D ground points, e.g.vegetation.Outliers are filtered by using the point to raster distance.Due to the ground point based raster, such outliers with a high point to raster distance are mainly non-ground points of the target point cloud.Therefore, the outlier removal also filters the ground points in the photogrammetric point cloud.
The calculation of the transformation parameters can be solved online, i.e. during the UAV flight, by implementing the least squares approach in an adaptive manner.Such registration can be used to enhance the 3D points derived from the aerial laser scanner acquisition with a higher resolution.Since this can be done online, the UAV pilot is able to see the preview result in the field and can decide if the point cloud derived from the captured images is of sufficient density and accuracy.Another advantage is, since the landscape itself is used for the registration, that the approach is independent of urban structures like buildings.This is shown on an example area of a river basin in southern Bavaria.It is a small side arm of the river Mangfall which can be considered as stagnant water due to the very low water velocity.Therefore, morphodynamic characteristics in the area of the river basin are not considered.Possible influences of morphodynamic characteristics are minimized by considering the whole landscape.The chosen area shows no buildings, only vegetation and bare ground.
The landscape shows a locally flat area but there are some hills and structures of the river basin which makes it possible to use the landscape for registration.The achieved accuracy of the registration is up to 0.4 m in translation and lower than 0.3 degrees in rotation, by using a raster size of the DEM of 2 m.Considering the ground sampling distance of the airborne data which is up to 0.4 m between the scan lines, this result is near to the result which can be achieved with an ICP algorithm but the proposed approach is able to solve the described registration while ICP is not.
The registration of UAV-based images and airborne laser scans should be used in the future for topobathymetric change detection.Since the approach is not limited to RGB-cameras, a NIRcamera can also be used which simplifies the detection of vegetation areas and water bodies for refraction correction.The point to grid distance can also be used to determine true geometric changes and differences caused by measurement, reconstruction or registration errors, this will be the key advantage of the proposed method for future change detection.
Section 2 shows a short overview of the state of the art in the context of point cloud registration.The proposed approach is described in detail in Section 3 and the used data is shown in Section 4. Results are shown in Section 5 and the summary is given in section 6.

RELATED WORK
Marker less registration for aligning different scans is a currently important research topic.The topic gets even more important by thinking of multi sensor registration.Methods in this field are generally based on geometric characteristics which are used for alignment.These geometric features used in the geometric approaches can either be point based features or features extracted from other primitives (e.g., planes).
The point-based approaches aim for finding corresponding point pairs in the different scans.The Iterative Closest Point (ICP) can be seen as the major algorithm in this field.ICP and its variants works in overlapping areas by minimizing point-to-point distances (Besl and McKay, 1992;Habib et al., 2010;Al-Durgham and Habib, 2013).This ICP-like methods have been proven to be effective in terms of accuracy but the matching of corresponding points needs good initial alignment and a time consuming matching process.To reduce the searching space for correspondences, there are also key point based methods which adopt classic key point descriptors of the photogrammetry in 3D.There are SIFT based (Böhm and Becker, 2007;Weinmann et al., 2011), DoG based (Theiler et al., 2014) or FPFH key point descriptors (Weber et al., 2015).Another strategy in the search for corresponding points is to look for intersection points which are considered as semantic points (Theiler and Schindler, 2012;Yang et al., 2016;Ge, 2017).In general all these methods are able to create a good alignment of different point clouds but get into troubles when dealing with different densities, noisy data or large scale datasets.Also a slightly different level of detail in the two point clouds can affect the alignment in a negative way.
There are also high level geometric features which can theoretically increase the robustness of the matching of feature pairs.For representative features lines (Habib et al., 2005), surfaces (Ge and Wunderlich, 2016), or planes (Xiao et al., 2013) can be found.Line features can be further splittet in several representations of the line which are represent in a large number of investigations for registration problems.Some examples of line representations are intersection lines of neighbouring planes (Stamos and Leordeanu, 2003), spatial curves (Yang and Zang, 2014) or 3D straight-lines (Habib et al., 2005).For the alignment of plane features, there are approaches searching for plane correspondences (Dold and Brenner, 2006;Von Hansen, 2006;Xiao et al., 2012) as well.But the extraction of planar surfaces by region growing or model-fitting algorithms can be time consuming and therefore limits the performance of the registration (Wang et al., 2016;Xu et al., 2017).In the field of segmentation plenty of studies shows the advantage of using voxel structures, by constructing a rasterized representation of volumetric elements (voxel) to simplify the dataset.These voxel structures can also overcome the uneven point density.Registration approaches which use such voxel structures can be found in the field of coarse registration by using EGI features (Wang et al., 2016) or planar surfaces extracted from the raw point cloud (Xu et al., 2017).
By comparing these two groups of registration approaches, it can be said that point based approaches are mostly located in the field of fine registration because of high accuracy (especially the ICP is used here).Beneath the point based approaches, the high level geometric primitives can overcome problems in different point densities but the majority of them focuses on residential or other urban areas.So the high level primitives will get into trouble when dealing with rural areas like in forests or bathymetric applications.Therefore, the proposed approach aims for a point based registration using the landscape itself as geometric primitive.This will combine the high accuracy of a point based approach and the stability of a geometric generalisation against noise and different density.

METHOD
The proposed method can be seen as an adoption of the ICP.But with the main difference, that the proposed method doesn't assume point correspondences, which are hard to find considering multi sensor data.Instead of point distances, the proposed method minimizes distances of 3D points to 2.5D grids.Therefore, only one point density takes account to the final registration accuracy, instead of the point density of each point cloud.The whole process is summarised in Algorithm 1 and explained in detail in the Sections 3.1 and 3.2.Beneath the iterative minimizing approach (Section 3.1 ) there is also a regularisation included (Section 3.2) which efficiently filters outliers in the target point cloud. Data

The minimizing problem
The proposed approach aims for the minimisation of the distances of a single 3D point to a source 2.5D grid.This grid can be described in a mathematical way as S(x, y) which should be a continuous function.One point on the grid is given by the 3D coordinates x; y; S(x, y) T .The continuity of the function S(x, y) is reached by a bilinear interpolation of the point inside a single grid cell.
with: s = x − xi : difference of the x coordinate of the current point to the x coordinate of the anchor edge(xi) on the grid inside the range [0,1] t = y − yj : analogue to s with the y coordinates of the current point z = heights of the grid cells edge points (i.e. the z coordinate of each edge of the current grid cell) i,j = coordinates of the anchor edge of the current grid cell, corresponding to the given x,y values The height of each grid cell are estimated like described in Förstner and Wrobel (2016) by the optimisation of with the regularisation terms where the transformed coordinates are generated by the use of Equation 6with the use of p n .The transformation parameters are estimated by the optimisation of where σn is the variance of the single observation.
The optimisation is solved by the least squares method which gives an analytical formula for the parameters: with: x = vector of the parameters with σn;j as the variance of the j-th coordinate of point n.The variance of the raster cells height to the point n is given as σz;n.This value can be calculated based on the stochastical output of the least squares reconstruction of the source grid.Therefore, the stochastical model for the observation n considers reconstructed grid cells with high variance in a lower weight for the observation itself.It is also noteworthy, that if the variances of 3D point coordinates are not taken into account, the observation weight fully consist of the reconstruction accuracy.
Since the shown functional model is a non-linear model, there are initial values needed to solve the derivatives.These initial values are fine tuned in an iterative way by calculating the parameter update in a single iteration.The iterative process is aborted if either a maximum number of iterations is reached or if the parameter update is small enough to be considered as convergence.

Outlier removal
The proposed approach considers a grid of ground points as source.
Therefore, an outlier removal can also be seen as a binary classification which returns the ground and non ground points in the target point cloud.This outlier removal is realised by a dynamic threshold during the iterative registration process.Therefore only points whose distance to the grid is below the dynamic threshold are considered as observations for the registration.The distance threshold is determined with a histogram based approach.
Each point in the target point cloud is considered to build a histogram of distances to the source grid within every iteration.To define the distance threshold which divides ground from non ground points, some conditions are assumed.First the initial values must be good enough to create the highest peak of the ground points in the histogram (e.g.there are no rotations of 90 degrees or This outlier elimination is important for the registration approach, because of the the fragility of the least squares method in case of outliers.Furthermore, it generates a binary classification to ground and non ground as a further output of the whole registration process.

DATA
For testing, an area of the river Mangfall in Bavaria, Germany was chosen.The chosen area is located in a side arm of the river which can be considered as an area of standing water.Here are some beaver dams which were the main focus of a past research project and were captured in 3D with some UAV flights in the past.These beaver dams are structures which are to small to be captured by airborne laser scanning and therefore, exemplary show the need of spatial enrichment of airborne data.
The airborne LiDAR Data was captured with an Riegl VQ880 Li-DAR sensor which uses a full waveform processing and a green wavelength for the LiDAR measurement.Therefore, the data can distinguish between ground and vegetation or water points.For the UAV flights a Falcon 8 UAV was used with a Sony Nex 7 camera which captures 24 MPix images.Figure 3 shows an overview of the 3D data in the chosen area.The beaver dam is located in the marked area of Figure 3b and shown in a closer look in Figure 4.
The data used for evaluating the proposed method have differences in the attributes which must be considered as border conditions of the matching approach.At the one hand the radiometrics can't be considered as similar since the one data is derived from RGB images and the other one from the measurements of a Li-DAR system.Therefore, the intensity values differs a lot.On the other hand, the geometric information is slightly different as well.Where are differences in the point density of the 3D data and also in the level of detail, as mentioned above.The ground sampling distance of the LiDAR data is 1 to 4 dm and the ground sampling distance of the potogrammetric data is 0.6 dm.Therefore, the photogrammetric point cloud captures more detail but also includes more noise.Furthermore, the landscape itself can be considered as a locally flat landscape in this region.Therefore, the test side shows several border conditions which the proposed method must handle.The method must be able to register multi sensor data by using a geometric approach, since there are no intensity similarities.The registration have to be robust against noise, differences in point densities as well as slightly different level of details.

RESULTS
The proposed approach is evaluated to get knowledge about two things: first the accuracy of transformation parameters in the shown area and second the influence of initial values.As mentioned in Section 3, the proposed approach needs good initial values to optimise the transformation.What are good initial values is therefore also tested in this section.The LiDAR point cloud is georeferenced by direct georeferencing using the GNSS and IMU sensor in the air plane.The photogrammetric point cloud from the UAV taken images is georeferenced using ground control points.Therefore, these two point clouds can be considered as registered which is used for comparison between the proposed automatic approach and a control point based approach.As a second comparison of the accuracy of the transformation parameters, the same LiDAR point cloud was registered to the source grid with an initial translation and rotation different to the identity transformation (transform one point to himself).Since it is the same data, the transformation parameters should be the one for the identity transformation.Differences of the resulting parameters can therefore be considered as accuracy differences.Furthermore, the evaluation of the influence of initial values is done by changing the initial values to bigger differences from the identity transformation.
For the evaluation of the transformation accuracy a slightly wrong position up to 5 m is considered, which can be considered as the accuracy of the UAVs GNSS sensor.The rotation error was considered as up to 3 degrees, which also can be considered as the UAVs IMU accuracy.This initial alignment is visualised by the overlapping point clouds in Figure 5.The misalignment can be seen in the marked area of the water body where radiometric break lines are visible in the LiDAR data.With these initial parameters the output of the proposed automatic approach can be compared to a control point based approach.Figure 7 shows visually results of the proposed approach with the given initial alignment.In comparison Figure 6 shows visual results of ICP (the implementation in cloud compare 1 was used) with the same initial parameters.A visual inspection of the radiometric break lines of the marked area in the LiDAR data shows that, the ICP has difficulties to handle this registration problem Figure 6.Result of ICP with given initial alignment (see Figure 5) the marked areas show the comparison between the initial alignment Figure 7. Result of the proposed approach with given initial alignment (see Figure 5) the marked areas show the comparison between the initial alignment while the proposed approach creates a suitable alignment.One possible reason for the difficulties of ICP can be the varying point density.By switching the target and source cloud to the source being the denser photogrammetric point cloud and the target being the LiDAR cloud, the alignment gets even worse which show an existing influence of point density in the ICP.Another possible reason can be the slightly different level of detail between the smooth LiDAR point cloud and the high resolution photogrammetric one.This difference makes it difficult to find one-to-one point correspondences.The proposed approach optimises distances to a continuous surface which can be calculated to each target point and enables an average between the two level of details.
A numerical evaluation results in the following numbers: the translation can be determined with an accuracy of about 4 dm and the rotation accuracy can be determined below 0.3 degrees.Considering a used grid size of 2 m and the ground sampling distance of the LiDAR cloud of about 5 dm the translation is inside the ground sampling distance.As a visual result of the outlier removal Figure 8 shows the raw input as well as the filtered output for the photogrammetric point cloud.The marked areas show successful filtered vegetation points which are determined as outliers.The final threshold is below 4 dm, which efficiently filters fallen trees out of the ground layer but includes little grass hills or structures in the river bed.
The numerical evaluation is also summarised in Table 1.A comparison between an ICP approach shows that the ICP gets into troubles with the given dataset.The difference of the transformation parameters between the ICP and the proposed approach reaches the maximal values of 22 m in translation and 3 degrees in rotation.As another visual result Figure 9 shows the aligned region of the beaver dam.The spatial enrichment of the photogrammetric point cloud can be seen and the scan lines can be located on the beaver dam.

Method
The evaluation of the initial parameters shows that suitable transformation results are reachable with a wrong position initialisation up to -5 to +5 m and a wrong rotation initialisation up to -8 to 8 degrees.These values are determined considering the opposite (translation or rotation) parameters as very close to the identity.If both parameters are considered as wrong, the borders are slightly different.In sum the translation can be inside -3 to +3 m and the rotation between -6 and 6 degrees simultaneously.These values also show, that initial values captured by good UAV sensors should be enough.If these sensors are not given or have a lower Numerical values for the registration error are given in Table 2 and compared to the ICP approach.These values are estimated using the same point cloud and initial parameters are given in the table.This evaluation is some kind of unfair because the ICP gets the exactly same point cloud which means that the ICP reaches 100 percent overlap while the proposed approach needs to handle a grid with 2 m cell size.But the evaluation shows that the proposed approach is able to handle the grid size and creates accuracies which are near to the ICP ones.Since the shown scene shows a locally flat area, there is one degree of freedom left in the translation parameters and also one degrees of freedom in the rotation parameters.Therefore this evaluation on the shown scene is to pessimistic for the actually approach.This can also be seen, on the x translation, which gives the dm accuracy, the other parameters have an error on cm level.So for future evaluations another scene which shows a better landscape for registration should be considered.
Compared to an ICP approach it can be said that the proposed approach creates similar results, but works also in strong border conditions.The approach handles different point density as well as differences in the level of detail of the landscape.

CONCLUSION
This paper presented a registration method for multi sensor point clouds and shows an example of a registration of airborne Li-DAR data and photogrammetric data from UAV taken images.
Beneath the multi sensor application and therefore the need for a geometric approach, there are other constraints which are taken into account and the proposed approach is able to handle these constraints.The main focus for applications goes to a rural area.Therefore it must be considered that no man made objects are in the data.Furthermore, the different sensor configurations also create differences in the geometry, which are different point densities or slightly different level of details.Another constraint can be seen by the further application for the registered data.To define significant changes of the geometry which does not rely on measurement or registration errors there is a need to have the variance of the transformed 3D points of one point cloud to the source.The proposed approach is able to consider the measurement accuracy of each point in the two point clouds and take this into account to determine the transformation parameters as well as the variance of this parameters.With this stochastical model, the variance of a transformed point from the target point cloud to the source can be calculated and taken into account for definition of a significant change.
The proposed approach uses a least squares minimisation of point to grid distances.The source grid is calculated from ground points.Therefore, an outlier removal creates a binary classification into ground and non ground points during the registration process.The achieved registration accuracy can be compared with control point based or ICP based methods.But compared to the ICP the proposed approach have no limitations due to different point densities or level of details.However, there is the need of good initial parameters, as by ICP-based methods.These initial parameters can be get from UAV mounted sensors or from the output of an coarse registration.
Future work aims for a change detection in the river bed.These change detection should divide geometric changes from changes caused by registration or sensor errors.Therefore the stochastical grid as well as the accuracy of the transformed point to this grid is needed.So for the change detection the same distances are used as in the proposed approach for registration.But the change detection will include a statistical test if the distance is significant or not.To exclude potential water points inside the registration, it is possible to mark the water cells of the grid as bad registration areas by assuming a high variance here and therefore ignore all observations which belongs to these cells.
and with: lm = z-coordinate of observation point m (S(x, y)) zm = z-coordinate of point m on the grid σ 2 n = variance of observation point σ 2 δ = variance of the regularisation term The parameters of the transformation to estimate are the ones of the 3D similarity transformation , z0 = translation along x,y,z axis x, y, z = coordinates of the point to be transformed m = scale R = rotation The observation of a single point p n = xn yn zn T is fn = S(txn, tyn) − tzn (7)

A
= functional matrix (derivatives of the observations with respect to the parameters pi) of conditional observations (derivatives of the observations with respect to the coordinates xi) matrix of observations f = vector of fn for each pointTo simplify the equation it is assumed that there are no correlations between the 3D points of the target cloud or between the cells of the source grid.This is truly not the case but simplifies the equations such that the parameters can be estimated in real time (meaning on the UAV flight) and no full matrix must be saved, which saves memory resources as well.Especially by considering big number of points in the target point clouds, the saving of working memory gets very interesting.In future works, it can be tested how this simplicity influenced the final registration accuracy compared to the estimation speed up.With no correlations, Q becomes a diagonal matrix and the multiplication result BQB T becomes a diagonal matrix with scalar values for a single point.For a single observation point p n it becomes

Figure 2 .
Figure 2. Sheme of threshold definition: The bins discretize the point to grid distances, values are the number of points corresponding to one bin, the threshold t is searched on the right side of the bin with maximal value.something similar).And second, ground points are the major count of points in the target data.With these conditions, the distance threshold can be found in the histogram of distances on the right side of the highest peak.Therefore, the dynamically determination of the distance threshold is to find the bin inside the histogram which is below a percentage threshold of the highest value and on the right side to the highest value.Figure 2 exemplary shows the determination of the distance threshold (t) in the histogram of distances.The values are the number of points corresponding to the distance which is discretized with the corresponding bin.

Figure 3 .
Figure 3. Data overview, a) LiDAR Data, b) photogrammetric data, the marked area shows the area of the beaver dam which can be seen in a closer look in Figure 4

Figure 4 .
Figure 4. Beaver dam a) in LiDAR data, b) in photogrammetric data

Figure 5 .
Figure 5.Initial alignment used for fine registration the marked areas show the visible misalignment

Figure 8 .
Figure 8. Results shown with eliminated outliers in marked areas: a) original data, b) data without outliers

Figure 9 .
Figure 9. Results shown on the exemplary beaver dam accuracy a coarse registration should be considered.

Table 1 .
numeric results for registration of LiDAR and photogrammetric point cloud