Comparison and Evaluation of Different Pit-Filling Methods for Generating High Resolution Canopy Height Model Using UAV Laser Scanning Data

: As a common form of light detection and ranging (LiDAR) in forestry applications, the canopy height model (CHM) provides the elevation distribution of aboveground vegetation. A CHM is traditionally generated by interpolating all the ﬁrst LiDAR echoes. However, the ﬁrst echo cannot accurately represent the canopy surface, and the resulting large amount of noise (data pits) also reduce the CHM quality. Although previous studies concentrate on many pit-ﬁlling methods, the applicability of these methods in high-resolution unmanned aerial vehicle laser scanning (UAVLS)-derived CHMs has not been revealed. This study selected eight widely used, recently developed, representative pit-ﬁlling methods, namely ﬁrst-echo interpolation, smooth ﬁltering (mean, medium and Gaussian), highest point interpolation, pit-free algorithm, spike-free algorithm and graph-based progressive morphological ﬁltering (GPMF). A comprehensive evaluation framework was implemented, including a quantitative evaluation using simulation data and an additional application evaluation using UAVLS data. The results indicated that the spike-free algorithm and GPMF had excellent visual performances and were closest to the real canopy surface (root mean square error (RMSE) of simulated data were 0.1578 m and 0.1093 m, respectively; RMSE of UAVLS data were 0.3179 m and 0.4379 m, respectively). Compared with the ﬁrst-echo method, the accuracies of the spike-free algorithm and GPMF improved by approximately 23% and 22%, respectively. The pit-free algorithm and highest point interpolation method also have advantages in high-resolution CHM generation. The global smooth ﬁlter method based on the ﬁrst-echo CHM reduced the average canopy height by approximately 7.73%. Coniferous forests require more pit-ﬁlling than broad-leaved forests and mixed forests. Although the results of individual tree applications indicated that there was no signiﬁcant difference between these methods except the median ﬁlter method, pit-ﬁlling is still of great signiﬁcance for generating high-resolution CHMs. This study provides guidance for using high-resolution UAVLS in forestry applications.


Introduction
As an active remote sensing technology, light detection and ranging (LiDAR) can penetrate the canopy to obtain the vertical structure of a forest, and it has been widely used in forest inventory analyses [1,2]. LiDAR has two common data formats: one is the original discrete point cloud, and the other is the spatially continuous raster surface obtained from raw laser points, such as the digital elevation model (DEM) and digital surface model (DSM) [3]. The canopy height model (CHM) is usually constructed by subtracting the DEM from the DSM or by interpolating the normalized point cloud, which is a direct manifestation of the absolute height distribution of the vegetation canopy above the ground [4]. Compared with the original point cloud data, the storage volume of the canopy surfaces [27]. Khosravipour et al. [17] proposed a spike-free method considering all laser returns for generating DSMs, and they found that compared with using only the first return DSM, the accuracy of treetop detection, especially for small trees, was significantly improved. Hao et al. [27] developed a canopy surface point filtering algorithm called graph-based progressive morphological filtering (GPMF), and they found that the CHM generated with the GPMF method produced few pits while retaining canopy details. These methods usually require several parameters to drive the model, and the effect of applying them to specific data sets requires further evaluation.
In recent years, as an emerging remote sensing dataset, unmanned aerial vehicle laser scanning (UAVLS) data have been widely used in the estimation of forest canopy structure and forest parameters in small-scale forest inventories [28,29]. UAVLS can generate data with point densities of 100-300 points per square metre, or even up to 1000 points per square metre, representing a significant increase in point density compared with the data provided by airborne laser scanning (ALS) [30]. High-density point clouds create high-resolution CHMs because the CHM resolution is usually determined by considering the average point cloud density [31]. A high-resolution (~10 cm) CHM can describe the morphological structure of tree crowns in detail [17,32], and more small trees can be seen compared with lower resolution data [7]. However, the high-resolution CHM has more complex height variations [9] and will create more data noise [33,34]. Although a considerable amount of literature focuses on the pit-filling method, to our knowledge, most UAVLS studies ignore this point, and few studies have used the first return [35] and highest point interpolation techniques [33,36,37] to generate CHMs. Moreover, the applicability of the abovementioned pit-free methods in a CHM derived from high-resolution UAVLS requires discussion.
The main purpose of this study is to compare and evaluate the performances of eight different pit-filling methods to generate high-resolution CHMs derived from UAVLS. The most widely used, recently developed and representative CHM generation methods are selected, including the first-echo interpolation method, smooth filtering (mean, medium and Gaussian) method, highest points interpolation method [12], pit-free algorithm [10], spike-free algorithm [17] and GPMF algorithm [27]. Furthermore, this study aimed to select the most representative raster to describe the upper canopy surface to further improve the accuracy of using high-resolution CHM data in individual tree detection and tree height estimation. Finally, we provide guidance for improving the quality of UAV-LiDAR-derived CHMs in forest inventory applications.

