Robust and Parameter-Free Algorithm for Constructing Pit-Free Canopy Height Models

Data pits commonly appear in lidar-derived canopy height models (CHMs) owing to the penetration ability of airborne light detection and ranging (lidar) into tree crowns. They have a seriously negative effect on the quality of tree detection and subsequent biophysical measurements. In this study, we propose an algorithm based on robust locally weighted regression and robust z-scores for the construction of a pit-free CHM. A significant advantage of the new algorithm is that it is parameter free, which makes it efficient and robust for practical applications. Simulated and airborne lidar-derived data sets are employed to assess the performance of the new method for CHM construction, and its results are compared to those of three classical methods, namely the natural neighbor (NN) interpolation of the highest point method (HPM), mean filter, and median filter. The results from the simulated data set demonstrate that our algorithm is more accurate compared to the three classical methods for generating pit-free CHMs in the presence of data pits. CHM construction using the lidar-derived data set shows that, compared to the classical methods, the new method has a better ability to remove data pits as well as preserving the edges, shapes, and structures of canopy gaps and crowns. Moreover, the proposed method performs better compared to the classical methods in deriving plot-level maximum tree heights from CHMs. Thus, the new method shows high potential for pit-free CHM construction.


Introduction
In addition to characterizing ground topography, airborne light detection and ranging (lidar) data have been widely used in forestry in recent years [1,2].Lidar-derived canopy height models (CHMs), which represent the absolute height distribution of the vegetation canopy above the ground [3,4], have a significant influence on the extraction of forest inventory information, such as individual tree crown delineation [5][6][7], tree height estimation [8][9][10], and biomass estimation [11,12].
A CHM is typically constructed by interpolating the first return lidar points and determining their heights above a digital elevation model (DEM) [13].Data pits with heights significantly lower than those of their neighbors commonly appear in raster CHMs.According to Leckie et al. [14], data pits appear due to the combination of different flight line data sets and the penetration of the laser beam into canopy branches and foliage before its first return.Data pits may even appear in individual flight lines at off-nadir scan angles.Ben-Arie et al. [15] identified that data pits may be caused by the combination of various factors such as data acquisition and post-processing.Vosselman [16] indicated that the planimetric error between overlapping flight strips caused by the global positioning system (GPS) and inertial measurement unit (IMU) measurements can also generate data pits, particularly for small-footprint lidar.Many algorithms have been developed to generate pit-free CHMs.The performance of the present algorithms is often subject to different parameters that should be optimally tuned.However, the process of parameter tuning is not easy in practice.Thus, parameter-free algorithms are more desirable.
In a statistical sense, data pits in a CHM can be considered outliers because the heights of data pits are significantly lower than those of their neighbors.Thus, the problems of data pits can be remedied by outlier detection techniques.We propose a robust statistical method based on this idea to detect data pits.The new method uses robust locally weighted regression to perform interpolations and robust z-scores to detect data pits.After the removal of data pits, the remaining lidar points are interpolated using the natural neighbor (NN) method to construct a pit-free CHM.NN was selected because it is parameter free, easy to apply, and highly accurate for interpolating lidar points [17].Compared to the classical pit-free algorithms, the main novelties of the new method are as follows: (i) it is parameter-free; and (ii) it is applied on the height-normalized lidar point clouds rather than the raster CHMs, thus avoiding errors in the transformation from the randomly distributed point clouds to the rasterized CHM.
The remainder of this paper is organized as follows.Section 2 presents related work.Section 3 elaborates the principle of the new method.In Section 4, a numerical test and a real-world example are employed to assess the performance of our algorithm, and its results are compared with the results of three classical methods, that is, NN interpolation of the highest point method (HPM), mean filter, and median filter.It should be noted that the mean and median filters use a 3 × 3 kernel in the tests.Discussion and conclusions are presented in Section 5.

