BARE-EARTH EXTRACTION AND DTM GENERATION FROM PHOTOGRAMMETRIC POINT CLOUDS WITH A PARTIAL USE OF AN EXISTING LOWER RESOLUTION DTM

A method of extracting bare-earth points from photogrammetric point clouds by partially using an existing lower resolution digital terrain model (DTM) is presented. The bare-earth points are extracted based on a threshold defined by local slope. The local slope is estimated from the lower resolution DTM. A gridded DTM is then interpolated from the extracted bare-earth points. Five different interpolation algorithms are implemented and evaluated to identify the most suitable interpolation method for such non-uniformly scattered data. The algorithm is tested on four test sites with varying topographic and ground cover characteristics. The results are evaluated against a reference DTM created using aerial laser scanning. The deviations of the extracted bare-earth points, and the interpolated DTM, from the reference DTM increases with increasing forest canopy density and terrain roughness. The DTM created by the method is significantly closer to the reference DTM than the lower resolution national DTM. The ANUDEM (Australian National University Digital Elevation Modelling) interpolation method is found to be the best performing interpolation method in terms of reducing the deviations and in terms of modelling the terrain realistically with minimum artefacts, although the differences among the interpolation methods are not considerably large. * Corresponding author


INTRODUCTION
Photogrammetry has long been limited by technological difficulties in acquiring dense 3D topographic data.However, recent innovations in imaging sensors and positioning and orientation tools have created the possibilities of acquiring images at high spatial and spectral resolutions.Simultaneous improvements in computer hardware performances and image matching software have enabled dense matching of stereo images with faster speed and higher precision making it a competing alternative to ALS (Leberl et al. 2010).
Dense 3D point clouds are created by computing the crossing points of all homologous pixels of two or more stereo images (i.e.dense matching).The crossing points can be either on bare ground surfaces or on top of any off-ground objects.Ground surfaces are represented only in open areas where imaging rays from at least two camera positions can reach.The chances of two imaging rays reaching the ground in partially open forests get better with higher overlap of the stereo images.The point clouds are therefore digital representations of the top surfaces, i.e. digital surface model (DSM).
The use of DSM created through dense matching of stereo images, for example in computation of the canopy height model (CHM) in forestry, necessitates normalization of the DSM.The DSM is normalized by subtracting an accurate digital terrain model (DTM) of the area which is not often available.In the absence of an accurate DTM, the choices are either to use low quality DTM that often exists at national level or to compute a new DTM from the photogrammetric DSM itself.The low resolution DTMs that often exist at national level have error margins often not significantly lower than the average CHM, making the computation of CHM and other forest attributes less reliable.Generating a new DTM from the DSM requires identification of bare-ground points from which the DTM could be generated using a robust interpolation method.The generation of new such DTM is only possible if satisfactory number of ground points can be identified which is not possible in areas with complete canopy cover.
The technique of identifying ground points from 3D point clouds is often referred to as ground filtering or bare-earth extraction (Wang et al. 2009;Meng et al. 2009).A number of ground filtering methods have been developed over the years especially for ALS point cloud data.The ground filtering methods commonly use some kind of discriminant function to classify ground and off-ground points.Behind each discriminant function lies an assumption or a model of some attributes of the terrain such as elevation, slope, curvature, etc.Such assumption is often made based on the point clouds themselves making the problem somewhat ill-posed.In ALS, there are often adequate local ground points as ALS can penetrate forest canopies, shadows and light clouds.Photogrammetry is typically constrained under those conditions.Such constraints make ground filtering algorithms used in ALS inapplicable of filtering photogrammetric point clouds, at least without modifications.A priori information about the local terrain might help improve the ground filtering in photogrammetry.One interesting such information is low resolution DTM that often exists at national level in many countries with varying qualities.Although the absolute elevation values of these DTMs may not be accurate, they may be useful in modelling some attributes of the terrain such as slope.Identification of ground points from photogrammetric point clouds with a partial input from such data needs to be thoroughly investigated.
Once ground points are extracted, a gridded DTM can be created by interpolating the elevation values of the areas lacking ground points.Theoretical and empirical reviews of the interpolation methods used for scattered data show that each method performs best under specific conditions (Zhang and Gruen 2006;Tan and Xu 2014;Li, Cheng, and Lu 2000;Yang et al. 2004;Babak and Deutsch 2009).The choice of the right interpolation method depends on the dataset, the type of terrain and the application field.
In this study, an algorithm that extracts ground points from point clouds made by dense matching of aerial stereo images is implemented and evaluated.The algorithm uses a lower resolution national DTM to model the local terrain on which the extraction of the ground points is based.The algorithm shares some elements with other ground filtering methods used on ALS data such as the slope based filtering and the directional scanning methods (Vosselman 2000;Meng, Currit, and Zhao 2010;Meng et al. 2009).The study also investigates the performances of some interpolation methods in creating gridded DTM from the extracted ground points.

