Estimating Stem Volume in Eucalyptus Plantations Using Airborne LiDAR: A Comparison of Area- and Individual Tree-Based Approaches

: Forest plantations are globally important for the economy and are signiﬁcant for carbon sequestration. Properly managing plantations requires accurate information about stand timber stocks. In this study, we used the area (ABA) and individual tree (ITD) based approaches for estimating stem volume in fast-growing Eucalyptus spp forest plantations. Herein, we propose a new method to improve individual tree detection (ITD) in dense canopy homogeneous forests and assess the e ﬀ ects of stand age, slope and scan angle on ITD accuracy. Field and Light Detection and Ranging (LiDAR) data were collected in Eucalyptus urophylla x Eucalyptus grandis even-aged forest stands located in the mountainous region of the Rio Doce Valley, southeastern Brazil. We tested ﬁve methods to estimate volume from LiDAR-derived metrics using ABA: Artiﬁcial Neural Network (ANN), Random Forest (RF), Support Vector Machine (SVM), and linear and Gompertz models. LiDAR-derived canopy metrics were selected using the Recursive Feature Elimination algorithm and Spearman’s correlation, for nonparametric and parametric methods, respectively. For the ITD, we tested three ITD methods: two local maxima ﬁlters and the watershed method. All methods were tested adding our proposed procedure of Tree Bu ﬀ er Exclusion (TBE), resulting in 35 possibilities for treetop detection. Stem volume for this approach was estimated using the Schumacher and Hall model. Estimated volumes in both ABA and ITD approaches were compared to the ﬁeld observed values using the F-test. Overall, the ABA with ANN was found to be better for stand volume estimation ( r y ˆ y = 0.95 and RMSE = 14.4%). Although the ITD results showed similar precision ( r y ˆ y = 0.94 and RMSE = 16.4%) to the ABA, the results underestimated stem volume in younger stands and in gently sloping terrain ( < 25%). Stem volume maps also ﬀ ered between the approaches; ITD represented the stand variability the stands and methods.

Abstract: Forest plantations are globally important for the economy and are significant for carbon sequestration. Properly managing plantations requires accurate information about stand timber stocks. In this study, we used the area (ABA) and individual tree (ITD) based approaches for estimating stem volume in fast-growing Eucalyptus spp forest plantations. Herein, we propose a new method to improve individual tree detection (ITD) in dense canopy homogeneous forests and assess the effects of stand age, slope and scan angle on ITD accuracy. Field and Light Detection and Ranging (LiDAR) data were collected in Eucalyptus urophylla x Eucalyptus grandis even-aged forest stands located in the mountainous region of the Rio Doce Valley, southeastern Brazil. We tested five methods to estimate volume from LiDAR-derived metrics using ABA: Artificial Neural Network (ANN), Random Forest (RF), Support Vector Machine (SVM), and linear and Gompertz models. LiDAR-derived canopy metrics were selected using the Recursive Feature Elimination algorithm and Spearman's correlation, for nonparametric and parametric methods, respectively. For the ITD, we tested three ITD methods: two local maxima filters and the watershed method. All methods were tested adding our proposed procedure of Tree Buffer Exclusion (TBE), resulting in 35 possibilities for treetop detection. Stem volume for this approach was estimated using the Schumacher and Hall model. Estimated volumes in both ABA and ITD approaches were compared to the field observed values using the F-test. Overall, the ABA with ANN was found to be better for stand volume estimation (r yŷ = 0.95 and RMSE = 14.4%). Although the ITD results showed similar precision (r yŷ = 0.94 and RMSE = 16.4%) to the ABA, the results underestimated stem volume in younger stands and in gently sloping terrain (<25%). Stem volume maps also differed between the approaches; ITD represented the stand variability better. In addition, we discuss the importance of LiDAR metrics as input variables for stem volume estimation methods and the possible issues related to the ABA and ITD performance.

Study Site
The study site is located in the Rio Doce Valley, Minas Gerais, southeastern Brazil ( Figure 1). This State is the first in terms of area of Eucalyptus plantations in Brazil [15]. The main characteristic of the region is the high variation in slope values and frequency of steep terrains. The total study area is 471.48 ha divided into 19 Eucalyptus urophylla x Eucalyptus grandis clone stands planted for pulpwood production. We used forest stands covering three age classes: young (<4 years), intermediate (between 4 and 6 years), and mature (>6 years) (see Table 1).

