AN ITERATIVE TERRAIN RECOVERY APPROACH TO AUTOMATED DTM GENERATION FROM AIRBORNE LIDAR POINT CLOUDS

This paper presents a hierarchical recovery method to generate DTMs from airborne LiDAR point clouds based one an idea of layering. The developed method first registers the last return points, and then layering them. The layering is done by dividing the points into different height layers and assigning layer numbers to each point. The layer numbers are comparing references in later identification process. Then a series of rasterized pyramid levels which consisted of lowest points are generated. Since the outliers have been removed after the layering, the cells in top level are considered to be terrain points and used as reference to identify cells in the following level. After the identification of the second level, an interpolation will occur in the cells which identified as offterrain. And the interpolated level will be used as reference in its following level and the same process is repeated at each level. Once this process of the bottom level finished, the proposed method adjusts the results based on the first return feedback and followed by the final interpolation. As a result, this produces the final DTM. The developed method is data driven, and does not assume a prior knowledge about the scene complexity. The proposed method was tested with the ISPRS WG III/3 LiDAR datasets covering different terrain types and filtering difficulties. The results show that the proposed method can perform well in flat terrain or gentle slope area.


INTRODUCTION
Compared to conventional methods such as aerial photogrammetry and field surveys, the generation of Digital Terrain Models (DTM) from airborne LiDAR point clouds is fast and cost-effective over a large area, especially in vegetation covered areas since laser pulses can penetrate some of the canopy.However ， developing an automated and robust approach to terrain point identification and DTM generation is challenging.As a preliminary task of DTM generation using airborne LiDAR point cloud data, filtering terrain and offterrain is critical and fundamental to feature extraction and classification (Briese, 2010).The identified terrain points are the input of further interpolation process in many developed algorithms.The inappropriate identification will cause the deviation in the following interpolation, which leads to further error and low accuracy of DTM products (Guo et al., 2010).Besides, filtering is usually very challenging and time consuming because of the algorithms have to processing a large amount of data.Therefore, an efficient and effective filter algorithm is important for DTM generation. 1 However, the current filtering algorithms are facing difficulties in handling complex circumstances such as outliers (points lie far above or below the most points), complex objects, steep slopes, attached objects, uncertainty of the terrain definition (such as the ramp of a bridge), vegetation (such as shrub), discontinuities of the terrain, low objects like curb and railway tracks, as well as the combined complex scene (Sithole and Vosselman, 2004;Meng et al., 2010).Some of these problems are critical.The outliers, especially low outliers, can affect the reference point selection in the algorithms which adopt the lowest points as reference terrain points.The scenarios with * Corresponding author.Phone: 1 519-888-4567, ext.34504.different sizes of the buildings will face a dilemma in choosing a filtering window size.Applying a small window size, the algorithm may identify a point on a big building point as terrain; and using a large window size, small terrain relief variations may be ignored.Low elevation objects are hard to remove because their heights are very close to terrain.A lot of research has been dedicated to DTM generation, especially to filtering, during the last decade (Meng et al., 2010).However the aforementioned difficulties have always been a barrier in developing a fast, robust, and reliable automatic filter and it is the major obstacle in DTM generation from airborne LiDAR data.
The rest of the paper is organized as follows.Section 2 describes the multi-scale terrain filtering method.Section 3 discusses the results obtained using the ISPRS WG III/3 test LiDAR datasets.Finally, Section 4 concludes the study.

METHOD
The proposed multi-scale terrain filtering (MTF) method identify terrain points by iteratively recover terrain models from multi-scale pseudo-grid images in a coarse-to-fine way.As shown in Figure 1, the method has three steps: point cloud preprocessing, multi-scale terrain filtering, and DTM refinement.
In the point cloud pre-processing step, all laser scanning points must be pre-processed to retain last-return points of multiple returns (laser echoes), and then are layered with regard to the statistical height histogram of the whole data set.Two objectives of height histogram based layering are to assign the layer numbers to each point for the following MTF implementation, and to remove lower outliers and noises.The second step is the multi-scale terrain filtering, which includes the rasterized pyramid level generation, iterative point identification and interpolation.Several rasterized pyramid levels are generated at first, and lowest points of every grid in every level are marked as representative points.The highest level is referred to an initial digital terrain model, from which the proposed MTF is employed as a reference.Then the identification and interpolation is iteratively processed in every level from the second highest level to the lowest level in the pyramid.The identification is based on comparing the layer numbers generated in layering and a slope-based threshold.This is followed by the interpolation at identified off-terrain points.The points from a processed level then become reference points in the identification of next level.Iteratively, DTMs are recovered and densified from coarse scales to fine scales.
In the third step, the terrain results are adjusted based on the normalized Digital Surface Model (nDSM).The original LiDAR point cloud data are filtered based on the refined DTM.The separated terrain points are applied to generate the final DTM through the IDW interpolation.As a result, this produces the final and complete DTM. Figure 2 presents the pseudo-code of the proposed method.

