Generating Tree-Lists by Fusing Individual Tree Detection and Nearest Neighbor Imputation Using Airborne LiDAR Data

Joonghoon Shin, Hailemariam Temesgen

Department of Forest Engineering, Resources and Management, Oregon State University, Peavy Hall, Corvallis, OR, USA.

DOI: 10.4236/ojf.2018.84032                

Abstract

Individual tree detection (ITD) and the area-based approach (ABA) are combined to generate tree-lists using airborne LiDAR data. ITD based on the Canopy Height Model (CHM) was applied for overstory trees, while ABA based on nearest neighbor (NN) imputation was applied for understory trees. Our approach is intended to compensate for the weakness of LiDAR data and ITD in estimating understory trees, keeping the strength of ITD in estimating overstory trees in tree-level. We investigated the effects of three parameters on the performance of our proposed approach: smoothing of CHM, resolution of CHM, and height cutoff (a specific height that classifies trees into overstory and understory). There was no single combination of those parameters that produced the best performance for estimating stems per ha, mean tree height, basal area, diameter distribution and height distribution. The trees in the lowest LiDAR height class yielded the largest relative bias and relative root mean squared error. Although ITD and ABA showed limited explanatory powers to estimate stems per hectare and basal area, there could be improvements from methods such as using LiDAR data with higher density, applying better algorithms for ITD and decreasing distortion of the structure of LiDAR data. Automating the procedure of finding optimal combinations of those parameters is essential to expedite forest management decisions across forest landscapes using remote sensing data.

Keywords

Tree-List Generation, Individual Tree Detection, Nearest Neighbor Imputation, Parameter Sensitivity, Airborne LiDAR

1. Introduction

A tree-list provides detailed data foresters often desire for management and planning such as tree species, diameter at breast height (DBH), tree height (HT), basal area (BA) and stem volume. Field cruising has been commonly used to obtain such data. Field cruising is costly, however, and remote sensing data can be used as auxiliary information to improve the accuracy and precision of estimates in forest inventory.

Among various remote sensing techniques, airborne light detection and ranging (LiDAR) has been increasingly used in forestry applications during the last decade. LiDAR has performed well in estimating forest attributes such as biomass (Næsset & Gobakken, 2008) , diameter distribution (Gobakken & Næsset, 2004) , volume and BA (Lindberg & Hollaus, 2012) . Tree-lists have also been estimated by LiDAR (Lindberg, Holmgren, Olofsson, Wallerman, & Olsson, 2010, 2013) or aerial photographs (Temesgen, LeMay, Froese, & Marshall, 2003) .

In general, there are mainly two approaches using LiDAR data in forestry, the area-based approach (ABA) and the individual tree detection (ITD) approach (Vauhkonen, Maltamo, McRoberts, & Næsset, 2014) . ABA assumes that the vertical height distribution of laser point clouds is related to variables of interest in an area. A host of summary statistics derived from the point cloud are used to predict many forest inventory attributes. Information on the LiDAR point cloud is not fully utilized in ABA, i.e., most of the studies have focused on vertical height distribution in a sample plot and only a few studies using horizontal information obtained from the LiDAR point cloud. Pippuri, Kallio, Maltamo, Peltola, and Packalén (2012) found horizontal texture metrics from a canopy height model (CHM) could be used to predict the spatial pattern of trees, and horizontal landscape metrics from a CHM used to predict the need for first thinning.

In contrast, ITD identifies individual trees and provides estimates of forest attributes based on the identified individual trees. Although many variations exist, ITD commonly uses a rasterized CHM to segment individual trees with horizontal location of treetop and height across the CHM area. Thus, ITD has apparent advantages over ABA regarding utilization of horizontal information in LiDAR point clouds and can be more suitable for tree-level forest inventories than ABA. However, information on understory vegetation is likely to be missed when using ITD (Koch, Kattenborn, Straub, & Vauhkonen, 2014) . This is because rasterizing LiDAR point clouds into CHM means that there is a rounding effect of summarizing all the point clouds within a range of cells into one cell height value mainly focusing on higher point clouds making it difficult to detect or estimate understory vegetation. Additionally, it is well known that LiDAR has weaknesses for detecting or estimating understory vegetation regardless of the approach used because LiDAR data lack information on understory vegetation (lower proportion of point clouds in understory) (Takahashi, Yamamoto, Miyachi, Senda, & Tsuzuku, 2006) .

Many approaches have been proposed to overcome the limitations above. Maltamo, Eerikäinen, Pitkänen, Hyyppä, and Vehmas (2004) combined a theoretical probability distribution function with the tree height distribution estimated from ITD to detect small and suppressed trees. ITD first estimated the height distribution and the number of large trees. For small trees, two approaches have been used including―the complete Weibull distribution with the parameter prediction method and the left-truncated Weibull distribution with estimation of parameters from the estimated height distribution by ITD. These approaches were tested for the estimation of the height distribution and the number of trees. DBHs for large and small trees were then predicted using the relationship between DBH and LiDAR metrics. Total timber volume and stem density were finally determined by summing the estimates from the two approaches for large and small trees. Lindberg et al. (2010) proposed a methodology to generate a tree-list combining a CHM-based ITD and ABA estimation. To better detect trees that are close to each other or small: 1) the number of trees per segment was estimated using a training dataset in which the number of field-measured trees for each tree crown segment was known, and 2) a candidate tree-list from the ITD was calibrated using the target distributions of HT and DBH estimated by a k-Nearest Neighbor (NN) approach. The combined approach improved the estimation of distributions for DBH and HT, and produced unbiased estimates of forest attributes. In addition to ITD based on CHM, Lindberg et al. (2013) utilized a 3D clustering method to model a tree crown using a priori information on the shape and proportions of tree crowns. The 3D clustering method identified more trees below the tallest canopy layer and with a DBH < 20 cm than ITD based on CHM. Hamraz, Contreras, and Zhang (2017) proposed the use of vertical stratification of point clouds and LiDAR data with high point cloud density (50 points/m²), which would have more information on understory vegetation than the one with low density, to detect understory trees. The proposed approach improved detecting understory trees without affecting the overall quality of segmentation for overstory trees.

Many parameters affect the performance of tree segmentation by ITD; these can be classified into two parameters, biological and technical. For the biological parameter, Vauhkonen et al. (2012) claimed that the performance of ITD methods depends more on forest structure, stand density, and tree clustering than on detection techniques. For example, an estimated tree segment by ITD could have no, one, or several trees in it (Breidenbach, Næsset, Lien, Gobakken, & Solberg, 2010) , and trees in an understory under a dense upper canopy are hard to detect with LiDAR (Maltamo et al., 2004) . On the other hand, the methods for ITD were reported as the primary parameter affecting the performance of ITD by Kaartinen et al. (2012) . Substantial differences in the percentage of matched and missed trees, and commission error were found among the ITD methods. Also, the accuracy of determining tree location, tree height, and crown delineation changed according to the ITD methods. In contrast, pulse density showed less impact on ITD.

A typical ITD method consists of the following two steps: 1) generating a rasterized CHM with appropriate smoothing and resolution using normalized LiDAR point cloud data, and 2) tree segmentation using a segmentation technique on the rasterized CHM (finding local maxima as treetops and delineating tree crowns) (Yu, Hyyppä, Holopainen, & Vastaranta, 2010) . Therefore, the performance of ITD is affected by the parameters (smoothing and resolution for CHM, and the algorithm used for tree segmentation). In addition to these parameters, Wiggins (2017) reported that excluding trees below a specific height (minimum height cutoff) improved ITD’s accuracy for overstory trees. Maltamo, Tokola, and Lehikoinen (2003) noted that a proper value of the truncation parameter of Weibull for DBH distribution, which can be considered the same as a height cutoff, should be further studied. According to McGaughey (2016) and Wiggins (2017) , there might be an optimal parameterization that balances the smoothing of the CHM, resolution of the CHM, and the height cutoff to best identify individual trees, although Koch et al. (2014) and McGaughey (2016) pointed out that the optimal parameterization can vary over large forest areas with diverse and complicated structure. To offset the variation of the optimal parameters, Koch, Heyder, and Weinacker (2006) proposed applying different intensities of smoothing according to HT. This method would prevent under- and over-representation of local HT maxima.