Materials
To quantitatively and accurately evaluate the differences between different CHM generation methods and the real canopy surface, both simulated and real-world UAVLS point clouds were employed for testing.

Simulated Data
We used two geometric models to generate simulated point clouds: cone and hemisphere [38]. In the range of 50 m × 50 m, 50 cones and 50 hemispheres are randomly distributed. The geometric surface points are randomly generated according to Equations (1) and (2), and the average point space is 0.08 m.
Cone : z = l − x 2 + y 2 r/l , x 2 + y 2 ≤ r 2 (1) Hemisphere : z = |r 2 − x 2 − y 2 | , x 2 + y 2 ≤ r 2 (2) where x, y and z are the spatial coordinates of geometric surface points, r is the radius of the cone and hemisphere and l is the height of a cone. The r and l values are set to 1.5-3.5 and 2-6 m, respectively, and a basic height of 2-6 m is added to the simulated geometry. The average height of the cone is 8.6 m, while that of the hemisphere is 8.8 m. If the simulated geometry overlaps horizontally, the highest point in the vertical direction of the overlapped part is retained to ensure that the sunlit part is the crown surface. In addition, ground points with the same density are added to the non-projection area of the simulated geometry, and the z-value of these ground points is set to 0. Thus, two simulated surface points are created (see Figure 1a,b). Moreover, the surface points were interpolated to the reference CHM (CHM Reference ) for subsequent comparison. Then, different proportions (10-60%) of noise points are randomly added below the simulated surface to create pits with different sizes. The final simulated point clouds are shown in Figure 1c,d. It is worth noting that the simulated points do not follow the laser scanning sensor geometry, so the distributions of these under-canopy points will not be distributed as the simulated along vertical strata in real data.
Hemisphere: z = | − − | , + ≤ where , and are the spatial coordinates of geometric surface points, is the radius of the cone and hemisphere and is the height of a cone. The and values are set to 1.5-3.5 and 2-6 m, respectively, and a basic height of 2-6 m is added to the simulated geometry. The average height of the cone is 8.6 m, while that of the hemisphere is 8.8 m.
If the simulated geometry overlaps horizontally, the highest point in the vertical direction of the overlapped part is retained to ensure that the sunlit part is the crown surface. In addition, ground points with the same density are added to the non-projection area of the simulated geometry, and the z-value of these ground points is set to 0. Thus, two simulated surface points are created (see Figure 1a,b). Moreover, the surface points were interpolated to the reference CHM (CHMReference) for subsequent comparison. Then, different proportions (10-60%) of noise points are randomly added below the simulated surface to create pits with different sizes. The final simulated point clouds are shown in Figure 1c,d.
It is worth noting that the simulated points do not follow the laser scanning sensor geometry, so the distributions of these under-canopy points will not be distributed as the simulated along vertical strata in real data.

Study Area
The study area is the Maoershan Forest Farm, Shangzhi, Heilongjiang Province, Northeast China, ranging from 127°18′0″ to 127°41′6″E and 45°2′20″ to 45°18′16″N ( Figure  2). The slope ranges from 5° to 25°, the terrain is high in the south and low in the north and the average altitude above mean sea level is approximately 400 m. The forest type is a typical natural secondary forest of Northeast China that is mainly composed of precious broad-leaved forest, poplar birch forest and oak forest, as well as a small number of coniferous plantations, such as Korean pine, larch and Scotch pine. In this study, a total of 3   (Figure 2). The slope ranges from 5 • to 25 • , the terrain is high in the south and low in the north and the average altitude above mean sea level is approximately 400 m. The forest type is a typical natural secondary forest of Northeast China that is mainly composed of precious broad-leaved forest, poplar birch forest and oak forest, as well as a small number of coniferous plantations, such as Korean pine, larch and Scotch pine. In this study, a total of 3 sites, which represent three forest stand types (broadleaf forest, coniferous forest and coniferous and broad-leaved mixed forest), were established.

UAV-LiDAR Data
The UAV-borne LiDAR equipment used in this study was an extremely lightweight RIEGL mini VUX-1UAV LiDAR scanner (www.riegl.com/products/unmanned-scanning/ riegl-minivux-1uav, accessed on 8 June 2021) carried by a Feima D200 UAV platform (www.feimarobotics.com/en/productDetailD200, accessed on 8 June 2021). The laser scanner operated at a pulse repetition rate of 100 kHz with scan speeds up to 100 scans per second. The maximum measurement range is 250 m and the range measurement accuracy is 15 mm. The footprint size is 160 mm × 50 mm at 100 m and the beam divergence is 1.6 × 0.5 mrad. The field of view is up to 360 • and the angle measurement resolution is 0.001 • . The scanner has multiple target capabilities and can generate up to 5 target echoes per laser shot. In addition to the laser sensor, the UAV platform is also composed of a GNSS antenna, a high-precision inertial measurement unit (IMU) and a high-speed storage control unit (see Figure 2). Three parallel batteries ensure the safe landing although two out of three are off, which could support a 48 min hover.
UAVLS data were acquired in August 2019 at three sites. All flights were designed as crossing transects with 80 m swath overlaps at 80 m altitude and 5.0 m/s speed. The average point density for each site ranged from 150 to 250 pt/m 2 . The point clouds of these three sites are shown in Figure 2