Datasets and Test Areas
The Photogrammetric point clouds: The aerial images were acquired with radiometric resolution of 12 bit and ground sampling distance (GSD) of 35 cm with 60% along-and 20% across-track overlap.The stereo images were matched using the software called MATCH-T version 5.4.1.The extreme terrain option with the 2.5D filtering was used for the matching.Hierarchical matching starting with feature-based matching and ending with block-based matching was used.The matching resulted in one point per pixel, i.e. mean point density of about 9 points per m 2 .
The lower resolution DTM: The national DTM (DTED10 henceforth), produced by the Norwegian Mapping Authority, has spatial resolution of 10 m and accuracy varying from region to region.The data sources used for the DTED10 are contours from the national thematic map at scales ranging from 1:500 to 1:100000 depending on the region.The Authority states that the mean square error of the DTM varies from 2 m to 6 m (Krtverket, 2015).Slope is computed from the DTED10 as the steepest downhill descent in a 3 by 3 window (Gallant and Wilson 2000).

The reference DTM:
A DTM created from ALS data is used as a reference against which the newly computed DTM is evaluated.The reference DTM is created at spatial resolution of 1 m.The ALS data is available only for a part of the country.
Test sites: Four test sites are selected based on their topographic characteristics and availability of the reference DTM.The sites get more complex in terms of topography and above ground cover from the first to the fourth (Figure 1

The New Ground Filtering Algorithm
The new algorithm is based on the principle that the value of terrain slope is not much affected by elevation biases of a DTM as long as the slope is not highly undulating and the spatial resolution of the DTM is not too low.Assume a case where the terrain is perfectly flat.Regardless of the DTM resolution and elevation bias, the slope remains zero as long as there is no random error and the distance over which the slope is computed is kept within the bounds of the flat area.The same is true for a terrain with a uniform (constant) slope.The slope starts to change with scale (DTM resolution) if there is local undulation.The odds of getting local undulations increase with ground distance or DTM pixel size.Hence, as long as the distance over which the slope is computed is kept fairly short, the effects of local undulation is expected to be minimal.If the resolution of a DTM is not very low and the terrain is not highly undulating, the slope of the terrain that is computed using low resolution DTM can be used to estimate the slope of the same area at a higher resolution.This could be formulated mathematically in the following paragraph.
Take a point P in the point cloud and create a window of, for example, 10 m by 10 m around it, the minimum elevation within this window is identified as, say, point M. Assume for now that point M is actually a ground point.Z P and Z M are the true elevation values of points P and M, while Ź P and Ź M are their respective values in the DTED10.The difference between Z P and Ẑ P is the error ξ of the DTED10 (Equation (1)), i.e.
bias.The slope (θ) between points P and M is computed based on the distance (d) and the elevation difference (∆Z) between the two points (Equations (2-4)), assuming that the direction of steepest descent actually passes through point M. The elevation difference can be computed from the true elevation values of points P and M (Equation ( 4)) which is related to their elevation values in DTED10, Ẑ P and Ẑ M respectively in Equations (5 and 6).Equation 5shows that the bias of DTED10 (ξ) does not a b C d have an effect on the slope as it is cancelled during the subtraction.Only the random error expressed in terms of error standard deviation (σ ξ ) remains.Therefore, the slope computed from the DTED10 can substitute θ to model the local terrain slope as long as σ ξ is not too high (Equation ( 6)).The slope can thus define the expected slope-based elevation difference between the minimum point and any point in the window, i.e. the threshold elevation difference (H τ ).The threshold elevation difference (H τ ) can be derived from Equation (2) by adding an adjustment term (h) to account for the random error σ ξ and/or the local undulations (Equation ( 7)).
Ground points have an elevation difference (∆Z) less than or equal to the H τ as in the case of point P while off-ground points, for example point C in Figure 2, have ∆Z greater than H τ .
For the computation to be valid, the minimum point M should necessarily be a ground point.First, point M is checked against the ground points in the neighbourhood.Neighbourhood is defined based on the distance over which the DTED10 slope is computed.A distance over 3 DTED10 pixels is 30 m; that is, the maximum distance used to search for the nearest ground point.The nearest ground points are identified if available, and their elevation differences with that of point M are evaluated in relation to the threshold set by the slope.If the slope is flat or nearly flat, i.e. below 5 degrees, all the elevation differences within the neighbourhood are constrained to be below the threshold.If the slope is steeper, the constraint is limited to only the nearest point as elevation changes are more common in such areas.Second, if there is no ground point in the defined neighbourhood, an approach similar to that of the slope-based method of Vosselman (2000) is used.The gradient at the point is evaluated against the threshold set based on the DTED10 slope for all the nearest points in all directions.If both point P and the minimum point M fulfil the requirements for ground point, both points are considered as ground points.The process continues scanning through the points line by line in a similar fashion to the directional scanning methods of Shan and Aparjithan (2005).To improve the computational efficiency, a certain interval is defined for scanning through the data instead of visiting all points.

DTM Interpolation
To create a gridded DTM from the identified ground points, five interpolation algorithms known for their performances in scattered data interpolation are implemented and their results are evaluated.These are the inverse distance weighting method (Shepard 1968), the natural neighbour interpolation algorithm (Sibson 1981), the regularized spline (Franke 1981(Franke , 1982)), the triangle-based piece-wise cubic interpolation (Ramos 2001;Alfeld 1984) and the ANUDEM (Australian National University Digital Elevation Modelling) method (Hutchinson, Xu, and Stein 2011;Hutchinson 1989).The interpolation methods have continuously been adapted over the years.The DTMs are created at 1 m spatial resolution.

Evaluation
The interpolated DTM are evaluated using the ALS DTM as reference.The deviations of the interpolated DTM (Ź i ) from that of the ALS DTM (Z i ) is computed (Equation ( 8)).The mean (Equation ( 9)), the standard deviation (Equation ( 10)), and the root mean square (Equation ( 11)) of the deviations are computed.Besides, the Median Absolute Deviation (M MAD , Equation ( 12)) is used to compute a robust estimate of the spread of the deviations.M Mad is computed as the median (M) of the absolute deviation of each error value from the median of the errors multiplied by the consistency constant (k).The consistency constant is a distribution dependent constant and is approximately equal to 1.28 for a normally distributed data.M Mad estimates the population standard deviation of the error (Rousseeuw, 1993).

RESULTS AND DISCUSSION
The qualities of the new DTMs are evaluated in relation to the ALS DTM by computing the deviations between the two.The statistics of the deviation (Table 1) is limited to 99% of the data to avoid the influence of large outliers.The deviations of the DTMs created using the five methods do not show large differences However, for all of the four test sites the least deviations are consistently recorded by the ANUDEM method, followed by the Natural Neighbour and IDW closely.The capability of ANUDEM to combine local and global approaches (Hutchinson, Xu, and Stein 2011;Hutchinson 1989) makes it more robust to a highly scattered data such as test site 4. The lack of large difference between the algorithms is in agreement with other studies that empirically compared most of the interpolation algorithms (Qi-uan, Wan-chang, and Jun-hui 2004; Yang et al. 2004).
For the open areas (site 1), the RMSE is down to 0.5 m as registered by the ANUDEM.The ground points are dense in comparison to the other sites as the surface is bare.This contributes to the better accuracy of the interpolation.Such areas have therefore double advantages: that is, dense ground points and better interpolation accuracy due to the dense points.The choice between interpolation methods is less important for this test site.
The relatively low accuracies in test sites 2 and 4 are basically related to the inability of photogrammetry to penetrate forest canopies.The view shed of site 2 is shown to visually compare the computed DTM with the original DSM and the ALS DTM (Figure 3).The view shed of the new DTM (d) is more similar to that of the ALS DTM (c) than that of the DTED10 (b).Closer look shows that the new DTM has some artefacts in the forested areas.The artefacts are created due to the sparseness of the ground points and the erroneous inclusion of some off-ground points.The artefacts are more observable in local interpolators than in ANUDEM.For test site 4 there is additional challenge that comes from the undulating terrain with narrow valleys, sharp ridges and shadows.The narrow valleys are not accurately captured by photogrammetry due to occlusion effect.Both the sharp ridges and narrow valleys are not well represented in the DTED10 slope and cannot be accurately accounted for by the undulation factor in the algorithm.Shadow areas have either no or erroneous values in the photogrammetric point clouds as matching often fails in these areas.Therefore, it is not unexpected that test site 4 shows the highest deviation.However, given those challenges to the data and to the algorithm, the deviations are not unexpectedly high.
One of the interesting results is the algorithms performance in the built-up areas.The buildings are successfully excluded.This shows the algorithm is robust against high above ground objects such as buildings.Besides, the terrain in the built-up areas is often not as complicated as other areas.The slope on which the modelling of the local terrain is based is also more accurate in these areas as the DTED10 has better quality in the built-up areas.As seen in the view shed (Figure 4 (d)), some details are lost in the photogrammetric DTM due to the sparse distribution of the ground points used in the DTM interpolation.The deviations of the DTED10 from the reference DTM are also included in the analysis to show the relative importance of creating the new DTM.In all of the four test sites, the deviation of DTED10 from the ALS DTM is the highest.In fact, the newly created DTMs reduce the RMSE of the DTED10 by 27% to 83% (table 1).The highest reduction of error is registered for site 1; whereas, the lowest is recorded for site 3, i.e. the built-up areas where the DTED10 itself is already close to ALS DTM.The error of the newly interpolated DTM is the lowest in test site 1 as the area is bare and the topographic roughness is the lowest.The RMSE of the new DTM are below 2 m indicating that they are well suited to be used for the normalization of the DSM.Although, the error is larger in the dense forest areas, leading to larger error in the CHM, it successfully reduces the CHM errors on bare and sparsely vegetated areas.The study also investigated the performances of different interpolation methods in creating gridded DTM from the extracted ground points.The method is implemented on data from four test sites that vary topographically and in above ground cover.

Interpolatio
The results show that the extracted ground points do not deviate from ALS DTM systematically except that they are influenced by dense forest canopies and highly undulating terrain.The algorithm, not unexpectedly, identifies very sparse ground points in dense forest areas which lead to larger DTM error for these areas.The DTM created by the algorithm has much lower deviation from the ALS DTM than the national DTM indicating the absolute importance of creating such DTM rather than using the national DTM for purposes such as CHM computation and forest structure investigations.
Although the difference between the five interpolation algorithms is not large, ANUDEM is found to be best performing in terms of reducing the RMSE and in terms of realistically modelling the terrain.However, the difference between the interpolation algorithms is not so large that one can chose any of them depending on the purpose of application.If the DTM is going to be used for hydrological and terrain analysis purposes, ANUDEM is most suitable.For purposes such as computation of CHM, where just local accuracy is important and computation time is an issue, any of the other algorithms specially the IDW and the natural neighbour maybe used.
The method proves to be useful in cases where there is lack of accurate DTM and the DSM is to be normalized for further applications such as forest attribute computation.Low resolution and poor quality DTM with error margins not significantly lower than the average CHM are not useful for such purposes.The computation of DTM thus becomes paramount.
The study presented here used a national DTM with spatial resolution of 10 m for the estimation of local slopes.The effect of changes in the spatial resolution on the computed DTM is not investigated and requires separate inquiry.Nonetheless, the possibility of using an existing lower resolution DTM in the process of extracting ground points and thereby creation of new higher resolution DTM from photogrammetric point clouds has wider application potentials.
a to d).Test site 1 (Figure 1 a) is a bare area with a slope varying gradually from flat to steep.Test site 2 (b) is a steep area partly covered by dense vegetation.Test site 3 (c) is a built up urban area and its surroundings.Test site 4 (d) is an area with a complex topography including steep and undulating terrain, dense forests and open areas.The area covered by each test site ranges between 4 km 2 and 17 km 2 .

Figure 1 .
Figure 1.Aerial images of the four test sites: The letter a-d correspond to the test site numbers 1-4 used in the text, respectively.

Figure 2
Figure 2 Sketches of ground and off-ground points together with the threshold height and the actual elevation difference in relation to the reference surface

Figure 3 .
Figure 3.The DSM of the test site 2 (a) together with the view shed of the DTED10 (b), the ALSDTM (c) and the newly computed DTM (d) of the same area

Figure 4 .
Figure 4.The DSM of the test site 3 (a) together with the view shed of the DTED10 (b), the ALSDTM (c) and the newly computed DTM (d) of the same area

Table 1 .
The error variables of the DTMs interpolated by the five different methods for the four test sites