Related Work
It has been demonstrated that data pits have a seriously negative effect on the estimation of forest parameters.For example, Persson et al. [9] found that data pits in a CHM complicate the detection of tree crowns.Gaveau and Hill [18] demonstrated that the penetration of laser beams into a tree crown leads to the underestimation of canopy height.Khosravipour et al. [19] indicated that in the context of treetop detection, data pits in CHMs lead to large commission and omission errors.Thus, a variety of algorithms have been developed to construct pit-free CHMs.For example, Leckie et al. [14] presented an algorithm that first assigns all lidar points to 25-cm grid cells, and then employs the finite difference method to interpolate the highest point in each grid.Ben-Arie et al. [15] presented a semiautomatic pit-filling algorithm to detect and replace data pits using Laplacian and median filters.This method was further developed by Zhao et al. [20].They used morphological crown control to recover crown coverage.Shamsoddini et al. [21] employed an adaptive mean filter with 3 × 3 and 5 × 5 kernels to detect and fill data pits.Liu and Dong [22] used a selecting and sorting scheme to select data points during CHM construction.
Although data pits can be removed using the aforementioned image smoothing methods, the heights of all pixels in the CHMs are altered.This can cause the underestimation of tree heights and crown diameters due to the omission of treetops and the reduction of crown shoulders [19].Thus, Khosravipour et al. [19] proposed a pit-free CHM algorithm by generating a stack of CHMs that were partially pit free-each representing a different height interval of the canopy-and finally merged into a pit-free CHM.However, the performance of the pit-free CHM algorithm is largely determined by a rasterization threshold and a series of height thresholds.This makes it difficult to apply.Thus, a parameter-free algorithm is more practical.

Principle of the Proposed Algorithm
The flowchart of the proposed algorithm for pit-free CHM construction is shown in Figure 1.As can be seen, the main steps are robust locally weighted regression (RLWR) and pit removal using robust z-scores.As in the case of outliers in statistics, data pits can be defined as the observations having large z-scores with respect to interpolation errors.In the presence of data pits, a robust interpolator should be used to interpolate lidar points.Here, RLWR is adopted due to its simplicity, robustness, and high accuracy with densely distributed data sets [23].Then, observations with robust z-scores exceeding a tolerance (i.e., 2.5) are flagged as data pits.
ISPRS Int.J. Geo-Inf.2017, 6, 219 3 of 13 observations having large z-scores with respect to interpolation errors.In the presence of data pits, a robust interpolator should be used to interpolate lidar points.Here, RLWR is adopted due to its simplicity, robustness, and high accuracy with densely distributed data sets [23].Then, observations with robust z-scores exceeding a tolerance (i.e., 2.5) are flagged as data pits.

RLWR
Locally weighted regression (LWR), originally proposed by Cleveland [24] and further developed by Cleveland and Devlin [25], has been widely used for function estimation.Compared to global polynomial fitting, LWR has many desirable statistical properties, such as numerical stability, minimax efficiency, and absence of boundary effects [26].Moreover, in contrast to geostatistical interpolation methods, LWR does not require the specification of a semivariogram, which should be determined beforehand.
Given a group of 3D points ( ) , the data is supposed to be fitted with the model as follows: ( ) where i e represents sample errors with zero mean and unknown variance and ( ) function to be estimated.In this study, x and z represent the horizontal coordinate and normalized height, respectively, of a lidar point.
When performing LWR, a low-degree polynomial is first fitted to the neighbors of xi using a weighted least square algorithm, where the weight of each point is obtained from a weight function.