UAVLS Data Pre-Processing
The raw UAVLS data are processed by a series of noise removal processes to remove many air points, points below the ground and isolated points in the canopy [29]. Air points and low points were removed manually, and isolated points were determined by the number of points inside the search neighbourhood for a given search radius (5 m). Subsequently, the remaining points were classified into ground and nonground points using the progressive triangulated irregular network (TIN) densification method developed by Axelsson [39]. Then, the ground points were interpolated into a digital terrain model (DTM) using kriging interpolation with a 1 m pixel size [40]. According to the DTM, the average slopes of the three sites were approximately 9 • , 11 • and 7 • . The normalized height of point clouds was obtained by subtracting the DTM value from the elevations of all points [41].

Description of CHM Generation Algorithms
The most widely used, recently developed and representative eight pit-filling methods are selected to generate CHMs, including the first-echo (FE) interpolation method, smooth filtering (mean, medium and Gaussian) method, highest point (HP) interpolation method [12], pit-free algorithm [10], spike-free algorithm [17] and GPMF algorithm [27]. The smooth filtering method and pit-free algorithm are post-processing methods, while the other four methods are pre-processing methods. Figure 3 shows the process of generating CHMs by these eight methods. FE interpolation is a common CHM generation method that is generated by interpolating all the first echoes in the normalized point cloud. Then, the mean filter, median filter and Gaussian filter are used to smooth the FE CHM. HP interpolation first generates 0.1 m grid cells and then extracts the highest point of each cell for interpolation. The pit-free method, spike-free method and GPMF method are described below.

Pit-Free Algorithm
The pit-free method is a semiautomated pit-filling algorithm that uses a user-defined threshold to detect pits and automatically fill them. Specifically, a Laplacian operator is first applied to the FE CHM and a percentage of the cumulative histogram of the Laplacian image is defined to determine the percent of pits. Then, the data pits are marked as one and the non-data pits are marked as zero to generate a binary mask. At the same time, a median filter (kernel 3 × 3) is used to smooth the FE CHM to generate a pit-filled CHM. Finally, the pixel marked data pits (one) were replaced with the corresponding 3 × 3 median filter CHM value, while the non-data pits (zero) were retained as the original value of the FE CHM.

Spike-Free Algorithm
This algorithm uses all laser returns to generate a TIN and then rasterizes it. The TIN generation starting from the HP and the Delaunay triangle increment is constrained by the lengths of triangle edges and vertical spaces of returns, which suggests that triangles with small edges are more relevant for accurately representing the canopy surface than those with long edges in TIN. Hence, the triangles satisfying the three edge constraints ("freeze distance") will be frozen in the Delaunay triangle increment. If the elevation of the next candidate point is smaller than the elevation of the triangle, it will be ignored to prevent the appearance of pits. To prevent the triangle from being frozen in advance and causing the pits to be overfilled, a height threshold called the "insertion buffer" is introduced. Before the triangle freezes, the insertion buffer defines a vertical zone to ensure that all triangles that still have points in the buffer will not be frozen. Thus, freeze distance and insertion buffer jointly determine the generation of a spike-free TIN.

Graph-Based Progressive Morphological Filtering
This algorithm uses an adaptive morphological operation to filter surface points from all normalized points in progressive filtering. First, a graph is constructed from all normalized points to model the relationships within these points. Second, adaptive morphological filtering is used to "pull up" the depression in the graph. In general, the more depressed the connected neighbourhoods of points in the graph are, the more likely there are data pits. A sharpness index (SI) is introduced to quantify concave tapering and construct the alterable height threshold adaptive to the spacing value between a certain point and its neighbourhood. The adaptive height threshold in graph-based morphological filtering can filter out non-surface points. Then, a progressive process is conducted to improve the adaptive morphological filter for eliminating spikes of various sizes until no points are excluded. The remaining points are regarded as canopy surface points. Finally, the pit-filled CHM is generated by interpolating the surface points.
To avoid errors caused by spatial interpolation methods, all methods used linear interpolation to generate CHMs. According to previous studies, the spatial resolution of a CHM is finer than one-fourth of the crown diameter [33], and it is determined by considering the mean point space [31,42,43]. In our study, the minimum pulse density was approximately 150 pulses/m 2 , which corresponds to a point spacing of 0.08 m. Therefore, the resolution of the CHM was set to 0.1 m. All subsequent processing was performed on the basis of CHM with 0.1 m resolution. The performance of the pit-free method was analysed by the graphical user interface (GUI) developed by Ben-Arie et al. [10], and the other methods were coded in MATLAB R2019b.