Other than ITD, detailed information on forest resources, such as a tree list or stand table, has been estimated by several methods that can be mainly classified into two categories: 1) diameter distribution modeling, and 2) imputation. In diameter distribution modeling, parameters of some theoretical distributions are estimated to describe the distribution of tree diameters. Three approaches commonly used are the parameter prediction method, parameter recovery method, and quantile prediction method (Temesgen et al., 2003) . Imputation methods directly substitute measured values from sample locations (references) for locations for which a prediction is desired (targets). The distance metric used to identify suitable references and the number of references used in a single imputation (k) are the key considerations to classify the imputation methods such as most similar neighbor, gradient nearest neighbor, or Random Forest NN (RF NN hereafter) (Eskelson et al., 2009) . Temesgen et al. (2003) used a set of proxy variables to represent a tree-list in NN imputations because there is no single variable to represent the tree-list. On the other hand, Strunk et al. (2017) used plot identities as a response variable in NN imputations in evaluating NN strategies to impute a tree-list.

In our study, we combined ABA and ITD to estimate tree-list using LiDAR data inspired by the ideas from Maltamo et al. (2003) , Maltamo et al. (2004) and Wiggins (2017) . This was for overcoming the weakness of LiDAR data and the ITD method in identifying understory trees, and utilizing the strength of ITD over ABA. Maltamo et al. (2003) combined pattern recognition of single trees with the truncated Weibull distribution to estimate forest characteristics using digital video imagery. Trees were grouped into large (DBH > 17 cm) and small (DBH ≤ 17 cm) trees. The cutoff DBH value (17 cm) was the minimum size of trees that could be detected by the pattern recognition method. The value of 17 cm in DBH was used as a truncation parameter of the left-truncated Weibull. Pattern recognition was applied to large trees (DBH > 17 cm), and the diameter distribution modeling to small trees (DBH ≤ 17 cm), respectively. This idea was improved upon by Maltamo et al. (2004 ) , who combined ITD based on CHM for large trees and diameter distribution modeling for small trees. HT distribution was modeled using LiDAR metrics as auxiliary variables. Wiggins (2017) examined the effect of height cutoff on the accuracy of LiDAR data for estimating forest structure of taller trees and found that a 12 m height cutoff produced better results in estimating forest structure and spatial pattern.

For ITD, we used watershed segmentation (Vincent & Soille, 1991) for overstory trees (trees taller than a height cutoff) and ABA by NN (k = 1) imputation for understory trees (trees shorter than the height cutoff). While the performances of diameter distribution modeling depended on the results from large tree estimation by the single tree pattern recognition in Maltamo et al. (2003) or the ITD based on CHM in Maltamo et al. (2004) , in this study, we used ITD and ABA independently. They were only linked by a height cutoff when generating a complete tree-list. Whereas Lindberg et al. (2010) estimated a tree-list for all trees by an ITD method and calibrated it, our approach separated a forest stand into overstory and understory trees, then applied different methods to the overstory and understory trees, respectively. We examined the effects of the combination of the three parameters, smoothing of CHM, resolution of CHM and the height cutoff, as well as LiDAR height classification of field plots on estimating tree-lists via ITD. The explanatory power of our approach was also investigated. We evaluated the performance of generating tree-lists in terms of BA, mean HT, stems per hectare (SPH), and distributions of DBH and HT.

2. Methods

2.1. Study Area

The study area (43.02435˚N, 124.056˚W) is located in southwestern Oregon with the extent of 647,951 hectares (Figure 1). The elevation of the area ranges approximately from 20 m to 1000 m above sea level in elevation. The range of slopes in the area is 0˚ to 89.97˚. Douglas-fir (Pseudotsuga menziesii) is the dominant tree species in the study area, and other important species are western hemlock (Tsuga heterophylla), red alder (Alnus rubra), Oregon myrtle (Umbellularia californica), bigleaf maple (Acer macrophyllum), tanoak (Notholithocarpus densiflorus), western redcedar (Thuja plicata), and grand fir (Abies grandis).

2.2. Airborne LiDAR

Airborne LiDAR data were collected between April 27^th, 2008 and April 5^th, 2009 using Leica ALS50 Phase II instrumentation. The collection was acquired as logistical constraints and weather allowed. The average pulse density (the average

Figure 1. Map of study area and plots.

number of pulses returned from surfaces) was 8.10/m² for the study area. Table 1 shows the specifications for the LiDAR survey. Laser points with elevations above ground level lower than 1 m and higher than 91.44 m (300 feet) were excluded from the computation because they did not likely represent vegetation of interest (the maximum tree height measured in the field data was 88.4 m).

2.3. Field Data

Stratified sampling based on the LiDAR metrics (Hawbaker et al., 2009) was used for field data collection. Only the lands owned by the BLM or the Coquille Tribe in the study area were considered. Then, the non-forested areas were removed. Within this pre-selected area, a set of LiDAR grid metrics (22.86 m by 22.86 m) were calculated from the LiDAR point clouds. Using the principal component analysis, the 80^th percentile and standard deviation of the LiDAR height were selected as describing best the variation in forest structure in the pre-selected area. Two thousand cells were randomly selected from the cells with the pre-selected area. Based on these random samples, the range of 80^th percentile heights was subdivided into ten classes with a length of 6.10 m, and the range of standard deviations within each height class into three equal-width classes. The maximum height of the uppermost 80^th percentile class was increased to 83.52 m to cover the values of the grid cells in the full dataset. A total of 30 bins (10 × 3) were created.

Table 1. LiDAR survey specifications.

*Point on the ground vertically beneath the laser sensor on the aircraft.

Every grid cell in the pre-defined area was assigned to the bins. Then, 30 primary and 20 alternate plot locations from each bin were randomly selected from the grid cells. 30 plot locations from each bin were measured by field crews from those 50 locations using the primary plot locations unless inconsistencies were found between the LiDAR measured structure and the actual state of the forest. Such inconsistencies were caused by disturbances, such as timber harvesting, fires, or wind throw that occurred after the LiDAR data acquisition. In that case, the next available alternate plot would replace the primary plot. Plot locations overlapping roads, and in tall shrub vegetation near the coast were discarded.

Field sampling was conducted between May 25, 2010 and May 10, 2011. Nested plots with two plot sizes (12.68 m and 5.09 m) were used to measure large (both live trees with a DBH larger than 14 cm and dead trees with a height of 3.05 m or greater and a DBH of 14 cm or greater) and small (only live trees with a height taller than 1.37 m and a DBH less than 14 cm) trees, respectively. Note that only the large tree data were used for this analysis. There was one missing plot, resulting in a total of 899 plots. Table 2 and Table 3 provide a plot-level and tree-level summary of the field measurements. The ten 80^th percentile classes for the stratification sampling were used as LiDAR height classes in the current study (from “1” to “10” as height increases) to investigate the effect of LiDAR height classification of field plots on the performance of our proposed approach.

2.4. Generating Tree-Lists

The general steps of our approach are shown in Figure 2. Trees taller than a specified height (a height cutoff) were estimated by ITD using LiDAR data yielding the number and HT of the taller trees. DBHs for the taller trees were predicted based on the estimated HT using the relationship between DBH and HT from field data. For estimating the trees shorter than the height cutoff, tree-lists for target plots were first imputed with the tree-list from reference plots by RF NN imputation using both LiDAR and field data. Then, the shorter trees

Figure 2. Flowchart of the approach.

Table 2. Plot-level summary statistics of attributes from the field measurements.

*Standard deviation.

Table 3. Tree-level summary statistics from the field measurements.

*Standard deviation.

were selected from the imputed tree-lists. A complete tree-list can be generated by combining those estimated taller and shorter trees. The variables in the complete tree-list were the tree ID, HT, and DBH.

2.5. Individual Tree Detection

