Hypsometric approaches to Eucalyptus experiments in Brazil

The objective of this study was to determine the best method to estimate tree height by different hypsometric approaches. This study is a provenance trial with Eucalyptus spp. in a 35-year-old pioneering plantation in the municipality of Lavras, Brazil. The census was taken by measuring the diameter at 1.30 m above ground and total height. All x and y coordinates of the individual trees were obtained and these trees were classified according to their stem quality. The ordinary kriging, cokriging and ordinary kriging of the residuals generated from the regression models were fitted by an exponential model. The Curtis model was selected for calculating the regression. The best method to estimate the height-diameter relationship was based on the statistical and graphical criteria. The spatial prediction model did not adequately represent the height–diameter relationship. The regression technique was more precise and accurate by inserting variables that helped capture different development standards of the trees. The regression-kriging produced more accurate total height estimates. Inclusion of the stochastic effect in the general spatial behavior of the total height variable helped to capture specific details about the stand. Therefore, this method is recommended for experimental areas.


INTRODUCTION
In Brazil, the forestry sector is important in generating products, taxes, foreign investment, jobs, and incomes, which represents approximately 5% of the national Gross Domestic Product (GDP).The importance of this sector in the Brazilian economy is mainly related to the adaptation of Eucalyptus spp. to different environments in the country as well as the high level of technology used in the silviculture and management of tropical plantations.Between 2006 and 2012, Eucalyptus plantations increased by 26.5% in Brazil (ABRAF, 2013).
Eucalyptus plantations in Brazil are primarily used for the production of cellulose and energy, and on a smaller scale, for sawnwood and lumber.Plantations also offer the additional benefit of minimizing the pressure on the commercial use of native forests.
In the context of development of the Brazilian forestry sector, the art of managing a Eucalyptus forest stand involves updating the growth and yield parameters of the forest at regular intervals, to optimize plans and future decisions that involve planting, harvesting, and planning the supply of the enterprise, in accordance with environmental and marketing variables.A key process for monitoring forest growth is the forest inventory.This facilitates the generation of benefits based on the information obtained; thus enabling, the quantification of the economic return of the forest.This involves continuous investment in the creation of maps, designing and installing forest inventory plots and measuring variables, such as diameter, height, and sometimes the canopies of individual trees (Avery and Burkhart, 2002).
Usually, in forest inventories the diameter of all the trees in the plot is measured at 1.30 m above ground.On the other hand, because of the difficulty of measuring height, only a few trees are measured.The hypsometric relationship is defined as the relationship between the height and diameter of a tree.Knowing this relationship is of great significance in the data collection system.By measuring only a portion of the heights and all the diameters in the plot, a mathematical relationship can be established using the height-diameter pairs measured, thus, enabling the heights of unmeasured trees in the plot to be estimated.This relationship is particularly important when it is more difficult to measure the height of the tree, as in older plantations, and it represents a financial saving of forestry inventory resources without much loss of precision (Silva et al., 2007).
Although, the hypsometric relationship is an excellent alternative for estimating unmeasured heights and lowering inventories cost, there are situations where deterioration of the height-diameter relationship makes it necessary to find other options.For example, in malformed stands or in locations with high variability, a weak correlation between diameter and height is expected, given that for the same diameter there will be different heights, and for the same height, different diameters (Scolforo, 2006).
Furthermore in clonal stands, the deterioration of the height-diameter relationship occurs.Although the genetic constitution of the plants is similar, there are interactions with the environment.This usually results in stands with similar heights but with a great variability in diameters, because this variable is very heavily influenced by management and the environment.This leads to a weak or nonexistent relationship between the two variables, that is, it is natural in these cases for the coefficient of determination to be close to zero (Barros et al., 2002).
Stands subjected to thinning present another case in which the deterioration of the height-diameter relationship occurs, because when thinning is performed, mainly smaller trees are removed, resulting in a trend towards height standardization.However, with wider spacing, the trees respond with an increased diameter growth, leading to greater differences in this variable.This occurs because the post-thinning or remaining trees have different diameters.Thus, the higher the number of thinnings, the lower the correlation between the height and diameter, and in some instances the correlation will be zero (Barros et al., 2002).Thus, factors impacting the hypsometric relationship should be observed.The relationship is affected by location, density, crown class, canopy size, species, and age (Curtis, 1967;Loetsch et al., 1973;Gómez-García et al., 2015).The development of geographic information systems and specific forestry inventory software enable the control of all these factors.Since there are several factors affecting the hypsometric relationship it is usually fitted at the plot level.The reason for this is because of the small plot size, similar location and canopy pattern, and the same density, species, and age (Tome et al., 2007).Therefore, consistent relationships are obtained that contribute to good forest inventory estimates.
Other options for large plantation areas are to make model fitting by adding plot covariables that help reduce variation, or combining the traditional hypsometric relationship model with the variables that affect this relationship, particularly age, density, and productivity measurements of the locations, as these are the variables that greatly impact this relationship.This method is the generic hypsometric relation model (Soares et al., 2004).
However, ignoring the spatial dependency that is often present in observations of forest data can cause incorrect inferences regarding the parameters of the model (Finley et al., 2013).Thus, geostatistics emerges as an alternative in forest inventory procedures, in addition to being associated with spatial dependency models for estimating height (Zhang et al., 2009;Lu and Zhang, 2011).The data obtained from specific collections, such as dendrometric variables from forest inventory, require the use of spatial prediction methods for their mapping, of which kriging is one such method.Besides the use of ordinary kriging, co-kriging, which is treated as a multivariate extension of kriging (Yamamoto and Landim, 2013), and regression-kriging, which comprises a combination of regression models and stochastic kriging behavior (Hengl et al., 2007), are alternatives with the potential to provide this spatialization.
In tropical countries that work with fast-growing species in rotations that often do not exceed eight years of age, the effect of environmental variables has much more impact than those species that grow in temperate climates (Almeida et al., 2004).Therefore, modeling studies in tropical countries are significantly more complex (Castro and Morrot, 1996).In a global context of Where: SQ = stem quality, N = number of individuals, CV (%) = coefficient of variation, HT and DBH measured in meters and centimeters, respectively.
more accurate height estimates in Eucalyptus stands, there is a gap in the scientific research on methods that evaluate the performance of methods that generate height estimates, which have exclusively focused on ordinary regression using traditional or generic models.This study proposes a novel form of comparative evaluation in the height estimation of Eucalyptus spp., a species that is synonymous of the forestry sector development in Brazil.Thus, the objective of this study was to define a method to estimate height, considering approaches that use the hypsometric relationship in spatial prediction models, ordinary regression model, regression models with addition of covariables, and a combination of regression models and kriging for experiment with Eucalyptus spp.

