DIGITAL TERRAIN MODELS FROM MOBILE LASER SCANNING DATA IN MORAVIAN KARST

During the last ten years, mobile laser scanning (MLS) systems have become a very popular and efficient technology for capturing reality in 3D. A 3D laser scanner mounted on the top of a moving vehicle (e.g. car) allows the high precision capturing of the environment in a fast way. Mostly this technology is used in cities for capturing roads and buildings facades to create 3D city models. In our work, we used an MLS system in Moravian Karst, which is a protected nature reserve in the Eastern Part of the Czech Republic, with a steep rocky terrain covered by forests. For the 3D data collection, the Riegl VMX 450, mounted on a car, was used with integrated IMU/GNSS equipment, which provides low noise, rich and very dense 3D point clouds. The aim of this work is to create a digital terrain model (DTM) from several MLS data sets acquired in the neighbourhood of a road. The total length of two covered areas is 3.9 and 6.1 km respectively, with an average width of 100 m. For the DTM generation, a fully automatic, robust, hierarchic approach was applied. The derivation of the DTM is based on combinations of hierarchical interpolation and robust filtering for different resolution levels. For the generation of the final DTMs, different interpolation algorithms are applied to the classified terrain points. The used parameters were determined by explorative analysis. All MLS data sets were processed with one parameter set. As a result, a high precise DTM was derived with high spatial resolution of 0.25 x 0.25 m. The quality of the DTMs was checked by geodetic measurements and visual comparison with raw point clouds. The high quality of the derived DTM can be used for analysing terrain changes and morphological structures. Finally, the derived DTM was compared with the DTM of the Czech Republic (DMR 4G) with a resolution of 5 x 5 m, which was created from airborne laser scanning data. The vertical accuracy of the derived DTMs is around 0.10 m. * Corresponding author


INTRODUCTION
Lidar technologies are commonly used in many applications during the last fifteen years.The derived three-dimensional data are regularly used for digital terrain and surface modeling.It is often the first task in the processing chain of Lidar data.Airborne laser scanning (ALS) has become a standard method for gathering topographic data used for the generation of digital terrain models (DTMs).However, processing of scanned data can be challenging in a forest, rock, and other hard rich places, due to occlusions, especially from airborne laser scanners.In our work, we decided to use mobile laser scanning system (MLS) to capture the road environment of the nature reserve called the Moravian Karst (Czech Republic).It is one of the most important karst areas of Central Europe, with many unique geological features including caves and gorges.The terrain is also covered by forest.
Mobile laser scanning systems lie somewhere in a middle between terrestrial and airborne laser scanners in the sense of point density and achievable 3D accuracies (Lehtomäki et al., 2016).MLS allows the high precision capturing of the reality in a fast way because of its close position to the objects.The newest MLS systems consist often of two laser scanners mounted on one platform allowing an almost 360 o coverage of the environment.Therefore, MLS have become popular not only for 3D city capturing but also for scanning nature objects.
For example, Bitenc et al. (2011) evaluate the use of MLS for monitoring sandy coast; Gorte et al. (2015) and Lindenbergh et al. (2015) present tree separation and extraction from mobile laser scanning data.To model water flow over the surface, Wang et al. (2014) compute road side environment from MLS data.Lehtomäki et al. (2016) present a study about the automatic machine-learning based object recognition of MLS point clouds in a road and street environment.Jaakkola et al. (2008) model the road surface as TIN network and propose automatic methods for road marking and kerbstone classification.None of these processes would be possible without detection of ground points.Accurate extraction of the terrain points helps to avoid mistakes in the mentioned further processing steps.There are many filter algorithms for terrain points filtering.Sithole and Vosselman (2004) with ISPRS Working Group WGIII/3 present eight solutions and their comparison for extraction ground points from ALS data.They state, that each data set needs its unique approach and still it cannot be fully automated.Due to the geometry of an MLS scanner, generation of DTM needs improved approach.The extraction of ground surface scanned by MLS (in urban areas) can be done by different ways: by robust local weighted regression (Nurunnabi et al., 2013); by image based extraction of terrain points based on a global height estimation and local planarity priors (Vallet and Papelard, 2015); analysis of extracted large planes at a certain distance below the 3D trajectory of the MLS (Pu et al., 2011).
In our work, we created a DTM from several MLS data sets in the neighbourhood of the road in the Moravian Karst.Classification of the terrain points is based on fully automatic, hierarchical interpolation and robust filtering (Kraus and Pfeifer, 1998).Robust filtering was originally designed for steep, wooded terrain (Pfeifer et al., 1999), which is typical for the Moravian Karst landscape.We have used three automatic interpolation methods for interpolation of the DTMrobust moving planes, delaunay (triangulation) and moving paraboloid.All algorithms are implemented in the software OPALS (Mandlburger et al., 2009).
Study area and data description are presented in Section 2. The proposed methodology for classification of the terrain points and interpolation methods for final DTM generation are described in Section 3. The results and the validation of the DTMs and its comparison with DMR 4G (national DTM of Czech Republic) and geodetic measurements are demonstrated in Section 4. Section 5 shows the conclusions.

