Multi-Level Interpolation-Based Filter for Airborne LiDAR Point Clouds in Forested Areas

Over the past decades, plenty of filtering algorithms have been presented to distinguish ground and non-ground points from airborne LiDAR point clouds. However, with the existing methods, it is difficult to derive satisfactory filtering results on rugged terrains with dense vegetation due to the low-level penetration ability of laser pulses. Therefore, a multi-level interpolation-based filter is developed in this paper. The novelty of the algorithm lies within its usage of multi-scale morphological operations and robust z-score to correctly select ground seeds as more as possible, and a terrain-adaptive residual threshold to adapt to various terrain characteristics. Rural samples provided by International Society for Photogrammetry and Remote Sensing (ISPRS) were employed to assess the performance of the proposed method and its results were compared with 15 filtering algorithms developed in recent 5 years (2015-2019). Results show that the proposed method with the optimized parameters produces the best accuracy with the average total error and kappa coefficient of 1.89% and 87.88%, respectively. We further filtered high-density point clouds in six forested areas with different vegetation covers and terrain slopes. Results demonstrate that the proposed algorithm is more accurate than the well-known filtering methods including morphological-based filter, progressive TIN densification filter (PTD), improved PTD and cloth simulation filter, with the average total error decreased by 26.2%, 19.9%, 3.8% and 40.4%, respectively. Moreover, the DEMs of the proposed method have lower average root mean square errors than the four classical filters. Therefore, the proposed method can be considered as an effective ground filtering algorithm for airborne LiDAR point clouds in forested areas.


I. INTRODUCTION
Nowadays, airborne light detection and ranging (LiDAR) with the ability of tree canopy penetration and little influence by shadows and surrounding light conditions has become a valuable and effective tool in forest inventory, such as individual-tree based height determination and volume estimation, and tree-species classification [1]- [3]. Generally, the extraction of ground points from raw airborne LiDAR point clouds is the pre-requisite and most important step for almost all LiDAR applications. This process is termed as filtering [4]. In the last two decades, plenty of filtering The associate editor coordinating the review of this manuscript and approving it for publication was Jenny Mahoney.
algorithms have been proposed. According to their working principles, the existing filtering algorithms can be divided into six broad categories: slope-based, morphology-based, interpolation-based, statistical-based, machine learningbased and segmentation-based approaches [5], [6].
Generally, these filters achieve a good performance in relatively flat terrains with simple landscapes, but show poor results in mountainous terrains, cliffs, ridges and discontinuities with complex environments [7]. Compared to the other filtering methods, the interpolation-based filters commonly obtain better results since they can use more contexts [5]; nevertheless, they cannot perform well in dense forested areas, because the canopy prevents most laser pulses from reaching the ground surface, resulting in few ground points beneath the canopy [8]. Although some specific methods [9]- [12] have been designed for forested areas, their performances are still unsatisfactory and require to be improved [13].
To solve the aforementioned problems, we propose a multilevel interpolation-based filter for airborne LiDAR data in complex forested areas in this paper. Compared with the existing interpolation-based filters, the proposed method has three main contributions: • Multi-scale morphological operations are used to select potential ground seed points. These points are more than those derived by the local minimum method with predefined grids, which is commonly used for the detection of initial ground seeds [4], [14]. Thus, detailed terrain features can be captured with the so many ground seeds, especially in dense forested areas.
• The non-ground points falsely recognized as the potential ground seeds are removed by robust z-score, which guarantees the accuracy of the initial reference DEM.
• An adaptive residual threshold with the consideration of terrain slope is introduced, which assures the adaptability of the filtering in forested areas with rugged terrains. The remainder of this paper is organized as follows. Related works are reviewed in Section II. The proposed method is detailed in Section III. In Section IV, experiments are designed to evaluate the proposed method. Section V comparatively analyzes the results of the proposed method. Finally, discussion and conclusions are given in Sections VI and VII, respectively.