Description of the study area and data collection
The data used in this study is from a provenance trial of Eucalyptus spp. in a 35-year-old seminal plantation of 1.2 ha located on the campus of the Universidade Federal de Lavras (Federal University of Lavras, UFLA) in the municipality of Lavras, Brazil.The municipality of Lavras, according to Köppen's classification, has a humid subtropical climate, with hot and rainy summers, and cold and dry winters (Cwa climate), an average annual temperature of 19.3°C, an average annual rainfall of 1530 mm, and an average annual relative humidity of 76% (Sa Junior et al., 2012).A census of the area was conducted, measuring the diameter at 1.30 m above ground (DBH) and the total height (HT) of all the individuals in the stand.All individuals were also georeferenced by adopting a system of arbitrary coordinates (x,y) in meters, in which the first tree in the stand was designated as the reference tree, for which the values of the arbitrary coordinates were zero (0,0).
A qualitative evaluation of the trees in the stand was also conducted in relation to stem quality (SQ).This variable was used as a covariable in the regression model.It is expected that this will be a determining factor in differentiating growth standards.
SQ 1: Trees of large dimensions with straight trunks and with a potential use in sawmills for the manufacture of furniture and structures; SQ 2: Trees of medium to large dimensions, with a potential use for charcoal, fence posts, or stakes; SQ 3: Trees of small dimensions, suppressed or dominated, with twisted trunks that can be used in the production of charcoal, firewood, and in some cases as fence posts or stakes.
Table 1 presents descriptive statistics separately for each SQ and also for the stand.Based on an analysis of the data, it can be seen that the number of individuals for SQ 1 is greater than those of SQ 2 and SQ 3, that is, there is a greater tendency for the individuals of these stands to be used as lumber.Furthermore, the table shows that the variability of the data on HT and DBH is also increasing from SQ 1 to SQ 3, as would be expected by the characterization of quality of each trunk type, described above.

