Introduction
Current methods of forest inventory are based on the direct surveying of trees and sampling of ground plots (^{Campos and Leite, 2013}). Statistical models generated to predict average estimates of stand attributes and their spatial variability are assessed only sporadically (^{Zhou et al., 2013}). Data interpolation techniques can be used to assess stand spatial structure, yet they are limited by the sampling intensity of ground plots, which is usually insufficient to yield precise estimates (^{Bouvier et al., 2015}; ^{Viana et al., 2012}). On the other hand, integrating remotely sensed data is an actual possibility, as it can be acquired to a great extent, within a short time and at reasonable prices (^{Hummel et al., 2011}). Airborne laser scanning (ALS) technology is a potential remote sensing tool, as it can provide data on multiple strata of the canopy, while passive sensors cannot (^{Magnussen and Boudewyn, 1998}; ^{Næsset, 1997}). The laser sensor records tridimensional coordinates over the surveyed area by collecting the time in which a laser pulse takes to go back to the aircraft after being reflected by the target. The result is a georeferenced 3D point cloud with high spatial resolution (^{Baltsavias, 1999}; ^{Reutebuch et al., 2005}).
One of its promising applications has consisted of characterizing the vertical structure of forests with apparent canopy height profiles (CHPs) (^{Coops et al., 2007}; ^{Dean et al., 2009}; ^{Lefsky et al., 1999}; ^{Lovell et al., 2003}; ^{Maltamo et al., 2004}). The apparent CHP is a vertical distribution of laser points, which can be represented by their empirical distribution or by probabilistic functions (^{Harding et al., 2001}; ^{Magnussen et al., 1999}). A CHP represented by a probability function is a convenient way of summarizing and retrieving the canopy vertical form (^{Coops et al., 2007}). The Weibull probability density function (pdf) has been used to model vertical profiles thanks to its flexibility in characterizing different types of vegetation, and its parameters of scale and shape were successfully correlated with above ground attributes, such as height of trees, density of stems, and diameter at breast height (^{Coops et al., 2007}; ^{Mori and Hagihara, 1991}). In Brazil, ALS technology has been used for terrain modeling, but applications involving assessment and mapping of biomass in eucalypt stands are still incipient (^{Packalén et al., 2011}; ^{Silva et al., 2014}; ^{Vauhkonen et al., 2011}). This study aimed to map the stem biomass stock based on canopy height profiles statistics using ALS data of an even-aged eucalypt plantation, in southeastern Brazil, and compare the results with alternative maps generated by ordinary kriging interpolation from field-derived measurements.
Materials and Methods
Study site
The study site is located in the state of São Paulo, southeastern Brazil (22º58’04” S; 48º43’40” W) (Figure 1). The area is approximately 200 ha and is composed of a 6.5-year-old Eucalyptus grandis (W. Hill ex Maiden) plantation with seedlings coming from a fourth-generation seed orchard (^{Campoe et al., 2012}). The stands were planted in December 2002, with an approximate density of 1600 trees ha^{−1} (3.75 m × 1.6 m) following minimum site preparation (^{Gonçalves et al., 2013}). The Köppen climate type is Cfa; with mean annual temperature equal to 19.4 ºC and mean annual precipitation totaling 1300 mm (^{Alvares et al., 2013}). The site relief is mainly flat to soft wavy, and the altitude approximately 750 m (^{Campoe et al., 2012}).
Field measurements and plot summaries
Figure 1 shows the network of inventory ground plots in the study site. They were divided into 22 training plots and 21 validation plots. Plots 1 to 12 (457 m^{2} to 574 m^{2}) were established for research related to ecophysiology. Their location covered the site’s productivity range after a census of diameter at breast height was carried out (^{Campoe et al., 2012}). In addition, plots 13 to 22 (811 m^{2} to 881 m^{2}) were designed specifically for ALS investigations. Plots 23 to 43 (240 m^{2}) are permanent ground plots of continuous inventory owned by the forest producer.
Field measurements were conducted in the training plots in July 2009, and July 2008, in the validation plots (Table 1). Diameter at breast height (DBH, 1.3 m above ground level) and total tree height were measured on all plots. Stem dry biomass was estimated using the following local-specific allometric equation derived from destructive sampling in Sep 2008 (^{Campoe et al., 2012}) (Eq. 1).
Training plots (n = 22) |
Validation plots (n = 21) |
||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|
Stat. | Surv. | DBH | G | H | B | Stat. | Surv. | DBH | G | H | B |
- | % | cm | m^{2} ha^{−1} | m | Mg ha^{−1} | - | % | cm | m^{2} ha^{−1} | m | Mg ha^{−1} |
min. | 82 | 13.1 | 25 | 21.5 | 115 | min | 72 | 12.7 | 22 | 20.7 | 94 |
avg. | 90 | 14.9 | 30 | 24.6 | 155 | avg. | 88 | 14.2 | 27 | 22.6 | 125 |
max. | 97 | 16.6 | 34 | 27.5 | 195 | max | 98 | 16.4 | 35 | 24.9 | 167 |
st. dev. | 3.9 | 1 | 3.5 | 1.5 | 27 | st. dev. | 6.8 | 0.9 | 3.4 | 1 | 19 |
n = number of sample plots; Stat. = statistics; Surv. = survival of trees; DBH = diameter at breast height; G = basal area; H = total height; B = dry stem biomass; min. = minimum observed value among plots; avg. = average value among plots; max. = maximum observed value among plots; st. dev. = standard deviation of mean values.
where: = dry stem biomass per tree (kg per tree), D_{i} = diameter at breast height (cm), H_{i} = total height (m).
Airborne laser scanning (ALS) data
The ALS survey was undertaken during Apr 2009, a period of the year with maximum stable leaf area index (LAI) (^{le Maire et al., 2011}). The flight mission was conducted by a twin-engined light aircraft equipped with a discrete-return small footprint laser scanner (laser wavelength = 1064 nm). The parameters from the flight mission and the scanning process were as follows: flight height = 900 m; flight speed = 132 km h^{−1}; swath width = 235 m; swath overlap = 30 %; scanning angle = 15º; scanning frequency = 74 Hz; frequency of pulse emission = 110 kHz; footprint = 21 cm; point density = 6.5 pts m^{−2}; standard deviation of point density = 2 pts m^{−2}.
GPS observations from ground and aircraft were processed in a method to obtain a unique and adjusted cinematic solution to a well-known coordinate system, using the Waypoint GraphNav software. A digital elevation model (DEM) was generated with a one meter resolution, after ground points were labeled with the Multiscale Curvature Classification algorithm (MCC-LiDAR) (^{Evans and Hudak, 2007}). The DEM was subtracted from all point elevations to remove topographic variations (normalization process). Then, the normalized point cloud was clipped to the same locations as the inventory field plots. The Fusion LiDAR Toolkit software was used to normalize the point cloud (clipdata), to clip plots using shapefiles (polyclipdata) and to extract ALS metrics (cloudmetrics).
We tried to select metrics based on a priori knowledge in an attempt to improve the capacity of the model generalization (^{Bouvier et al., 2015}). Among the metrics extracted from each ALS plot, we focused on using height percentiles. A percentile x is the height z in which x % of points in the ALS cloud are beneath z. We tested the percentiles which predicted better basal area and volume in ^{Zonete et al. (2010)}: first return percentiles 10, 30, 70, 90. Moreover, we selected metrics such as mean height, standard deviation and variance of heights, all of them also being derived from first returns. These metrics were calculated disregarding laser points with a height less than 2 m (^{Næsset, 2002}; ^{Zonete et al., 2010}). We also parameterized a metric of density, which we named p_understory consisting of the proportion of laser points between 0.3 and 15 meters in height and the laser points between ground, (0 meters) and 15-meters in height. The 15-meter threshold was set based on visual inspection of the point cloud and on field experience. In this case, all return echoes were used to increase the sample of points in the understory layer.
Apparent canopy height profiles (CHPs)
From a set of possible theoretical distributions, we chose the Weibull pdf due to its flexibility in characterizing foliage distributions of different types of vegetation and its potential for fitting skewed data, which is a common feature of forest-derived ALS data (^{Coops et al., 2007}; ^{Dean et al., 2009}; ^{Lovell et al., 2003}; ^{Magnussen et al., 1999}).
The apparent canopy height profiles (CHPs) were obtained by curve-fitting all ALS returns with a height greater than 5 m. The threshold aimed to exclude laser points on the ground, in shrubs, and dominated trees, to eliminate the bimodal effect on the vertical profile (^{Coops et al., 2007}; ^{Lovell et al., 2003}). Curve-fitting analysis was carried out applying the maximum likelihood estimation technique (^{Cohen, 1965}) using fitdistr in R (R Foundation for Statistical Computing; MASS package). The Weibull pdf with two parameters (scale and shape) is shown in Eq. (2).
where: f (x / α,β) = probability density of x; α = scale parameter; β = shape parameter.
A Weibull pdf with a shape parameter equal to 1 reduces to an exponential distribution, while shape values between 3 and 4 will approximate a normal curve. The scale parameter is at the 63.2nd percentile of the distribution (^{McCool, 2012}). The Weibull scale and shape parameters were obtained from the CHPs at each training plot, and they were used as candidate predictors for stem biomass modeling.
Regression modeling of stem biomass with ALS data
The linear relationships between the ALS predictors and stem biomass were explored by carrying out a paired sample correlation t-test. We suspended ALS metrics from further analysis when the Pearson’s correlation coefficient (ρ) was not significant at 99 % confidence level. Regression models were fit by the ordinary least squares method, and the ALS predictors were selected using the best subset approach (^{Lumley, 2009}). We used the variance inflation factor (VIF) to detect multicollinearity of explanatory variables (^{Fox and Monette, 1992}). Models with VIF greater than 5 were excluded from further analysis (^{d’Oliveira et al., 2012}). Graphical analyses of residuals and hypothesis testing were performed to check the assumptions underlying the linear regression theory.
The regression models were submitted to leave-one-out cross validation (^{Picard and Cook, 1984}). This method uses the training data set without one of its observations (n-1) to predict the value removed from the sample (n is the sample size). This process occurs n times, so all observations are excluded once. The cross validation output was assessed by the relative root mean square error (rRMSE) statistic (^{Meng et al., 2009}; ^{Nyström et al., 2012}) (Eq. 3).
where: rRMSE = relative root mean square error (%); y_{i} = observed stem biomass in plot i (Mg ha^{−1}); = predicted stem biomass in plot i (Mg ha^{−1}); n = number of observations; mean of observed stem biomass (Mg ha^{−1}).
Furthermore, we compared observed values of stem biomass in the validation dataset with predictions by the best regression models. To make such a comparison, we used Pearson’s correlation coefficient statistic (Eq. 4).
where: = sample Pearson’s correlation coefficient; X_{i} and Y_{i} = observed values of variables X and Y; and = mean of variables X and Y.
Stem biomass interpolation with ordinary kriging
We used ordinary kriging to interpolate stem biomass from the training dataset as an alternative method (^{Viana et al., 2012}). The experimental semivariograms were obtained according to Eq. (5):
where: = semivariance estimate of class k in distance = mean distance of class k; N_{k} = number of pairs observed in class k;= residual (random error) observed in x_{i}; x_{i} = position i with coordinates x and y.
The following types of theoretical models were tested: spherical, exponential, linear, and Gaussian, and we picked the one with the least residual squares summation (^{Hiemstra et al., 2009}). Ordinary kriging interpolation was carried out according to Eq. (6):
where:= stem biomass estimate at position x_{0}; n = number of observations; Z(x_{i}) = observed value of stem biomass at position x_{i}; λ_{i} = weight assigned to observation Z(x_{i}), in which the sum of weights is 1 (^{Viana et al., 2012}).
The fitting of theoretical semivariograms and the ordinary kriging were conducted using autofitVariogram and autoKrige in R (automap package). As for the regression models, we used leave-one-out cross validation and the validation dataset to assess the ordinary kriging performance.
Spatial representation of stand attributes
The prediction maps were generated with a spatial resolution of 16 m, an approximate size of the validation plots. We also generated two extra basal area maps (linear regression and ordinary kriging) to visually compare results with one former basal area map obtained from a census of DBH on Feb 2008 (^{Campoe et al., 2012}). To our knowledge, no study has yet compared the quality of stem biomass maps generated from regression with ALS data and generated from ordinary kriging of field-derived data, in even-aged eucalypt plantations.
Results
All plots presented negatively skewed data, and considerable variability between plots (Figure 2). On average, 88 % of the laser points were observed above a height of 20 meters. The quantile-quantile plot showed the best fit for the Weibull pdf in the canopy upper layer (Figure 3), partially due to crowns occluding the penetration of laser beams at lower parts (^{Harding et al., 2001}). The integration of terrestrial laser scanning (TLS) data with ALS-derived CHPs was presented as an alternative for improving canopy modeling at lower layers (^{Zhao et al., 2013}).
The Weibull shape parameter (β_CHP) ranged from 5.7 to 11.3 confirming the negative asymmetry observed in the histograms of Figure 2 (Table 2). All the ALS metrics except β_CHP were positively correlated with the stem biomass. The metric P10 did not present strong correlation with the stem biomass, and it was excluded from further analysis (Table 3).
Stat. | α_CHP | β_CHP | P10 | P30 | P70 | P90 | Mean | St. dev. | Variance | p_understory |
---|---|---|---|---|---|---|---|---|---|---|
- | m | - | -------------- m -------------- | m^{2} | % | |||||
min | 23.4 | 5.7 | 18.6 | 21.9 | 23.8 | 25.0 | 22.4 | 2.8 | 8.0 | 4.0 |
avg. | 26.6 | 8.4 | 20.8 | 25.1 | 27.6 | 29.0 | 25.4 | 4.4 | 20.1 | 10.7 |
max | 28.9 | 11.3 | 24.4 | 28.0 | 30.7 | 32.2 | 27.7 | 6.0 | 36.4 | 16.5 |
st. dev. | 1.7 | 1.6 | 1.4 | 1.8 | 2.1 | 2.1 | 1.6 | 0.9 | 8.3 | 3.7 |
Stat. = statistics; α_CHP = Weibull distribution scale parameter; β_CHP = Weibull distribution shape parameter; P## = percentiles in height ##; Mean = mean height of first returns; St. dev. = standard deviation of heights from first returns; Variance = variance of heights from first returns; p_understory = proportion of points in the canopy understory layer; min. = minimum observed value among plots; avg. = average value among plots; max. = maximum observed value among plots; st. dev. = standard deviation of mean values.
# | Model | adj.R^{2} | AIC | VIF |
---|---|---|---|---|
1 | B = -287.3 + 16.2*P30 + 4.3*β_CHP | 0.93 | 151.1 | 3.43 |
(49.4) (1.5) (1.7) | ||||
2 | B = -207.2 + 13.7*α_CHP | 0.89 | 159.8 | - |
(27.63) (1.03) | ||||
3 | B = 254.0 - 11.6*β_CHP | 0.5 | 193.4 | - |
(20.9) (2.5) | ||||
4 | B = 104.9 + 4.9*p_understory | 0.5 | 193.5 | - |
(11.8) (1.0) | ||||
5 | ln(B) = -3.9 + 2.6*ln(P30) + 0.16*sqrt(β_CHP) | 0.93 | -71.6 | 3.33 |
(0.88) (0.23) (0.06) | ||||
6* | ln(B) = 2.9 + 0.09*P30 | 0.91 | -65.6 | - |
(0.15) (0.006) | ||||
7* | ln(B) = 4.6 + 0.1*P30 - 0.7*ln(P90)^{ns} | 0.91 | -64.8 | 22.4 |
(1.6) (0.03) (0.67) | ||||
8* | ln(B) = 2.6 + 0.09*Mean + 0.0002* sqrt(P70)^{ns} | 0.89 | -60.9 | 40.4 |
(0.86) (0.04) (0.38) | ||||
9 | B = -304.1 + 14.5*α_CHP + 0.05*stem.density | 0.92 | 153.4 | 1.1 |
(40.1) (0.9) (0.02) |
adj.R^{2} = adjusted coefficient of determination; AIC = Akaike information criterion; VIF = variance inflation factor; ns = not significant (p > 0.01); *= adapted from ^{Zonete et al. (2010)}.
The best model resulted from using the explanatory P30 variable combined with β_CHP (adj. R^{2} = 0.93), since using the P30 alone had already yielded a good fit (adj. R^{2} = 0.91). It was already expected that height metrics would have greater explanatory power than metrics of height dispersion (^{Bouvier et al., 2015}). The metric p_understory had only moderate correlation with the stem biomass (adj. R^{2} = 0.5) and no interaction effect with other ALS metrics was observed. On the other hand, there was an improvement issue when α_CHP was combined with the field-derived statistic density of stems (adj. R^{2} = 0.89 to adj. R^{2} = 0.92). No improvements were noted through using the natural logarithm of the stem biomass (Table 3).
The observed and predicted values for the stem biomass regressed on P30 and β_CHP (model 1) and on α_CHP (model 2) are shown in Figure 4. The observations close to the 1:1 diagonal indicate a good model fit. This was also corroborated by the leave-one-out cross validation, which demonstrated rRMSE to be equal to 5 % for both models.
Figure 5 shows the Gaussian semivariogram for the stem biomass constructed from the training dataset. The model presented a nugget effect equal to 128 (Mg ha^{−1})^{2} and a sill equal to 1129 (Mg ha^{−1})^{2} at a range of 574 m. The leave-one-out-cross validation from the ordinary kriging resulted in an rRMSE equal to 13 %, which was 160 % greater than the regression models.
The predictions from the regression models showed slightly higher correlation with the validation dataset than the ordinary kriging (ρ = 0.8, ρ = 0.82, and ρ = 0.71, respectively, Figure 6). The result is coherent with the work of ^{Meng et al. (2009)}, who studied different methods of kriging to estimate the basal area in pine forests in the state of Georgia, USA. They found that applying regression kriging using Landsat ETM+ data as the auxiliary variable improved the results compared to ordinary kriging interpolation (R^{2} = 0.9 and R^{2} = 0.75, respectively).
The map of the stem biomass generated from the metrics P30 and β_CHP were shown to be sensitive to the overlapping effect along adjacent flight lines (pixels scattered in the northwest-southeast direction) (Figure 7).
All maps show the gradient of productivity similar to that described by ^{Campoe et al. (2012)}; i.e., lower elevations in the terrain had higher stem biomass stock than the highest locations. However, the maps differed considerably in relation to local spatial patterns. A visual validation with the collinear variable of the stem biomass, basal area, is shown in Figure 8.
Discussion
There was great variability between the observed apparent canopy height profiles (CHPs), even though all trees in the stand are about the same age. The CHP is a signature of the forest structure, and is useful in a variety of contexts, like monitoring spatial and temporal changes (^{Coops et al., 2007}), mapping homogeneous strata (^{Nelson et al., 2003}), and identifying vegetation types (^{Harding et al., 2001}; ^{Jaskierniak et al., 2011}). ^{Dean et al. (2009)} were able to retrieve DBH values in a 36-year-old even-aged loblolly pine (Pinus taeda L.) stand from variable height to the base of live crown and height to crown median, which were retrieved from ALS-derived CHPs (n = 17, R^{2} = 0.97).
The ALS metric P30 presented the greatest correlation with stem biomass. ^{d’Oliveira et al. (2012)} observed a similar correlation between the height percentile P25 (all returns) and above-ground biomass in the Amazon rainforest, Brazil. In planted forests, similar results were observed by ^{Stephens et al. (2012)} and ^{Zonete et al. (2010)}. According to ^{Stephens et al. (2012)}, the lower percentiles can combine information of tree height and crown density, in which denser stands will present greater values for such metric. The significant correlation obtained between the scale parameter of the Weibull distribution (α_CHP) and stem biomass was consistent with the work of ^{Mori and Hagihara (1991)}. They succeeded in correlating the scale parameter of the Weibull distribution with DBH, when studying the vertical foliage profile from a hinoki stand (Chamaecyparis obtusa (Sieb. et Zucc.) Endl.) in central Japan. The improvement from adding the field-derived variable density of trees together with α_CHP seems to be promising owing to the potential of ALS data to quantify trees at the stand level (^{Görgens et al., 2015b}; ^{Popescu et al., 2003}; ^{Oliveira et al., 2012}).
The shape parameter of the Weibull distribution (β_CHP) was the only ALS metric negatively correlated with stem biomass. This result was also observed in ^{Coops et al. (2007)} in relation to the DBH variable in different mixed stands of Vancouver Island, Canada. The Weibull pdf shape parameter is related to the degree of data dispersion in the distribution. With the scale parameter fixed, the greater the value of the shape parameter, the smaller the curve width at the mode. In models where β_CHP was used together with P30, there was a mediation effect; i.e., β_CHP correlated positively with stem biomass. The hypothesis is that for similar canopy height layers, a smaller vertical dispersion of heights is indicative of homogeneity, which would lead to more productive stands (^{Stape et al., 2010}).
The ALS metric p_understory showed only moderate correlation with stem biomass and did not capture the variation in the forest horizontal structure as expected. This is because the most productive plots also had more ALS points intercepted by the crown, thereby underestimating the density of trees in the understory layer. Additionally, observations were influenced by the point density heterogeneity within the ALS data, partially created by the overlapping of swaths during the flight (^{Bater et al., 2011}; ^{Görgens et al., 2015a}). The selected regression models had at most two explanatory variables, and they were able to explain 93 % of the stem biomass variation. ^{Stephens et al. (2012)} explained 70 % of total carbon stock in forests of New Zealand with just one ALS height metric, observing that the addition of a density metric slightly improved the quality of the model (2 %), but a third variable did not bring any improvement.
The stem biomass maps constructed from the linear models with ALS metrics and from the ordinary kriging interpolation were consistent with the existing gradient of productivity shown by ^{Campoe et al. (2012)}. However, we observed from the validation statistics and visual inspection that the regression models fitted from the ALS metrics (P30 and α_CHP) generated more realistic maps, corroborating the initial hypothesis.