II. RELATED WORKS
Slope-based filters classify the point as the non-ground point when its slope is greater than a threshold [15]. In general, this method obtains good results in gentle slope areas but cannot perform well in steep slopes. To overcome this problem, Wang and Tseng [16] proposed an adaptive dual-directional slope filter, where the filtering results of different directions are complementary to each other, so that over-filtering problem can be avoided. Susaki [17] proposed an adaptive slope filter, where the slope parameter adopted in the algorithm is updated after an initial production of the digital terrain model (DTM). Although the slope-based filters are theoretically simple, they are very sensitive to the slope threshold.
Morphology-based filters are based on the assumption that elevation differences of the non-ground points before and after the opening operations are greater than a threshold. However, their success significantly depends on the selection of an appropriate structure element (SE). Thus, Zhang et al. [18] proposed a progressive multi-scale morphological filter by gradually increasing SE. Pingel et al. [19] presented an improved simple morphological filtering algorithm (SMRF). It uses a linearly increasing window and simple slope thresholding with the combination of an image inpainting algorithm. Li et al. [20] developed a filtering algorithm based on geodesic transformations of mathematical morphology, where the transformations only use the elementary SE and converge after a finite number of iterations. Overall, since morphology-based filters work on rasterized data, they have a highly computational efficiency. Yet, the transformation from raw point clouds into the gridded data inevitably causes accuracy loss.
The interpolation-based filters first produce a reference surface using certain interpolation method [21], [22] and then gradually select more and more ground points. For example, Axelsson [23] proposed a progressive triangulated irregular network (TIN) densification (PTD) filtering method, where a point with the angle and distance to the corresponding triangle less than the thresholds is labeled as the ground point. Evans and Hudak [12] proposed a multi-scale curvature classification method, which used regularized spline with tension method instead of TIN interpolation to construct raster surfaces. Chen, et al. [14] developed a multi-resolution hierarchical classification algorithm (MHC), where the reference ground surface is produced with thin plate spline [24]. Hu et al. [25] presented an adaptive surface filter, which used bending energy to realize adaptive filtering threshold. However, the interpolation-based filters suffer from the problem of filtering out terrain features, such as discontinuities.
The statistical-based filters are based on the assumption that natural ground points follow a normal distribution and the raw point clouds can be seen as a mixture of Gaussian models. Considering the fact that the mixture of ground and non-ground points makes the skewness of the point clouds positive, Bartels et al. [26] proposed a skewness balancing filter (SBF), which was achieved by iteratively removing the highest point until the skewness of the remaining points is equal to zero. Bao et al. [27] updated SBF by the additional use of a statistical measure, i.e. kurtosis. Bartels and Wei [28] further improved this method to adapt to slope terrains. Crosilla et al. [29] employed the sequential skewness and kurtosis analysis of elevation and intensity point distribution values to filter and classify point clouds. Taking the filtering as a separation of mixed Gaussian models, Hui et al. [30] proposed a threshold-free algorithm based on expectationmaximization. However, the statistical-based filters are prone to classify the non-ground points as the ground points when the number of the former are greater than that of the latter [28].
To improve the robustness of the filtering algorithms, many machine learning-based filters have been proposed. For example, Lu et al. [31] classified the LiDAR points as being either ground or non-ground using supervised learning techniques with a variety of features. Jahromi et al. [32] proposed a novel filtering algorithm based on artificial neural networks to select ground points. Hu and Yuan [33] developed a filtering method based on deep learning using deep convolutional neural networks. Rizaldy et al. [34] used fully convolutional networks to classify point clouds.
Hui et al. [35] proposed an active learning filtering method, where the initial training samples are automatically prepared by a morphological-based filter and then a support vector machine is employed to filter the point clouds. Nevertheless, Segmentation-based filters first divide the point cloud into a set of segments, and then classify the segments according to certain algorithms. For example, Tovari and Pfeifer [36] employed a weighted moving least squares method to classify segments. Lin and Zhang [37] classified the segments using the PTD method. Chen, et al. [4] used the MHC to filter the segments. Yang et al. [38] filtered the segments according to the elevation differences between boundary points of the segment and the provisional DEM. Nevertheless, the performance of segmentation-based filters heavily depends on the quality of segmentation.
In addition, to incorporate advantages of different filtering methods, some researchers suggested the combination of different filters. For example, Podobnikar and Vrečko [39] first divided the study area into different subareas according to the characteristics of the terrain, and then the point clouds in the subareas were classified based on different filters. Deng and Shi [40] integrated PTD and hierarchical robust interpolation. Yang, et al. [38] employed segment-based and multi-scale morphological filters to classify segmented points and one set of individual points, respectively. Zhao, et al. [13] proposed an improved PTD (IPTD) filter for dealing with a variety of forested landscapes, where the initial ground seeds are first obtained by the morphological method and then densified by the progressive TIN. Bigdeli, et al. [3] integrated the results of slope-based and morphological-based filters for DEM extraction in dense forested areas. Cai et al. [41] developed a novel filtering algorithm that combines cloth simulation (CS) and PTD, where an initial DEM is obtained by CS and the parameter thresholds of the PTD are derived from the initial DEM. Finally, the PTD with the adaptive parameter thresholds is used to update the initial DEM. Although the fusion of different filtering methods could improve the results, the overall accuracy significantly relies on the performance of each algorithm. Firstly, potential ground seeds are selected with morphological opening operations. Then, non-ground points mixed in the potential ground seeds are removed by robust z-score under a statistical manner. Hence, true ground seeds are obtained. Finally, an adaptive multi-level interpolation-based filter is developed to densify the ground points, where the initial reference DEM is interpolated by thin plate spline (TPS) on the true ground seeds.

