Introduction 

Intensively managed forest plantations have an important role, both nationally and internationally, due to their ability to efficiently produce large quantities of biomass. This biomass can be utilized in forms suitable for bioenergy or, after processing into cellulose-derivatives, in materials such as paper. The myriad uses of fibers derived from intensively managed plantation forests can aid in the reduction of anthropogenic pressures on natural forests, while simultaneously generating millions of jobs ([22]).

In Brazil, 4.4 million jobs are generated directly and indirectly by forest plantations, and the total reforested area exceeds 7.6 million hectares ([9]). Of this area, 72% is composed of Eucalyptus spp., a key genus in Brazilian forest plantations. The reasons for the success of Eucalyptus spp. plantations in Brazil include Eucalyptus’ rapid growth, genetic improvement, and its tolerance of a wide range of climatic conditions ([32]).

In order to manage the growth and yield of these plantations, annual forest inventories are required. In general, the traditional approaches to conducting plantation inventories follow the precepts of sampling theory, which is based on ground measurements collected from simple random sampling plots. The forest attributes of interest are then inferred for an entire stand based on observations collected at each sample plot ([21]). However, considering the difficulties in field surveys such as sampling error, bias, and the intensive use of time and labor, remote sensing represents an important additional tool for monitoring large areas accurately and efficiently.

Both active and passive remote sensors have been presented as feasible alternatives for estimating forest variables. High spatial resolution optical data is useful at providing spectral information on species and condition ([35], [18]), while active remote sensing technologies, such as Airborne Laser Scanning (ALS), can provide highly accurate forest height data ([37]). When predicting forest inventory attributes using ALS data or aerial images, two main approaches can be applied: the Area-based Approach (ABA - [19], [36]) and the Individual Tree Detection (ITD) approach ([27], [11]).

To estimate volume in Eucalyptus spp. plantations Packalén et al. ([20]) applied the ABA approach. In the ABA, the response variable (e.g., wood volume per unit area in m³ ha^-1) is aggregated over a grid cell of predetermined size or over a plot size ([3]). This response variable is estimated using predictor variables, which are derived from the ALS metrics directly or indirectly, such as the mean height, height percentiles, densities and proportions, or from image, such as the reflectance ([19], [39]).

The ABA approach is most suited to the estimation of a number of crucial forest stand attributes over large areas ([19], [36]). The prediction of stand attribute variables from ALS data, such as biomass, volume, basal area, mean diameter, or average height have been shown to have similar or improved accuracy compared to traditional inventory methods ([36]). However, when tree-level information such as stem number or species is desired, the ABA approach is less suitable ([3]).

In this case, the ITD approach is the most appropriate to estimate individual tree characteristics. However, algorithms normally used in the ITD approach detect, on average, 75% of trees ([26], [18]). Therefore, the use of automatic ITD may lead to the non-detection of some trees, and can interfere with forest attributes prediction ([3]).

Vastaranta et al. ([34]) applied an alternative method to overcome the shortcomings of ITD detection, involving the aggregation of the tree numbers in a grid cell, according to ABA. By aggregating trees into grid cells it is possible to know how many trees were detected or not in each cell, so long as detection rates are compared with ground truth data and the most appropriate methodology to correct deviations is used. This approach, also used in this study, can lead to a reduction in the number of systematic errors, and to an increased accuracy in the estimation of the number of trees per hectare, a variable extremely important to forest managers ([33], [3]).

According to Packalén et al. ([20]), the advantage of using ITD approaches is that the correct detection of individual trees could reduce the number of plots needed for an accurate and unbiased inventory. Breindenbach & Astrup ([3]) also note that the ITD approach is inherently more intuitive than ABA because the response variables are related to each tree (the smallest unit in a management system), providing coordinates that can be tied to individual tree growth, as well as explicit links to tree growth and yield models. In the case of ALS data, ITD approaches require high density point clouds, which can be impractical because of the cost of data collection ([17]).