RESULTS AND DISCUSSION
The dataset used in this study was released from the ISPRS WG III/3, have been made available through the web site (www.commission3.isprs.org/wg3/).A total of 15 sites were selected to test the performance of our multi-scale terrain filtering algorithm and compare the results with other methods evaluated by ISPRS (Sithole and Vosselman, 2004).The LiDAR point cloud data were captured by an Optech ALTM scanner and the reference data were generated manually.Those data are located along seven study sites over the Vaihingen test field and the centre of Stuttgart City, Germany.The study cites have varied terrain characteristics and diverse feature content (e.g., open field, vegetation, building, road, railroad, river, bridge, powerline, water surface, etc.).Those sites are listed in Table 1.This dataset covers many different land features and filtering difficulties.However, it does not contain small woods and residence in urban area.And the reference data is only available for the 15 samples.The reference data for entire area is not available, which means will limit the algorithm testing for a large area.Sithole and Vosselman (2003b) reviewed and compared eight filtering algorithms, and their comparing method and data are frequently cited and applied in many researches of the laser scanning data filtering (Pfeifer and Mandlburger, 2008;Liu, 2008;Briese, 2010;Meng et al., 2010).This paper adopts part of their assessment method and tests the performance of the MTF method.According to Sithole and Vosselman (2003), the cross-matrices are applied in this study to quantitatively analyze the Type I, Type II errors and their relationship.Type I errors are the errors which wrongly identified terrain points as offterrain points, and Type II errors are the errors which wrongly identified off-terrain points as terrain points.(Sithole and Vosselman, 2003) Figure 3. Type I errors, Type II errors, Accuracy Rates and Kappa coefficients of the 15 sample sites from ISPRS tested by the proposed MTF method

Quantitative Analysis
The 15 sample sites from ISPRS WG III/3 are selected on city and forest areas, the ground features such as slope gradient, vegetation density are various.Therefore, Meng et al. (2010) divided the 15 ISPRS study sites into three groups.The sites in the first group (Samples 11, 24, 41, 54) have rough slope and dense vegetation; the sites in the second group (Samples 12,21,22,23,31,42) are relatively flat urban area; and the sites in the third group (Samples 51, 52, 53, 61, 71) contain rough terrain and discontinuous (e.g., river banks and mining fields).The following discussion will refer to these three groups.
Figure 3 shows the quantitative assessment results of the 15 sites, while Table 2 lists the parameters used to generate these results.The overall accuracy and Kappa coefficient for one site may be required from tests with different combination of parameters.For example, "samp11k" in Table 2 refers to the parameter combination for the Kappa coefficient value of Sample 11 showed in Figure 3, the numbers of Type I, Type II errors and accuracy for Sample 11 are generated by the parameters listed as "samp11".
As shown in Figure 3, the average, best, worst values of the accuracy rate are 85%, 96% and 70% respectively.The standard deviation of the accuracy rate and Kappa coefficient in fifteen sites are 8% and 25%, which means the overall accuracy is relatively stable while the Kappa coefficient varies depending on the study sites.However, the parameters in the proposed MTF method have to be tweaked to obtain the best results during the finite number of experiments, and the optimal result is not guaranteed in these experiments.In order to analysis the performance of the proposed MTF method on different situations, a series of comparisons are carried out as follows.Since the 15 sample sites are divided into three groups, a comparison is shown in Figure 4. Group 2 sites shows the lowest errors and highest accuracy rate and kappa coefficient, which means that the MTF method can handle Group 1 sites better than the other two groups, this number is also good enough to compare with filters compared by ISPRS (Sithole and Vosselman, 2003).The performance of the MTF method on Group 2 is average.However, the performance on Group 3 shows a very low Kappa coefficient because of the high Type II errors.Group 3 sites contain features like steep slope and high percentage of terrain points.Therefore, the MTF method probably has flaw on process this type of areas.The ISPRS dataset is originally sorted as city sites and forest sites, the performance of the MTF method on city sites and forest site are shown in Figure 5.It is obvious that the performance on city site is better since it has lower Type I, II errors and higher accuracy and kappa.The steep slope and the dense vegetation coverage might be the reason why the MTF method has an unsatisfactory result on forest sites.The forest sites are basically overlapping with the Group 3 sites.The high Type II error is probably from the buildings on the slope which is a difficulty mentioned by Sithole and Vosselman (2003).It is also the key to improve the value of the Kappa coefficient.The MTF is based on layering, which is a global analysis of the data height value.Therefore, the terrain points' portion of all points can affect the result.Figure 6 demonstrates the MTF performance based on different terrain point percentage.It seems along with the growth of the percentage, the errors especially Type II error become higher, while the accuracy and kappa become lower.But it needs to be noticed that there is only one sample for the terrain point percentage smaller than 40%, and the sample which have higher than 80% terrain points are all in the Group 3. Therefore, the uncertainty of this feature still requires further discussion.In order to know which feature of the data has more influence to the result, the average standard deviations of the previous three types of sortation are calculated as shown in Figure 7.The chart shows that they are all in the same range of each characteristic; however, the group sortation has the lowest average standard deviation among the three types of sortation.To compare with the algorithms analyzed by ISPRS WG III/3, the error rate and average of Kappa coefficients are shown in Figures 8 and 9, respectively.Unfortunately, the error rates are the worst of the four comparing method in nine of fifteen sites.However, it has better results than Roggero's method in rest six sites.Similarly, the average of Kappa coefficients chart shows a 61.2% of the proposed MTF method, which is the 6th of all 9 methods, only higher than Elmqvist, Brovelli and Sithole's methods.Therefore, a further improvement of the MTF method is required.
Figure 9. Average of Kappa Coefficients in 15 sites of MTF method and eight method tested by ISPRS (Meng et al., 2009)