A. SELECTING GROUND SEEDS
To select potential ground seeds, the raw point cloud is first gridded to form the minimum surface (z min ). Then, the morphological opening operation is applied to the z min . The grids are flagged as the non-ground grids if their elevation differences before and after the opening operation exceed a threshold. The specific steps are as follows: Step1: Gridding point cloud The point cloud is first covered with grid cells, and the elevation of the lowest point in each cell is taken as the grid value (Fig. 2). This is because the lowest points are more likely to be ground points. Note that the size of the grid cell depends on the data density. Specifically, the denser the point cloud, the smaller the grid cell size. In general, the size of the grid cell equals to the average point space. For the empty cells where no points are located, their values are filled by the spring-metaphor inpainting algorithm [42]. After the above operations, the minimum surface (z min ) is obtained, which is taken as the input of the following morphological opening operation.
Step2: Morphological opening operation The opening operation, which consists of erosion followed by dilation, is sensitive to the size of the structuring element (SE), since a large SE fails to remove low vegetation points whereas a small SE wrongly removes terrain features. To this end, a multi-scale morphological opening operation developed by Pingel, et al. [19] is used here, where SE and residual threshold gradually increase. Note that the surface after each opening operation is taken as the input of the next operation (Fig. 3). For each opening operation, the points with the elevation differences exceeding the residual threshold are flagged as non-ground points. Thus, after all the operations, the points without being flagged as non-ground points are taken as the potential ground seeds.
Step3: Removal of non-ground points from the potential ground seeds It is inevitable that the potential ground seeds contain nonground points (here, low outlier is also considered as a nonground point), which must be removed before the production of initial ground surface. In statistics, the non-ground point can be considered as an outlier, which can be detected by the z-score (z). It is expressed as where x i is the elevation of ith point,x is the mean value, and σ is the standard deviation. It is well known that in a normal distribution, 99.4% of the observations fall within 2.5 standard deviations from the mean, and the z-scores of the observations follow the standard normal distribution with the mean of 0 and standard deviation of 1. Thus, it is reasonable that x i is flagged as an outlier if |z i | ≥ 2.5. In robust statistics [43], breakdown point is the percentage of outliers that an estimator can resist before giving wrong results. The mean estimator has a breakdown point of 0%, since a single outlier can destroy it. Likewise, standard deviation also has a breakdown point of 0%. Thus, the traditional z-score has zero breakdown point, thereby tending to eliminate true ground points due to the swapping effect [44]. The median estimator, in turn, has a breakdown point of 50%, since it can resist 50% of outlier observations to disturb it. Thus, the median is a robust alternative to the mean. Similarly, there is a robust alternative to the standard deviation estimator, termed median absolute deviation (MAD). Like the median, MAD also has a breakdown point of 50%. Based on the two estimators, a robust z-score (rz) is developed, where the mean and standard deviation estimators in the traditional z-score are respectively replaced with the corresponding robust estimators, namely, where MAD is expressed as In the context of non-ground point detection, the robust z-score of each point in the candidate ground seeds is estimated based on its k nearest neighbors (e.g. k = 12). Thus, x i is flagged as a non-ground point and then removed from the ground seeds if |rz i | ≥ 2.5. After the removal, the true ground points are obtained. Note that some abrupt terrain points may be wrongly removed. However, these removed points could be restored based on the slope-adaptive threshold, which will be discussed in following Section.