Field Measurements
A total of 33 rectangular field plots of~280 m 2 were measured in these stands as part of previous forest inventory. Tree dbh was measured with tree calipers and heights with Haglof digital clinometer (Haglof Inc., Madison, MS, USA). The outside-bark volumes (m 3 ) were calculated using the Schumacher and Hall [45] equation fitted using historical tree scaling data in the area of study (see Supplementary Material). The field measurements were used as observed values to compare to the LiDAR volume estimates (Table 1). In order to make this comparison, one geographic position from each rectangular plot was collected with a dual-frequency GNSS receiver TOPCON HiPer II (Topcon Corporation, Tokyo, JPN) with estimated uncertainty <2 m in stand-alone mode. We created a 9.44 m-buffer around the geographic point to represent the total area of each plot. The single point in each plot was the only positional reference the stand owners had for the plots, and previous stand harvest hampered us in obtaining more geographic positions for this study. This procedure, although not ideal, was usable because of the nature of our target. Even-aged clone trees have considerably low structural variability, especially at this spatial scale (Table S2), and plantation spacing is standardized within stands.

Lidar Data Acquisition and Processing
LiDAR data was acquired in an airborne mission with the aircraft Cessna model 206, which had an average height above ground of 618 m, speed of 55 m/s and pulse frequency of 300 kHz, yielding a point cloud with average 5 points per square meter. The laser beam divergence resulted in a footprint of 0.31 m on the ground. The maximum scan angle of the sensor was 30 • covering an average of 713 m per swath. The survey was accomplished in February 2014 near to the dates when field data was collected (June-August 2014).
Prior to the application of the two approaches, we derived the Digital Terrain Model (DTM), normalized point cloud, and Canopy Height Model (CHM). We classified the point cloud into ground and non-ground returns with the adapted Kraus and Pfeifer [45] morphological filter, in the software FUSION [46]. This filter works by classifying points as ground (1) and non-ground (0) in an iterative process following a weight function. Detailed information on the algorithm operation is available in [47]. We tested combinations of parameters g (−4, −3, −2, −1, 0 and 1) and w (1, 2, 3 and 4) from the morphological filter equation and held the best result (smoother DTM) with the values −3 and 3, respectively. Parameters a and b were kept constant as 1 and 4 since these values have shown adequate performance in other applications [46,47]. We created the DTM by filtering the terrain points and rasterizing it to a 0.45 × 0.45 m cell size [48]. The average elevation of points within each cell is defined as the cell value. We then normalized the whole point cloud (all ground elevations equal to 0) and created the CHM. The CHM represents the canopy heights and we used it to detect trees in the ITD approach, while the normalized point cloud was used to calculate LiDAR metrics for the ABA.

Modeling Forest Stand Volume
Under ABA, we tested two parametric and three nonparametric methods: Linear regression, Gompertz model [49,50], Artificial Neural network (ANN) [51,52], Random Forest (RF) [53], and Support Vector Machine (SVM) [54,55]. Feature selection, tuning and accuracy assessment of the models described in Section 2.5.2. In the ITD, we first assessed tree detection accuracy, proposing a Remote Sens. 2020, 12, 1513 5 of 25 new method to improve it, and estimated stand volume by summing up tree volumes calculated with the Schumacher and Hall [45] model (Equation (1)).
where: V = Tree volume, dbh = tree diameter at breast height (1.3 m above ground) and h = tree height.

LiDAR Metrics and Feature Selection
From all the elevation metrics available within the software FUSION [46], we defined as candidate features those related to forest structure and that had been used in other studies of forest applications of ALS [20,23,56] (Table 2). The lower threshold for the metrics was 1.3 m (i.e., only points above this height were considered in the metric calculation). We selected the LiDAR metrics for the nonparametric methods considering two aspects: collinearity and importance. First, we excluded metrics that were highly correlated using a ± 0.9 correlation coefficient threshold [20,22]. Subsequently, we used a method based on Recursive Feature Elimination (RFE) [57], considering only the features not excluded in the first step to account for feature importance in the models. This method is a backward selection algorithm that computes feature importance at each iteration by ranking them from the most important to the least, removing a user-defined subset iteratively at each stage [58]. Although feature collinearity may not severely affect nonparametric methods, the exclusion of highly correlated ones was important for making RFE iterations more constant, as features could be interchangeable within the models. We ran the algorithm 100 times with each tested method and selected the most frequently held features considering their position in the importance rank defined by the algorithm. The RFE algorithm is present at the "Caret" package [59] and was run through the "recursive_feature_elimination" function of the Labgeo package [60], both in the R environment [61]. For the parametric methods, we computed the Spearman's correlation coefficient to rank the relationships between the candidate LiDAR metrics and volume. We selected the metric with the highest correlation coefficient with respect to volume as the explanatory variable in the parametric models.