RLWR
Locally weighted regression (LWR), originally proposed by Cleveland [24] and further developed by Cleveland and Devlin [25], has been widely used for function estimation.Compared to global polynomial fitting, LWR has many desirable statistical properties, such as numerical stability, minimax efficiency, and absence of boundary effects [26].Moreover, in contrast to geostatistical interpolation methods, LWR does not require the specification of a semivariogram, which should be determined beforehand.
Given a group of 3D points , the data is supposed to be fitted with the model as follows: where ε i represents sample errors with zero mean and unknown variance and g(x) is a smooth function to be estimated.In this study, x and z represent the horizontal coordinate and normalized height, respectively, of a lidar point.When performing LWR, a low-degree polynomial is first fitted to the neighbors of x i using a weighted least square algorithm, where the weight of each point is obtained from a weight function.Generally, the weight function of LWR assigns high weights to points near the point of interest and low weights to points far away.This is based on the idea that points near each other are more likely to be related than points that are further apart.Cleveland [24] suggested using the tri-cube weight function, where N(x i ) is the set of the k nearest neighbors of x i in 2D space, d ij is the distance between x i and its j-th neighbor, and The objective function of LWR is defined as min ∑ Thus, the function value at x i can be obtained by evaluating the local polynomial with the polynomial coefficients estimated from the solution of Equation ( 3).The commonly used polynomial is either linear or quadric.This is based on the idea that any function can be accurately and easily approximated in a small neighborhood by a low-order polynomial.High-degree polynomials are prone to overfitting the data in each subset and are numerically unstable [26].
LWR is highly sensitive to outliers due to the least squares -based objective function.Thus, if data pits exist in the neighborhood of the point of interest, LWR may lead to a biased estimation.Therefore, to reduce the influence of outliers, a robust LWR (RLWR) is adopted that assigns a robust weight to each point in the neighborhood.Cleveland [24] employed the bisquare weight function to assure the robustness of RLWR.The bisquare weight function is defined as follows where r j is the regression residual of the j-th nearest neighbor of the point of interest, r j = z j − g x j , and s is the median of |r i |.Thus, the objective function of RLWR is expressed as min ∑ Therefore, the estimation of g(x i ) at x i by RLWR involves the following steps: (i) Finding the k nearest neighbors N(x i ) of x i from the point cloud.Here, we set k = 12.
(vi) Computing the robust weights of x j , j = 1, 2, • • • , k based on the robust weight function, i.e., Equation ( 4).(vii) Estimating the polynomial coefficients by minimizing the objective function of RLWR, i.e., Equation ( 5).(viii) Repeating (v)-(vii) until the polynomial coefficients are stable, or the maximum number of iteration is reached.Cleveland [24] indicated that two iterations are sufficient to obtain a good fit.(ix) Estimating g(x i ) at x i with the computed polynomial coefficients.The interpolation error of x i is expressed as, e i = z i − g(x i ).

Pit Detection Using Robust Z-Scores
As in the case of outlier detection in statistics, the well-known z-score with respect to interpolation error is employed to detect data pits in this study.Z-score is a distance-based measure that is defined as the standardized residual [27], where e and σ e denote the mean and standard deviation, respectively, of e.Both e and σ e have a breakdown point (BP) of zero and are sensitive to outliers.Their robust alternatives are the median and median absolute deviation, respectively, both of which have a BP of 50% [27].The robust z-score is expressed as ti where e median and σ MAD represent the median and median absolute deviation, respectively, of e.
In statistics, observations with ti exceeding 2.5 are regarded as outliers.The threshold value of 2.5 is derived from the fact that the proportion of a random variable with a standard normal distribution exceeding 2.5 is less than 1.24%.The threshold value of 2.5 for a criterion to detect the data-pits as outliers is derived from a statistical viewpoint rather than the forest types.Therefore, this threshold is suitable for any forest type.Because the heights of data pits are always lower than those of their neighbors, we define data pits as the observations with ti less than −2.5.

Experiments
We use a numerical test and a real-world example to assess the robustness of our method for CHM construction.Furthermore, plot-level maximum tree heights derived from CHMs are used to indirectly evaluate the performance of the proposed method in the real-world example.

Numerical Test
Two simple geometric models using cones and hemispheres [28] are employed to simulate tree crowns.The mathematical formulations of the two models are expressed as: Hemisphere : The procedure for comparing the performance of our method with that of the present methods is as follows: (i) 1000 points subject to the condition of x 2 + y 2 ≤ 1 are randomly sampled from the model, e.g., cone.The average point space is 0.056.(ii) n data pits are randomly selected and their heights are artificially changed, i.e., their elevations are reduced by 0.3.Here, n = 100 or 200; the contaminating proportion (α) is 10% or 20%.
(iii) Our method and the three other pit-free methods (mean filter, median filter, and HPM) are used to construct CHMs with the contaminated data set.(iv) Taking the original sample points as check points, the CHMs are assessed in terms of root mean square error (RMSE) and mean error (ME).These are expressed as: where z i and ẑi are the true and simulated values, respectively, at the i-th check point and m is the number of check points, i.e., m = 1000.
For mean and median filters, a raw CHM was first constructed using NN interpolation of the contaminated data set.Then, the CHM was smoothed by the two filters.For HPM, the highest point in each grid cell with a resolution of 0.056 was extracted, and then, NN was used to construct a CHM with the highest points.
The results (Table 1) indicate that, irrespective of accuracy measures, NN produces the poorest CHM with the contaminated data sets, demonstrating that data pits have a significantly negative effect on the quality of CHMs.The performance of all pit-free algorithms decreases as the contaminating proportion increases.Our method yields the best results on all simulated data sets.On average, the new method is ~2.7,1.7, 1.8, and 2.1 times as accurate as the raw CHM, median filter, mean filter, and HPM, respectively, in terms of RMSE.In terms of ME, the corresponding figures are 31.3,25.6, 31, and 16.6.Figure 2 shows CHMs of the cone model constructed by different methods with the original points and contaminated points under α = 10%.The raw CHM has many data pits randomly distributed in the image (Figure 2b).Among the pit-free algorithms, HPM yields the poorest result, because many data pits are clearly visible (Figure 2c).This is because some of the highest points still contain data pits.Even though the mean filter removes data pits, it produces smoothed square artifacts (Figure 2d).Both the median filter and our method are successful in removing data pits (Figure 2e,f).This is expected, as the two methods have high breakdown values, especially the median filter with a value of 50%.However, both filters alter the values of all CHM cells; this results in a slight distortion of the crown shape in the right-bottom (Figure 2e).The proposed method (Figure 2f) yields the best result, which closely approximates the original (Figure 2a).