Variographic study
The variographic study was established from the generation of experimental semivariogram and its modeling for the application in ordinary kriging, co-kriging, and ordinary residual kriging; the latter for the application of regression-kriging.
In this study, the exponential model, described by Journel and Huijbregts (1978), was fitted to obtain the set of parameters to be used in the estimation of ordinary kriging, co-kriging, and regression-kriging.This model was chosen because it demonstrated good performance in previous studies in forestry, such as that of Galeana-Pizana et al. (2014).
The exponential model fitted for each situation was conducted by the method of weighted least squares.The fittings were made by using R software (R development core team 2013), using the geoR package (Ribeiro Jr. and Diggle, 2001) and gstat (Pebesma, 2004).

Spatial prediction models -Ordinary kriging and co-kriging
The prediction of unsampled points based on ordinary kriging was first described by Journel and Huijbregts (1978), and more recently by Meusburger et al. (2012) and Hengl et al. (2004).In this approach, weight was obtained with the parameters of the semivariogram fitted for the characteristic being evaluated (HT).The ordinary kriging estimator (Equation 1) is: Where: The estimation of unsampled points made from co-kriging was described by Yamamoto and Landim (2013), among other authors, and is an extension of kriging, in which one or more easy-tomeasure variables helps explain a difficult-to-measure primary variable (Wu et al., 2009).For this approach, the identification of weights was obtained with the parameters of crossed semivariance fitted for the variable HT, where the secondary variable was the DBH, which was obtained for all the individuals in the stand.The linear correlation of 0.79 between the variables HT and DBH was statistically significant.This calculation is relevant because a good correlation between the auxiliary and principal variables make cokriging more effective (Basaran et al., 2011).The co-kriging estimator (equation 2) is: Where:

Regression model
The estimation of HT was done by fitting the Curtis model (Equation 3), which has desirable efficiency describing height-diameter relationships, as the studies of Ribeiro et al. (2010) and Filho et al. (2010) have shown: Where: HT and DBH were previously defined; ln is the natural logarithm; and 0  and 1  are the parameters of the model;  is the error term.First, the fit was performed by the traditional method, which is based on ordinary least squares regression and the model used HT only as a function of DBH.The second form of fit was performed by adding categorical covariables into the model.For the second approach, additional characteristics that differ from the standard distribution of the variable of interest (HT) were incorporated.This addition being understood for this study as a categorical covariable.
The covariable considered in this study was SQ, that is, this variable was considered to be an important variable in changing growth patterns.Therefore, to fit the ordinary regression models by adding categorical covariables, two distinct situations were considered for the Curtis model.In the first, the three SQ's were inserted as dummy variables in the model intercept (ß0), which means that the fitted model have different intercepts for each SQ, but same slope.In the second situation, these three SQ's were inserted as dummy variables in the intercept and slope (ß0 e ß1) parameters of the fitted model, which means different ß0 and ß1 according to each SQ.

Regression-kriging
This is a hybrid form of fitting for the estimation of HT, that is, it is a combination of a regression and a spatial prediction model that lends a stochastic aspect to the estimates.The regression identifies the spatial behavior of the interest variable throughout the area, although without specific details of more specific locations (Mello et al., 2013;Scolforo et al., 2015).For the final result of the estimates to have more details about the specific points, a correction of the estimates developed by the regression model is necessary.This correction is made by applying ordinary kriging to the residuals generated by this model, with the main objective of correcting trends and detailing the spatial behavior of the principal variable (Meusburger et al., 2012).
Thus, after fitting the regression models, the next step was to study the behavior of the residuals derived from each regression model in the space.For each of the three situations (traditional regression model and two regression models with the addition of a covariable), the study of its spatial autocorrelation as a function of distance was performed, fitting the theoretical semivariogram models to the experimental semivariogram models by means of the R software (R Development core team, 2013) using the geoR package (Ribeiro Jr. and Diggle, 2001).Then the residual for each individual was estimated by ordinary kriging.
Given that all individuals in the area were georeferenced, continuous cells with dimensions of 3 × 3 m (regular planting spacing) were developed throughout the stand using the ArcMap program (ESRI, 2010).In each of these cells, the independent variable DBH was characterized and the regression model was applied to generate the HT estimates.Afterwards, the values based on the residuals kriging were extracted for each cell, enabling correction of the values produced by the fitted regression model.