Estimation Methods Tuning and Fitting
We trained the nonparametric methods (i.e., ANN, RF, and SVM) varying their tuning parameters selecting the best tuning based on correlation coefficient (r yŷ ; Equation (2)) and absolute and relative Root Mean Square Error (RMSE and RMSE%; Equations (3) and (4), respectively). The ANN is a multi-layer perceptron with backpropagation learning function and logistic activation function for the hidden layers. The tuning parameter is the number of units in the hidden layer (hereafter referred to as "size"). For the RF, we set the number of decision trees to 500 and varied the number of candidate features randomly sampled at each tree-node (hereafter referred to as "mtry"). The SVM was with the Radial Basis Function Kernel tuned varying the parameters "Cost" (C) and sigma. All algorithms were fitted in R [61] using the package Labgeo [60], with dependencies of the packages RSNNS [62], random Forest [63] and kernlab [64] for ANN, RF and SVM, respectively. Before model fitting, we randomly split the dataset into training and test data with 60% and 40% of the data, respectively. The training set was used to train and define the best non-parametric methods using leave-one-out cross validation and to fit the parametric models. The test set was held out from this process and used to apply and assess the methods' performance (external validation; see Section 2.7).
whereŶ i = estimated stem volume; Y i = observed stem volume; Y = mean observed stem volume; Y m = estimated mean stem volume and n = number of observations.

Tree Detection and Crown Delineation
In this study, we tested three methods to detect the treetops: two Local Maxima Filter (LMF) algorithms [27,29,31], varying window size; and the watershed method [65] for crown delineation (WTS), using the ArcGIS hydrology tools (ESRI-Redlands, CA, USA) [66]. The latter comes from watershed creation using the inverted CHM, defining its boundaries. This resembles classic watershed algorithms [28]. Treetops are the maximum value within each delineated crown in this method. To improve the tree detection performance, we proposed a new algorithm we call Tree Buffer Exclusion (TBE), to be applied after the tree detection methods.
To understand how it works, consider N as the set of all detected trees in a stand. This exclusion consists of iterations through N, selecting each individual (N i ) and creating a buffer of radius r. If the analyzed tree (N i ) is the highest among those within the buffered zone, it is considered as a real treetop, if not, it is eliminated ( Figure 2). Therefore, it allows the removal of points too close to each other, which likely belong to the same crown. This is based on the a priori consideration that trees on even-aged planted stands, like those in this study, have similar canopy structures. Such homogeneity, together with the high spatial resolution of the CHM might induce the algorithms to count more than one "treetop" per tree.
Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 25 watershed creation using the inverted CHM, defining its boundaries. This resembles classic watershed algorithms [28]. Treetops are the maximum value within each delineated crown in this method. To improve the tree detection performance, we proposed a new algorithm we call Tree Buffer Exclusion (TBE), to be applied after the tree detection methods.
To understand how it works, consider N as the set of all detected trees in a stand. This exclusion consists of iterations through N, selecting each individual ( ) and creating a buffer of radius r. If the analyzed tree ( ) is the highest among those within the buffered zone, it is considered as a real treetop, if not, it is eliminated ( Figure 2). Therefore, it allows the removal of points too close to each other, which likely belong to the same crown. This is based on the a priori consideration that trees on even-aged planted stands, like those in this study, have similar canopy structures. Such homogeneity, together with the high spatial resolution of the CHM might induce the algorithms to count more than one "treetop" per tree. The combination of method, window size, and TBE exclusion radius yielded 35 possibilities for treetop detection (Table 3). We defined the best method validating tree count to the estimated stem density in the stand, which was calculated using the sample plots data. We found that to be more reliable than plot validation since we did not have geolocation for trees in the plots, geolocation accuracy was low for this purpose and a slight change in plot position might severely affect the presence/absence of a tree in a plot. Table 3. Summary of parameters used in each tree detection and delineation algorithm. NA = not applicable. TBE stands for our proposed method of Tree Buffer Exclusion.