Study area
The study area is located in the eastern part of Czech Republic to the north of Brno (Figure 1).The scanning process was performed in two hours, by scanner Riegl VMX -450 in May 2015.The total length of two covered areas is 3.9 km and 6.1 km respectively, with an average width of 100 m.The two Riegl VQ -450 laser scanners, the IMU and the GNSS equipment are integrated into measuring head, which was mounted on a roof of a car.The measurement rate of each scanner is 550,000 meas./sec, with a scan rate about 200 profiles/sec.Each scanner provides 360 o gapless profile, with high accuracy (8 mm) and precision (5 mm).Тhe first data set (Road 1) contains 8 LAS files (4 from scanner_1 and 4 from scanner_2) with different size from 0.3 to 4.5 gigabytes; the second data set (Road 2) contains 12 LAS files (6 from each scanner) from 0.5 to 3.5 gigabytes.

DMR 4G
The Digital Terrain Model of the Czech Republic of the 4th generation (DMR 4G) represents relief of the full Czech Republic in a regular grid of 5 x 5 m.DMR 4G is a product created by The Czech Land Survey Office, in cooperation with the Military Geographical and Hydrometeorological Office from airborne laser scanning data.Laser scanning was performed with LiteMapper 6800 system during 2010 -2013 years, on an average height 1200 -1400 m above the terrain in individual blocks.For data filtering and interpolation the software SCOP ++ 5.4 was used.DMR 4G is presented like shaded relief in S-JTSK coordinate system and Height in the Baltic Vertical Datum.The total standard error in bare terrain is near 0.3 m and in forest areas around 1 m.All other information about DMR 4G can be found in Technical Report for DMR 4G in Czech language (Brázdil et al., 2015).In Figure 2, the shaded representation of the DMR 4G is shown for the study area in the Moravian Karst.

Algorithms
The proposed methodology for deriving the digital terrain model (DTM) from the MLS data is summarized in the following steps.Each step corresponds to a specific OPALS module and is applied for filtering the terrain points and for interpolating the terrain surface model through the terrain points.The average point density of the original point clouds is about 1000 [pts/m 2 ].In some areas on the road or where the roadside trees are very close to the MLS, this number is even higher.The big capacity of data sets makes it hard to handle and process them.For this purpose, we imported files, one by one, to OPALS and applied morphological filtering (steps 1, 2).In each 0.15 x 0.15 m cell we chose the point with smallest Z-value, including all echoes.However, this filter is not enough to get rid of all non-terrain points, which is shown in Figure 3.The thinned data sets can now be joined (step 3) for further processing.Using the combination of triangulated and robust moving planes grids, we filled gaps in the second one, which were still appearing at trees positions.Then, attribute "normalized Z" can be added to the point cloud by subtracting the preliminary DTM from the z-values.The second interpolation was implemented in the same way, but with other limitations (steps 9 -13).Now, for each (0.25 x 0.25) m cell we chose the second minimum Zvalue and height threshold from (-1 to 0.3 m).This asymmetric height threshold was used to ensure, that all points which are below the fitted plane will be considered.Then, we again interpolated two grids with cell size (0.25 x 0.25) m and renewed our attribute "normalized Z".Finally, extracted terrain points are used for digital terrain generation (step 15).For this purpose, we used three interpolation techniques: robust moving planes, moving paraboloid and triangulation.Each method has its own parameters.For moving planes method, the number of points was set to 30 and maximum search radius 2 m.We tried different combinations of these two parameters.With other settings, the resulted DTM became much smoother or output can be zero.For moving paraboloid method, the number of points was set to 100 and search radius 2 m.For triangulation, we also applied 2 m search radius.Each DTM is generated in (0.25 x 0.25) m grid.