ITD was implemented by the function “TreeSeg” in the FUSION software (McGaughey, 2016) with the argument “ht_threshold” to estimate the tree-list for large trees. This function applies a generalized watershed segmentation algorithm by Vincent and Soille (1991) to a CHM. It should be noted that over-segmentation, known as one of the disadvantages of the watershed algorithm, may be produced with noisy imagery (Romero-Zaliz & Reinoso-Gordo, 2018) . Conceptually, the CHM is inverted, so tree crowns appear as basins. Water fills the basins from local height minima in the CHM by the algorithm, and the basins fill and join with adjacent basins, then watershed edges are established (McGaughey, 2016) . This also can be explained at the pixel level on the CHM. In every CHM pixel above a height threshold, a path is placed by iteratively moving to the neighboring pixel with the largest height value until a local height maximum is reached. A tree crown segment is defined by cells that reach the same local height maximum (Lindberg & Holmgren, 2017) . The “ht_threshold” sets minimum height (height cutoff) for tree segmentation. Fractions of CHM below this height cutoff were excluded in the segmentation process. The other two parameters, the amount of smoothing and the resolution of the CHM, were applied in generating the CHM implemented by the function “CanopyModel” in FUSION. We set levels of those three parameters as follows: 1) 3 levels of smoothing of CHM―no smoothing, median filter using a 3 by 3 neighbor window and median filter using a 5 by 5 neighbor window, 2) 24 resolutions of CHM―0.2, 0.3, …, 2.4, and 2.5 m, 3) 9 percentile height cutoffs on the LiDAR height for each plot―10^th, 20^th, …, 90^th. Because the range of HT is extensive, the LiDAR height percentiles were used as height cutoffs instead of absolute heights as in Wiggins (2017) .

After implementing ITD, we obtained a tree-list above a height cutoff including information on individual tree count, HT, a location of tree, and a number of CHM cells within a tree crown at a combination of smoothing, resolution of CHM, and height cutoff. To predict the DBHs of trees in the estimated tree-lists, an RF regression model for DBH was fitted with the HTs from the field data (16,200 trees). With this model, the DBHs of trees in the estimated tree-lists were predicted using the HTs of those trees. Then, those predicted DBHs were added to the estimated tree-lists. The model was fitted in R version 3.3.3 (R Core Team, 2017) using the R package “randomForest” (Liaw & Wiener, 2002) .

2.6. Nearest Neighbor Imputation

To estimate tree-lists for understory trees, we used RF NN imputation instead of diameter distribution modeling because there were many sample plots with multimodal or irregular shapes in diameter distribution and some plots had a small number of trees. NN imputation directly substitutes measured values from references for targets. The type of NN imputation is determined mainly by the distance metric and number of neighbors (k) (Eskelson et al., 2009) . The distance metric measures the similarity between target and reference observations, and the k indicates how many reference observations are used in a single imputation (prediction). Four distance metrics, Euclidean, Mahalanobis, most similar neighbor and RF (Breiman, 2001) , were tested. RF appeared the best for BA, SPH and error index (EI; will be defined in the following section), and Euclidean showed the best for HT (this result is not presented in this manuscript). Thus, we selected the RF algorithm as the distance metric and chose k = 1. RF builds multiple classification (or regression) trees, called forests, with bootstrap samples of training data, while selecting predictors randomly for the best split at each node in the trees. Distance in RF NN is computed as one minus the proportion of classification trees where a target observation is in the same terminal node as a reference observation (Crookston & Finley, 2008) . To estimate tree-lists by RF NN imputation, we imputed plot identities as in Strunk et al. (2017) .

To fit an NN model, it is necessary to define response and predictor variables. Predictor variables were derived from LiDAR point clouds at each filed plot location. It is not clear which a single response variable or multiple response variables should be used for estimating tree-lists because many attributes can be extracted from a tree-list. For example, Temesgen et al. (2003) used a set of 22 proxy variables to represent a tree-list. We considered several forest inventory attributes (basal area, stem volume, Lorey’s height, quadratic mean diameter, stems per ha) simultaneously to select appropriate predictor variables for estimating tree-lists via RF NN imputation. “Best subsets” was used as a variable selection method producing the best three predictors for each forest inventory attribute.

From the best predictors for each forest inventory attribute, we obtained a total of 11 predictors after removing duplicates. The selected predictors are shown in Table 4. Like leave-one-out validation, the target plot was excluded from training data when modeling. Nine different tree-lists for each height cutoff were generated from the estimated tree-lists by subtracting trees above the corresponding height cutoff. The variable selection and imputation modeling were implemented in R version 3.3.3 (R Core Team, 2017) using R packages “yaImpute” (Crookston & Finley, 2008) and “randomForest” (Liaw & Wiener, 2002) .

2.7. Performance Measures

Bias and root mean squared error (RMSE) (Walther & Moore, 2005) for mean HT, BA and SPH were computed as follows:

$\text{bias}=\frac{{\displaystyle {\sum }_{i=1}^n\left({\hat{y}}_i-y_i\right)}}{n}$ (1)

Table 4. Selected predictor variables for RF NN imputation.

$RMSE=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^n{\left({\hat{y}}_i-y_i\right)}^2}}{n}}$ (2)

where ${\hat{y}}_i$ is the prediction at the i^th plot, $y_i$ is the field-measured value at the i^th plot, and n is the number of total sample plots.

Large trees would produce greater uncertainty in estimation than small ones because the larger trees have greater values of HT, DBH, etc. To see the effect of several parameters on tree-list estimation free from the influence of greater value, relative bias (RBias) and relative RMSE (RRMSE) were also calculated for each LiDAR height class by the equations below:

$RBias\left(\%\right)=\frac{{\displaystyle {\sum }_{i=1}^n\left({\hat{y}}_{ih}-y_{ih}\right)}}{n_h}\times \frac{100}{{\bar{y}}_h}$ (3)

$RRMSE\left(\%\right)=\sqrt{{{\displaystyle {\sum }_{i=1}^n\left({\hat{y}}_{ih}-y_{ih}\right)}}^2}\times \frac{100}{{\bar{y}}_h}$ (4)

where ${\hat{y}}_{ih}$ is the prediction at the i^th plot in the h^th LiDAR height class, $y_{ih}$ is the field-measured value at the i^th plot in the h^th LiDAR height class, ${\bar{y}}_h$ is the average of field-measured values at in the h^th LiDAR height class, h is the number of LiDAR height classes, and $n_h$ is the number of sample plots in the h^th LiDAR height class.

The error index (EI) (Reynolds, Burk, & Huang, 1988) was used to evaluate the size distributions of DBH and HT, respectively. EI measures the proportions of absolute deviation between the predicted and field-measured number of trees to the total number of field-measured trees over the entire distribution. EI for a plot was computed as:

$EI\left(\%\right)=\frac{{\displaystyle {\sum }_{i=1}^k\left|n_{pi}-n_{oi}\right|}}{N}\times 100$ (5)

where $n_{pi}$ and $n_{oi}$ are the predicted and observed numbers of trees, respectively, in DBH or HT class i. k is the number of DBH or HT classes. N is the total number of field-measured trees. The bin widths for classifying DBH and HT were 10 cm and 5 m, respectively.

The coefficient of determination measures (R²) the proportion of variance in a response variable that is explained by predictor variables. It shows that how well a model’s predictions fit the observed values of the response variable, which means the actual explanatory power of the model. The R² is calculated as:

$R^2=1-\frac{{\displaystyle {\sum }_{i=1}^n{\left({\hat{y}}_i-y_i\right)}^2}}{{\displaystyle {\sum }_{i=1}^n{\left(y_i-\bar{y}\right)}^2}}$ (6)

where ${\hat{y}}_i$ is the prediction at the i^th plot, $y_i$ is the field-measured value at the i^th plot, $\bar{y}$ is the average of field-measured values of the total sample plots, and n is the number of total sample plots.

3. Results

3.1. Effects of Smoothing, Resolution, and Height Cutoff on Tree-List Estimation

All the resolutions with pixel size less than 1 m produced too large of estimates of SPH and yielded unreasonable estimates of other attributes regardless of the amount of smoothing and the height cutoff. Hence, resolutions with pixel sizes less than 1 m were dropped from the analysis. The amount of smoothing in CHM had a relatively small effect on tree-list estimation compared to the other parameters. The smoothing generally decreased the variability of estimation among the resolutions at a given height cutoff or the height cutoffs at a given resolution. For this reason, we show the performance only from the smoothing of 3 by 3 neighbor window.

Most cases of the combinations of resolution and height cutoff resulted in the underestimation of SPH (Figure 3(a)). Unbiased SPH estimations were found around 1.1 m to 2.0 m in CHM resolution with the various height cutoffs. Generally, a higher cutoff had a smaller absolute bias compared with the absolute bias from a lower cutoff. In terms of precision, Figure 3(b) shows that a higher cutoff had a relatively consistent RMSE along with resolutions in CHM, which means that higher cutoffs were less affected by resolution for SPH estimation than lower cutoffs as also shown in Figure 3(a). The combinations of the finer resolutions (1.2 ~ 1.3 m) and the lower height cutoffs (p20 and p30) provided the lowest RMSEs. For overstory trees, the patterns of performance measures were similar to the patterns from the combined approach, but the best RMSEs were always found at height cutoff p90. For understory trees, bias and RMSE increased as height cutoff increased except for the bias at height cutoff p90.