Method
Window We assessed tree detection accuracy by age, slope and scan angle classes. Age classes were the same as defined in Section 2.1. Slope classes were: gentle (<25%), moderate (25-45%) and steep slope (>45%). A relatively higher threshold (45%) is used to define mountainous slopes in Brazil [67]. The classes for scan angles were: nadir (0-7 • ), off-nadir (7-23 • ) and large off-nadir (>23 • ) [41]. Scan angles are likely to overlap, mainly between flight tracks. Thus, it was not possible to obtain a reasonable number of plots with a single scan angle class. To assess its effects, we calculated the scan angle class mode using values retrieved for each returned point in a plot and defined it as the single scan-angle class for that plot. Trees in the border of plots were kept for the analysis only if the treetop was within the plot.

Individual Tree Dbh and Volume Estimates
We tested four methods to estimate dbh (see Supplementary Files Figure S1) and found ANN as the best unbiased method to estimate dbh for trees in the stand (r yŷ = 0.94, RMSE% = 5.52, bias% = 0.0001; Figure S1). The input variables were height, soil type, stem spacing and regeneration method. With the dbh and LiDAR retrieved height for all trees, we estimated volume with the Schumacher and Hall [55] equation. The parameters for this equation were provided by the stand owners and were obtained with historical tree scaling data in the area of study (Table S4). The performance assessment for this approach was conducted comparing the field observed values to the sum of the detected tree volumes in the plot and stand.

Model Accuracy, Statistic Tests and Comparison of the Approaches
Assessment of all models' performance was carried out computing the correlation coefficient (r yŷ ) absolute and relative RMSE, and absolute and relative bias (Equation (2) to Equation (6)). We also tested all estimates with the statistic F modified by Graybill [68]. This statistic calculates the significance of a simple regression testing the null hypothesis that the parameter b 0 (intercept) is equal to zero and b 1 (slope) to one (H 0 : The comparisons were for the pairwise combination of estimated (Ŷ i ) and observed (Y i ) values [69].
whereŶ i = estimated stem volume; Y i = observed stem volume; Y = mean observed stem volume and n = number of observations.

Volume Maps
Volume maps (i.e., wall-to-wall estimates) were created by sectioning the stand into tiles and then estimating the stem volume for each one. The area of each tile was equal to those of the inventory plots (280 m 2 ). In the ABA, we computed the selected LiDAR metrics for each tile and then applied the estimation methods. For the ITD, we summed up all estimated tree volumes within a tile. Hence, the volume for the whole stand was calculated and can be assessed (Figure 3).
Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 25 where = estimated stem volume; = observed stem volume; = mean observed stem volume and n = number of observations.

Volume Maps
Volume maps (i.e., wall-to-wall estimates) were created by sectioning the stand into tiles and then estimating the stem volume for each one. The area of each tile was equal to those of the inventory plots (280 m 2 ). In the ABA, we computed the selected LiDAR metrics for each tile and then applied the estimation methods. For the ITD, we summed up all estimated tree volumes within a tile. Hence, the volume for the whole stand was calculated and can be assessed (Figure 3).

Area-Based Approach (ABA)
In the first feature selection step for the nonparametric methods (ANN, RF, SVM), we removed 15 from the 28 candidate features, based on a correlation threshold of 0.90 ( Figure 4). Feature importance ranks were then built using the RFE algorithm with the remaining 13 features.
In the first feature selection step for the nonparametric methods (ANN, RF, SVM), we removed 15 from the 28 candidate features, based on a correlation threshold of 0.90 ( Figure 4). Feature importance ranks were then built using the RFE algorithm with the remaining 13 features. The most frequently ranked features did not significantly change among estimation methods, especially until the fourth most important metric (Table S4). The four most important LiDAR metrics to estimate volume were: H.P95, Hmod, H.P50, and H.P10. Therefore, these metrics were selected as input features for the nonparametric methods. From the fifth position forward, the ranked features changed over the iterations, indicating that the importance of such features for model fitting is not as clear as the first four ranked metrics.
For the parametric methods (linear regression and Gompertz), the variable H.P95 (height at the 95th percentile) was selected based on Pearson's correlation coefficients. Its correlation with volume was the highest (ρ = 0.8361) amongst LiDAR metrics. The estimation methods were then fitted using the selected metrics (Table 4).