The merging of ALS and optical remote sensing data can assist in the accurate estimation of forest attributes, with algorithms based on height or brightness ([8], [19], [38], [35], [29], [18], [20]). Saarela et al. ([30]) integrated ALS and Landsat image data to improve accuracy of timber volume estimation as compared with estimates based on ALS and images variables alone. This comparison is also employed in the present study, though we used high spatial resolution images. For instance, the ABA and ITD approaches can be combined to estimate wood biomass or volume in forest plantations ([19], [2]).

In this study, we use ABA and ITD to estimate stand volume in Eucalyptus spp. plantations using ALS and high spatial resolution multispectral image data. To this goal, we first modeled the stand volume using standard ABA approaches, using different combinations of ALS-based and image-based metrics as independent model variables. Second, we used an ITD-based approach to estimate tree count, applying a correction method to improve single tree detection, which is integrated with the stand volume predicted using ABA. Finally, we merge the two approaches (ITD and ABA), that yields volume predictions at the single tree level. We conclude by discussing the possibilities of using ABA for individual tree volume predictions.

  Materials and methods 

Study area

The study was performed in a forest plantation owned by FIBRIA S.A. company, located in São Luiz do Paraitinga in São Paulo State, Brazil (23° 33′ S, 45° 34′ W and 23° 30′ S, 45° 31′ W - Fig. 1). The dominant species in the study area is Eucalyptus urograndis, which is a hybrid between Eucalyptus urophylla and Eucalyptus grandis ([7]). According to Köppen & Geiger ([12]), the climate is humid subtropical (Cwa) characterized by cold winters and rainy summers, with average annual precipitation of 1210 mm and a mean monthly temperature of 20.2 °C ([5]). The topography of the study area is complex, with elevations varying from 974 to 1266 m a.s.l. and the average slope is 24.3 degrees.

Fig. 1 - The study area near São Luíz do Paraitinga (São Paulo, Brazil), including the 6.5 year-old Eucalyptus urograndis forest plantation analyzed in this investigation.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

The forest stands in the study area were planted in December 2006 and harvested between April and May 2013, following a standard 6.5 year rotation period. All harvested wood was taken to a mill. The trees were planted in rows with fixed tree spacing (2 × 3 m), resulting in a density of 1667 stems per hectare.

Field measurements

Reference data were based on field measurements collected at 48 circular plots (r = 11.28 m, area of 400 m²) between April and May of 2013 prior to harvesting. In each plot, all diameter at breast height (DBH) were recorded, and the height of 10 trees were measured. Heights for all remaining trees in the plot were then estimated by applying the Curtis ([6]) hypsometric model, adjusted based on the 10 measured trees in each plot. Tree volume was estimated with locally-developed volume models based on the model proposed by Schöepfer ([31] - Tab. 1). The geographical location of the plots were obtained using high-precision GPS receivers (ProXRT^®, Trimble Inc., Sunnyvale, CA, USA) at the center of each plot.

After harvesting, all wood arriving at the mill was weighed. The wood density was obtained according to Archimedes’ principle and, therefore, the total volume of wood at the mill could be defined by weight-to-volume (W-V) relationships obtained in the density measurement process. A complete description can be found in Batista ([1]).

Remote sensing data

ALS data were acquired by a RIEGL LMS - Q680I^® system (Riegl Gmbh, Horn, Austria) simultaneously on 15 January 2012 with aerial photographs under free-cloud conditions. The average ALS point density was 5 pulses per m². An orthophoto was generated from aerial images with a pixel size of 0.15 cm. The 4 spectral bands were blue (429-514 nm), green (514-600 nm), red (600-676 nm) and near-infrared (695-831 nm).

Initially we separated ALS ground returns from the vegetation returns, to generate a 15 cm resolution DTM and DSM, which is the minimum size resolution based on point density ([15]). The DTM was used to normalize the point cloud, followed by the generation of the most usual metrics (such as tree height, canopy cover, topography, image bands and ratios) from each type of remote sensing data at a 5-meter resolution. This resolution was chosen following the ABA concept that pixel size should be larger than a tree crown ([36]), with Eucalyptus spp. crown diameter varying from 2 to 3 meters in the study area.

Predictor variables