For BA estimation, bias decreased as resolution decreased as shown in Figure 4(a). Unbiased BA estimation was achieved for the combination of several cutoffs from p10 to p60 and resolutions with pixel sizes larger than 1.7 m. RMSE in

Figure 3. (a) Bias and (b) RMSE in SPH estimation: the left graph is for overstory and understory trees via the combined approach by resolution of CHM and height cutoff; the middle graph is for overstory trees via ITD by amount of smoothing, resolution of CHM and height cutoff; and the right graph is for understory trees via RF NN by height cutoff (the horizontal dashed line indicates unbiased estimates).

BA estimation also decreased as resolution decreased (Figure 4(b)). Lower cutoffs yielded lower RMSE. The lowest RMSEs appeared for resolutions around 1.8 ~ 2.0 m. For overstory trees, the differences in RMSE between height cutoffs at a given resolution were smaller than the differences for the combined approach except for 1.0 m resolution. For understory trees, bias and RMSE increased as height cutoff increased, and all height cutoffs overestimated SPH.

HT estimation had better performance than the other attributes. The pattern for HT estimation was different from the other attributes. The best accuracy in HT estimation was found with the cutoff at p50 or p60 for any resolution. The poorest accuracy in HT estimation appeared only for the cutoff p10, which had a worse bias for HT estimation as resolution decreased. HT estimation became unbiased as resolution decreased except with cutoffs p10 and p20 (Figure 5(a)). Height cutoffs showing better RMSEs were p50 and p60 with the middle and higher resolutions, and p80 in the lower resolutions at any smoothing level. RMSE increased as resolution decreased especially for cutoffs p10, p20, and p30

Figure 4. (a) Bias and (b) RMSE in BA estimation: the left graph is for overstory and understory trees via the combined approach by resolution of CHM and height cutoff; the middle graph is for overstory trees via ITD by amount of smoothing, resolution of CHM and height cutoff; and the right graph is for understory trees via RF NN by height cutoff (the horizontal dashed line indicates unbiased estimates).

(Figure 5(b)). Bias and RMSE of HT estimation for overstory trees only by ITD increased as the resolution decreased. For understory trees, Bias and RMSE for HT estimation also increased as height cutoff increased.

For the lower resolutions, the lower cutoffs showed better DBH distribution estimation than the higher cutoffs, while it was the opposite with the higher resolutions (Figure 6(a)). This pattern was also observed in HT distribution estimation. The best DBH distribution was found with cutoffs p30 and p40 and lower resolutions while cutoff p90 had the best DBH distribution for the higher resolutions. The HT distribution estimation, in most cases, had the better result with the lower cutoffs than the higher cutoffs (Figure 6(b)). The cutoff p50 had the best performance in most cases, except p90 for 1 and 1.1 m resolutions, and p30 for 1.3 ~ 1.5 m resolutions. The resolutions with medium pixel sizes were better for estimating the HT distribution. For overstory trees, EI for DBH decreased as resolution decreased, and the lowest height cutoff p10 always yielded the best DBH distribution estimation at every resolution. DBH distribution

Figure 5. (a) Bias and (b) RMSE in HT estimation: the left graph is for overstory and understory trees via the combined approach by resolution of CHM and height cutoff; the middle graph is for overstory trees via ITD by amount of smoothing, resolution of CHM and height cutoff; and the right graph is for understory trees via RF NN by height cutoff (the horizontal dashed line indicates unbiased estimates).

estimation for overstory trees was poorer than for both overstory and understory trees. The best height cutoff for HT distribution of overstory trees estimation increased as resolution decreased. For understory trees, EIs for DBH and HT were reduced as height cutoff increased except for cutoff p80. Contrary to HT estimation for both overstory and understory trees by the combined approach, the best height cutoffs in the estimation of the HT distribution for overstory trees by ITD was for higher cutoffs from p60 to p80 except for resolutions higher than 1.4 m.

Compared to the combined approach for all trees or the ITD for overstory trees, NN imputation produced much lower biases for understory trees’ SPH, BA, and HT (Figures 3-5). The smallest biases for understory trees for SPH, BA and HT estimation were found at cutoffs p10, p20, and p40, respectively. The smallest RMSEs in the three attributes were observed only at cutoff p10.

3.2. Effects of Classification of Field Plots by LiDAR Height on Tree-List Estimation

The absolute and relative performance measures separated by LiDAR height class were calculated for each forest attribute estimated. The smallest group,

Figure 6. (a) EI for DBH and (b) EI for HT: the left graph is for overstory and understory trees via the combined approach by resolution of CHM and height cutoff; the middle graph is for overstory trees via ITD by amount of smoothing, resolution of CHM and height cutoff; and the right graph is for understory trees via RF NN by height cutoff.

class 1, showed distinct properties in those performances. For the absolute measures, such as bias and RMSE, lower LiDAR height classes, especially the lowest class, generally yielded comparable or better performances for BA and SPH than the higher classes. However, based on the relative measures, the lowest class had much poorer results. Similar patterns were found in EIs for DBH and HT as well. The effect of the amount of smoothing in CHM by LiDAR height class was relatively small. The performances by height cutoff in a given resolution were averaged for this section because it is better to show the general effect of height class on tree-list estimation performance. For SPH estimation (Figure 7), bias decreased as resolution decreased for every height class, but the resolutions showing unbiasedness varied among height classes. Lower height classes had larger variability in bias among resolutions than higher height class. Height class 1 had much larger RBias at higher resolutions than the other height classes. Larger RMSE occurred in height classes 1 through 6, and the largest RMSE was found in height class 3. RRMSE in height class 1 was largest at every resolution. Relatively larger RRMSEs at higher resolutions were observed in the taller height classes.

Figure 7. Bias, RBias, RMSE, and RRMSE for SPH estimation via the combined approach by LiDAR height class and resolution of CHM: the values of each performance by height cutoff in a given resolution are averaged.

In BA estimation (Figure 8), biases in the taller height classes were generally larger than biases in the shorter height classes. This pattern was similar for RBias except for height class 1. RBias in height class 1 was larger than the other height classes at resolutions less than or equal to 2.3 m. Lower height classes generally had smaller RMSE than higher height classes, but height class 1 had a much larger RRMSE than the other height classes. Figure 9 shows the performance measures for HT estimation by height class. The pattern of HT estimation among height classes was different from the pattern of SPH and BA estimation. Height class 1 had comparable or better performance in bias, RBias, and RMSE. The primary difference in bias and RBias between class 1 and the other classes was that class 1 mainly underestimated HT while the other classes overestimated. RRMSE for HT in height class 1 had slightly larger values than RRMSE from other height classes. Estimated distributions of DBH and HT for height class 1 were much poorer than the distributions for the other classes. Except class 1, lower height classes showed better performance in EIs for both DBH and HT than higher height classes. Lower resolution generally had lower EIs (Figure 10).

Figure 8. Bias, RBias, RMSE, and RRMSE for BA estimation via the combined approach by LiDAR height class and resolution of CHM: the values of each performance by height cutoff in a given resolution are averaged.

3.3. Explanatory Power of Individual Tree Detection for Overstory Trees and Random Forest Nearest Neighbor Imputation for Understory Trees

Tables 5-7 show R²s for SPH, BA and HT estimation for trees over a given height cutoff (overstory trees) via ITD by resolution of CHM and height cutoff with smoothing using a 3 by 3 window. For SPH estimation (Table 5), the best R² was found at resolutions between 1.2 m and 1.7 m for each height cutoff. Height cutoff p90 yielded the largest R², 0.501, and the best R² decreased as the height cutoff decreased. The lowest height cutoff p10 had negative R² at all the resolutions. BA estimation by ITD showed poor explanatory power for overstory (Table 6). Most combinations of resolutions and height cutoffs had negative R²s, and the best R² was 0.338 with the resolution 2.0 m and the height cutoff p10. Larger height cutoffs, from p70 to p90 provided negative R² at every resolution. In HT estimation (Table 7), the 1.0 m resolution yielded the best R² at every height cutoff except p90. The middle height cutoffs, p50 or p60, had better R² than the other height cutoffs. Inferior explanatory power was found at height cutoffs p10 and p90. The explanatory power for HT estimation generally decreased as resolution increased.