CONCLUSIONS
This paper has presented an automatic method called multiscale terrain filtering to generate DTM from last return points of high resolution airborne LiDAR point clouds.The developed method consists of three main steps.In the first step, outliers and noise points are eliminated.The method separates the point clouds into several layers based on the distribution of the elevation value of the points, and layer numbers are assigned to the points.In the second step, rasterized pyramid levels are generated from the lowest points in each cell.Then a series of iterative identifications and interpolations are processed to generate a rough DTM.The identification is comparing the layer numbers of the points with the reference points.The interpolation replaces the layer number and height value of offterrain cells by the average value of their neighbors.The last step is to refine the DTM.The terrain points are adjusted by comparing with the nDSM and then separated from the original data by the generated rough DTM.By using these identified terrain points, an IDW interpolation is processed to produce a final DTM.To verify the effectiveness of the developed method, ISPRS data sets with eight study sites and fifteen samples were used in this study.ISPRS also provided the results of eight existing algorithms for this data set.The result of the proposed multi-scale terrain filtering method indicates that it works as well as other filters in the flat terrain or terrain with gentle slopes.The proposed method can also overcome the difficulties like bridges, complex scenes, outliers and vegetation.However, the performance of the proposed method drops very much when handling the steep slope or discontinuities of the terrain.The total accuracy of the proposed method can be higher than 90% in some samples; however, it can be as low as around 65% in one study site.And the average Kappa coefficient in all fifteen study sites is 61.2%, which is low than average performance of all tested algorithms.

Figure 1 .
Figure 1.Flowchart of the MTF method

Figure 4 .
Figure 4. Average values of Type I, Type II errors, Accuracy and Kappa sorted by three groups

Figure 5 .
Figure 5. Average values of Type I, Type II errors, Accuracy and Kappa sorted by City Sites and Forest Sites

Figure 6 .
Figure 6.Average values of Type I, Type II errors, Accuracy and Kappa sorted by percentage of terrain point

Figure 7 .
Figure 7. Average values of Standard Deviations of Type I, Type II errors, Accuracy and Kappa in three types of sortation

Table 1 .
The proposed method was developed on C++ by Visual Studio 2008.The morphological filter and adaptive TIN filter compared in this research are included in the ALDPAT Version 1.0, which was developed by the International Hurricane Research Center, Florida International University in 2007.The final IDW interpolation and accuracy evaluation are processed on ArcGIS 10.The processor of the computer is equipped with Intel Core2 Duo CPU T5800 @ 2.00 GHz and 4 GB RAM.International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XXXIX-B4, 2012 XXII ISPRS Congress, 25 August -01 September 2012, Melbourne, Australia Features of the ISPRS dataset