We extracted LiDAR metrics from the ALS data ([16], [36], [28], [30]) to extract image metrics (bands, ratios between bands and Normalized Difference Vegetation Index), as listed in Tab. 2. We generated the most common metrics (such as tree height, canopy cover, topography, image bands and ratios) from each type of remote sensing data by describing the properties of the point clouds and spectral properties of the images. For the generation of these metrics, only first returns of the ALS point cloud were used and the reflectance value of each pixel in the aerial image was aggregated into a 5 × 5 m grid cell. These metrics were then used in subsequent statistical modeling as explanatory variables. In total, 39 metrics were generated (Tab. 2).

Modeling approach to estimate total wood volume

The first step was the estimation of wood volume using the ABA. Three different methods were applied, depending on the input data source. In the first approach, only variables derived from ALS data were used, in the second approach only variables derived from the imagery were used, and in the final one we combined both sets of variables.

One of the most important considerations when working with a large number of predictor variables in a multiple linear regression is variable selection. To that end, we applied the Best Subset Selection (BSS), which matches an appropriate model to each model size ([10], [14], [28]). To select a model among those suggested by BSS, we applied a k-fold cross-validation (k=5). The application of a cross-validation, similar to the leave-one-out method, allows for a better understanding of model behavior, especially when using a database that was not utilized for model adjustment ([10]). After applying the cross-validation we chose the model according to the RMSE. In this study, we used a confidence interval with a p-value equal to 0.05 generated from 1000 simulations.

The total volume of the forest stands was obtained by the mean ABA estimation multiplied by the total area.

Variable selection and model adjustment

The k-fold estimate of RMSE is computed by averaging the values obtained by eqn. 1:

where CV_k is the k-fold cross validation and MSE is the mean square error. The linear model assumptions were verified according to ([23]), where the assumptions were verified for each of the chosen models with statistical tests and graphical analysis ([24], [28]). The approach was then chosen based on the lowest RMSE and the linear model assumptions analysis, generating a predicted wood volume map for the entire study area. This prediction was validated with stand inventory data and information about tree volume received by the mills, after harvesting.

Individual tree volume

Individual Tree Detection based on aerial images

For the second step we used high resolution multispectral images for detecting individual treetops, applying the local maximum filtering method, where the central location of the crown is assigned the highest spectral value ([35], [13]). We applied fixed and varying sizes of moving window, based on the methodology proposed by Kumar ([13]) and Pereira ([25]). The window sizes varied from 0×3 m, 0×5 m, 0×7 m to 0×9 m ([13], [25]). The individual tree detection procedure was performed for each of the 4 spectral bands, and the first principal component (first PCA) of the bands ([35]). According to Wang et al. ([35]), the PCA effectively aggregates a large quantity of information and shows good results for ITD.

To generate the final ITD product, all five intermediate outputs were used. First, the layer with the highest number of detected trees, with their locations visually corrected, was taken as the basis for ITD. From this layer, we removed all the treetops in a 0.7 m buffer around each tree to avoid multiple peaks from single tree. Treetops detected in the second layer, outside this buffer were added creating a new layer of additional ITD points. This layer was then taken as the new basis for ITD. This procedure was repeated until all layers had contributed to the generation of a single ITD layer.

Merging the ITD approach with ABA

Based on the final ITD layer, stem counts were aggregated into 5 × 5 m grid cells, according to ABA. Following this, the number of trees in each grid cell was identified by local maxima derived from the ITD approach, resulting in a graphical representation of the number of trees per hectare for the entire study area.

Detection and correction of the mean difference (MD) in total area

The mean difference (MD) detection of the ITD approach was done by comparing the number of detected trees and the number of observed trees in the inventory plots. The difference between the two was considered the response variable, while the ALS and image variables were considered the predictor variables, leading to the modeling of the ITD mean difference, and its prediction for the entire study area. In our case this MD was modeled with 3 variables, maximum height, slope, and ITD detection, as explanatory variables. These variables were selected by applying the Best Subset Selection method and cross-validation, as described earlier. The MD detection model was applied to the entire study area, representing the MD per hectare, to correct the initial estimates of tree count per pixel.

Estimating individual tree volume

In order to generate the volume per tree in each 5 × 5 m cell, the previously-estimated stand volume for the entire study area (m³ ha^-1), obtained by the integration of ALS data and high-resolution imagery, was divided by the number of trees per hectare detected by the ITD approach.

  Results 