Comparison with DMR 4G and quality validation
To compare the different data sources (step 16) it has been ensured that the coordinate system are the same.Mobile mapping data and products have ellipsoid Height and UTM 33 coordinates.The Czech Land Survey Office provides data only in S-JTSK/Krovak (East North) -ESPG:102067 coordinate system and Baltic Height.Because DMR 4G has less grid points, we decided to transform it to UTM 33.For this purpose, we used software QGIS.To compute ellipsoid heights for DMR 4G we used quasigeoid CR-2005.To accomplish the comparison, raster images should also have the same grid size.
Before the coordinate transformation, we resampled DMR 4G to grid size (1 x 1) m, and interpolated new DTMs from classified terrain points with grid size (1 x 1) m. Results of the comparison are presented in Table 2, Figure 6 and Appendix.
For the quality check of interpolated DTM, 87 points were measured by a total station on the road area of data set 1 (step 17).The measured profile has a total length about 1 km.For comparison, we chose DTM interpolated by robust moving planes with grid size (0.25 x 0.25) m.The DTM is subtracted from control points.The results are presented in Figure 7.

RESULTS
Shadings of the generated DTMs can be seen in Figures 4 and  5.In Figure 4a-c and 5a-c the shadings of DTMs derived from the different interpolation methods are shown.
As we can see in Figure 4b and Figure 5b triangulation interpolation better fill gaps.Robust moving planes method claims to be better for DTM generation in forest areas, during fitting tilted plane, it also detects outliers, which do not take part in interpolation.That is why another solution for DTM generation in our landscape is the combination of robust moving planes and triangulation methods (main images in Figure 4 and Figure 5).
Table 2 presents the height differences between DMR 4G and DTMs interpolated with different interpolation techniques.These differences are summarized from the histograms in Аppendix.As can be seen in Table 2 and in Appendix histograms C and F, the height difference between DMR 4G and triangulated DTM shows the worst results.The average height difference between mentioned grids is more than 3 m when between others it is near 1 m.On the histograms C and F it is shown, that even height differences more than 10 m take almost 15 % in both cases.Triangulation interpolation is not suitable for DTM generation from MLS data.It fills gaps better, but also connects outlier and isolated points that affect the whole grid accuracy.
The average height difference between DMR 4G and robust moving planes DTM is -0.877 m; -1.415 m and moving paraboloid DTM is -0.861 m; -1.343 m.However, for moving paraboloid DTM, the minimum and maximum height differences with DMR 4G are much bigger than from robust moving planes DTM.These results indicate that gross errors are still appearing in generated DTMs and some outlier points take place.
In the histogram A (DMR 4G and robust moving planes DTM comparison) it is shown, that mostly 70% of compared data, has height differences in range from -1 to 1 m.As was mentioned before, the total standard error of DMR 4G in forest areas is 1 m.Let us have a look in Figure 6.It shows the visual profile comparison of DMR 4G (red line) and DTM interpolated by robust moving planes (green line).Axis X is the profile length in meters and axis Y is the height of the terrain.In Figure 6a, it is shown that differences between data are increasing on borders.From 10 to 100 m range the planes are overlapping or difference is less than 1 m, but on borders it is near 2 m or more.Commonly these areas are covered by rocks and are therefore badly reachable for airborne laser scanners especially in dense forests.Histogram D shows the differences between DMR 4G and DTM interpolated by robust moving planes from data set 2. It is shown, that only 50% of the compared data has differences between -1 and 1 m.The area on data set 2 is more steep and rocky than for data set 1.In Figure 6b, it is shown that differences are appearing along all profile, even on a road area.Lastly, the height differences between the reference points and the derived DTM are shown in the histogram in Figure 7.The average height difference is 0.095 m, the maximum is near 0.20 m and minimum 0.02 m.In the histogram it is shown, that all values are negative, this means that reference heights are lower than DTM heights.This can be impacted by different measurement equipment (total station vs MLS), their accuracy, weather conditions and placement of satellites, etc.Also, this height shift can be explained by the quality of GNSS signal.For good results we need at least 4 satellites, however in such forest areas, this requirement cannot always be observed.These results show the high quality of derived DTM.However, other areas should also be explored.Firstly, MLS data sets are processed and filtered using hierarchical interpolation and a robust filtering.Filtered terrain points take in average only 2% of the data set.After extraction of ground points, different interpolation methods are applied for DTM generation.The grid size of interpolated DTMs is (0.25 x 0.25) m.The average vertical accuracy is 0.10 m, checked by geodetic measurements of 87 points on the 1 km long road area.The derived DTM can be used for analysing terrain changes and morphological structures, modelling of water flows and exploring geological features.We have compared generated DTMs with DMR 4G, which has been generated from airborne laser scanning data.The biggest height difference occurs in rock places, which are difficult to reach, and which are often affected by occlusions.
In future work, we will have to consider more filtering algorithms or their combination, to improve a better quality of terrain extraction from MLS data in rock areas to eliminate outlier, non-terrain and other classification errors that have a big influence on resulting products.

