ASSESSMENT OF THE QUALITY OF DIGITAL TERRAIN MODEL PRODUCED FROM UNMANNED AERIAL SYSTEM IMAGERY

Production of digital terrain model (DTM) is one of the most usual tasks when processing photogrammetric point cloud generated from Unmanned Aerial System (UAS) imagery. The quality of the DTM produced in this way depends on different factors: the quality of imagery, image orientation and camera calibration, point cloud filtering, interpolation methods etc. However, the assessment of the real quality of DTM is very important for its further use and applications. In this paper we first describe the main steps of UAS imagery acquisition and processing based on practical test field survey and data. The main focus of this paper is to present the approach to DTM quality assessment and to give a practical example on the test field data. For data processing and DTM quality assessment presented in this paper mainly the in-house developed computer programs have been used. The quality of DTM comprises its accuracy, density, and completeness. Different accuracy measures like RMSE, median, normalized median absolute deviation and their confidence interval, quantiles are computed. The completeness of the DTM is very often overlooked quality parameter, but when DTM is produced from the point cloud this should not be neglected as some areas might be very sparsely covered by points. The original density is presented with density plot or map. The completeness is presented by the map of point density and the map of distances between grid points and terrain points. The results in the test area show great potential of the DTM produced from UAS imagery, in the sense of detailed representation of the terrain as well as good height accuracy. * Corresponding author


INTRODUCTION
Digital terrain model (DTM) is one of the main topographic classes and gives the digital representation of the bare terrain.It is quite difficult to represent it in all details in a digital form as it is a complex surface.Collected data of the terrain are usually in a form of points in regular or irregular pattern that can be additionally augmented with vector lines as break and structural forms of the terrain.
Different methods for terrain data acquisition are available.Traditional aerial photogrammetry has long been the only economical and accurate approach to mass topographic data acquisition.In the last decade, aerial laser scanning or lidar became a complementary approach that has significant advantages in forest and dense vegetation areas.In the recent years, small unmanned aerial systems (UAS) equipped with consumer-grade imaging sensors offer low-cost and automated production of point clouds, although usually limited to smaller terrain coverage.
Digital terrain model that is filtered from the automatically generated point cloud from UAS imagery, usually resulting from the Structure from Motion (SfM) approach, and the consequent orthophoto are the main products from the UAS imagery.As the UASs are financially affordable also for small enterprises, we are globally-wide witnessing an increasing use of the UAS for surveying purposes.This is very positive, but there is a concern about the quality and reliability of the processed data.The more the process is automated more strictly the quality assessment of the products should be implemented.
Different approaches to assess the capability and quality of UAS-based photogrammetric data collection and products can be found in publications.Dayamit et al. (2015) compared UAS elevation model with lidar data and found out that lidar technology is still more accurate method, but UASs are more flexible and bring good solutions for many applications.The capability of UAV-based data collection was thoroughly investigated by Haala et al. (2011) with comparison to conventional aerial survey with state-of-the-art digital airborne camera systems.The research of Ruiz et al. (2013) clearly proved that UAV imagery is strongly affected by positioning error due to the use of low-cost GPS/INS equipment.There is correlation between the position errors and the final DEM accuracy.An approach to assess the accuracy of georeferenced point cloud of a natural coastal site with total station survey is presented in Harwin et al. (2012).
The aim of this paper is to present a practical example of the quality assessment of the DTM produced from the UAS imagery.First we describe the study area, the aerial surveying mission, and terrain field survey in Chapter 2, all the necessary steps of data processing are then described in Chapter 3. The main results of the quality assessment are presented in Chapter 4, where methodology and the outcomes are presented.We would like to emphasise that the computer programs and algorithms used to accomplish this research project have been in-house developed.

Study Area
The study area of approximately 17 ha is located at the edge of Ljubljana basin in Slovenia (village Brezje).This is partly open and modestly undulating karstic landscape with some buildings and forest areas (Figure 1).

Surveying Mission with UAS and Field Survey
For aerial surveying mission a Microdrone md4-1000 (Figure 2) has been used, equipped with position and navigation system, communication data link, battery and consumer-grade camera.The UAS can fly in autonomous or manual mode.Prior to the surveying mission, flight and data acquisition parameters have to be planned with dedicated software.We start from the area of interest, required ground sample distance (GSD) and the intrinsic parameters of the used digital camera.In our case, the study area has been divided in two photogrammetric image blocks (Kerin, 2014).For each block the average flying height above the ground was 88 m, the overlapping of images in both directions was 66 %, GSD of the colour (RGB) images was 2.2 cm.Flight time for image acquisition of one mission was 22 minutes.
As indirect method of image orientation was applied, sufficient homogenously distributed control points and additional check points for quality assessment had to be planned.We defined the position of 15 control points over the entire area; this is 9 control points in each image block, and additional 18 check points.In Figure 1 the location of control points (red triangles) and check points (blue squares) is presented.The control and check points were signalized by black circular targets on a white background plate (Figure 3).

DERIVATION OF DIGITAL TERRAIN MODEL
Digital terrain model of the study area was derived from filtered photogrammetric point cloud and was further gridded into resolution size of 0.20 m.