Modeling approach to estimate total wood volume

The best fitting models obtained by cross-validation contained 4 parameters, including the intercept, for all the three approaches tested (Tab. 3). In the first approach (using ALS data variables only), the chosen model was based on height percentiles (h_p20, h_p40, h_p50) and elevation (A_mean). In the second approach, where only image predictors were used, the selected variables were B3, NDVI, B1/B3 and B1/B4. In the third approach, where both remote sensing data sources were used as predictors, two ALS-based and two image-based variables were selected (h_mode, A_mean, B4, B2/B3). The RMSE for each approach was 6.84%, 8.45% and 5.23%, respectively. Among the three analyzed approach, we chose the model that combined ALS data and image information, according to the corresponding adjusted R² value (Tab. 3).

Fig. 2 displays the relationship between the observed and predicted values of stand volume obtained for each approach. Although the homoscedasticity of residuals was not rejected, their distribution was asymmetric in the first and second approach, suggesting non-normality, while in approach 3 the distribution of residuals showed to match a normal distribution. Therefore, the third approach was selected for wall-to-wall stand volume prediction and subsequent comparisons with the inventory and mill-based volumes as reference values.

Fig. 2 - Scatterplot of the observed and predicted values of stand volume for each of the approaches analyzed in this study. (1) ALS data; (2) aerial images; (3) both datasets.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Across the whole study area, the mean estimated volume per hectare based on the ABA method was 319.9 m³, totaling 44.466.1 m³. When this result was compared with the volume measured at the mill (considered the true volume by the company), we found a difference of 2.98%, while when comparing the forest inventory data with the mill volume, it overestimated by 3.29%. In addition, it was observed that the combination of both ALS and imagery data could provide more information regarding the occurrence of drought and canopy gaps in the forest.

Semi-ITD approach to estimate wood volume per tree

Overall, using the ABA method 151.048 trees were detected in the study area, with an average of 994 trees ha^-1, while the average number of trees per hectare estimated by the forest inventory was 1708 trees ha^-1. Therefore, approximately 42% of the trees were not detected using the above approach, reinforcing the need to correct the remotely sensed tree count by modeling the tree count mean difference (MD).

After identifying the predictive mean difference model, we performed a linear adjustment using the least squares method. The model selected used three variables: maximum height (h_max); the slope (S) estimated from ALS data; and the number of trees detected in the aerial images found using the semi-ITD approach (Tab. 3), which had an RMSE of 12.66%.

The average number of trees per hectare detected by the conversion of ITD into ABA (semi-ITD) for the study area after mean difference (MD) correction was 1707.7, while the inventory data provided a figure of 1708.3 trees per hectare. Therefore, the mean volume per tree estimated from study data was 0.188 m³, while the mean volume of an individual tree derived from forest inventory was 0.195 m³. Fig. 3 shows the wall-to-wall tree volume maps generated in this work.

Fig. 3 - Maps of ABA-derived predictions for stand volume (top left), corrected number of trees (top right) and mean tree volume (bottom).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

  Discussion 

Accurate stand attribute estimates are critical for effective forest plantation management, partly due to their markedly shorter rotation ages and faster growth rates relative to natural forests ([20]). Moreover, for proper management of forest plantations, information on individual trees is of high value ([33]). Remote sensing has proven to be an excellent data source for supporting sustainable forest management ([37]). However, obtaining accurate estimates of forest attributes is often still elusive ([30]). Our results indicate that the combination of different remote sensing data sources to estimate stand volume can lead to improved accuracy of predictions (RMSE = 5.23%).

In this research, we used two approaches to estimate stand volume in Eucalyptus spp. plantations using ALS as well as high spatial resolution multispectral image data. We first modeled the stand volume using ABA with different combinations of ALS- and image-derived metrics as input data, in three different models. Despite ABA slightly underestimated stand volume (2.98%) compared to the reference volume (total wood volumed measured at the mill) the results were close to those obtained by the forest inventory. However this approach could be enhanced by predicting the number of trees using the ITD approach, this introducing a correction method to improve tree count predictions. We then integrated that information with stand volume previously predicted using the ABA method, allowing for the calculation of mean individual tree volume.