Accuracy Assessment
First, the CHM quality can be judged by its visual effect. A high-quality CHM not only has no or few unnatural black pixels within the crown area but also has clear edges and sufficient crown surface details, and this CHM can even reflect the crowns of small trees [27]. Second, the difference between the optimized CHM and the real canopy surface height is also an indicator of CHM quality [7]. Finally, the CHM application performance should be evaluated, such as the accuracy of individual tree detection (ITD) and tree height estimation [4,7]. Herein, we evaluate the CHM generated by the eight methods in a simulated dataset and a real-world dataset.

Accuracy Assessment of Simulated CHMs
The main purpose of constructing a simulated point cloud is to obtain an exact canopy surface height (the upper surface points of the cone and hemisphere, see Figure 1a,b), which can be used as a reference, and it can accurately quantify the height difference from these processed CHMs. We introduced root mean square error (RMSE), mean bias error (Bias), relative RMSE (RMSE%) and relative Bias (Bias%) as evaluation indicators [7,44] (Equations (3)-(6)). A paired t-test was also used to examine whether there were significant differences between the means of the eight CHMs and between them and the reference CHM.
where n is the number of observations (pixels), y i is the reference height value of the i-th observation,ŷ i is the i-th pixel value in the processed CHM and y i is the mean of the observations.

Accuracy Assessment of UAVLS-Derived CHMs
In addition to visual interpretation of the CHM generated by the UAVLS, a cross comparison was carried out in pairs. The RMSE and Bias were calculated among these CHMs derived from the eight methods. Unlike the simulated point cloud, the laser scanning point cloud has no exact actual surface points [27]. Thus, a surface point selection process was conducted to validate the accuracies of various CHMs. We interpreted canopy surface points from the normalized point clouds in each subplot of three sites, including two types: one is the peak of the canopy, and the other is the valley of the canopy. These two types represent the ability of CHMs to express crown heights and crown edges, respectively. A total of 454 points were eventually selected as reference canopy surface points. The coefficient of determination (R 2 ), RMSE and Bias were calculated to evaluate the performances of different methods to express the canopy surface. Additionally, the application of CHMs was evaluated by the accuracy of ITD and tree height estimation from the CHMs generated by different methods. For ITD, a traditional local maximum (LM) algorithm was introduced to automatically detect trees, which was the most commonly used ITD method [6]. The detected trees were matched with field measurement trees based on matching criteria developed by Hao et al. [28]. Both the horizontal distance constraint (less than the reference crown radius) and the height constraint (less than 20% of the top height of the plot) are used to determine the detected tree that matches the reference tree. All correctly matched trees are true positives (TP). Detected trees without a link to references were commission errors (false positives, FP), while reference trees that were not matched to any detected trees were classified as omission errors (false negative, FN). Accuracy assessments in terms of precision (Pr), recall (Re) and Fscore (F score ) were calculated as Equations (7)- (9). For the tree height estimation, RMSE% was used to evaluate the difference between matched trees.