Creation of Photogrammetric Point Cloud
A georeferenced photogrammetric point cloud is a result of image block triangulation and image matching algorithms.The input data for the image block adjustment are measured image coordinates of tie and control points.Different commercial products are available in the market for this purpose.However in our project we used in-house developed software 3Dsurvey by Modri planet.3Dsurvey is a software solution using image data processing for surveying of land and built environment.The basic input for data acquisition are aerial images taken from the UAS, terrestrial images or combination of both.Image matching algorithms and photogrammetric bundle block adjustment are implemented to compute the image orientation parameters and point clouds of different density level from high to low (Peterman, 2015).Classification of point cloud to terrain points and other classes can be accomplished interactively or automatically.The software can then produce digital terrain model and orthophoto.Additionally, some applications are implemented, for example calculation of areas and volumes, measurement of distances, derivation of contour lines of chosen equidistance, automatic production of longitudinal and transverse profile plots at selected locations.

Filtering of Point Cloud and Interpolation of Digital Terrain Model
The raw photogrammetric point cloud of the study area (Figure 4, top), containing terrain, low and high vegetation, and buildings, has 11.45 million points, all containing RGB values.The raw point cloud has been filtered to extract only points of the terrain and has 11.17 million points.Terrain points were then interpolated into regular grid of heights with the resolution of 0.20 m (Figure 5).Visualisation of the DTM surface was done by triangulation method.In Figure 5 we can see slightly undulating landscape due to the karstic ground.We accomplished all these tasks (filtering, interpolation and triangulation) using the 3Dsurvey software.

ASSESSMENT OF THE QUALITY OF DTM
As described in previous chapter, digital terrain model of the study area has been derived from filtered photogrammetric point cloud and has been further gridded into resolution of 0.20 m.
For the assessment of the quality of this DTM we applied the methodology provided and recommended by EuroSDR (Höhle and Potuckova, 2011).The theoretical background, described in chapters 4.1 and 4.2, is extracted exclusively from this reference.Although this methodology has not been specifically developed for the assessment of the quality of point clouds and DTMs produced from UAS imagery, but rather for "nationwide DTMs supporting orthoimage production of larger area project" (citation from Höhle and Potuckova, 2011), we found the potential to use this approach also in UAS imagery originated digital terrain models.However, to implement this approach and to assess the quality of terrain data in our study case we used our own computer programs and developed the needed computer algorithms.
We assessed two accuracy parameters: positional accuracy and completeness.The positional accuracy of DTM was estimated from point cloud.The completeness was assessed by density of point cloud and calculated distances between each grid point of DTM and the nearest terrain point.

Positional Accuracy of DTM
Photogrammetric point cloud is the result of image block adjustment and image matching algorithms.When using UAS for image acquisition, the flying height is usually low and image resolution is high, thus very dense point cloud could be achieved.
A common approach to assess the positional accuracy is to compare the DTM data with reference values for a sample.The accuracy of reference data must be higher at least for factor 3. In the case of photogrammetric point cloud derived from UAS images and georeferenced with ground control points, the expected positional accuracy is in a range of a few centimetres, thus the reasonable method for measuring check points is land surveying, most often applied surveying method is GNSS (Global Navigation Satellite System).
The positional accuracy of DTM can be divided into vertical and horizontal or planimetric accuracy, which are usually assessed separately.The problem with the point cloud or gridded DTM is that it is difficult or even impossible to define the exact position of the corresponding reference point.Thus, in our study case we used targeted check points.We compared the GNSS measured check points with the points manually extracted from the point cloud in the positions identical as much as possible.We computed the difference between DTM value and the reference value for each component (X, Y, Z), which is defined as an 'error', for the sample of check points n.From these values we compute some accuracy measures under assumption that errors are normally distributed.These measures are: Root Mean Square Error (RMSE), mean error, standard deviation.Because gross errors may exist in the DTM and we must eliminate them, we set the threshold for outliers as being higher or equal to 3 times value of RMSE.The complete number of outliers is denoted as N.According to Höhle and Potuckova, 2011 (page 37) and using the mathematical notation for object coordinates (X, Y, Z), the accuracy measures for Z coordinate are computed using the equations in The computed accuracy measures for our study case are shown in Table 3.We generated also Q-Q plot that is presented in Figure 7, which ascertain some deviations from the normal distribution (straight line).
Figure 7. Q-Q plot for the distribution of ΔZ Also the Q-Q plots for X (Figure 8) and Y (Figure 9) coordinates show some deviations from the normal distribution.We can see that the robust accuracy measures (Table 4) give slightly different assessment of the accuracy as the calculated standard measures (Table 3).From the RMSE value, which is usually the standard measure for positional accuracy, we would conclude that in our case the accuracy is: 0.034 m for X coordinates, 0.031 m for Y coordinates and 0.025 m for Z coordinates.These values are similar to the 68.3 % Q computed values: 0.039 m for X coordinates, 0.032 m for Y coordinates and 0.026 m for Z coordinates.However, more reliable measure is the 95 % Q, where we can assess the positional accuracy with the probability of 95 %, in our case: -accuracy of X coordinates is 0.064 m, -accuracy of Y coordinates is 0.051 m, -accuracy of Z coordinates is 0.042 m.