Figure 9. Bias, RBias, RMSE, and RRMSE for HT estimation via the combined approach by LiDAR height class and resolution of CHM: the values of each performance by height cutoff in a given resolution are averaged.

Figure 10. EIs for DBH and HT estimation via the combined approach by LiDAR height class and resolution of CHM: the values of each performance by height cutoff in a given resolution are averaged.

Table 8 shows the explanatory power of RF NN imputation for trees under a given height cutoff (understory trees). For HT estimation, the R²s were around 0.5. However, the R²s for BA and SPH estimation were much poorer than the R²s for HT estimation or even had negative values. For understory trees for each

Table 5. Explanatory power (R²) of SPH via ITD for trees taller than given height cutoffs with the smoothing of 3 by 3 neighbor window.

Table 6. Explanatory power (R²) of BA via ITD for trees taller than given height cutoffs with the smoothing of 3 by 3 neighbor window.

forest inventory attribute, the scatter plots of observed vs. predicted via RF NN imputation did not show any anomaly. The lower height cutoff we used, the more observations with zero values we had. The prediction results for those observations with zero values were inferior for every height cutoff.

4. Discussion

No single combination of smoothing, resolution and height cutoff was found to produce the best results for all performance measures (Table 9). Koch et al. (2014) and McGaughey (2016) also reported similar findings. Similarly, ITD’s performance varied depending on the algorithm used to delineate trees in the CHM (Kaartinen et al., 2012) . Differences in performance between the lowest LiDAR height class and the other classes were found based on both absolute and

Table 7. Explanatory power (R²) of HT via ITD for trees taller than given height cutoffs with the smoothing of 3 by 3 neighbor window.

Table 8. Performance measures of RF NN imputation by inventory attributes for trees shorter than given height cutoffs.

*Standard deviation of field-measured inventory attribute under given height cutoffs.

relative performance measures. Kaartinen et al. (2012) reported that the HT class did not generally impact the accuracy of HT estimation, but greater uncertainty was observed for ITD methods capable of finding small trees. According to Hopkinson et al. (2005) , vegetation classes with short height, such as low shrub and aquatic vegetation, yielded the largest relative errors in canopy height estimation, whereas tall vegetation classes showed the largest absolute errors. The low level of penetration of LiDAR returns into the sub-canopy surface might be an essential reason for the high relative bias for low shrub and aquatic vegetation. For aquatic vegetation, it was also believed that the weak laser backscatter from the saturated ground caused the high relative bias. These results were very similar to ours although the smallest height class in our research almost exclusively consisted of trees.

As we reported above, resolutions with pixel sizes less than 1.0 m were dropped in the analysis because it yielded unreasonably large SPH estimations.

Table 9. Best performance for each assessment by estimation method.

*The first argument indicates the amount of smoothing, the second resolution in CHM, and the third percentile height cutoff for the combined method. †This represents percentile height cutoff.

Pouliot, King, Bell, and Pitt (2002) claimed that in high-resolution imagery, tree detection and crown delineation became more complicated. This is because high-resolution imagery can display very detailed objects such as branches causing tree crowns to deviate from the conic shape. Thus, more tree crowns could be estimated at higher image resolutions. Conversely, in low-resolution imagery, it is more challenging to identify crown boundaries because they become less distinct. Another reason for our large SPH estimation might be data pits, which are height irregularities in a CHM. The function ‘CanopyModel’ in FUSION used for generating CHMs in our study fills pixels without LiDAR point clouds using an eight-way search and a distance-weighted average (McGaughey, 2016) . However, it might be difficult to avoid irregularities in height on a CHM if laser pulses used for our LiDAR data acquisition penetrated deeply into tree crowns causing large height variations within individual tree crowns (Persson, Holmgren, & Soderman, 2002) . Image smoothing with various filters using mean, median, or Gaussian approaches have been applied to reduce data pits (Persson et al., 2002; Yu, Hyyppä, Vastaranta, Holopainen, & Viitala, 2011) . In our study, the smoothing did not work well at resolutions with pixel sizes larger than 1.0 m although the smoothing using a 5 by 5 window showed smaller SPH estimation than no smoothing and the smoothing with a 3 by 3 window. A pit-free CHM proposed by Khosravipour, Skidmore, Isenburg, Wang, and Hussin (2014) was found to improve the accuracy of tree detection based on either high or low-density LiDAR data; however, this approach could help solve our large SPH estimation at the finer resolutions.

The ratio of average crown diameter to image pixel size was proposed as a guide to determine an optimal image resolution for tree detection and crown delineation using digital camera imagery (Pouliot et al., 2002) . With a small crown diameter to pixel size ratio, it is hard to have distinct crown boundaries in an image, resulting in under-segmentation. However, a large crown diameter to pixel size ratio might cause high variability within a crown in an image resulting in over-segmentation. Although our data from field surveys do not have information on crown diameter, there might be significant variations in the tree crowns considering the diversity of forest stands in our study area. This might be one of the reasons why the high CHM resolutions overestimated SPH in our study. Barnes et al. (2017) found that no single CHM resolution produced the best performance of ITD for both healthy and diseased larch trees, claiming that not only the tree crown size but also the maximum tree height governed an optimal size of CHM resolution. The performance of ITD with high-resolution CHMs (0.15 m) was best for plots with low maximum height (<20 m), and the performance with low-resolution CHMs (0.5 m) was best for plots with high maximum height (>30 m).

LiDAR point cloud density might be related to the optimal CHM resolution as with the tree crown diameter. With LiDAR data of high point cloud density, high CHM resolution could yield high within-crown variations on a CHM. Inversely, with LiDAR data of low point cloud density, low CHM resolution could produce less distinct crown boundaries making it difficult to identify tree crowns. The high CHM resolutions should have yielded good performance in that the LiDAR data used for this study had low point cloud density. However, the high diversity of forest stands in the study area might add more within-crown variations. Even though an optimal resolution of CHM was set based on the crown diameter to CHM resolution ratio, it should be noted that the performance of ITD was still affected by LiDAR point cloud density for trees with small DBH (<20 cm) as reported in Khosravipour et al. (2014) .

The results of large SPH estimation are quite different from previous studies. Stere?czak, Będkowski, and Weinacker (2008) reported that the 0.25 and 0.5 m resolutions in CHMs were better than the 1.0 m resolution for estimating SPH through individual tree delineation based on a similar method to Heurich and Weinacker (2004) . It was found that the number of detected trees decreased as the resolution of CHM decreased (Stere?czak et al., 2008) , and this was also observed in our work, excluding height cutoffs p10 and p20. Smreček et al. (2018) showed very similar results to ours for SPH estimation based on ITD. At the highest resolution (0.5 m), the number of trees identified was hugely overestimated; the number of trees identified decreased as the CHM resolution decreased from 0.5 m to 2.0 m, as was the case in our study. The optimal resolutions for tree identification were 1.0 and 1.5 m depending on the sample plot. Smreček et al. (2018) claimed that this was because the CHM with 0.5 m resolution was too detailed. We observed many estimated trees from ITD with extreme small areas compared to their estimated heights. Those trees should have been removed from the estimated tree-list using an appropriate criterion. With this filtering process, overestimation at high resolutions would be decreased.

Most combinations of parameters resulted in underestimating SPH. According to Lindberg et al. (2010) , ITD underestimates SPH because ITD often misses trees below dominant trees or recognizes trees close to each other as one tree. It was expected that there would be more underestimation as pixel size increased. The larger pixel size we have, the more aggregated information we would get, so lower resolution also would result in underestimating SPH. For this reason, estimates of BA decreased as resolution decreased. The approach of Lindberg et al. (2010) could give an improvement for estimating overstory trees for our study. Considering that most of the combinations for overstory trees by ITD produced negative biases in SPH estimation in our study (Figure 3(a)), estimating of the number of trees per segment would improve the negative biases in SPH estimation by increasing the number of detected trees.

Performance measures for HT estimation were better than measures for the other variables tested. This might be because LiDAR directly measures heights of target objects, so there is less uncertainty in height estimation than other attribute estimation. According to Stere?czak et al. (2008) , there was no difference between the three resolutions (0.25, 0.5, and 1.0 m) in CHM for HT estimation. For understory trees, biases in HT estimation less than 0.15 m in absolute value were produced by RF NN at every height cutoff. The higher the height cutoff applied, the larger the RMSE obtained. This is attributed to the fact that RF NN will have more and larger trees to estimate with higher height cutoffs.