Sensitivity Analysis
Among the eight CHM generation methods, only the FE method and the HPs methods are parameter-free, and the other six methods require at least one user-defined parameter. These parameters will indeed impact the performance of the methods. To obtain the optimal result of each method and enable a fair comparison, a sensitivity analysis of the parameters was carried out using simulated point clouds with pits of different proportions.
For the global filtering method, mean, medium and Gaussian filtering were performed on CHM FE with 3 × 3, 5 × 5 and 7 × 7 window sizes, respectively. The results of the comparison of the CHM Reference against filtered CHMs are shown in Supplementary S1 in the Supplementary Materials, Table S1. The results showed that the median filter and Gaussian filter have the highest accuracies in the 5 × 5 window size, while the mean filter has the highest accuracy in the 3 × 3 window size. The accuracies of the mean filter and median filter were not affected by the proportion of pits or the geometric shape (cone or hemisphere). For a hemisphere with 10% pits, a smaller window is more suitable for Gaussian filtering, and large pits use a large window for Gaussian filtering with higher accuracy.
A user-defined threshold of the pit-free method was controlled from 1% to 30%. The results showed that the accuracy of CHM pit-free was progressively increased by gradually increasing the threshold (starting from 1 at 5% intervals) (Supplementary S1 in the Supplementary Materials, Table S2). The threshold corresponding to a small proportion of pits reached the optimal solution earlier than the large proportion of pits. When the threshold was 20%, the accuracy reached the highest and remained unchanged.
For the spike-free method, the "freeze distance" and the "insertion buffer" constrained the construction of the TIN. When the freeze distance is increasing, filtered pits also increase, and when the insertion buffer is increased, the constraints loosen. We tested the two parameters through the control variables, and the results showed that although the accuracy was nearly the highest when the freezing distance and the insertion buffer were 0.3 and 0.3 (Supplementary S1 in the Supplementary Materials, Table S3), the visual effect of CHM Spike-free was the best when the freezing distance and the insertion buffer were 0.4 and 0.5, respectively (Supplementary S1 in the Supplementary Materials, Figure S1).
As the core parameter of the GPMF method, an SI was used to express the depression degree between a point and its neighbourhoods. Here, we adjusted the SI from 0 to 10 with a step of 1 to test the changes in RMSE between the CHM Reference and the CHM GPMF . The results showed that the accuracy of the cone gradually decreased from SI = 4, while that of the hemispherical sphere first increased and then decreased with increasing SI (Supplementary S1 in the Supplementary Materials, Table S4). Since the accuracy difference between SI = 3 and 4 was only one in ten thousand, we chose 4 as the optimal parameter of SI.
The comparison of the accuracy of CHMs generated by various methods with optimal parameters and reference CHMs with different proportions of pits is shown in Table 1. Notably, the accuracy of all methods decreases with the increase in the proportion of pits. The trend of different methods in different proportions of pits is consistent, which indicates that the data quality has no significant influence on the performance of different methods. Overall, the GPMF method had the lowest RMSE among all the methods, and the spike-free method was the second lowest. The FE method was undoubtedly the worst. The optimal parameters of these methods were used for subsequent applications. Taking the data with 20% pits as an example, the 0.1 m resolution CHMs generated by the simulated surface points and the eight processing methods are shown in Figure 4. A visual comparison of these CHMs was implemented. In general, CHM FE has the largest number of randomly distributed dark pixels, while other CHMs have different degrees of pit-filling. Among them, the visual performance of CHM Spike-free and CHM GPMF are closest to that of CHM Reference . The smooth filtered CHM (CHM Mean , CHM Median and CHM Maussian ) made the image blurry, and the pits did not disappear. The performance of CHM HP and CHM Pit-free were also not satisfactory. The small pits in the CHM Pit-free sample were removed, but clustered pits still occurred. A further comparison of CHM Spike-free and CHM GPMF found that CHM GPMF has a clearer crown edge, and the connection between adjacent canopies is also better processed, while there is adhesion at the adjacent canopy of CHM Spike-free (see the red circle in Figure 4). In addition, we found that there are many pits in the overlapping crown of the hemispherical model, and only the spike-free and GPMF methods can filter these pits well.

Quantitative Analysis
A further quantitative comparison of the differences between these methods is shown in Figure 5. The results indicated that only CHM Spike-free and CHM GPMF have no significant difference from CHM Reference (RMSEs are 0.1602 m and 0.1096 m, and RMSE% values are 1.86% and 1.27% of cone, respectively). The difference between CHM Reference and CHM FE was the largest (RMSE = 0.9114 m, RMSE% = 10.6%), which was approximately 9% lower than CHM Spike-free and CHM GPMF . There was no significant difference between the mean heights of CHM Mean and CHM Gaussian and that of CHM FE . Additionally, the mean heights of CHM HP and CHM Pit-free showed no significant difference. In terms of Bias, we found that the spike-free and GPMF methods slightly overestimated the crown surface height, and the other six methods underestimated it. The results also showed that the median filtering method has the highest accuracy among the three smooth filtering methods. The application of the eight processing methods was not different in the cone and hemisphere.