Predictive validation of the different hypsometric approaches
A comparison of the estimates generated from the hypsometric approaches was made.The statistic of standard residual error (Syx) in meters (Equation 4), and residual average were analyzed.In addition, the predicted values around a 1:1 line was analyzed, to assess the variance behavior in terms of its homogeneity.
he dataset was separated into data to perform the fits and data for the predictive validation of the models, corresponding to 322 (70%) and 138 (30%) of the individual trees, respectively.The predictive validation means that 30% of the data was not used in fitting the models.Thus, the predictive validation enabled the evaluation of the best technique by statistical criteria (Crescente-Campo et al., 2010).
The data for fitting and predictive validation were separated randomly, with weighting proportional to the number of individuals per SQ: Where: () MS E is the mean square error; n the number of observations; and p the number of parameters.

Variographic modeling for ordinary kriging and cokriging of HT -Spatial prediction models
Figure 1 shows the experimental and modeled semivariograms for HT (a) and crossed experimental and modeled semivariograms for this same variable with the aid of the secondary variable DBH (b).The exponential model was appropriately fitted in both cases, allowing spatial continuity modeling of the variable of interest.
For the univariate (semivariogram) and multivariate (crossed semivariogram) approaches, nugget effect and sill values enabled similar spatial dependency degree (around 35%), which characterizes a moderate spatially dependent structure for both cases according to Cambardella et al. (1994).
The practical ranges of both modeled semivariograms that represent the maximum distance, within which the characteristic is spatially correlated, were 37 and 46 m for the univariate and multivariate approaches, respectively.
Thus, the HT variation was captured better by using the crossed semivariogram, which indicates the gain from the use of the auxiliary variable DBH.

Regression models for estimating HT
Table 2 contains the fitted parameters of Curtis models (via traditional regression and regression considering SQ as a covariable), in addition to the residual standard error (S yx ) in meters.The fitted regression models using the addition of covariables had better performance than the ones without, according to values of S yx (m).This improvement in results when compared with the traditional fitted regression model was expected, given that SQ captures the different growth standards, reducing the variability effect of the data and minimizing the estimation errors.

Variographic modeling for regression-kriging
Figure 2 shows the experimental and modeled semivariograms for the residuals of the HT variable derived for each of the three situations of the Curtis fitted regression model.An evaluation of Figure 2 showed that in each of the three situations the exponential model was fitted adequately, enabling modeling of the spatial continuity of the variable of interest.The nugget effect and sill values in all three situations allowed characterization of a moderate spatially dependent structure, according to the concept of Cambardella et al. (1994).