While RF NN imputation showed better performance in estimating SPH than the combined approach and tree segmentation (Figures 3-6), this does not mean that RF NN imputation is better than the combined approach or ITD. It is because the target trees for those two methods are different from each other (tall trees above a height cutoff for ITD and short trees below the height cutoff for RF NN imputation). Therefore, the values dealt with in RF NN imputation were smaller than ITD. Based on relative measures not included on this manuscript, RF NN imputation was generally better in RBias, comparable in EIs, and worse in RRMSE.

The errors for BA and mean HT estimation in taller height classes were larger than in shorter height classes contradicting the fact that airborne LiDAR has difficulty in detecting understory vegetation. This might be because large trees have larger DBH and HT than small trees. To offset this potential issue, relative performance measures such as RBias and RRMSE were calculated. These relative measures revealed that the performance of estimation in shorter height classes was poorer than for the trees in taller height classes. Stere?czak et al. (2008) found a similar phenomenon for young stands.

There was also no single combination of the three parameters tested for explanatory powers that proved best overall. While HT estimation was good, estimation of BA and SPH were poor. Especially, BA estimation was very poor. The negative R² indicates (Tables 5-8) our results were worse than the mean value of the data. However, the combinations of parameters for each forest attribute could be a partial guide of generating tree-lists for the forest attributes. First, for overstory trees, R² for SPH had larger values with the finer resolution and the higher height cutoff. On the contrary to SPH, R² for BA had larger values with the coarser resolution and the lower height cutoff. R² for HT had larger values with the finer resolution and the middle height cutoff. For understory trees, R² for SPH had larger values with lower height cutoff, and R² for BA with higher height cutoff. BA estimation for overstory trees via ITD had more uncertainty sources than the other attributes, including SPH estimation and subsequent prediction of DBH for each detected individual tree (estimated HT used to predict DBH provided additional uncertainty source to the DBH prediction). These uncertainty sources might partially explain the poor performance in BA estimation. Utilizing the limited information in LiDAR data might affect the poor performance for the explanatory powers. We used CHM-based ITD; this method has limitation summarizing LiDAR point clouds within a range of cell into one cell height value regardless of generating a pit-free CHM. Instead, 3D ITD methods have been recently studied using information in LiDAR as much as possible (Kandare, Ørka, Chan, & Dalponte, 2016) . However, the 3D ITD methods required more complex algorithms to implement, and also processing time could be a new parameter to consider (Pirotti, Kobal, & Roussel, 2017) .

It is well known that it is difficult to estimate characteristics of understory vegetation. Eskelson, Madsen, Hagar, and Temesgen (2011) used beta regression to estimate percent shrub cover, and it yielded poor explanatory power. Rahman and Gorte (2008) developed a tree filtering technique to separate dominant tree and undergrowth vegetation, but it was found difficult to separate undergrowth vegetation very close to a tree using the filtering. Liu, Shen, Zhao, and Xu (2013) suggested a method to extract individual tree crowns from airborne LiDAR in residential areas showing promising applications, but also reported that small trees were omitted if there were an only small number of points representing them in the dataset. Our results for understory trees via RF NN were not good (Table 8). To improve NN estimation with LiDAR data having low point cloud density, we investigated many LiDAR metrics such as metrics from LiDAR point clouds under several height cutoffs as Wing et al. (2012) proposed to estimate understory vegetation cover with airborne LiDAR. Some of the metrics from understory point clouds were selected for NN imputation (Table 4). However, it did not greatly improve the performance of NN imputation compared to NN imputation without those metrics (not presented here). This might fundamentally be because our LiDAR data lacked information on understory vegetation.

In NN imputation, one of the critical parameters is the selection of a number of neighbors for imputation modeling or distance metrics used to measure the similarity between the reference and target plot using auxiliary variables (Eskelson et al., 2009) . While their result varied among different forest types, Strunk et al. (2017) reported that k = 3 and Mahalanobis distance metric produced better performance over other NN strategies in estimating tree-lists. In this study, we only used k = 1 and RF as a distance metric in NN modeling. Combination of the two parameters needs to be examined for understory vegetation. In addition to these two factors, implementing variable selection procedure for NN imputation to each LiDAR height class could have the potential to improve NN modeling performance.

Compared to the results of HT estimation, results related to DBH estimation such as BA and EI for DBH showed poorer performance. It is known that predicting tree-level DBH from height-derived metrics has considerable variability (Matti Maltamo & Gobakken, 2014) . Kaartinen et al. (2012) reported that estimation of DBH based on HT and crown size would have considerable uncertainty because allometric equations used for estimating DBH are sensitive to errors in input data such as the size of tree crown or HT. Another potential reason is the dead trees in the field data. The Pearson’s correlation coefficients between the field-measured HT and DBH for live and dead trees are 0.771 and 0.212, respectively. Even though the dead trees account for only 8.8% of a total number of field-measured trees, appropriate handling for dead trees would give opportunities to improve estimating tree-lists.

The scanning angle is another parameter to consider for LiDAR projects (Gatziolis & Andersen, 2008) . If the scanning angle increases, it facilitates changes in pulse propagation direction and increases the distance the pulse moves through the canopy. The change in pulse direction and the increased distance are related to LiDAR data artifacts such as returns below the ground. Therefore, with a wide scanning angle, LiDAR data might have more data artifacts than with a narrow-angle. Additionally, these data artifacts could increase when data acquisition is carried out on a slope, as an off-nadir scanning angle increases on the slope (Gatziolis & Andersen, 2008) .

39.4% of our field plots had slopes more than 30˚ based on digital terrain models from the study site. Khosravipour, Skidmore, Wang, Isenburg, and Khoshelham (2015) showed that normalized LiDAR point clouds could distort tree locations detected from CHM and height estimation depending on the steepness of slope and crown shape. For the slope of more than 30˚ 44.6% of correctly detected trees with wider and irregular crown shapes were affected by the horizontal and vertical displacements. They suggested using a non-normalized CHM to avoid the adverse effect of the distortion by steep slopes, especially in a heterogeneous forest with multiple species. The slope was also found to affect the ABA approach by distorting heights of LiDAR point clouds (Hansen et al., 2017) . They proposed two methods, Procrustean transformation and histogram matching, to counter the distortion of LiDAR point clouds on slope terrain for extracting LiDAR metrics. These point cloud distortions by slope terrain could worsen our results for both overstory and understory estimations.

Another issue is that there was the time lag between LiDAR acquisition and field surveys. This might have the potential source of error, particularly for younger fast-growing stands. Also, there were seasonal differences in the LiDAR acquisition dates (e.g., April through June in the spring, June through August in the summer, and September and October in the fall). According to Gatziolis and Andersen (2008) , the seasonal differences can induce considerable variability in canopy penetrability by LiDAR pulses especially for deciduous forests (e.g., leaf-on and leaf-off conditions) and weather-related limitations. The variability in canopy penetrability might increase uncertainty in modeling forest attributes, and the weather-related limitations could make it difficult to keep the quality of LiDAR data consistent over our whole study area. Time windows, part of LiDAR data acquisition considerations in Gatziolis and Andersen (2008) , should be carefully planned according to project objectives.

5. Conclusion

We proposed an approach to combine ITD and ABA to generate a tree-list using airborne LiDAR data and field measured data. The approach aimed to compensate for the disadvantage of LiDAR data and ITD in estimating understory trees, and to keep the strength of ITD in estimating overstory trees in tree-level. The selected parameters, smoothing, resolution and height cutoff, were examined to determine how they affected the performance of the proposed approach. There was no single combination of the three parameters that provides the best estimation results for all the forest attributes in this study. For each attribute, the best results depended on different combinations of those parameters. This is concurrent with what Koch et al. (2014) and McGaughey (2016) reported. However, our study provided the ranges and patterns of the selected parameters that yielded better performance results for each forest attribute, which could be a partial guide of estimating tree-lists using airborne LiDAR. It would be practical and useful to determine how to automatically find the optimal combinations of those parameters across the forest landscape using remote sensing data. In addition to the three parameters tested in the present study, the automation for the optimal combinations would require considering additional parameters such as forest types, tree species, tree-size parameters (tree crown width or maximum tree height) and topography.