Visual Performance
For clarity, the CHMs of four subplot sizes in Plot 3 were selected as examples and are shown in Figure 6. The median filter method eliminates almost all the pits, but the image is excessively smooth and the crown shape is blurred. Although the HP method makes the crown clearly visible, it produces more clumps and black pits than the FE method. The Gaussian filter is more blurred than the mean filter, and at the same time, it does not improve the pit problem. Overall, the spike-free and GPMF methods achieved the best visual performances of all the methods, followed by the pit-free method. Figure 7 provides an example of a 0.1 m wide canopy profile highlighting the difference between laser points and the CHMs generated by eight methods. The FE points are not only distributed on the surface of the canopy but also distributed in the canopy or even on the ground, which directly leads to CHM FE with many depressions. Among all the methods, we first noticed the HP profile because it produced sharper and deeper peaks than FE. The pit-free method only fills in the local pits, and the other parts are consistent with the FE profile. In general, the CHM Spike-free and CHM GPMF were closest to the surface of the point cloud, while the smoothed CHMs were seriously underestimated. Notably, the GPMF pulled up some valleys between the two crowns, but this phenomenon is likely due to the influence of the Y-axis neighbourhood and not just the X-axis neighbourhood; thus, the GPMF accuracy needs to be further evaluated.  Table 2 shows the RMSE and Bias values between the eight CHM generation methods. The cross-comparison results showed that there was no significant difference between CHM Mean and CHM Gaussian . The mean heights of the CHM Spike-free and CHM GPMF were not significantly different in broad-leaved forests. The difference between FE and other CHMs was higher than 1 m, of which the largest was HP (2.4499 m of Plot 3). The Bias between mean filtering and median filtering and FE is 0, indicating that these two filtering methods have not changed the average height of the CHM. Among the eight CHMs, CHM Spike-free had the highest average height in broadleaf forests, while CHM GPMF had the highest average height in coniferous and mixed forests. The average value of CHM HP is the lowest, which may be because too many extremely deep pits reduce the average HP based on the abovementioned canopy profile (Figure 7). In coniferous forests, the difference between the eight methods was the largest.  When the reference surface points were introduced to evaluate the performance of the eight methods, the results are shown in Table 3. Overall, the spike-free method showed the best performance in terms of overall accuracy, and the GPMF was suboptimal. The HP method showed good performance in the peak canopy but poor performance in the valley. In contrast, the median performed well in the valley of the canopy but poorly in the peak. The mean filter performed the worst not only in the peak of the canopy but also in the valley. Compared with the canopy valley, the accuracy loss of the canopy peak was the main disadvantage of FE. The spike-free and GPMF methods have a strong ability to explain canopy surface height variation (R 2 close to 1, and RMSE% less than 3%). Overall, all the methods underestimated the reference surface (Bias > 0), of which the mean filter method underestimated the surface the most (Bias% = 7.73%) and the spike-free method underestimated the surface the least (Bias% = 1.32%).

Individual Tree Application Evaluation
The local maxima algorithm was applied to all the processed CHMs. Figure 8 illustrates the accuracy of ITD and tree height estimation plotted against window size (WS) of the local maxima algorithm. Overall, with the increase in WS, the Pr of all the methods increases, the Re and RMSE% decrease and F score presents a "mountain" shape, which first increases to a maximum value and then decreases. That is, the commission error always increases with increasing WS, while the omission error is the opposite. When the commission and omission errors become balanced, the overall accuracy reaches the maximum value. The accuracy of tree height estimation decreased with increasing WS. Further analysis showed that the sensitivity of different CHMs to the LM of different WSs was slightly different. In terms of the accuracy of ITD, the median method was the main outlier in the eight methods. Basically, the median filtering CHM performed best in controlling omission errors, but mass commission errors led to the lowest Pr and F score values. In addition, almost all methods achieved the highest overall accuracy (F score ) when WS was 17 × 17 (1.7 m × 1.7 m). Therefore, the small window in Figure 8 is partially enlarged when WS is equal to 17. The dotted line and the solid line represent four post-processing methods and four pre-processing methods, respectively. As there is little difference between these methods, we evaluated them by 50% of the rankings. Specifically, for broadleaf forest, the HP, GPMF and spike-free methods were in the top 50% due to their excellent performances in controlling the commission error. For coniferous and mixed forests, the mean and Gaussian filter showed remarkable performances. For the spike-free method, even if there is an enormous omission error, the extremely high commission error control gives it a high overall accuracy in each stand type. In addition, we found that the post-processing method (mean, median, Gaussian and pit-free method) had a small omission error (Re is usually high), which may benefit from excessive detection.
In terms of tree height estimation (RMSE%), we found that coniferous forest has the highest accuracy compared with broad-leaved forest and mixed forest, and the WS has the greatest influence on the estimation accuracy of mixed forest. The GPMF and spike-free method always achieved high accuracy regardless of the forest stand. The HP method performed well in broadleaf and coniferous forests, and the mean filter performed well in mixed forests. In coniferous forests, when WS gradually increases, the tree height estimation accuracies of the three smooth filtering methods become increasingly lower, and the difference compared with other methods increases significantly. The accuracy of the median filter method was significantly lower than that of the other methods in mixed forests. In addition, the difference from other methods was not obvious.