B. ADAPTIVE MULTI-LEVEL INTERPOLATION-BASED FILTER (MIF)
After the selection of true ground seeds, the raw point cloud is classified into two groups: ground points and candidate points. Thus, in the following steps, the other ground points mixed in the candidate points are to be selected by the adaptive MIF.
The adaptive MIF includes m levels (e.g. m = 3), where the reference DEM resolution and the residual threshold gradually increase from the low to the high levels. In the ith level VOLUME 8, 2020 (i = 1, 2, . . . , m), the steps for the selection of ground points is as follows: Step1: Produce a reference DEM with the resolution of h i by TPS on the ground seeds. Step2: Compute the height residuals of each candidate point to the 9 surrounding grid cells of the reference DEM. If 4 of 9 residuals are less than the threshold, the point is flagged as the ground point. In the classical MHC method [14], the residual threshold only considers the elevation difference, so it is prone to misclassify points on steep slopes [19]. To overcome this problem, a slope-adaptive threshold (t) is adopted, which is expressed as t = t i + s f × slope, where t i is a fixed threshold value, s f is a slope scale and slope is the terrain slope derived from the interpolated DEM. As shown in Fig. 4, with a fixed threshold, many ground points on steep slopes are misclassified as the non-ground points. However, this misclassification is avoided by the slope-adaptive threshold. Step3: Update the ground points and the ground seeds based on the selected points in Step2. The ground seeds are the local lowest points in the grid cells with the resolution of h i . Step4: Repeat steps1-3 until no ground points are found.
Step5: Go to the (i + 1)th level and repeat steps1-4 until the highest level is reached, where h i+1 = h i /2 and t i+1 = t i + 0.1. It can be found that the adaptive MIF has three key userdefined parameters, namely, initial DEM resolution (h 1 ), initial fixed threshold value (t 1 ) and slope scale (s f ). Note that the code of our method was written using MATLAB R2018b software on a Windows 7 operating system with the hardware configuration of Intel Core i5-3470 CPU @ 3.2 GHz and 4.0 GB RAM.