There are several topics for further study to improve the combined approach. A denser point cloud data would have more information on both overstory and understory vegetation in a forest, thus could increase the combined approach’s performance. The algorithm used to generate a CHM and to delineate trees on the CHM is another critical parameter in ITD. Comparison of different algorithms for processing the CHM is an active area of research. Estimating the number of trees per crown segment would help obtain unbiased SPH estimation. A point cloud based ITD method could lead to improvement by utilizing more information in LiDAR data. A minimum crown area by ITD should be examined so that tiny crown would not degrade the quality of the predicted tree-lists. The effect of slope on CHM generation and LiDAR metrics extraction need to be considered for better estimation. Fusing ITD and ABA to predict overstory and understory vegetation shown in this research indicates that forest analysts can benefit from the predictive abilities of the imputation approach and the quality information provided by LiDAR. In that, the approach presented herein can be sufficient for strategic inventory purposes.

Acronyms

ABA: Area-based approach

BA: Basal area

CHM: Canopy height model

DBH: Diameter at breast height

SPH: Stems per hectare

EI: Error index

HT: Tree height

ITD: Individual tree detection

NN: Nearest neighbor

LiDAR: Light detection and ranging

R²: Coefficient of determination

RBias: Relative bias

RF: Random forest

RMSE: Root mean squared error