Study Site and Raw Data Sets
The study site chosen was the Tianlaochi catchment (38°23′55″-38°26′57″ N, 99°53′45″-99°57′12″ E), Gansu Province, China.It is characterized by a cold semi-arid climate with annual average temperature and precipitation of 0.6 °C and 437.2 mm, respectively.The site was mainly covered by coniferous plantation forests and shrubs.The forest vegetation types included picea crassifolia and sabina przewalskii.A mean elevation of 3292 m, with a minimum of 2596 m and a maximum of 4411 m, and a mean slope of 29.6° characterized the topography of the study area.
A total of 30 rectangular plots of area 20 × 10 m 2 were chosen in the study site in August 2013 (Figure 3).Among the 30 plots, 26 were located in picea crassifolia stands, two in sabina przewalskii stands and two in stands of a mixture of the two species.All tree heights in each plot were measured by a handheld laser rangefinder (Haglof Vertex IV Ultrasonic Hypsometer made by HAGLOF INC), and the maximum tree height was recorded.
Raw lidar data for the study site were acquired using Leica ALS70 with a laser wavelength of 1064 nm in July 2012.During data capture, the absolute flying height was about 4800 m with a mean data density of 1 point/m 2 .The lidar system recorded up to four returns for each laser pulse, depending on the ground cover.The TerraScan (http://terrasolid.fi)software, produced by Terrasolid Ltd., Helsinki, Finland, was employed to differentiate ground and non-ground points.
The resolution of raster CHMs plays an important role in estimating individual tree attributes [29].The CHM resolution should not be larger than half the minimum crown size, and overly smaller than the mean point space [19,30].Thus, a CHM resolution of 0.5 m was determined by recognizing an average crown diameter of 1.8 m and an average point space of 1 m.

Study Site and Raw Data Sets
The study site chosen was the Tianlaochi catchment (38 • 23 55"-38 • 26 57" N, 99 • 53 45"-99 • 57 12" E), Gansu Province, China.It is characterized by a cold semi-arid climate with annual average temperature and precipitation of 0.6 • C and 437.2 mm, respectively.The site was mainly covered by coniferous plantation forests and shrubs.The forest vegetation types included picea crassifolia and sabina przewalskii.A mean elevation of 3292 m, with a minimum of 2596 m and a maximum of 4411 m, and a mean slope of 29.6 • characterized the topography of the study area.
A total of 30 rectangular plots of area 20 × 10 m 2 were chosen in the study site in August 2013 (Figure 3).Among the 30 plots, 26 were located in picea crassifolia stands, two in sabina przewalskii stands and two in stands of a mixture of the two species.All tree heights in each plot were measured by a handheld laser rangefinder (Haglof Vertex IV Ultrasonic Hypsometer made by HAGLOF INC), and the maximum tree height was recorded.
Raw lidar data for the study site were acquired using Leica ALS70 with a laser wavelength of 1064 nm in July 2012.During data capture, the absolute flying height was about 4800 m with a mean data density of 1 point/m 2 .The lidar system recorded up to four returns for each laser pulse, depending on the ground cover.The TerraScan (http://terrasolid.fi)software, produced by Terrasolid Ltd., Helsinki, Finland, was employed to differentiate ground and non-ground points.
The resolution of raster CHMs plays an important role in estimating individual tree attributes [29].The CHM resolution should not be larger than half the minimum crown size, and overly smaller than the mean point space [19,30].Thus, a CHM resolution of 0.5 m was determined by recognizing an average crown diameter of 1.8 m and an average point space of 1 m.

Data Processing
The non-ground points were first height-normalized by replacing the elevation of each point with its height above the ground.Then, the four pit-free methods were employed to construct CHMs with the normalized lidar data.Since tree height underestimation is one of the major concerns for CHMs in forest areas [22], plot-level maximum tree heights were derived from the CHMs and compared with field measurements.This provides a comparison between our method and the other algorithms [3].

Results
Taking one plot as an example, Figure 4 shows the 0.5-m CHMs produced by different methods.The raw CHMs made by the NN interpolation are seriously influenced by data pits (Figure 4a).The CHM of the HPM appears visually similar to the original, indicating the poor performance of HPM in terms of the removal of data pits (Figure 4b).Mean and median filters (Figure 4c,d) eliminate the data pits; however, their CHMs are heavily blurred compared to the original.Our method performs better than the others for removing data pits, while preserving the edges, shapes, and structures of canopy gaps and crowns (Figure 4e).
Figure 5 shows the relationship between the measured and simulated plot-level maximum tree heights for the 30 plots.For the mean and median filters (Figure 5a,b), the linear models explain ~66.3% of the variation of plot-level maximum tree heights, whereas HPM (Figure 5c) and our method (Figure 5d) explain ~71% and 73%, respectively.The new method is more accurate compared to the classical methods, because the points of our method approximate the line y = x more closely than those of the other methods.
Figure 6 shows the relationship between the measured plot-level maximum tree heights and simulated residuals, in other words, the difference between measured and simulated plot-level maximum tree heights, for the 30 plots.It can be seen that mean and median filters (Figure 6a,b) underestimate the maximum tree heights in almost all plots.HPM (Figure 6c) and our method (Figure 6d) seem significantly better than median and mean filters.In comparison, the simulated residuals of our method are more prone to zero than those of the other methods.
Table 2 quantitatively shows a comparison between the four methods of the accuracy in the estimation of the plot-level maximum tree heights in the 30 plots.All methods have the defect of tree height underestimation, and our method has a lower systematic error in terms of mean error

Data Processing
The non-ground points were first height-normalized by replacing the elevation of each point with its height above the ground.Then, the four pit-free methods were employed to construct CHMs with the normalized lidar data.Since tree height underestimation is one of the major concerns for CHMs in forest areas [22], plot-level maximum tree heights were derived from the CHMs and compared with field measurements.This provides a comparison between our method and the other algorithms [3].

Results
Taking one plot as an example, Figure 4 shows the 0.5-m CHMs produced by different methods.The raw CHMs made by the NN interpolation are seriously influenced by data pits (Figure 4a).The CHM of the HPM appears visually similar to the original, indicating the poor performance of HPM in terms of the removal of data pits (Figure 4b).Mean and median filters (Figure 4c,d) eliminate the data pits; however, their CHMs are heavily blurred compared to the original.Our method performs better than the others for removing data pits, while preserving the edges, shapes, and structures of canopy gaps and crowns (Figure 4e).
Figure 5 shows the relationship between the measured and simulated plot-level maximum tree heights for the 30 plots.For the mean and median filters (Figure 5a,b), the linear models explain ~66.3% of the variation of plot-level maximum tree heights, whereas HPM (Figure 5c) and our method (Figure 5d) explain ~71% and 73%, respectively.The new method is more accurate compared to the classical methods, because the points of our method approximate the line y = x more closely than those of the other methods.
Figure 6 shows the relationship between the measured plot-level maximum tree heights and simulated residuals, in other words, the difference between measured and simulated plot-level maximum tree heights, for the 30 plots.It can be seen that mean and median filters (Figure 6a,b) underestimate the maximum tree heights in almost all plots.HPM (Figure 6c) and our method (Figure 6d) seem significantly better than median and mean filters.In comparison, the simulated residuals of our method are more prone to zero than those of the other methods.
(ME) than that of the classical methods.Additionally, the new method has the lowest-level variation in terms of standard deviation, and its maximum error is significantly smaller than that of the other methods.(ME) than that of the classical methods.Additionally, the new method has the lowest-level variation in terms of standard deviation, and its maximum error is significantly smaller than that of the other methods.Table 2 quantitatively shows a comparison between the four methods of the accuracy in the estimation of the plot-level maximum tree heights in the 30 plots.All methods have the defect of tree height underestimation, and our method has a lower systematic error in terms of mean error (ME) than that of the classical methods.Additionally, the new method has the lowest-level variation in terms of standard deviation, and its maximum error significantly smaller than that of the other methods.

Discussion
Data pits and tree height underestimation are two major problems of lidar-derived CHMs.To overcome the two problems, we proposed a robust algorithm based on locally weighted regression (LWR) and z-scores to remove data pits.Rather than removing data pits from CHM rasters, our method operated on the height-normalized lidar points in 3D space.Thus, it avoided the alteration of CHM values, and accurately preserved CHM contrast.Both the numerical test and the real-world example validated the advantages of the proposed method.
In recent years, many algorithms have been developed to generate pit-free CHMs.However, they commonly suffer from the necessity of selecting appropriate parameters.For example, the method of Khosravipour et al. [19] required a series of height thresholds and a rasterization threshold for constructing a number of partial CHMs.Liu and Dong [22] showed that a percentage threshold should be set for selecting non-ground points for CHM construction.Rizaev et al. [31] suggested choosing a proper moving window size for local minimum identification.However, it is non-trivial to tune the optimal parameters in practice.In contrast, our method is parameter free for lidar point interpolation and data pit detection, making it useful for practical applications.

Conclusions
A robust algorithm based on robust locally weighted regression and robust z-scores was proposed to reduce the negative effect of data pits on CHM derivativeness.A numerical test indicated that with the existence of data pits, the proposed method is more accurate than HPM, mean filter, and median filter methods for CHM construction.A real-world example demonstrated that the proposed method performed better than the classical methods in removing data pits and preserving the edges, structures, and shapes of canopy crowns and gaps.Moreover, the proposed method reduced the underestimation of plot-level maximum tree height.In conclusions, the newly developed method shows high potential for pit-free CHM construction.

Figure 1 .
Figure 1.Flowchart of the proposed algorithm for pit-free canopy height model (CHM) construction.

Figure 1 .
Figure 1.Flowchart of the proposed algorithm for pit-free canopy height model (CHM) construction.

13 Figure 2 .
Figure 2. CHMs produced by (a) natural neighbor (NN) interpolation of original sample points; (b) NN interpolation of the contaminated data set with the contaminating proportion of 10%; (c) highest point method (HPM); (d) mean filter; (e) median filter; and (f) our method.

Figure 2 .
Figure 2. CHMs produced by (a) natural neighbor (NN) interpolation of original sample points; (b) NN interpolation of the contaminated data set with the contaminating proportion of 10%; (c) highest point method (HPM); (d) mean filter; (e) median filter; and (f) our method.

Figure 5 .
Figure 5. Relationship between the measured and simulated plot-level maximum tree heights for the 30 plots.(a) Mean filter; (b) median filter; (c) HPM; and (d) our method.

Figure 6 .
Figure 6.Relationship between the measured plot-level maximum tree heights and simulated residuals, i.e., difference between measured and simulated plot-level maximum tree height, for the 30 plots.(a) Mean filter; (b) median filter; (c) HPM; and (d) our method.

Figure 5 .
Figure 5. Relationship between the measured and simulated plot-level maximum tree heights for the 30 plots.(a) Mean filter; (b) median filter; (c) HPM; and (d) our method.

Figure 5 .
Figure 5. Relationship between the measured and simulated plot-level maximum tree heights for the 30 plots.(a) Mean filter; (b) median filter; (c) HPM; and (d) our method.

Figure 6 .
Figure 6.Relationship between the measured plot-level maximum tree heights and simulated residuals, i.e., difference between measured and simulated plot-level maximum tree height, for the 30 plots.(a) Mean filter; (b) median filter; (c) HPM; and (d) our method.

Figure 6 .
Figure 6.Relationship between the measured plot-level maximum tree heights and simulated residuals, i.e., difference between measured and simulated plot-level maximum tree height, for the 30 plots.(a) Mean filter; (b) median filter; (c) HPM; and (d) our method.

Table 1 .
Accuracy comparison between our method, raw CHM, and three other methods for removing data pits under different contaminating proportions (α) in the numerical test.

Table 2 .
Accuracy comparison between the four methods (unit: m).ME, Std, and Max represent mean error, standard deviation, and maximum error, respectively, between the measured and simulated plot-level maximum tree heights in the 30 plots.