Density of Point Cloud and Completeness of the DTM
Checking the completeness quality parameter for digital terrain model depends of the form of the model.If we are dealing with the original data, the recommended quality measure is the density of data.On the other hand, if we have DTM in a form of regular grid, which is basically the result of interpolation, it is recommended to calculate the distances between each grid point of DTM and the nearest terrain point.As in the original data areas with no or sparse points can exist (e.g. after filtering under vegetation and buildings), these areas are usually filled by interpolation using surrounding points.We have to set a maximum gap distance as a threshold (for example 3 times the grid width) in order to avoid unreal values, what can lead to void areas in the DTM.To visualize the results, a map can be generated (map of distances), which shows areas of distances within selected intervals.
In our study case we produced the density map and the map of distances.Density map of filtered point cloud (Figure 10) shows the density span from 0 to 465 points per square meter.From both, the density plot and histogram of density could be seen that there are a lot of void areas, which is mainly due to elimination of vegetation and building points from the point cloud in the filtering process, and also other reasons (e.g.not enough texture for successful image matching).However, the most of the areas have point density between 100 and 150 points/m 2 .
The map of distances is presented in Figure 12.Points are coloured according to the computed class intervals that are presented in Table 5.Table 5. Classes of distances presented by colours in figure 12 The areas in Figure 12 where orthophoto is shown from the background are void areas.This means that DTM has extrapolated values of component Z that exceed the specified maximum gap distance (i.e. 3 times the grid width).
Additional useful tool for inspection is a histogram that shows the distribution of distances.This histogram in our case (Figure 13) shows point cloud points classified into first 5 classes.From this histogram we can see that the distance of 0.045 m has the highest frequency.Furthermore, 94.35 % of DTM grid points are less than 0.60 m far from the original point.
Figure 13.Histogram of distribution of distances

CONCLUSIONS
The results from our research show that using UAS as the platform for image acquisition and photogrammetric approach to point cloud production can derive a quality terrain data or DTM in larger scales.However, we have to check the quality of the data and this should be done in a proper way.
In the practice, usually the RMSE or standard deviation values are computed for the sample points to check the quality of the terrain data (point clouds or DTMs).However, especially with the UASs the derived point clouds and DTMs can be affected by known or unknown systematic influences, thus it is better to use robust accuracy measures.In addition to this, the filtered point cloud or DTM in regular grid are usually the input data in the orthophoto generation process.The influence of the quality of terrain data on the quality of produced orthophoto is direct, but is not visually expressed in the final product, i.e. orthophoto.The users do not know in which areas the orthophoto might be vague in quality due to the influence of DTM.We are highly convinced that a density plot or a map of distances between interpolated grid points and original terrain points should become obligatory supplement to the products (terrain data, orthophoto).The users can thus easily inspect the quality and can avoid using the data in areas of bad accuracy from the objective reasons.

Figure
Figure 1.Study area

Figure 2 .
Figure 2. Microcopter md4-1000 with remote control The images have been taken by the Olympus PEN E-P2, and lens with 17.0 mm focal length.The sensor has 12.3 megapixels of the resolution.

Figure 3 .
Figure 3. Example of targets Control and check points were measured in the field with GNSS fast static RTK method.

Figure 4 .
Figure 4. Part of RGB coloured raw photogrammetric point cloud (top), classified point cloud (middle; terrain points in red) and filtered terrain points (bottom)

Figure 5 .
Figure 5. Detailed view of coloured DTM with 0.20 m resolution

ΔX
that errors are normally distributed is not always true (different systematic influences can yield to systematic errors in the data set), we have to check the results.For this purpose we generate a histogram of error distribution (by coordinate axes) and / or quantile-quantile (Q-Q) plot.The quantiles of the empirical distribution function are plotted against the theoretical quantiles of the normal distribution.If the actual distribution is normal, the Q-Q plot should yield a straight line.If the histogram of the errors or the Q-Q plot shows deviation from the normality, the robust accuracy measures should be applied, for example the Median of the distribution or Normalized Median Absolute Deviation (NMAD).From the histogram of the error ΔZ in our study case (Figure6) we can see some systematic influences in the results, as the differences in height are mostly positive (photogrammetrically defined points are lower than by GNSS measured points).

Figure 6 .
Figure 6.Histogram of the errors ΔZ

Figure 8 .
Figure 8. Q-Q plot for the distribution of ΔX

Figure 10 .
Figure 10.Density plot of filtered point cloud

Figure 11 .
Figure 11.Histogram of density (points/m 2 ) for filtered point cloud

Figure
Figure 12.Map of distances

Table 1 .
The accuracy measures for X and Y coordinates are computed in the same way.

Table 1
(Höhle and Potuckova, 2011;oordinate(Höhle and Potuckova, 2011; page 37)In our study case the computed difference between DTM value and the reference value (check points) for each component are presented in Table2.

Table 2 .
Positional errors of checkpoints

Table 4 .
Robust accuracy measures