Modelling stand biomass fractions in Galician Eucalyptus globulus plantations by use of different LiDAR pulse densities

Introduction Aboveground biomass is the total amount of biological material (usually oven-dried to remove water) present above the soil surface in a specified area. Because almost 50% of plant biomass is carbon, estimates of the total aboveground biomass in forest ecosystems are critical for carbon dynamics studies at multiple scales (Drake et al., 2003). Estimation of the biomass and carbon stock of trees has gained importance in recent years, especially since the Kyoto Protocol of the United Nations Framework Convention on Climate Change (UNFCCC) entered into force on 16 February 2005. Signatory countries must estimate carbon stocks in 1990 and any changes since 1990 from all afforestation, reforestation and deforestation activities (UNFCCC, 1997). In addition, the continued growth in energy demands in technologically developed societies and the requirement to reduce the substantial use of fossil fuels have made it necessary to diversify means of energy production. In this sense, the forest biomass in established plantations is attracting great interest as a renewable resource for biomass energy production, because of the following advantages: reduction in net warming and environmental Modelling stand biomass fractions in Galician Eucalyptus globulus plantations by use of different LiDAR pulse densities E. Gonzalez-Ferreiro1*, D. Miranda1, L. Barreiro-Fernandez1, S. Bujan1, J. Garcia-Gutierrez2and U. Dieguez-Aranda3 1 Laboratorio del Territorio. Departamento de Ingeniería Agroforestal. Escuela Politécnica Superior. Universidad de Santiago de Compostela. C/ Benigno Ledo. Campus Universitario, s/n. 27002 Lugo, Spain 2 Departamento de Ciencias de la Computación, Lenguajes y Sistemas. Escuela Técnica Superior de Ingeniería Informática. Universidad de Sevilla. C/ Reina Mercedes, s/n. 41012 Sevilla, Spain 3 Unidad de Gestión Forestal Sostenible. Departamento de Ingeniería Agroforestal. Escuela Politécnica Superior. Universidad de Santiago de Compostela. C/ Benigno Ledo. Campus Universitario, s/n. 27002 Lugo, Spain Abstract Aims of study: To evaluate the potential use of canopy height and intensity distributions, determined by airborne LiDAR, for the estimation of crown, stem and aboveground biomass fractions. To assess the effects of a reduction in LiDAR pulse densities on model precision. Area of study: The study area is located in Galicia, NW Spain. The forests are representative of Eucalyptus globulus stands in NW Spain, characterized by low-intensity silvicultural treatments and by the presence of tall shrub. Material and methods: Linear, multiplicative power and exponential models were used to establish empirical relationships between field measurements and LiDAR metrics. A random selection of LiDAR returns and a comparison of the prediction errors by LiDAR pulse density factor were performed to study a possible loss of fit in these models. Main results: Models showed similar goodness-of-fit statistics to those reported in the international literature. R2 ranged from 0.52 to 0.75 for stand crown biomass, from 0.64 to 0.87 for stand stem biomass, and from 0.63 to 0.86 for stand aboveground biomass. The RMSE/MEAN · 100 of the set of fitted models ranged from 17.4% to 28.4%. Models precision was essentially maintained when 87.5% of the original point cloud was reduced, i.e. a reduction from 4 pulses m–2 to 0.5 pulses m–2. Research highlights: Considering the results of this study, the low-density LiDAR data that are released by the Spanish National Geographic Institute will be an excellent source of information for reducing the cost of forest inventories. Key words: Eucalypt plantations; airborne laser scanning; aboveground biomass; carbon stocks; remote sensing; forest inventory. * Corresponding author: edu.g.fer[email protected] Received: 19-12-12. Accepted: 18-07-13. Instituto Nacional de Investigación y Tecnología Agraria y Alimentaria (INIA) Forest Systems 2013 22(3): 510-525 Available online at www.inia.es/forestsystems ISSN: 2171-5068 http://dx.doi.org/10.5424/fs/2013223-03878 eISSN: 2171-9845 pollution, decreased risk of sacrificing natural areas, enhanced regional welfare, and avoidance of competition with food production (Freppaz et al., 2004; Field et al., 2008). Therefore, the measurement of forest biomass provides an indication of carbon sequestration in trees and also an estimate of cellulosic material as a potential source of renewable energy (Popescu, 2007). Eucalyptus spp. are the most productive tree species in Galicia (NW Spain). Eucalypt plantations, which occupy 215,899 ha in the region, accumulated approximately 52% of the harvested volume in the period 2005-2010 (3,574,500 m3harvested in 2010, with an increment of 15% in the period 2009-2010) (Confemadera, 2010). The most widespread management regime for these plantations is short rotation forestry (SRF), with initial densities of 1,000 to 2,400 stems ha–1 in single stem stands at first rotation and subsequent replanting or coppicing (Pérez-Cruzado et al., 2011). The plantations are generally managed, without thinning, for wood pulp production (and to a lesser extent, chip-board), and they provide small wood that is not used by the sawnwood industry. Bioenergy production is one potential use for these plantations (Pérez-Cruzado et al., 2011). Despite the high productivity and potential of eucalyptus forest as carbon sinks, there are important deficiencies in the characterization and quantification of the biomass resource available in Galicia. This has arisen because forest information is usually acquired through conventional forest inventory methods, which often use detailed measurements taken from a small set of sampling plots distributed randomly or systematically over the forest area. The quality of the estimates is thus limited by the cost of establishing sufficient sample plots to measure the existing variability (Lovell et al., 2005; Musk and Osborn, 2007; Rombouts et al., 2008). Such inventories are extremely labour intensive and expensive (Hall et al., 2005), especially in small areas of privately owned land, and in practice they do not allow full inventory coverage of large areas on the ground or can extend over several years (Boudreau et al., 2008). Airborne Light Detection and Ranging (LiDAR) is currently considered the most promising remote sensing technology for forest inventory and biomass estimation (Boudreau et al., 2008). It is capable of offering detailed tri-dimensional information regarding the size and structure of the forest canopy (Reitberger et al., 2008; Wagner et al., 2008), which is closely related to aboveground biomass (Lim et al., 2003b). This has potential for timely and accurate measurement of tree biomass components and carbon stored over time (Naesset and Gobakken, 2008; Rosenqvist et al., 2003) for the whole area, with an equal or better accuracy than other remote sensing techniques (Bortolot and Wynne, 2005). Two main approaches are used for biomass estimation with commercial small-footprint discrete return laser scanning data. The “single-tree” approach used with spatially dense LiDAR data is based on laser detection and on the measurement of individual tree parameters, usually tree height and crown diameter (Bortolot and Wynne, 2005; Popescu et al., 2003; Popescu, 2007). The “stand-level” approach used with low resolution LiDAR data establishes empirical relationships between plot-level stand measurements and the height distribution obtained from laser returns or other laser measurements related to canopy density. This approach has been used in regression analysis to estimate forest biomass and carbon stock levels across a range of forest types (Lim et al., 2003a; Lim and Treitz, 2004a, 2004b; Patenaude et al., 2004; Andersen et al., 2005; Hall et al., 2005; Thomas et al., 2006; Naesset and Gobakken, 2008; Sherrill et al., 2008; García et al., 2010; Treitz et al., 2010; González-Ferreiro et al., 2012) as well as shrub biomass (Estornell et al., 2011). Metrics derived from LiDAR data are highly dependent on the tree species involved (Heurich and Thoma, 2008) and are also correlated with field assessment of various aspects of vegetation structure, which may influence the derived relationships (Goodwin et al., 2006). Small footprint discrete return LiDAR has been used in eucalyptus forests to estimate a variety of stand variables (Wack et al., 2003; Tesfamichael et al., 2010; Gonçalves-Seco et al., 2011), and to explore the effects of airborne LiDAR acquisition parameters on vegetation structural assessment (Goodwin et al., 2006). However, there are no reports of biomass estimation of Atlantic eucalyptus plantations of medium to high density, characterized by very small crowns with sparse leaves that may facilitate laser penetration. Laser pulse density and vegetation structure are the factors with the greatest effects on the height accuracy of the laser-derived Digital Terrain Model (DTM). If the ground elevation or the uppermost portion of the forest canopy is not well detected, the normalised heights of the trees and the Digital Crown Model (DCM) obtained will be underestimated (Hyyppä et al., 2008). Goodwin et al. (2006) argued that point density is even more important than footprint size or flight altitude in determining certain forest variables such as crown area Modelling stand biomass fractions with LiDAR 511 and volume. Magnusson (2006) observed a significant increase in the root mean square error for mean height and stand volume estimates when the laser pulse density was greatly decreased, while other authors (Maltamo et al., 2006; Thomas et al., 2006; Gobakken and Næsset, 2007; García et al., 2010; Treitz et al., 2010; González- Ferreiro et al., 2012) argued that the estimation of most forest stand variables may not be affected. This should be verified in different forests and regions to determine whether LiDAR data can be collected at a lower average point spacing for a more cost-effective ground coverage, while maintaining the accuracy of the estimates (Thomas et al., 2006). The aim of this study was to estimate biomass fractions in Eucalyptus globulus Labill. plantations from height and intensity data gathered with a small footprint discrete return LiDAR system. The effect of a reduction in pulse density on model precision was examined, in light of the availability of low-density LiDAR data (0.5 pulses m–2) for most of the Spanish territory. Materials and methods Study area The study area is a 1 ×4 km rectangle located in the municipality of Vilapena (Galicia, NW Spain) (see Fig. 1). This area is characterized by a wide variety of landform types and an elevation range of 150 to 530 m. The forests in this area are representative of Eucalyptus globulus stands in NW Spain, characterized by lowintensity silvicultural treatments and by the presence of tall shrub. The canopy of mature Eucalyptus globulus plantations is usually high and sparse, thus enabling sunlight to penetrate to the ground, which is often aggravated by severe defoliation caused by the eucalyptus weevil Gonipterus scutellatus Gyll. Such low resistance of the canopy to light penetration allows growth of herbaceous sub-shrub cover dominated by shade-intolerant fruticose species and both shadeintolerant and shade-tolerant herbs below the canopy, which often produce a large amount of biomass (about 2.5 to 3.2 Mg ha–1 total dry matter) (Silva-Pando et al., 1993). The already well-known loss of precision in DTMs in forest areas, together with steep slopes and the high variability and quantity of understory, poses a serious challenge to the accuracy of LiDAR technology (Raber et al., 2002; Hodgson et al., 2003, 2005). LiDAR data The LiDAR data were acquired in November 2004 with an Optech Airborne Laser Terrain Mapper (ALTM) 2033 sensor (www.optech.ca) operated at a laser wavelength of 1,064 nm from a flight altitude of 1,500 m above sea level. The beam divergence was 0.3 mrad, the pulsing frequency 33 kHz, the scan frequency 50 Hz, and the maximum scan angle ±10∫. The first and last return pulses were registered. The whole study area was flown in 18 strips and each strip was flown three 512 E. Gonzalez-Ferreiro et al. / Forest Systems (2013) 22(3): 510-525 Figure 1. Inventory plots. 500 250 0 500 1,000 m 1:25,000 644500 645000 645500 646000 644500 645000 645500 646000 486500 487000 487500 488000 488500 489000 489500 4810000 4810500 4811000 486500 487000 487500 488000 488500 489000 489500 4810000 4810500 4811000 N Study area Centroids of square plots times, which gave an average measurement density of about 4 pulses m–2. Field data A total of 39 square plots of 225 m2were located and measured in the Eucalyptus globulus plantations in the study area, between February and March 2005. The plots were subjectively selected to represent the existing range of ages, stand densities and sites in the regions. Topographical surveys were carried out using total stations and GPS to determine the location of the four corners and the position of every tree within the plots. First, a Trimble®5800 GPS (Trimble, Sunnyvale, CA, USA, www.trimble.com) (dual-frequency realtime kinematic receiver with a planimmetric precision of ±5 mm+0.5 ppm and a altimetric precision of ±5 mm +1 ppm) was used to obtain the coordinates of a densified geodetic network for the study area by applying realtime kinematics. Based on the network established with the GPS, a topographic survey of the plots was conducted using a Trimble®5603 Robotic Total Station (Trimble, Sunnyvale, CA, USA, www.trimble.com) with a precision in distances measurement of ±2 mm+2 ppm and a precision in angles measurement of 3 to 5”. Observations of the four corners of each plot and of the position of each tree within the plot were made during the survey. For all the trees in each sample plot, two measurements of diameter at breast height (1.3 m above ground level) were made at right angles, with a tree calliper. Measurements were made to the nearest millimetre, and the arithmetic mean of the two measurements was calculated. Total tree height was measured to the nearest decimetre with a Vertex III hypsometer (Haglöf Sweden AB, Långsele, Sweden, www.haglof.se). The dry weight of the biomass fractions of each tree were estimated from the following equations for Eucalyptus globulus in Galicia reported by Diéguez-Aranda et al. (2009): [1] [2] [3] [4] [5] where ww is stem wood biomass (kg), wb7is wood and bark biomass on branches with 7 cm minimum top diameter (kg), wbis bark biomass on stem (kg), wb2–7 is wood and bark biomass on branches with 7 cm maximum butt diameter and 2 cm minimum top diameter (kg), wb0.5–2 is wood and bark biomass on branches with 2 cm maximum butt diameter and 0.5 cm minimum top diameter (kg), wb0.5 is wood and bark biomass on branches with 0.5 cm maximum butt diameter (kg), wl is needles biomass (kg), d is diameter at breast height outside bark (1.3 m above the ground level, cm), and h is total tree height (m). Finally, crown biomass (wcr), stem biomass (wst) and aboveground biomass (wabg) were calculated from the sum of the biomass fractions included: [6] [7] [8] The field measurements (heights and diameters) and the estimated dry weight of the biomass fractions were used to estimate the following stand variables in each plot, on a per hectare basis: stand crown biomass (Wcr), stand stem biomass (Wst), and stand aboveground biomass (Wabg). The estimates were used to develop models to derive these stand variables from LiDAR data. Preparation of LiDAR data The LiDAR data provided by the contractor only provided information, for each laser pulse emitted by wst =ww+wb7+wb wst =ww+wb7+wb wcr =wb2−7+wb0.5−2+wb0.5 +wl wl=0.02949 d1.917 wb0.5 =0.01258 d1.705 wb2−7+wb0.5−2=0.003685 d2.654 wb=0.01010 d2.484 ww+wb7=0.01308 d1.870h1.172 Modelling stand biomass fractions with LiDAR 513 Table 1. Forestry stand parameters of the sample plots (n = 39) Variable Average Minimum Maximum Standard deviation N1,663 622 3,378 770 d 13.5 9.2 19.3 2.7 Hm17.4 11.8 22.3 2.4 Hd23.5 13.3 36.5 4.1 G 25.8 7.2 47.0 9.6 V 229.4 39.9 511.1 105.1 Wcr 17.4 4.5 35.6 7.0 Wst 114.6 19.9 276.9 53.9 Wabg 132.0 24.4 312.4 60.7 N: number of stems per hectare. d: diameter at breast height outside bark (1.3 m above ground, cm). Hm:mean height (m). Hd:dominant height (m). G: stand basal area (m2ha–1). V: stand volume over bark (m3ha–1). Wcr:stand crown biomass (Mg ha–1). Wst:stand stem biomass (Mg ha–1). Wabg:stand above ground biomass (Mg ha–1). the sensor, about return type (first and last), X, Y and Z coordinates, and intensity values. Single returns, i.e. those that hit a solid surface, were recorded twice in the original dataset with the same information but different return type code. The theoretical nominal density of the LiDAR data was 4 pulses m–2, which implied collecting 8 returns m–2. Original LiDAR data with different resolutions, i.e. obtained from flights carried out in the same area, at the same time, and with the same flight parameters but at different LiDAR pulse densities, are the most desirable for assessing the influence of density on the extraction of forest information. However, the lack of such data can be overcome by artificially reducing the original LiDAR point cloud. Although this approach does not fully mimic the actual reduction that would be obtained with different flights, it allows study of the most important factor in canopy height modelling. For this purpose, several alternatives have been proposed, as follows: (i) a random reduction over the entire dataset (Anderson et al., 2006; Liu and Zhang, 2008; Puetz et al., 2009), (ii) a systematic reduction in each scan line (Gueudet, 2004; Raber et al., 2007; Magnusson et al., 2007; Treitz et al., 2010), and (iii) a reduction by randomly maintaining one return within a grid cell of a specific size (Magnusson, 2006; Gobakken and Næsset, 2007). The first alternative may not preserve the order and regularity of the original LiDAR data. The second requires information about the scan line and the order in which points were recorded. The third option overcomes the possible lack of regularity and allows for the use of data without scan line information. Because of the limited information provided by the data used in this study, random selection of LiDAR returns in a grid cell of 1 m2was carried out (note that this procedure does not allow for the simulation of other potentially important parameters such as flight height or scan angle). Examination of the original point cloud showed that many cells included more than 8 returns m–2. Therefore, to obtain a regular distribution of LiDAR returns, and to investigate the effect of the LiDAR point cloud density on the estimation of stand variables, two datasets were generated: one with a final density of 1 return m–2 and other with 8 returns m–2 (equivalent to 0.5 and 4 pulses m–2). Fig. 2a and 2b show the two LiDAR datasets generated for the plot number 3. Intensity is a radiometric constituent of LiDAR data (Singh et al. 2010), and it is recorded by the sensor as the amount of energy backscattered from objects or earth’s surface. The intensity values recorded by the sensor remain unchanged under different conditions of illumination caused by, e.g., shadows or occlusions (Donoghue et al., 2007; Höfle and Pfeifer, 2007), but they are affected by other factors, such as properties of the terrain, topography, flight and sensor characteristics, and atmospheric conditions (Donoghue et al., 2007; Höfle and Pfeifer, 2007; Mazzarini et al., 2007). Some authors have recommended considering range, incidence angle and atmospheric attenuation for intensity normalization (Höfle and Pfeifer, 2007; Gross et al., 2008; Jutzi and Gross, 2010). However, such data are not always available, and other authors recommend normalizing intensity values to a userdefined standard range in order to remove the range 514 E. Gonzalez-Ferreiro et al. / Forest Systems (2013) 22(3): 510-525 Figure 2. LiDAR cloud for the plot number 3 with (a) 4 pulses m–2 and (b) 0.5 pulses m–2. Scale bars show ellipsoidal height of the LiDAR returns for the plot number 3. Heigth 391.85 384.77 377.69 370.60 363.52 Heigth 390.94 384.19 377.44 370.68 363.93 Y+ X+ X+ Y+ LDV - V1.62 - USDA Forest Service Pacific Nortwest Research Station LDV - V1.62 - USDA Forest Service Pacific Nortwest Research Station a) b) dependency of the intensity signal (Donoghue et al., 2007; Mazzarini et al., 2007; García et al., 2010). In this approach, the normalized intensity values (I’) is obtained by multiplying the raw intensity value (I) by the quotient of the range of each return (rg), calculated as the difference between the average flying height and the height of each return, i.e. the difference between the sensor-target distance, and the standard range (rgs) (e.g. 1,000 m) (Eq. [9]): [9] This method eliminates the path length variations in the intensity recorded by the system, by providing values equivalent to the intensity that would have been recorded if all points were at the same range (García et al., 2010). In the present case, range and scan angle data were not available for every return. The terrain was very steep, with slopes above 40° in some areas and a topographic range of 150 to 530 m above mean sea level, and therefore the range was normalized to a userdefined standard range using Eq. [9]. For each return, the range (m) was estimated as the difference between the average altitude of the flight (1,500 m above sea level) and the ellipsoidal altitude of the return (m). This approach should not cause large errors in range computation because of the small scan angle (±10°) (García et al., 2010). Fig. 3a and 3b show the 0.5 pulses m–2 LiDAR dataset for the plot number 3, before and after the intensity normalization process. Extraction of LiDAR variables For the generated datasets (0.5 and 4 pulses m–2), FUSION software (McGaughey, 2009) was used to perform filtering, interpolation and DTM/DCM generation operations, as well as to compute the following variables related to the height and return intensity distributions metrics within the limits of the 39 field plots: mean, maximum and minimum values, mode, standard deviation, variance, interquartile distance, coefficients of skewness and kurtosis, average absolute deviation, and percentiles. The percentage of returns above a specific height threshold was also calculated. The following steps were carried out with several processing programmes implemented in the FUSION LIDAR Toolkit (McGaughey, 2009). First, ground returns were extracted from the LiDAR point cloud with the GroundFilter tool, which implements a filtering algorithm adapted from Kraus and Pfeifer (1998) and based on linear prediction (Kraus and Mikhail, 1972). Second, these returns were used to generate a DTM grid with the GridSurfaceCreate tool, which computes the elevation of each grid cell from the average elevation of all points within the cell; if the cell does not contain any points, it is filled by interpolation from the neighbouring cells; the cell size value was 1m 2. Third, the normalized LiDAR point cloud was obtained by subtraction of the ellipsoidal height of the DTM from the Z coordinate of each LiDAR return with the ClipData tool (the height switch in combination with the previously generated DTM were used); this ′ I=Irg2 rgs 2 Modelling stand biomass fractions with LiDAR 515 Figure 3. 05 pulses m–2 LiDAR dataset for the plot number 3. Scale bars show the intensity of the LiDAR returns: (a) before normalization and (b) after normalization. X+ X+ Y+ Y+ LDV - V1.62 - USDA Forest Service Pacific Nortwest Research Station LDV - V1.62 - USDA Forest Service Pacific Nortwest Research Station Intensity 140.00 105.25 70.50 35.75 1.00 Intensity 57.74 43.41 29.07 14.74 0.40 a) b) tool was used also to exclude returns below a normalised height of 2 m (zmin switch=2), which were considered as not belonging to tree crowns (e.g. hits on shrubs, rocks and logs). Fourth, the normalised LiDAR point cloud was clipped using the boundaries of the 39 field plots, which were previously stored as polygons in ESRITM shape files (ESRI, 1998). The PolyClipData tool allowed extracting the LiDAR cloud within the limit of each plot area, so that, an independent LiDAR cloud file was created, covering each plot area. Fifth, the metrics of heights and return intensity distributions of these 39 clipped and normalized point clouds were computed with the CloudMetrics tool. Fig. 4 shows the 4 pulses m–2 LiDAR dataset for the plot number 3 with the normalized height of the LiDAR returns and the generated DTM. Regression models Linear, (multiplicative) power function and exponential models were used to establish empirical relationships between field measurements and LiDAR variables. Their respective general expressions are as follows: [10] [11] [12] where Yare field values of Wcr (kg ha–1), Wst (kg ha–1), and Wabg (kg ha–1), and X1, X2…, Xnmay be variables related to the metrics of heights and return intensity distributions or measurements related to canopy closure. All variables were computed from the datasets with resolutions of 4 and 0.5 pulses m–2. The following variables can be related to height distribution: hmin, hmax, hmean, hmedian, hmode,hSD, hskw, hkurt, hID, hAAD, h05, h10, h20,…, h90, h95, h25 or h75, which are the minimum, maximum, mean, median, mode, standard deviation, coefficients of skewness and kurtosis, interquartile distance, average absolute deviation, percentiles, and first and third quartiles values of the height distribution of laser returns for each plot (m), respectively. The variables related to return intensity distribution are based on the same statistics as the variables related to height distribution, but in this case they are denoted with the letter iinstead of h. The variables related to canopy closure can be either r2, which is the number of returns above 2 m height for each plot, or c2-FP, which is the ratio of the number of laser hits above 2 m height to the number of first returns for each plot, expressed as a percentage. The additive error term ε is assumed to be normally, independent and identically distributed with zero mean. Model fitting and selection Linear models were fit by ordinary least squares, by applying the REG procedure of SAS/STAT®(SAS Institute Inc., 2004). Power function and exponential models were fit by nonlinear regression, by use of the Gauss-Newton method implemented in the NLIN procedure of the same statistical package. In a previous step, the latter two types of models were linearized by taking natural logarithms from both sides of Eqs. [11] and [12], in order to select the best subset of independent variables to be included in each and to obtain initial estimates of the parameters using the linear regression technique (Myers, 1990, p: 444). Once the complete linear form of the models was specified, the models were examined to determine whether all terms should be retained in the final regression equations. This involved fitting a number of subset models and comparing the relative performance of each (Clutter et al., 1983, p: 318). Although there are different approaches for selecting the subset models to be fit by linear regression (Draper and Smith, 1998, chapter 15), the Mallows’ Cpselection method of the Y=exp β0+β1X1+β2X2+…+βnXn () +ε Y=β0X1 β1X2 β2…Xn βn+ε Y=β0+β1X1+β2X2+…+βnXn+ε 516 E. Gonzalez-Ferreiro et al. / Forest Systems (2013) 22(3): 510-525 Figure 4. 4 pulses m–2 LiDAR dataset for the plot number 3 and the generated DTM. Scale bar shows the normalized height of the LiDAR returns. X+ Y+ LDV - V1.62 - USDA Forest Service Pacific Nortwest Research Station Heigth 22.53 16.90 11.27 5.63 0.00 REG procedure, which performs all possible subset regressions and lists the models in ascending order of Cp, was used. Heteroscedasticity was examined visually, by plotting residuals as a function of predicted values, but any discernible trend or heteroscedasticity was found. Multicollinearity among the explanatory variables was also checked with the condition index. Regressions with a condition index above 30 were disregarded, as recommended by Belsley (1991). Finally, only models in which all the parameter estimates were significant at 5% level were considered. Comparison of the estimates for the different models (linear, power function and exponential) was based on numerical and graphical analyses of the residuals. The following three statistics were calculated: the coefficient of determination (R2), the root mean square error (RMSE), and the Bayesian Information Criterion (BIC) proposed by Schwarz (1978). Although BIC was used as the final criterion for model selection (Peña, 2002, p: 570), it does not provide an intuitive idea of model precision. R2(also referred to as pseudo-R2when applied in nonlinear regression) indicates the proportion of the total variance of the dependent variable that is explained by the model. Although there are several shortcomings associated with use of R2in nonlinear regression, the general usefulness of some global measure of model adequacy appear to override some of those limitations (Ryan, 1997, p: 424). The RMSE provides an idea of the precision of the estimates in the same units as the dependent variable. Effects of density reduction Using a similar procedure to that proposed by Treitz et al. (2010), but considering only laser pulse density as the factor of interest, prediction errors from the best linear, power function, and exponential models developed for each dependent variable were calculated for each plot: [13] where eij is the prediction error of the ith plot (from 1 to 39) associated with the jth LiDAR pulse density (0.5 pulses m–2 and 4 pulses m–2), Yiis the stand variable obtained from field measurements for the ith plot, and yˆij is the corresponding predicted value for the ith plot and the jth LiDAR pulse density. The nonparametric Kruskal-Wallis one-way analysis of variance by ranks (Kruskal, 1952; Kruskal and Wallis, 1952) was used to compare the prediction errors distributions by LiDAR pulse density factor because of the non-normal distribution of the prediction errors for each model type and dependent variable, as suggested by the Shapiro-Wilk test for normality (Shapiro and Wilk, 1965, 1968) and the examination of Quantile- Quantile plots. If the computed value of the test suggests rejecting the null hypothesis, there is a high likelihood that the two samples represent populations with different median values (Sheskin, 2004, pp: 757-761). Results Regression models The parameter estimates and goodness-of-fit statistics of the best model developed by type (linear, power function and exponential), dependent stand variable (Wcr, Wst, and Wabg) and generated dataset (0.5 and 4 pulses m–2) are shown in Tables 2-4. For all the dependent variables and datasets, exponential models performed best (as shown by the BIC values), followed by linear and power function models. In Wcr modelling, exponential models provided R2values of 75.3% and 71.8% for the 4 and 0.5 pulses m–2 datasets, respectively (Table 2); the respective differences between the best and worst (power function) models, in terms of R2, were 11.7% and 19.6%. In Wst modelling, exponential models provided R2values of 86.6% and 84.1% for the 4 and 0.5 pulses m–2 datasets, respectively (Table 3); the respective differences between the best and worst (power function) models, in terms of R2, were 9.2% and 19.7%. In Wabg modelling, exponential models provided R2values of 85.7% and 83.1% for the 4 and 0.5 pulses m–2 datasets, respectively (Table 4); the respective differences between the best and worst (power function) models, in terms of R2, were 9.7% and 19.9. Only independent variables related to the metrics of height distribution proved to be reliable statistics for predicting the three dependent variables Wcr, Wst, and Wabg (Supplementary data includes plot information of the values of dependent and explanatory variables used in the final models presented in the Tables 2, 3 and 4). Further analyses were only carried out for exponential models. Figs. 5a, 5b, 6a, 6b, 7a and 7b show the field-mea- sured versus predicted values (using exponential models) of crown, stem and aboveground biomass fractions for the 0.5 and 4 pulses m–2 datasets, respectively. eij =Yi−ˆ Yij Modelling stand biomass fractions with LiDAR 517 Effects of density reduction For all the pairwise comparisons, the computed chisquare approximation of the Kruskal-Wallis test statistic (H) was below the tabled critical 0.05 chi-square value (χ1 05 =3.84) for 2-1 degrees of freedom (df). The alternative hypothesis (no equality of medians) was not supported at the 0.05 level in any of the cases (H<χ1 05, α=0.05 and p>0.05; see Table 5). Therefore, the results suggest that the reduction in the LiDAR point cloud has no effect on exponential model fit. Comparison in terms of R2can provide a more intuitive idea of the non loss of fit caused by reduction of LiDAR density. Variations in R2of 3.5, 2.5 and 2.6% were obtained for Wcr,Wst, and Wabg exponential models, respectively. The selected exponential models for Wcr, Wst, and Wabg include the same independent variable (h75) for the two datasets (4 and 0.5 pulses m–2) (Tables 2-4). Discussion A statistical approach based on regressors, which were calculated directly from the previously normalized laser-derived canopy height and intensity distribu- 518 E. Gonzalez-Ferreiro et al. / Forest Systems (2013) 22(3): 510-525 Table 2. Results of Wcr modelling for 0.5 and 4 pulses m–2 datasets Pulses Model Independent Parameter Standard t-p>|t|R2RMSE BIC m–2 variable estimate error1value (kg ha–1) 0.5 Linear β0parameter –20,897 4,086 –5.11 < 0.0001 0.708 3,812 650.5 h90 1,916 202.3 9.47 < 0.0001 Power β0parameter 108.2 101.4 1.07 0.2929 0.522 4,879 669.8 function h60 1.770 0.319 5.55 < 0.0001 Exponential β0parameter 7.544 0.264 28.56 < 0.0001 0.718 3,745 649.1 h75 0.1164 0.0129 9.02 < 0.0001 4 Lineal β0parameter –17,319 3,610 –4.80 < 0.0001 0.719 3,736 648.9 h75 1,835 188.4 9.74 < 0.0001 Power β0parameter 51.81 43.76 1.18 0.2421 0.636 4,256 659.1 function h60 2.009 0.285 7.05 < 0.0001 Exponential β0parameter 7.576 0.231 32.78 < 0.0001 0.753 3,504 643.9 h75 0.1126 0.011 10.24 < 0.0001 1For the power function and exponential models, the standard errors are approximates values. Table 3. Results of Wst modelling for 0.5 and 4 pulses m–2 datasets Pulses Model Independent Parameter Standard t-p>|t|R2RMSE BIC m–2 variable estimate error1value (kg ha–1) 0.5 Linear β0parameter –16,3950 23,804 5,424 < 0.0001 0.801 24,732 800.0 h75 15,802 1,319 11.98 < 0.0001 hskw 1,2920 5,425 2.38 0.0226 Power β0parameter 103.8 107.0 0.97 0.3383 0.644 32,594 817.9 function h60 2.434 0.349 6.97 < 0.0001 Exponential β0parameter 8.785 0.240 36.60 < 0.0001 0.841 21,773 786.4 h75 0.1493 0.0114 13.10 < 0.0001 4 Linear β0parameter –169,209 23,446 –7.22 < 0.0001 0.803 24,261 794.9 h75 15,017 1,223 12.28 < 0.0001 Power β0parameter 39.57 33.29 1.19 0.2422 0.774 25,991 800.2 function h60 2.745 0.282 9.73 < 0.0001 Exponential β0parameter 8.902 0.200 44.51 < 0.0001 0.866 19,985 779.7 h75 0.1406 0.0093 15.12 < 0.0001 1For the power function and exponential models, the standard errors are approximates values. Modelling stand biomass fractions with LiDAR 525 plantations using small-footprint, discrete return lidar data. Prog Phys Geog 34: 515-540. Thomas V, Treitz P, McCaughey JH, Morrison I, 2006. Mapping stand-level forest biophysical variables for a mixedwood boreal forest using lidar: an examination of scanning density. Can J For Res 36: 34-47. Treitz P, Kevin L, Murray W, Doug P, Nesbitt D, Etheridge D, 2010. LiDAR data acquisition and processing protocols for forest resource inventories in Ontario, Canada (Koch B, Kendlar G, eds). In: Silvilaser 2010: the 10th International Conference on LiDAR Applications for Assessing Forest Ecosystems Freiburg, Germany. pp: 1-10. UNFCCC, 1997. Kyoto Protocol to the United Nation Framework Convention on Climate Change [online]. Available in: http://unfccc.int/kyoto_protocol/items/ 2830.php [19 Dec 2012]. Wack R, Schardt M, Lohr U, Barrucho L, Oliveira T, 2003. Forest inventory for Eucalyptus plantations based on airborne laser scanner data. Int Arch Photogram Rem Sens Spatial Inform Sci 34: 40-46. Wagner W, Hollaus M, Briese C, Ducic V, 2008. 3D vegetation mapping using small-footprint full-wave- form airborne laser scanners. Int J Remote Sens 29: 1433- 1452.