Figure 2 .
Figure 2. Shading of the DMR 4G for the study area in the Moravian Karst (grid size 5 x 5 m), © Data Copyright 2016, ČÚZK

Figure 3 .
Figure 3. Profile views colored by echoes (blue (first echo)greenyellowred (last)): (a) row point cloud from 1 scanner (b) filtered point cloud with attributes -(0.15 x 0.15) m cell and smallest Z-value from 1 scannerIn steps 4 to 8 and 9 to 13, different combinations of hierarchical interpolations and morphologic filterings are applied to prepare data for robust filtering.Before the first interpolation, we filtered point cloud again.Now, for each (1 x 1) m cell, we chose the second smallest Z-value.The aim of this step is to minimize outliers and make points more homogenous.Then, two interpolation techniques (i.e.robust moving planes and triangulation) are applied with a cell size of 1 x 1 m.For robust moving planes, the number of points used for interpolation was set to 12 with a max.search radius of 2 m.Using the combination of triangulated and robust moving planes grids, we filled gaps in the second one, which were still appearing at trees positions.Then, attribute "normalized Z" can be added to the point cloud by subtracting the preliminary DTM from the z-values.The second interpolation was implemented in the same way, but with other limitations (steps 9 -13).Now, for each (0.25 x 0.25) m cell we chose the second minimum Zvalue and height threshold from (-1 to 0.3 m).This asymmetric height threshold was used to ensure, that all points which are below the fitted plane will be considered.Then, we again interpolated two grids with cell size (0.25 x 0.25) m and renewed our attribute "normalized Z".

Figure 4 .
Figure 4. Shading of the DTM (grid size 0.25 x 0.25 m): main imageroad 1 combination of robust moving planes and triangulation; (a) moving paraboloid; (b) triangulation; (c) robust moving planes The histograms A, B and C show height differences from comparison DMR 4G and data set 1, the histograms D, E, F -comparison DMR 4G and data set 2. For DTM interpolation new parameters were used: grid size (1 x 1) m, search radius 2 m for triangulation, search radius 1.5 m for moving planes and paraboloid.The number of points used for robust moving planes interpolation are 35 and for moving paraboloid interpolation -50.The new grid was generated from resampled points from DMR 4G using 50 m search radius and 200 points by robust moving planes interpolation.

Figure 6 .
Figure 6.Visual profile comparison of DMR 4G and robust moving planes DTM (a) data set 1 (b) data set 2

Figure 7 .
Figure 7. Height difference between control points and DTM robust moving planes 5. CONCLUSIONS In this work, we have presented a fully automatic algorithm for generation of Digital Terrain Model from several Mobile Laser Scanning data sets captured in the neighbourhood of a forest road.
In Table1we can see the number of captured points of an original point cloud, and points before and after robust filtering.Filtered terrain points take in average 2% of measured data set.98% of other points were filtered based on combinations of hierarchical interpolation and a robust filtering.The average point density of classified terrain points is about 33 [pts/m 2 ].

Table 2 .
Height difference between DRM 4G and DTMs