Discussion
In recent years, an increasing number of studies have directly used UAVLS as a measurement tool due to the high-density point cloud data, which requires data processing to minimize the loss of accuracy. As one of the inherent problems in the processing of LiDAR data, pits will seriously affect the quality of the CHM [9]. Although many studies focus on pit-filling of the CHMs, most of these algorithms were developed and applied to low-density ALS (0.7-8 pt/m 2 ) to generate CHMs from 0.2 m to 0.5 m [7,9,11,12,[25][26][27]. In the present study, we evaluated the performances of several pit-filling methods for generating 0.1 m resolution CHMs by using both high-density simulated data and UAVLS data. A comprehensive visual performance, quantitative analysis and application evaluation framework was used to compare these CHMs.
Compared with real laser scanning data, artificially constructed simulated PCs can define the exact canopy surface to intuitively and quantitatively compare the effects of different methods [11,26]. In terms of the simulated data results, we obtained a clear ranking of the advantages and disadvantages in the eight methods (the order is GPMF, spike-free, pit-free, median, HP, Gaussian, and mean). In contrast, there are many uncertainties in the real UAV point cloud, which makes the results more complicated. In particular, the UAVLS data for this study were collected from natural secondary forests with diverse tree species, high stand density and complex forest layers. The sources of data pits may vary, and it is difficult to find a true canopy surface as a reference to verify different CHMs. Previous studies have used different techniques to assess the accuracy of these pit-filling algorithms. Ben-Arie et al. [10] visually compared the efficacy of the pit-filling algorithm through CHM and its X-axis profile. Khosravipour et al. [4,17] evaluated the spike-free method using the accuracy of ITD. Hao et al. [27] first segmented individual trees and then selected convex hull vertices above the maximum crown diameter of each individual tree as reference surface points. Mielcarek et al. [7] tested and evaluated five CHM generation methods by the accuracy of tree height estimation. Chen et al. [11] and Zhang et al. [26] evaluated the accuracy of only plot-level maximum tree height estimation. By comparison, the evaluation framework we used can comprehensively and systematically compare the effects of different methods.
Regarding the results of the UAVLS-derived CHM evaluation, the GPMF and spikefree method showed optimal visual performances, similar to the simulated data, and they produced few pits while preserving details of the crown. The verification of canopy surface points indicated that the spike-free method performed slightly better than the GPMF method. The accuracy of ITD for broadleaf forest and tree height estimation of the GPMF method was slightly better than that of the spike-free method. For the FE method, since the first return not only exists on the canopy surface but is also distributed inside the canopy and on the surface (Figure 7), it produced the most data pits among all methods. This confirmation is consistent with the results from [9,10,28]. For the global filtering method, the accuracy of the median filter CHM construction was better than that of mean and Gaussian, but the accuracy of the median filter CHM application showed outliers. Further analysis found that the median filter produced many consecutive equal pixel values, which led to multiple continuous treetops in the LM filter and significantly increased the commission error (shown in Figure 9). The canopy height was significantly underestimated, and this method produced a mass commission error [27]. For the pit-free method, the ITD and tree height estimation results are almost consistent with those of the FE method because it only fills the filtered pits based on FE [7,27]. Thus, although the post-processing methods based on the original CHM have eliminated partial pits in the visual effect, they show limited improvement in the quality of the CHM and are even lower than the original CHM [27]. For the HP method, although it eliminates some small pits, this method produces deeper pits compared with the FE method, which will greatly reduce the average height. However, a small number of deep pits do not affect LM detection, so it is not restricted to individual tree applications. Hao et al. [27] claimed that the HP method is sensitive to the density of the point cloud; when the pixel is too large, the details of the crown are missed. Leckie et al. [12] set the grid pixel size at 25 cm according to the average point spacing. Our average point spacing was between 0.06 m and 0.08 m, so a 0.1 m resolution was appropriate. According to the description of Shamsoddini et al. [9], a robust pit-filling algorithm needs to satisfy the following criteria: identify pits using an adaptive threshold, perform oriented detection of pits without changing the maximum value, fill the pits with suitable values to represent the crown shape as closely as possible and preserve canopy gaps. Considering the driving parameters of the algorithm, the GPMF method and pit-free algorithm are parameter-driven algorithms, the spike-free method has two user-defined parameters, the filter method requires the WS be set, and the FE and HP methods have no parameters. The spike-free and GPMF differ fundamentally from the others; they used all laser returns to generate CHMs while improving the potential pits. Using all returns avoids the loss of crown information and eliminates some errors caused by using only the first return [7,45]. The difference between them is that the spike-free algorithm prevents spike formation during TIN construction [17], while GPMF excludes non-surface points from all returns in a progressive filtering process [27]. GPMF is similar to ground point filtering [11]; it screens surface points from laser scanning data and can be used for other canopy surface studies. Compared with other methods, the spike-free method was developed based on high-density PCs, which is beneficial to the application of this method in high-resolution CHM generation. Although the CHM filtering method underestimates the canopy height, it is simple to process and could prove useful for other applications. For example, Shamsoddini et al. [9] pointed out that CHM filtering combined with linear regression has an excellent performance in estimating tree height. The error of the postprocessing method mainly comes from the canopy peak, and the error of the pre-processing method mainly comes from the valley (Table 3), which proves that the pre-processing method can orient the detection of pits without changing the maximum value.
By comparing the results of different forest types, the average height difference in the CHM generated by different methods in coniferous forests is the largest ( Figure 5 and Table 2), indicating that the impact of pits is worth consideration when generating a CHM of coniferous forests. The FE and filter CHM severely underestimated the heights of coniferous forests, especially for the peak canopy (Table 3), because coniferous crowns are conical and easy to wear during filtration. The accuracy difference between the GPMF method and spike-free algorithm mixed forests is larger than that of the other forest types, which indicates that the adaptability of the GPMF method in complex forests is not as good as that in spike-free forests. Regarding the detection of individual trees, coniferous forests have higher overall accuracies than broad-leaved forests and mixed forests because of their lower commission errors ( Figure 8). Likewise, the tree height estimation accuracy of coniferous forest was the highest among the three forest types when the detection WS was the same. The accuracies of ITD and tree height estimation of mixed forest were the lowest.
This study demonstrated that the effect of the pit-filling method on individual tree applications of high-resolution CHMs is small. The outliers of the median filter in broadleaved forest, coniferous forest and mixed forest were 12%, 9% and 8%, respectively. The differences from other methods for different forest types are very small (the difference of F score was 1-3%), while the accuracies of different methods for estimating tree height in any forest was less than 0.1%. This result indicated that the effect of data pits on individual tree applications (ITD and tree height estimation) is not significant when using a highresolution CHM. Perhaps with the decrease in resolution, the difference will be more obvious. For example, the results of Mielcarek et al. [7] showed that the difference in tree height estimation accuracy between spike-free and Gaussian filters is 3.29% when using a 0.5 m resolution CHM. Hao et al. [27] showed that the difference in the overall accuracy of ITD between GPMF and FE is 17.43%. The results of Khosravipour et al. [4] show that the accuracy difference between the pit-free CHM developed by Khosravipour et al. [4] and the Gaussian smoothed CHM is 3.6% when using a 0.15 m resolution CHM, while the difference is 32% when using a 0.5 m resolution CHM.
This study aims to compare and evaluate the performances of different pit-filling methods to generate high-resolution CHMs. To reduce the uncertainty of complex parameter selection in the ITD algorithm, a single parameter and commonly used local maxima detection method are used. Although there is little difference in the application of individual trees, notably, CHMs without artefacts are more useful and reliable, with better reproductions of crown shapes [7]. Other CHM applications should be further evaluated, such as crown delineation and other tree metric estimations, and more data set from complicated natural forests would be recommended to take into account for evaluating the effectiveness of the framework on the diverse tree species, high stand density, and complex forest layers. Notably, this study area is a natural secondary forest with high-density stems, and it is difficult to measure individual tree position and height. Due to many uncertainties in field measurements [7], many studies suggest that LiDAR be directly used as a measurement tool [46,47]. Based on the topography of the study area (7 • to 11 • ), there is no significant difference between the CHM generated by normalized PCs and the CHM obtained by the DSM minus DEM, regardless of the method used.

Conclusions
The high sampling density near-ground UAVLS system can extract fine forest canopy structures. Quantifying the accuracy of high-resolution CHM construction derived from UAVLS is an important precursor to the use of UAVLS data to obtain other forest parameters. This study compared and evaluated several CHM generation methods from multiple aspects. The results demonstrated that different processing methods can affect the quality of the CHMs generated by high-resolution UAVLS data. The GPMF [27] and spike-free [17] methods have excellent performance in generating CHMs; they effectively remove data pits and preserve crown details, while minimizing the loss of canopy height under complex stand conditions. Using the first echo of laser scanning will produce massive amounts of data noise and result in the incorrect expression of the canopy surface. Although the global filtering method (mean, median and Gaussian filter) of the CHM decreases the noise visually, it significantly reduces the canopy surface height. A median filter is not suitable for individual tree applications based on LM detection when using a high-resolution CHM. The pit-free method [10] detects and fills pits locally in the original CHM, and it can improve limited accuracy, similar to other post-processing methods. The HP interpolation method [12] has a good performance in the generation of high-resolution CHMs, but a small number of deep points are introduced. In high-density and complex forests, the suppressed tree crown is also an important source of CHM pits. Compared with broadleaved and mixed forests, coniferous forests are more sensitive to data pits, so selection of the CHM generation method should be considered. Although these high-resolution CHM generation methods show little difference in the application of ITD and tree height estimation, a CHM without noise is also valuable. This study provides a comprehensive reference for CHM generation using high-resolution LiDAR data and demonstrates the potential of UAVLS as a forest measurement tool in obtaining forest inventories.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/rs13122239/s1. Supplementary S1: Table S1, the RMSE of CHM Reference and filtered CHMs with different window sizes in different proportion of pits; Table S2, the RMSE of CHM Reference and CHM Pit-free with different thresholds in different proportions of pits; Table S3, the RMSE of CHM Reference and CHM Spike-free with different parameters in different proportions of pits; Figure S1, the visual effects of CHM Spike-free with different parameters; Table S4, the RMSE of CHM Reference and CHM GPMF with different SI in different proportions of pits.  Data Availability Statement: Due to the nature of this research, participants of this study did not agree for their data to be shared publicly, so supporting data is not available.