Predictive validation of the hypsometric approaches
Figure 3 presents the predicted values around a 1:1 line, and Table 3 shows the statistical precision for each hypsometric approach evaluated.The hypsometric estimates from the regression models and the regression-kriging models presented an adequate balance of predicted values around a 1:1 line (Figure 3).This situation characterizes these methods as having good statistical performance, that is, indication of Where: SQ = quality of the trunk, ß0 and ß1 are the parameters of the models.homoscedasticity (Mello et al., 2013).The spatial models (essentially a geostatistical approach) showed an overall bias in the predicted values around a 1:1 line.The spatial models presented low prediction capacity due to high variability between individuals (Table 1), thus, concentrating the estimates on intermediate values (Figure 3).The regression models did not present any trends and the addition of SQ as a categorical covariable improved the results when compared to the traditional Curtis regression model (Figure 3).Also, Figure 3 highlights that the incorporation of the stochastic aspect to the fitted regression (regression-kriging) added more precision to the results, that is, the regression-kriging was the most efficient fit technique.Furthermore, this technique resulted in a tighter dispersion of HT estimates and a more balanced distribution of residuals between the under and overestimates around a 1:1 line plot (Figure 3).Nanos et al. (2004) conducted a study in central Spain in a stand of Pinus pinaster, where they observed the The closer the S yx statistics and the average residual are to zero, the more reliable and less prone to error the generated estimates become.Thus, Table 3 shows that the regression-kriging and regression modeling approaches are marginally superior to the geostatistical approach (spatial modeling).The fact that co-kriging produced better results than ordinary kriging is also pointed out, although these geostatistical approaches are more often indicated for mapping of areas than point estimates in situations where the vegetational variability is great in a reduced space.
It is inferred that the Curtis model with regressionkriging, considering different ß 0 and ß 1 (different growth standards of the trees classified according to each SQ), presented the best estimates, as shown in Table 3.Given the marked variability of the data due to the advanced age of the Eucalyptus stand and pointed out that analyses done here came from the predictive validation, these validation values from the S yx and the residual average are satisfactory.It should also be highlighted that the residual graphic tends to go towards a balance between the sub and superestimates, which in fact characterizes the Curtis model with regression-kriging procedure as an estimator with low bias.
A major advantage of the application of regressionkriging models is that in addition to prediction, the productivity map with spatial information about the stand is also produced (Scolforo et al., 2015).As a continuation of the developments regarding the hypsometric relationship approach presented in this study, it is believed that combining the results obtained here with the use of LiDAR (Light Detection and Ranging) technology will make it possible to model the height-diameter relationship using regression-kriging, achieving greater benefits.LiDAR enables an extraordinary understanding of the terrestrial biomass through direct threedimensional measurements of the biophysical properties of the vegetational profile (Jaskierniak et al., 2011), as well as by obtaining the heights of all the individuals in the stand and their locations in the space (Lingnau et al., 2008;Oliveira et al., 2012).Giongo et al. (2010) pointed out that LiDAR has great potential for forestry applications, suggesting the use of this technique for operational forestry measurements in the near future.

Conclusions
When spatial modeling was applied with the goal of generating the hypsometric relationship in a stand of mature Eucalyptus spp. with high variability, it was not possible to fit the model adequately, because of the difficulty in generating kriging weights, which affected the estimates generated by the spatial models.
The use of the regression technique by adding stem quality (SQ) as a covariable to predict height estimates was more precise and produced greater accuracy than estimates obtained through traditional regression.The insertion of variables that help capture the different standards of tree development make this technique more accurate when compared with traditional regression.Regression-kriging was more accurate in estimates of HT.Inclusion of the stochastic effect in the general spatial behavior of the total height variable, helped to capture specific details about the stand.In other words, it makes corrections for locations where the behavior of the regression was not followed and presented higher rates of estimation errors.From a statistical and accuracy point of view, kriging combined with regression is the recommended method for modeling the height-diameter relationship.
); λi is the i th weight attributed to each i th observation of the variable of interest in position X ; and Zxi, is the observed value.
the estimate of the value in position 0 X ; n is the number of sampled neighboring points used to predict the unsampled points 0 ˆX Z (kriging neighborhood); λ1i and λ2i are the i

Figure 1 .
Figure 1.Theoretical and experimental unidirectional semivariogram (a) and theoretical and experimental unidirectional cross semivariogram (b) modeled for the HT variable.

Figure 2 .
Figure 2. Experimental and modeled semivariograms of the residuals relating to the models: (a) Traditional Curtis regression; (b) Curtis with different ß0 for each SQ; (c) Curtis with different ß0 and ß1 according to each SQ.

Figure 3 .
Figure 3. Distribution of the predicted values around a 1:1 line of the HTs estimated by the different approaches -(a) Spatial modelordinary kriging; (b) Spatial modelco-kriging; (c) Ordinary regression model -Curtis; (d) Ordinary regression model -Curtis with different ß0 for each SQ; (e) Ordinary regression model -Curtis with different ß0 and ß1 for each SQ; (f) Regression-kriging model -Curtis; (g) Regression-kriging model -Curtis with different ß0; (h) Regression-kriging model -Curtis with different ß0 and ß1.

Table 1 .
Coefficient of variation (CV), minimum, average and maximum values for the HT and DBH variables according to each stem quality (SQ).

Table 3 .
Statistical precision Syx and residual average generated by the hypsometric approaches.