Estimation Method
Tuning parameter / equation The most frequently ranked features did not significantly change among estimation methods, especially until the fourth most important metric (Table S4). The four most important LiDAR metrics to estimate volume were: H.P95, Hmod, H.P50, and H.P10. Therefore, these metrics were selected as input features for the nonparametric methods. From the fifth position forward, the ranked features changed over the iterations, indicating that the importance of such features for model fitting is not as clear as the first four ranked metrics.
For the parametric methods (linear regression and Gompertz), the variable H.P95 (height at the 95th percentile) was selected based on Pearson's correlation coefficients. Its correlation with volume was the highest ( = 0.8361) amongst LiDAR metrics. The estimation methods were then fitted using the selected metrics (Table 4). Volume estimates in the parametric and nonparametric methods showed overall Pearson's correlation coefficient close to 0.80 in the test data (i.e., 40% of the data held out from model fitting) ( Table 5). The values ranged from 0.79, with linear regression to 0.83 with ANN. The ANN also showed the lowest values for RMSE% (19.39) while the SVM showed the highest (24.71). All estimates with parametric methods were different (p < 0.05) from observed values, despite the lowest overall bias with the Linear model.

Individual Tree Approach (ITD)
The best method to detect trees varied over the stands (Table 6). However, all tested methods had their performance improved by adding the TBE. The best overall results for the number of detected trees were observed by using LMF1 with a 2 m window and TBE radius equal to 1.5 m. This method had the lowest tree detection errors for 11/19 stands, followed by the WTS + 1.5 m buffer (7/19). Upon comparing LiDAR-and field-measured heights by age classes in the plots, the only difference (p < 0.05) found was on the young stands' plots (Figure 5a), where the trend was for underestimating tree mean height (bias% = −11.17). Slope classes also affected tree mean height for the LiDAR measurements. In gentle slopes (<25%) the differences from LiDAR to field-measured tree heights were significant showing an underestimation pattern (Figure 5d). In scan-angle classes, the LiDAR retrieved mean height was different (p < 0.05) from the field measures at nadir angles (0-7 • scan angle class) (Figure 5g). Although the differences were not significant on large off-nadir angles (23-30 • ), the number of observations in this class was too low to make inferences.

−7.39
Upon comparing LiDAR-and field-measured heights by age classes in the plots, the only difference (p<0.05) found was on the young stands' plots (Figure 5a), where the trend was for underestimating tree mean height (bias% = −11.17). Slope classes also affected tree mean height for the LiDAR measurements. In gentle slopes (<25%) the differences from LiDAR to field-measured tree heights were significant showing an underestimation pattern (Figure 5d). In scan-angle classes, the LiDAR retrieved mean height was different (p<0.05) from the field measures at nadir angles (0°-7° scan angle class) (Figure 5g). Although the differences were not significant on large off-nadir angles (23°-30°), the number of observations in this class was too low to make inferences.

Stand Volume Estimates with ABA and ITD
The stand volume estimations r yŷ and RMSE (%) were overall around 0.9% and 15%, respectively ( Figure 6). The best result was obtained with the ABA with ANN where r yŷ , RMSE (%) and bias (%) were equal to 0.95, 14.41 and −1.59, respectively. The estimates with the ITD approach, despite the relatively high r yŷ and low RMSE (%) (0.94 and 16.4, respectively), presented a bias (%) of −9.55 which was the highest absolute value. This tendency was observed mainly in lower volumes, underestimating the observed values (Figure 6f).

Stand Volume Estimates with ABA and ITD
The stand volume estimations and RMSE (%) were overall around 0.9% and 15%, respectively ( Figure 6). The best result was obtained with the ABA with ANN where , RMSE (%) and bias (%) were equal to 0.95, 14.41 and −1.59, respectively. The estimates with the ITD approach, despite the relatively high and low RMSE (%) (0.94 and 16.4, respectively), presented a bias (%) of −9.55 which was the highest absolute value. This tendency was observed mainly in lower volumes, underestimating the observed values (Figure 6f).

Stem Volume Mapping
The volume maps showed differences in ABA and ITD estimates (Figures 7-9). The main differences were in the volume distribution across the stand. Although both seemed to have normal distribution ITD stand volume distribution presented a wider range than ABA's.
In addition, ABA maps using ANN had the mean volume estimates closer to the ITD maps. In intermediate stands (Figure 8), the ANN had two peaks where estimates were more frequent, showing that stem volume estimates may have different behavior even within one age class.
Stem volume estimates using Gompertz differed from estimates with the ITD mainly in young stands (Figure 7). Values were more frequent in the classes around 140-160 m 3 ha −1 using Gompertz, while they were around 110-140 m 3 with ITD. This change made the distribution with Gompertz slightly skewed in relation to the estimated mean stem volume from the ITD approach. This pattern was attenuated in intermediate and mature stands (Figure 8, Figure 9).