RRMSE: Relative root mean squared error

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Barnes, C., Balzter, H., Barrett, K., Eddy, J., Milner, S., & Suárez, J. (2017). Individual Tree Crown Delineation from Airborne Laser Scanning for Diseased Larch Forest Stands. Remote Sensing, 9, 231. https://doi.org/10.3390/rs9030231
[2]	Breidenbach, J., Næsset, E., Lien, V., Gobakken, T., & Solberg, S. (2010). Prediction of Species Specific Forest Inventory Attributes Using a Nonparametric Semi-Individual Tree Crown Approach Based on Fused Airborne Laser Scanning and Multispectral Data. Remote Sensing of Environment, 114, 911-924. https://doi.org/10.1016/j.rse.2009.12.004
[3]	Breiman, L. (2001). Random Forests. Machine Learning, 45, 5-32. https://doi.org/10.1023/A:1010933404324
[4]	Crookston, N. L., & Finley, A. O. (2008). Yaimpute: An R Package for Knn Imputation. Journal of Statistical Software, 23, 16. https://doi.org/10.18637/jss.v023.i10
[5]	Eskelson, B. N. I., Madsen, L., Hagar, J. C., & Temesgen, H. (2011). Estimating Riparian Understory Vegetation Cover with Beta Regression and Copula Models. Forest Science, 57, 212-221.
[6]	Eskelson, B. N. I., Temesgen, H., Lemay, V., Barrett, T. M., Crookston, N. L., & Hudak, A. T. (2009). The Roles of Nearest Neighbor Methods in Imputing Missing Data in Forest Inventory and Monitoring Databases. Scandinavian Journal of Forest Research, 24, 235-246. https://doi.org/10.1080/02827580902870490
[7]	Gatziolis, D., & Andersen, H.-E. (2008). A Guide to Lidar Data Acquisition and Processing for the Forests of the Pacific Northwest (Gen. Tech. Rep. PNW-GTR-768) (32 p). Portland, OR: US Department of Agriculture, Forest Service, Pacific Northwest Research Station. https://doi.org/10.2737/PNW-GTR-768
[8]	Gobakken, T., & Næsset, E. (2004). Estimation of Diameter and Basal Area Distributions in Coniferous Forest by Means of Airborne Laser Scanner Data. Scandinavian Journal of Forest Research, 19, 529-542. https://doi.org/10.1080/02827580410019454
[9]	Hamraz, H., Contreras, M. A., & Zhang, J. (2017). Vertical Stratification of Forest Canopy for Segmentation of Understory Trees within Small-Footprint Airborne Lidar Point Clouds. ISPRS Journal of Photogrammetry and Remote Sensing, 130, 385-392. https://doi.org/10.1016/j.isprsjprs.2017.07.001
[10]	Hansen, E., Ene, L., Gobakken, T., Ørka, H., Bollandsås, O., & Næsset, E. (2017). Countering Negative Effects of Terrain Slope on Airborne Laser Scanner Data Using Procrustean Transformation and Histogram Matching. Forests, 8, 401. https://doi.org/10.3390/f8100401
[11]	Hawbaker, T. J., Keuler, N. S., Lesak, A. A., Gobakken, T., Contrucci, K., & Radeloff, V. C. (2009). Improved Estimates of Forest Vegetation Structure and Biomass with a Lidar-Optimized Sampling Design. Journal of Geophysical Research: Biogeosciences, 114, G00E04. https://doi.org/10.1029/2008JG000870
[12]	Heurich, M., & Weinacker, H. (2004). Automated Tree Detection and Measurement in Temperate Forests of Central Europe Using Laser Scanning Data. Freiburg: The ISPRS Working Group on Laser-Scanners for Forest and Landscape Assessment.
[13]	Hopkinson, C., Chasmer, L. E., Sass, G., Creed, I. F., Sitar, M., Kalbfleisch, W., & Treitz, P. (2005). Vegetation Class Dependent Errors in Lidar Ground Elevation and Canopy Height Estimates in a Boreal Wetland Environment. Canadian Journal of Remote Sensing, 31, 191-206. https://doi.org/10.5589/m05-007
[14]	Kaartinen, H., Hyyppä, J., Yu, X., Vastaranta, M., Hyyppä, H., Kukko, A., Wu, J.-C. et al. (2012). An International Comparison of Individual Tree Detection and Extraction Using Airborne Laser Scanning. Remote Sensing, 4, 950. https://doi.org/10.3390/rs4040950
[15]	Kandare, K., Ørka, H. O., Chan, J. C.-W., & Dalponte, M. (2016). Effects of Forest Structure and Airborne Laser Scanning Point Cloud Density on 3d Delineation of Individual Tree Crowns. European Journal of Remote Sensing, 49, 337-359. https://doi.org/10.5721/EuJRS20164919
[16]	Khosravipour, A., Skidmore, A. K., Isenburg, M., Wang, T., & Hussin, Y. A. (2014). Generating Pit-Free Canopy Height Models from Airborne Lidar. Photogrammetric Engineering & Remote Sensing, 80, 863-872. https://doi.org/10.14358/PERS.80.9.863
[17]	Khosravipour, A., Skidmore, A. K., Wang, T., Isenburg, M., & Khoshelham, K. (2015). Effect of Slope on Treetop Detection Using a Lidar Canopy Height Model. ISPRS Journal of Photogrammetry and Remote Sensing, 104, 44-52. https://doi.org/10.1016/j.isprsjprs.2015.02.013
[18]	Koch, B., Heyder, U., & Weinacker, H. (2006). Detection of Individual Tree Crowns in Airborne Lidar Data. Photogrammetric Engineering & Remote Sensing, 72, 357-363. https://doi.org/10.14358/PERS.72.4.357
[19]	Koch, B., Kattenborn, T., Straub, C., & Vauhkonen, J. (2014). Segmentation of Forest to Tree Objects. In M. Maltamo, E. Næsset, & J. Vauhkonen (Eds.), Forestry Applications of Airborne Laser Scanning: Concepts and Case Studies (pp. 89-112). Dordrecht: Springer Netherlands. https://doi.org/10.1007/978-94-017-8663-8_5
[20]	Liaw, A., & Wiener, M. (2002). Classification and Regression by Randomforest. R News, 2, 18-22.
[21]	Lindberg, E., & Hollaus, M. (2012). Comparison of Methods for Estimation of Stem Volume, Stem Number and Basal Area from Airborne Laser Scanning Data in a Hemi-Boreal Forest. Remote Sensing, 4, 1004-1023. https://doi.org/10.3390/rs4041004
[22]	Lindberg, E., & Holmgren, J. (2017). Individual Tree Crown Methods for 3D Data from Remote Sensing. Current Forestry Reports, 1-13. https://doi.org/10.1007/s40725-017-0051-6
[23]	Lindberg, E., Holmgren, J., Olofsson, K., Wallerman, J., & Olsson, H. (2010). Estimation of Tree Lists from Airborne Laser Scanning by Combining Single-Tree and Area-Based Methods. International Journal of Remote Sensing, 31, 1175-1192. https://doi.org/10.1080/01431160903380649
[24]	Lindberg, E., Holmgren, J., Olofsson, K., Wallerman, J., & Olsson, H. (2013). Estimation of Tree Lists from Airborne Laser Scanning Using Tree Model Clustering and K-Msn Imputation. Remote Sensing, 5, 1932-1955. https://doi.org/10.3390/rs5041932
[25]	Liu, J., Shen, J., Zhao, R., & Xu, S. (2013). Extraction of Individual Tree Crowns from Airborne Lidar Data in Human Settlements. Mathematical and Computer Modelling, 58, 524-535. https://doi.org/10.1016/j.mcm.2011.10.071
[26]	Maltamo, M., & Gobakken, T. (2014). Predicting Tree Diameter Distributions. In M. Maltamo, E. Næsset, & J. Vauhkonen (Eds.), Forestry Applications of Airborne Laser Scanning: Concepts and Case Studies (pp. 177-191). Dordrecht: Springer Netherlands. https://doi.org/10.1007/978-94-017-8663-8_9
[27]	Maltamo, M., Eerikäinen, K., Pitkänen, J., Hyyppä, J., & Vehmas, M. (2004). Estimation of Timber Volume and Stem Density Based on Scanning Laser Altimetry and Expected Tree Size Distribution Functions. Remote Sensing of Environment, 90, 319-330. https://doi.org/10.1016/j.rse.2004.01.006
[28]	Maltamo, M., Tokola, T., & Lehikoinen, M. (2003). Estimating Stand Characteristics by Combining Single Tree Pattern Recognition of Digital Video Imagery and a Theoretical Diameter Distribution Model. Forest Science, 49, 98-109.
[29]	McGaughey, R. J. (2016). Fusion/Ldv: Software for Lidar Data Analysis and Visualization, Forest Service. Pacific Northwest Research Station, United States Department of Agriculture.
[30]	Næsset, E., & Gobakken, T. (2008). Estimation of Above- and Below-Ground Biomass across Regions of the Boreal Forest Zone Using Airborne Laser. Remote Sensing of Environment, 112, 3079-3090. https://doi.org/10.1016/j.rse.2008.03.004
[31]	Persson, A., Holmgren, J., & Soderman, U. (2002). Detecting and Measuring Individual Trees Using an Airborne Laser Scanner. Photogrammetric Engineering and Remote Sensing, 68, 925-932.
[32]	Pippuri, I., Kallio, E., Maltamo, M., Peltola, H., & Packalén, P. (2012). Exploring Horizontal Area-Based Metrics to Discriminate the Spatial Pattern of Trees and Need for First Thinning Using Airborne Laser Scanning. Forestry, 85, 305-314. https://doi.org/10.1093/forestry/cps005
[33]	Pirotti, F., Kobal, M., & Roussel, J. R. (2017). A Comparison of Tree Segmentation Methods Using Very High Density Airborne Laser Scanner Data. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 42, 285-290. https://doi.org/10.5194/isprs-archives-XLII-2-W7-285-2017
[34]	Pouliot, D. A., King, D. J., Bell, F. W., & Pitt, D. G. (2002). Automated Tree Crown Detection and Delineation in High-Resolution Digital Camera Imagery of Coniferous Forest Regeneration. Remote Sensing of Environment, 82, 322-334. https://doi.org/10.1016/S0034-4257(02)00050-0
[35]	R Core Team (2017). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. http://www.R-project.org/
[36]	Rahman, M. Z. A., & Gorte, B. (2008). Tree Filtering for High Density Airborne Lidar Data. In Proceedings of Silvilaser 2008: 8th International Conference on LiDAR Applications in Forest Assessment and Inventory (pp. 544-553). Edinburgh: Heriot-Watt University.
[37]	Reynolds, M. R., Burk, T. E., & Huang, W.-C. (1988). Goodness-of-Fit Tests and Model Selection Procedures for Diameter Distribution Models. Forest Science, 34, 373-399.
[38]	Romero-Zaliz, R., & Reinoso-Gordo, J. F. (2018). An Updated Review on Watershed Algorithms. In C. Cruz Corona (Ed.), Soft Computing for Sustainability Science (pp. 235-258). Cham: Springer International Publishing. https://doi.org/10.1007/978-3-319-62359-7_12
[39]	Smrecek, R., Michnová, Z., Sackov, I., Danihelová, Z., Levická, M., & Tucek, J. (2018). Determining Basic Forest Stand Characteristics Using Airborne Laser Scanning in Mixed Forest Stands of Central Europe. iForest—Biogeosciences and Forestry, 11, 181-188.
[40]	Stereńczak, K., Bedkowski, K., & Weinacker, H. (2008). Accuracy of Crown Segmentation and Estimation of Selected Trees and Forest Stand Parameters in Order to Resolution of Used Dsm and Ndsm Models Generated from Dense Small Footprint Lidar Data. Beijing: The ISPRS Congress.
[41]	Strunk, J., Gould, P., Packalen, P., Poudel, K., Andersen, H.-E., & Temesgen, H. (2017). An Examination of Diameter Density Prediction with K-Nn and Airborne Lidar. Forests, 8, 444. https://doi.org/10.3390/f8110444
[42]	Takahashi, T., Yamamoto, K., Miyachi, Y., Senda, Y., & Tsuzuku, M. (2006). The Penetration Rate of Laser Pulses Transmitted from a Small-Footprint Airborne Lidar: A Case Study in Closed Canopy, Middle-Aged Pure Sugi (Cryptomeria Japonica D. Don) and Hinoki Cypress (Chamaecyparis Obtusa Sieb. Et Zucc.) Stands in Japan. Journal of Forest Research, 11, 117-123. https://doi.org/10.1007/s10310-005-0189-0
[43]	Temesgen, H., LeMay, V. M., Froese, K. L., & Marshall, P. L. (2003). Imputing Tree-Lists from Aerial Attributes for Complex Stands of South-Eastern British Columbia. Forest Ecology and Management, 177, 277-285. https://doi.org/10.1016/S0378-1127(02)00321-3
[44]	Vauhkonen, J., Ene, L., Gupta, S., Heinzel, J., Holmgren, J., Pitkänen, J., Maltamo, M. et al. (2012). Comparative Testing of Single-Tree Detection Algorithms under Different Types of Forest. Forestry: An International Journal of Forest Research, 85, 27-40. https://doi.org/10.1093/forestry/cpr051
[45]	Vauhkonen, J., Maltamo, M., McRoberts, R. E., & Næsset, E. (2014). Introduction to Forestry Applications of Airborne Laser Scanning. In M. Maltamo, E. Næsset, & J. Vauhkonen (Eds.), Forestry Applications of Airborne Laser Scanning: Concepts and Case Studies (pp. 1-16). Dordrecht: Springer Netherlands. https://doi.org/10.1007/978-94-017-8663-8_1
[46]	Vincent, L., & Soille, P. (1991). Watersheds in Digital Spaces: An Efficient Algorithm Based on Immersion Simulations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 13, 583-598. https://doi.org/10.1109/34.87344
[47]	Walther, B. A., & Moore, J. L. (2005). The Concepts of Bias, Precision and Accuracy, and Their Use in Testing the Performance of Species Richness Estimators, with a Literature Review of Estimator Performance. Ecography, 28, 815-829. https://doi.org/10.1111/j.2005.0906-7590.04112.x
[48]	Wiggins, H. L. (2017). The Influence of Tree Height on Lidar’s Ability to Accurately Characterize Forest Structure and Spatial Pattern across Reference Landscapes. Master of Science, Missoula, MT: University of Montana. http://scholarworks.umt.edu/etd/11025
[49]	Wing, B. M., Ritchie, M. W., Boston, K., Cohen, W. B., Gitelman, A., & Olsen, M. J. (2012). Prediction of Understory Vegetation Cover with Airborne Lidar in an Interior Ponderosa Pine Forest. Remote Sensing of Environment, 124, 730-741. https://doi.org/10.1016/j.rse.2012.06.024
[50]	Yu, X., Hyyppä, J., Holopainen, M., & Vastaranta, M. (2010). Comparison of Area-Based and Individual Tree-Based Methods for Predicting Plot-Level Forest Attributes. Remote Sensing, 2, 1481. https://doi.org/10.3390/rs2061481
[51]	Yu, X., Hyyppä, J., Vastaranta, M., Holopainen, M., & Viitala, R. (2011). Predicting Individual Tree Attributes from Airborne Laser Point Clouds Based on the Random Forests Technique. ISPRS Journal of Photogrammetry and Remote Sensing, 66, 28-37. https://doi.org/10.1016/j.isprsjprs.2010.08.003