IV. EXPERIMENTAL DESIGN A. DATASETS
The benchmark dataset provided by the International Society for Photogrammetry and Remote Sensing (ISPRS) are employed to assess the performance of the proposed method.  Table 1 presents detailed information about each sample regarding number of ground points (#G), number of non-ground points (#NG), standard deviation of the elevations (STD), and landscape features.
Since the ISRPS benchmark dataset was published about two decades ago, they have low sample density due to the limitation of hardware devices in that era. Thus, 6 forested sites with high sample density are further used to evaluate the performance of the proposed method. Considering the fact that vegetation cover and terrain slope have a serious effect on the success of filtering algorithms [45]- [48], we selected the 6 sites covered with different terrain characteristics and vegetation covers (Fig. 5). This can fully assess the performance of the proposed method under different landscapes. The area for each site is 25 ha. Detailed data collection summaries of the 6 sites based on airborne LiDAR systems are shown in Table 2.
In addition, some statistics of the study sites are shown in Table 3. It can be found that the 6 sites have the sample density ranging from 8.9 to 15.6 pts/m 2 , and include various terrain characteristics with the mean terrain slopes from 16 • to 37 • . Furthermore, the vegetation conditions change greatly with the mean canopy covers ranging from 41% to 98%, and mean tree heights from 2.6 to 12.8 m. The reference data for the 6 sites were generated using commercial LiDAR data processing software named TerraScan along with manual editing by the industry data processing experts. The reference DEMs were produced by interpolating the reference ground points using natural neighbor, as shown in Fig. 6.

B. COMPARATIVE METHODS
Since the publication of the ISPRS dataset, they have been widely employed by many researchers to assess the accuracy of their methods. Thus, we compare the performance   of the proposed method with those of 15 filtering algorithms developed in recent 5 years (2015-2019) [13], [20], [30], [38], [49]- [59], which also took the ISPRS dataset as the benchmark data. Moreover, because the proposed method can be considered as an improved version of the classical MHC [14], the results of MHC are also compared.
For the forested data, the proposed method is compared with four well-known filtering algorithms including VOLUME 8, 2020 morphological filter (MF), progressive TIN densification (PTD), improved PTD (IPTD) [13] and Cloth Simulation Filter (CSF) [60]. MF and PTD were respectively implemented in the publicly available software ALDPAT (http://lidar.ihrc.fiu.edu/lidartool.html) and Ter-raScan (http://www.terrasolid.com/home.php). IPTD was implemented with the private source code provided by the first author of the paper [13], and CSF was conducted with the publicly available MATLAB code (https: //www.mathworks.com/ matlabcentral/fileexchange/58139-c sf-cloth-simulation-filter /). The best results of the filtering methods are obtained by the cautious tuning of their parameters with the trial-and-error skill. Namely, the best result corresponds to the minimum total error.

C. ACCURACY MEASURES
Four accuracy measures including Type I errors (T.I), Type II errors (T.II), total errors (T.E) and kappa coefficient (κ) are used to quantitatively evaluate the filtering performance [5]. A cross-matrix of filtering result is shown in Table 4. Based on the cross-matrix, the four accuracy measures are respectively expressed as In addition, the DEMs of all the filtering methods on the forested datasets are produced by interpolating their filtered ground points using natural neighbor and then their root mean square errors (RMSEs) are computed to further assess the filtering performance. RMSE is expressed as where e i is the elevation difference between the interpolated and the reference DEMs (Fig. 6) at the ith grid cell and n is the total number of grid cells in the DEM.

V. RESULTS AND ANALYSIS A. ISPRS DATASET
To assess the average filtering performance, the 6 ISPRS samples were filtered by the proposed method with one single set of parameters. Through trial-and-error, the best parameter settings are obtained in terms of the minimum average kappa coefficient (Table 5). That is, the initial resolution (h 1 ) is 2 m, initial elevation threshold (t 1 ) is 0.29 m and slope scale (s f ) is 1. Furthermore, we also derive the highest kappa coefficient with the optimized parameters for each sample. Table 5 illustrates that with one single set of parameters, the average kappa coefficient and total error of the proposed method are 82.31% and 3.17%, respectively, while with the optimized parameters, the kappa coefficient and total error are 87.88% and 1.89%. Thus, the difference between the total errors for the single set of parameters and the optimized ones are less than 2%. This indicates that there is no obvious accuracy difference between the two groups of filtering results, which shows the insensitivity of the proposed method to the parameter settings.
Regardless of parameter settings (Table 5), samp51 and samp61 obtain the best results with respect to kappa coefficient and total error, respectively. The kappa coefficient of samp53 is much smaller than the others due to the steep slopes in this area (Table 1). In fact, this sample also puzzles the other filtering methods [19]. Yet, except for samp53, the kappa coefficients of the other samples with the optimized parameters are always greater than 84% and the total errors less than 3%. This demonstrates the robustness of the proposed method on the rural samples.
It can be found that the type II errors exceed the type I errors in all samples (Table 5). Thereinto, samp53 shows the largest difference between the type I and II errors while samp54 shows the smallest. For example, with the optimized parameters, the type II error of samp53 is 26.98 times larger than the type I error while the type II error of samp54 is 1.03 times larger than the type I error. On average, the type II error is about 773% larger than the type I error.
The computational cost of the proposed method mainly depends on the number of input points (Table 5). For example, samp54 with the least number of points (Table 1) has the lowest computing cost while samp61 with the largest number of points (Table 1) has the highest cost. On average, the proposed method requires 137.11 s to filter the point clouds. Thus, although the proposed method aims to improve the filtering quality of airborne LiDAR point clouds, its processing cost is acceptable due to the limited hardware configuration. Note that the processing cost of a filtering method also relies on the programming code (e.g. C++, C#, python and MATLAB). Hence, it is difficult to compare the computing efficiency of filtering methods under different running environment.   7 shows the distributions of the type I and II errors of the proposed method. It is obvious that the type I and II errors are mainly located on steep slopes, e.g. samp52 (Fig.7b) and terrain discontinuities, e.g. samp53 (Fig. 7c). In addition, parts of bridges with the attachment to the ground are incorrectly accepted as the ground points (Fig. 7f). Overall, the type I and II errors are relatively few and terrain features are well retained, especially for samp51, 54 and 61 (Figs. 6a, d and e). Table 6 gives an accuracy comparison between the proposed method, MHC and the recently developed 15 filtering methods with respect to total error. The proposed method VOLUME 8, 2020  produces the lowest errors in 4 of 6 samples (including samp51, 53, 54, and 61). However, the total errors of the remaining samples (including samp52 and 71) are approximate to the best results. Specifically, the proposed method in samp52 and 71 is only 11.6% and 22.6% less accurate than the corresponding best filters. On average, the proposed method has the total error of 1.89%, whereas the other filters have the errors ranging from 2.92% to 13.56%. Namely, the proposed method is at least 35.3% more accurate than the classical filters. Thus, the proposed method obtains promising and reliable results for the rural landscapes. Table 7 exhibits an accuracy comparison between the proposed method and the well-known filtering methods in the forested sites. We can see that for the proposed method, the type II errors are obviously larger than the type I errors in almost all samples, except for Data1. This conclusion is consistent with that derived from the ISPRS dataset. Compared with the classical filters, the proposed method achieves the lowest type I errors on 3 out of 6 samples (i.e. Data4, Data5 and Data6) and the lowest type II errors on 2 samples including Data1 and Data4. With respect to total error and kappa coefficient, all the methods produce more accurate results on Data4-6 than on Data1-3. This is because the latter are characterized by low vegetation on steep slopes (Table 3), which tend to puzzle the filtering methods. In terms of total error and kappa coefficients, the proposed method outperforms the other methods on almost all samples, except for Data2. On average, the proposed method with the total error of 6.80% and kappa coefficient of 85.61% achieves the best result, which is followed by IPTD, and CSF produces the worst results. More specifically, compared to MF, PTD, IPTD and CSF, the proposed method reduces the average total error by 26.2%, 19.9%, 3.8% and 40.4%, respectively. The DEM RMSEs of all the methods for each site are given in Table 8. The proposed method performs better than the other methods in 4 of 6 sites (i.e. Data1, Data3, Data5 and Data6), while IPTD has the best results on the other two sites (i.e. Data2 and Data4). On average, the proposed method is 46.9%, 39.7%, 35.3% and 46.5% more accurate than MF, PTD, IPTD and CSF, respectively. Fig. 8 shows the DEMs of the proposed method and the classical filtering methods on Data1. Results demonstrate that the proposed method (Fig. 8b) produces a similar DEM to the reference (Fig. 8a). In comparison, the classical methods including MF (Fig. 8c), PTD (Fig. 8d) and CSF (Fig. 8f) tend to preserve low vegetation on rugged terrain, thereby resulting in the production of unnatural coarse surfaces, such as denoted by the circles. The surface of IPTD (Fig. 8e) shows a good appearance to the reference, yet has a loss of terrain FIGURE 9. Error distributions of the proposed method (a) without and (b) with using morphology-based method for obtaining potential ground seeds on samp54 ( is the type I error and is the type II error). detail on the steep slope, as denoted by the ellipse. In contrast, CSF suffers from a serious loss of surface details on the rugged terrains, like those denoted by the ellipses (Fig. 8f).

VI. DISCUSSION
Theoretically, the proposed method belongs to the interpolation-based filters. However, compared to the classical interpolation-based algorithms, the new method has three main contributions, namely, multi-scale morphology-based method for obtaining potential ground seeds, robust z-scorebased method for removing non-ground points and slopeadaptive threshold for filtering point clouds. Thus, to further verify the performance of the proposed method with and without each contribution, three tests are conducted on the ISPRS dataset.
Taking samp54 as an example, we assess the effect of the morphological-based method on the detection of potential ground points. Results (Table 9) indicate that the proposed method with the improvement shows a clearly higher accuracy than that without the improvement. Specifically, compared to the latter, the former decreases the total error by 56.7% and increases the kappa coefficient by 5.5%. Fig. 9 shows that the proposed method without the improvement suffers from obvious type II errors on rugged terrain. Thus, more ground seeds can greatly improve the filtering accuracy. with using robust z-score for removing non-ground points on samp53 ( is the type I error and is the type II error).

FIGURE 11.
Error distributions of the proposed method (a) without and (b) with using slope-adaptive residual threshold on samp52 ( is the type I error and is the type II error).
We further assess the effect of the robust z-score-based method on the removal of non-ground points on samp53. Results (Table 9) demonstrate that the proposed method with the improvement performs better than that without the improvement. Nevertheless, the accuracy difference is slight. Specifically, the differences of total error and kappa coefficient are 0.07% and 2.97%, respectively. This indicates that the morphological-based method has a good ability of accurately detecting potential ground seeds. Fig. 10 demonstrates that the two error distributions are similar, except for that denoted by the ellipse.
Finally, the effect of the slope-adaptive residual threshold is evaluated on samp52. Compared to the method without the improvement, the proposed method with the improvement decreases the total error by 59.6% and increases the kappa coefficient by 22.9% (Table 9). Moreover, the type I errors are significantly reduced on the steep slopes (Fig. 11). Therefore, the slope-adaptive residual threshold has a great effect on the filtering performance on steep slope areas.
In conclusion, the three contributions make the filtering method adaptive, accurate and robust with respect to various terrain characteristics and vegetation covers.

VII. CONCLUSION
To improve the filtering accuracy of airborne LiDAR data in forested areas, we developed a multi-level interpolationbased filtering algorithm in this paper. The proposed method uses morphological operations and robust z-score to accurately detect ground seeds. Moreover, a terrain-adaptive residual threshold is developed to accommodate for different terrain characteristics. We demonstrated the effectiveness of our method by conducting a detailed comparison with the most recent and with some well-known approaches for filtering point clouds. The methods were respectively evaluated on ISPRS rural samples and 6 forested datasets with various terrain slopes and vegetation covers. The proposed algorithm achieved the average total error of 1.89% and average kappa coefficient of 87.88% on the ISPRS rural samples, while the average total error and kappa coefficient on the forest datasets are 6.80% and 85.61%, respectively. In comparison, our method performs better than the competitors, demonstrating its robustness and reliability. It should be noted that the proposed method have much larger type II errors than type I errors. Namely, some non-ground points are mixed in the ground points. Thus, our further researches will focus on the decrease of type II errors under the help of other simple filtering methods, such as slope-based approaches.