Stem Volume Mapping
The volume maps showed differences in ABA and ITD estimates (Figure 7, 8 and 9). The main differences were in the volume distribution across the stand. Although both seemed to have normal distribution ITD stand volume distribution presented a wider range than ABA's.
In addition, ABA maps using ANN had the mean volume estimates closer to the ITD maps. In intermediate stands (Figure 8), the ANN had two peaks where estimates were more frequent, showing that stem volume estimates may have different behavior even within one age class.
Stem volume estimates using Gompertz differed from estimates with the ITD mainly in young stands (Figure 7). Values were more frequent in the classes around 140-160 m 3 ha -1 using Gompertz, while they were around 110-140 m3 with ITD. This change made the distribution with Gompertz slightly skewed in relation to the estimated mean stem volume from the ITD approach. This pattern was attenuated in intermediate and mature stands (Figure 8, Figure 9).

Forest Stands Strata Characterization from LiDAR Metrics
The selected metrics used as input features in the estimating methods were H.P95, Hmode, HP50 and HP10. These metrics appear frequently in other studies that estimate stand volume using LiDAR data, and typically, at least one height percentile is used [9,11,17]. The selected metrics support the concept of using metrics that are related to the canopy height and its variance for volume estimation [71]. They represent specific heights in the point cloud and the same metric can be in a different forest stratum depending on the forest structure ( Figure 10). The highest height percentile used was H.P95, which corresponds roughly to the top of the canopy. However, using this metric alone might not account for the whole vertical variability in the stand, even in homogeneous forests as Eucalyptus even-aged stands. Thus, vertical changes in forest structure were detected when other height metrics were inserted into the models (Figure 10).
In this study, the most important metrics after the H.P95 were Hmode, HP50 and HP10. Figure  10 shows a representation of how these metrics account for the stand variability. If the canopy is uniform in a single layer, lower percentiles (e.g., H.P10) give us values near the highest percentiles (e.g., Figure 10a1 and Figure 10d1 where H.P10 = 24.23 m and H.P95 = 27.42 m, difference = 3.19 m). Therefore, in the case of a homogeneous forest, using more than one height metric will add meaning

Forest Stands Strata Characterization from LiDAR Metrics
The selected metrics used as input features in the estimating methods were H.P95, Hmode, HP50 and HP10. These metrics appear frequently in other studies that estimate stand volume using LiDAR data, and typically, at least one height percentile is used [9,11,17]. The selected metrics support the concept of using metrics that are related to the canopy height and its variance for volume estimation [71]. They represent specific heights in the point cloud and the same metric can be in a different forest stratum depending on the forest structure ( Figure 10). The highest height percentile used was H.P95, which corresponds roughly to the top of the canopy. However, using this metric alone might not account for the whole vertical variability in the stand, even in homogeneous forests as Eucalyptus even-aged stands. Thus, vertical changes in forest structure were detected when other height metrics were inserted into the models (Figure 10).

Individual Tree Detection and TREE Buffer Exclusion
Tree detection is the most sensitive part of the ITD approach [75][76][77]. The most common methods to detect trees using LiDAR data are raster-based (e.g., LMF and watershed algorithms) [32] and in homogenous forests, such as even-aged Eucalyptus stands, a simple local maxima filter is often selected to detect trees [35,78]. In this study, we used three raster-based methods and all of them benefited by adding the TBE. The best overall method was the 2-m window of [27] LMF with 1.5-m- In this study, the most important metrics after the H.P95 were Hmode, HP50 and HP10. Figure 10 shows a representation of how these metrics account for the stand variability. If the canopy is uniform in a single layer, lower percentiles (e.g., H.P10) give us values near the highest percentiles (e.g., Figures 10a 1 and 10d 1 where H.P10 = 24.23 m and H.P95 = 27.42 m, difference = 3.19 m). Therefore, in the case of a homogeneous forest, using more than one height metric will add meaning to the model if the crown height is a parameter of interest. Meanwhile, when there are more returns (i.e., intercepted objects) between the upper stratum and the ground, differences from lower and higher percentiles will be larger (Figures 10a 3 and 10d 3 where H.P10 = 9.38 m and H.P95 = 33.31 m, difference = 23.93 m). Points below the upper canopy are likely from suppressed trees or limbs not self-pruned during tree growth. Attention should be given to the fact that stands reaching the end of the rotation are frequently without recent silvicultural treatments and may contain undergrowth (e.g., herbs and shrubs) [20]. In addition, lower percentiles (below H.P50) may be unstable if the LiDAR data acquisition is taken at significantly high altitudes due to attenuation of the pulse signal [23].
Removing highly correlated metrics prior to feature selection is commonly applied when estimating forest parameters with LiDAR metrics, performed for example with PCA or the Pearson's coefficient of correlation [20,23,26,72]. Some regression models are more sensitive than others, but even for nonparametric methods, this step is important to remove redundancy and lead to more interpretability of the selected features [72]. The RFE has been applied in other fields of study to select the best subset of features to input in a model, mainly using a large number of candidate features [57,73]. This method benefits from refitting fewer models at each step due to the variable rank creation [58]. This approach directly outputs the most important features for the estimation method and is an option to make decisions when a large number of candidate features are available, which is likely to happen in studies with the integration of LiDAR data from other platforms (e.g., terrestrial and orbital) or with other sensors (e.g., optical and radar). We applied the RFE algorithm oriented by each tested nonparametric method (ANN, RF, and SVM). Although they led to similar results, computation labor was remarkably lower with RF. This algorithm has a well-known internal method for computing variable importance [74].