The model applied in the volumetric prediction used two variables derived from aerial imagery, the near infrared band (B4) and the green/red (B2/B3) ratio, and two ALS-derived variables of height (h_mode) and cover (A_mean). Using the ratio between the green and red bands is analogous to using the near infrared and red bands ratio. The reflectance in the near infrared band is linked to scattering in the mesophyll, which occurs due to the interaction of solar radiation with the internal leaf structure ([4]). Similarly, green wavelengths are reflected by mesophilic leaf tissues, while red light is absorbed by chlorophyll ([4]). Therefore, the ratio between these bands is indicative of vegetative vigor and vegetative abundance. Both vigor and abundance information were important complements to the ALS metrics.

Despite the use of high-resolution imagery, the second approach (which used the aerial image variables) showed the poorest performance among the three approaches tested (RMSE = 8.45%).

Our results show that use of information from a single remote sensing technology can meet the predictive needs of plantation managers, so long as expectations of accuracy and precision are attained accordingly. However, our results also highlight that the integration of data from multiple remote sensing technologies, such as the ALS and high-resolution passive imagery data used here, can lead to more accurate results than the use of a single technology.

This improvement in accuracy is achievable even in situations where the difference in time between data acquisition (2012) and field measurement (2013) is up to one year, as is the case in this study. Because ALS data acquisition can be costly, the field measurement-remote sensing hybrid approach can be of remarkable benefit to forest managers. The integration of data from different sources with different acquisition dates can increase the useful life of forest attribute data collected by any means, potentially decreasing the frequency of acquisition and therefore reducing costs.

While ABA provided estimates of total wood volume for the entire study area, individual tree delineation, detected from passive image data, allowed us to downscale the estimates to the individual tree level. Images were used for tree detection instead of ALS data because the trees were planted at a relatively high density (3 × 2 m) resulting in small canopy sizes. Their small size made it challenging for crowns to intercept a sufficient number of LiDAR points for crown extent to be accurately defined. Therefore, the performance of ALS was poorer when detecting individual trees compared to image-based approaches.

Despite improvements over ALS data, we still found the accuracy of single tree detection from images to be unsatisfactory (58% correctly detected trees) and so a correction factor was introduced. This correction, based primarily on slope, allowed us to increase the accuracy of single tree detection and achieve unbiased individual tree volume estimates. Slope proved to be a key variable to include in the modeling of tree detection error. With increasing slope more trees were located in the same projected area, leading to an aggregation of the canopies and an increase in the maximum tree height of the inventory plots. Increased maximum tree height, therefore, showed a higher correlation with the standard deviation (i.e., heterogeneity) of a plot, which made tree detection from the images more difficult.

In general, the correction procedure applied to our tree counting process performed well, showing good agreement with existing tree count estimates and leading to individual tree volumes with negligible mean difference and low RMSE (12.66%). Our corrected estimates of tree count were greater than uncorrected values identified from the image-based detection by 713.7 trees per hectare, very close to extant estimates based on the forest inventory. The final estimates of individual tree volume after mean difference correction were 0.188 m³ per tree, 3.5% less than the inventory estimate ([33], [20]).

We demonstrated that a traditional ITD approach can be modified by using ABA, in spite of losing the exact location of each tree, but enabling to obtain the desired tree density information ([3]).

  Conclusions 

Accurate information about wood volume at both the stand and individual tree levels is required to support the effective management of forest plantations. In this study we enhanced the traditional area-based approach using individual tree detection. The integration of geometrical and spectral information provided by ALS and aerial imagery, respectively, allows for the reliable predictions of volume at both spatial scales. Moreover, remote sensing data integration results in lower mean difference and RMSE. Our approach allowed for the estimation of individual tree volume with negligible bias.

  Acknowledgements 

The authors would like to acknowledge the Brazilian Research Foundation, CNPq, and the National Institute for Space Research (INPE) for supporting this research, FIBRIA S/A cellulose company for supplying data from their plantation, and the Integrated Remote Sensing Studio (IRSS) at the Faculty of Forestry, University of British Columbia, Canada, for its assistance, support and friendship. David Williams provided some editorial assistance.

  References