Individual Tree Detection and TREE Buffer Exclusion
Tree detection is the most sensitive part of the ITD approach [75][76][77]. The most common methods to detect trees using LiDAR data are raster-based (e.g., LMF and watershed algorithms) [32] and in homogenous forests, such as even-aged Eucalyptus stands, a simple local maxima filter is often selected to detect trees [35,78]. In this study, we used three raster-based methods and all of them benefited by adding the TBE. The best overall method was the 2-m window of [27] LMF with 1.5-m-buffer exclusion. This LMF, implemented in the software FUSION [46], has an optional enhancement to vary the window size based on an equation relating field measured height and crown diameter which we did not use due to the lack of required information. Thus, defining the best window size is done empirically by reasonably relating it to tree spacing. Wider windows than tree spacing cause omission errors (i.e., existing trees that are not detected), while narrower windows cause commission errors (i.e., detecting more than one maximum per tree). In our study, the commission errors were reduced using the TBE, based on a buffer around each treetop. Errors in tree detection here were likely caused by the dense canopy and round-shaped crown that makes treetops harder to define than in conifers, for example [38]. The watershed method also performed relatively well to detect trees and might be a good option in Eucalyptus stands. One should notice that the LMF and watershed method will not have perfect performance unless all trees are in the upper canopy (i.e., no presence of understory trees). Further, for LMF, the direction and start of the moving window is a factor not completely explored that can affect tree detection. The accurate information on stem density and tree height variation across the stand obtained by detecting the trees is important for precision forest and timber assortment planning [79,80].
Treetop detection and height retrieval are highly influenced by stand and sensor characteristics [3,39,40]. We observed significant differences between LiDAR and field measures when assessing the results in young stands (<3 years) and gentle slopes (<25%), underestimating tree heights. In the early stages of Eucalyptus growth, the competition still has not had its effects on the height of trees as much as in a mature forest. In addition, crown shape can affect LiDAR measurements of height [38,81]. We highlight that correctly identifying trees in young stands is important for growth and yield modeling nonetheless, for pre-harvest volume stock estimation which uses mainly mature forest data, this bias should not be an issue.
The differences were also significant in gentle slopes (<25%). At first, we expected that steep slopes would bias tree detection and height, which we did not observe. In gentle terrains, a dense upper canopy might act as a barrier to pulse penetration and return to the sensor, likely generating errors in the DTM. Height retrieval from ALS and the forest parameter estimates are sensitive to these errors [82].
Our results show that younger stands were more affected for tree height retrieval with LiDAR. Increasing point density should be considered for example by varying acquisition parameters such as pulse frequency and flight height. However, it implies a higher cost for LiDAR data acquisition and processing. An alternative should be to acquire a first (more expensive) dense point cloud, in order to model the terrain more accurately, and successive low-density point clouds for long-term volume estimation [44].

Area-and Individual Tree-Based Approach for Estimating Stem Volume
The ABA with ANN was the best overall method. ANNs are currently used in forestry in many applications, for instance, developing growth and yield models [83], predicting height [84] and stem volume [85]. As a nonparametric method, it benefits from its capacity to account for data variability and non-linear relationships. Alternatively, parametric models are simpler and of more widespread knowledge, besides being easier to share and explain. The issues arise when problem complexity increases, for instance by increasing the number of predictors, and in the presence of considerable variability and outliers.
The error associated with the stem volume estimation methods was around 15%. Although not as low as found in other studies in the same type of forest (e.g., [20]) it is reasonable given the reduced number of plots and the variability of age and regeneration methods in the stands in our study. Clearcut and coppice stands might have structural differences and can in practice be separated during modeling, for example. The recommended number of plots to fit estimation models with ALS data varies depending on the characteristics of the study area and, as in conventional inventory, the desired precision [86]. Maltamo et al. [87] demonstrated that plot selection to train stem volume models can be improved using ALS data as a priori information. Laranja et al.
[88] using a double sampling strategy with LiDAR, decreased the number of necessary plots from 35 to 10 while reducing the inventory error. Our study was carried out with few sample plots (n = 33), which can be sometimes referred to as an issue for training non-parametric methods. Nevertheless, this number is consistent with what is used to record homogeneous Eucalyptus plantation variability (one plot every 10-15 ha; [20,24]). Furthermore, the number of features (explanatory variables) was relatively low (four) and we were concerned about the tuning parameter levels to avoid overfitting.
The ITD approach had lower values for r yŷ (0.80), and higher for RMSE (23.96%) in the plot validation. ABA benefits in this type of validation for retrieving a mean value per area. Whereas, in the ITD, each miss-detected tree can increase the stem volume estimate error per plot. In our study, we used a buffer from a single georeferenced point in the sample plots to validate estimates and a slight change in its position affects the validation. When we validate the volume for the whole stand (i.e., field measurements from sample plots interpolated to the stand), ITD yields better precision and accuracy (r yŷ = 0.94 and RMSE = 16.4%), though biased in lower values (Figure 6f). The trend to underestimate lower volumes comes from the errors in tree height measurements in younger stands discussed in the previous section. Breidenbach et al.
[89] introduced a method for reducing bias from the ITD by calculating the volume by segments, integrating ABA and ITD. Some other studies have used calibration with linear regression to reduce biased predictions in ITD [76,77]. These are some of the ways to obtain unbiased estimates in the ITD approach.
In the volume maps we observe that the wall-to-wall estimates vary differently for both approaches. The main difference is the distribution of estimated volumes. The range is wider with ITD since ABA concentrate estimates only in a few volume classes. This is an indicator of how ITD can be more thorough in estimates since it is sensitive, for instance, to bigger and smaller trees, which are present in the stand. Extreme values (e.g., 1 or 640 m 3 ha −1 ) appear because one cell can concentrate more trees than others. If they appear too often it could represent estimation errors.
From an operational standpoint, ABA is related to lower costs because it is possible to create models even with low-density point clouds [16]. The drawbacks of this approach are that it is entirely dependent on field plot measurements, size and accuracy [75,76]. Because of that, reductions in fieldwork are limited. On the other hand, for the ITD, high-density point clouds are required to yield better results, implying higher costs [78]. This approach, however, gives a direct measurement of stem density and height distribution to all trees (nearly 520,000 in this study, Figure S2) in addition to dbh and volume estimates per tree [75,76]. Furthermore, fieldwork can be reduced to calibrate the models. Tavares Júnior et al. [85], for instance, found that tree volume could be predicted with four sample trees per diameter class using ANN.
The benefits of using LiDAR data to estimate volume is that we have precise spatially explicit estimates for the whole stand. Inventory plots have been used historically and are usually reliable but yield on a single mean value for the whole stand. They do not account for local growth differences affected by, for example, local site variations, diseases or pests spread. An enhanced forest inventory [86] for forest parameters estimation and monitoring can be reached by using the ABA and ITD approaches. Especially with the ITD, the results are promising for even-aged Eucalyptus stands. These forests are uniform, having most trees in the upper canopy favoring simple raster-based methods for tree detection (e.g., LMF and watershed algorithms), with minimized commission errors using algorithms such as the TBE presented here.

Conclusions
This study presented a comparison of area-(ABA) and individual tree-based (ITD) approaches for estimating stem volume in Eucalyptus forest plantations. The ABA was more accurate and precise especially with ANN for estimating stem volume as compared to the ITD approach. However, with ITD it was possible to obtain direct measures of tree density and height distribution as well as capture more variability when mapping volume, especially for homogeneous Eucalyptus stands. The Tree Buffer Exclusion (TBE), presented and tested here as our secondary objective, improved all raster-based tree detection algorithms used in this study. It contributes by reducing treetop detection commission errors, commonly present in studies associated with even-aged stands containing trees exhibiting irregular crown shapes. We expect this method to contribute to ITD implementation in Eucalyptus forest stands in the future.