Use of LiDAR to estimate stand characteristics for thinning operations in young Douglas-fir plantations

Light Detection and Ranging (LiDAR) has been successfully used to describe a wide range of forest metrics at local, regional and national scales. However, little research has used this technology in young Douglas-fir stands to describe key stand characteristics used as criterion for operational thinning. The objective of this research was to develop models of Douglas-fir mean top height, basal area, volume, mean diameter (at breast height), green crown height and stand density from LiDAR and stand information. Data for this study were obtained from four widely separated young (age range of 9 to 17 years) Douglas-fir plantations in the South Island, New Zealand. LiDAR was acquired for the entire area and stand metrics were measured within 122 plots established across the study area. Spatially synchronous stand and LiDAR metrics were extracted from the plots. Using this dataset, multiple regression models were developed for each of the six stand metrics. The final models constructed for mean top height, green crown height, total stem volume, mean diameter, basal area, and stand density had R2 values of 0.85, 0.79, 0.86, 0.86, 0.84 and 0.55, respectively, with root mean square errors of 1.02 m, 0.427 m, 20.2 m3 ha-1, 13.9 mm, 3.81 m2 ha-1 and 355 stems ha-1, respectively. With the exception of stand density, all relationships were relatively unbiased. Variables with the greatest contribution (with the partial R2 in brackets) to models of mean top height, green crown height, volume, mean diameter and basal area included the 75th (0.85), 1st (0.76), 10th (0.83), 95th (0.74), and 10th (0.72) LiDAR height percentiles. The LiDAR height interquartile distance was the most important contributor (partial R2 = 0.33) to the model of stand density. With the exception of stand density, the final models for stand metrics were sufficiently precise to be used for scheduling thinning operations. This study demonstrates the utility of LiDAR to accurately estimate key structural attributes of young Douglas-fir and to assist with forest management over a widely dispersed resource.


Background
High productivity and superior wood properties have made Douglas-fir (Pseudotsuga menziesii [Mirb.] Franco) one of the premier and most widely planted species throughout temperate forest areas. Globally, there are ca. 15 million hectares of Douglas-fir plantations with substantial areas occurring in Europe, South America, New Zealand, Australia and western North America (Hermann and Lavender 1999). Douglas-fir forest plantations are some of the most productive in the world (McMurtrie 1993).
Thinning is an important management operation and has a major impact on Douglas-fir stand development. Douglas-fir responds well to thinning and this operation is necessary to avoid stand stagnation and achieve merchantable sized logs over economically viable time frames. Typically, one to two thinnings are undertaken during a rotation and depending on the size of felled trees these are classed as either pre-commercial or commercial thinnings. A number of relatively simple (Reukema 1975;Emmingham and Green 2003) and more complex growth modelling methods (MacLaren and Knowles 2005) have been developed that simulate the impact of different thinning regimes on future growth and rotation end stand dimension. Important metrics describing the state of the current stand that affect thinning decisions include height, stand density (stocking), tree diameter at breast height of 1.4 m, volume, basal area and green crown height (Reukema 1975;Emmingham and Green 2003).
Stand-level inventory is the most widely used method for determining current characteristics and the need for thinning in Douglas-fir plantations. However, over the last two decades, traditional field based assessment has been progressively replaced or supplemented with Light Detection and Ranging (LiDAR) in a range of plantation species. LiDAR is a laser-based remote sensing technology that has been widely used to accurately predict stand height Coops et al. 2007;Means et al. 2000;Naesset 2002;Means et al. 1999) and volume Coops et al. 2007;Means et al. 2000;Naesset 2002;Means et al. 1999). Correlations of moderate to high strength have been found between LiDAR metrics and basal area (Naesset 2004b(Naesset , 2005(Naesset , 2002Nord-Larsen and Schumacher 2012;Means et al. 1999;Means et al. 2000), stem diameter (Naesset 2002) and green crown height (Naesset and Okland 2002). However, stand density is typically only predicted with a moderate degree of precision from LiDAR cloud metrics (Naesset and Bjerknes 2001;Hall et al. 2005;Naesset 2002). In contrast to traditional stand level inventory, LiDAR is often less expensive and more accurate (Naesset et al. 2004). Predictions from LiDAR metrics can be used to characterise spatial variation in tree characteristics allowing stratification of within-stand thinning operations. Despite the relatively wide adoption of LiDAR within the forestry sector, little research has used the technology to characterise a wide variety of stand characteristics in young Douglas-fir plantations.
The objective of this research was to investigate the utility of discrete LiDAR for predicting variation in important stand characteristics in young Douglas-fir plantations.

LiDAR data and stand measurements
The data used in this study were obtained from four Douglas-fir plantation forests in the South Island, New Zealand ( Figure 1). These four forests named Doon, Pentland, Saddle Peak and Shag River had areas of 1393, 1526, 561 and 2685 ha, respectively. Each forest was divided into a number of stands, based on age, with stand ages respectively ranging from 9-14, 14-17, 9-12 and 14-16 years old at Doon, Pentland, Saddle Peak and Shag River.
The LiDAR survey was conducted using a fixed wing aircraft between May and August 2012 using a small footprint (~0.20 m) Optech ALTM 3100EA system and a swath overlap of 50%. The LiDAR and flight parameters used to achieve first-return densities that averaged 15.6 points m -2 (range 5.7 to 28.2 points m -2 ) are summarised in Table 1. The system also utilised an Applanix 510 Position and Orientation System (POS) that uses the GPS and inertial measurement unit (IMU) sensors, and a GPS-based computer controlled navigation system.
Measurements of stand characteristics were made from July to October 2012 across the resource in plots ranging in area from 0.015 to 0.1 ha. High grade differential GPS units were used to locate the plots. After exclusion of plots with species other than Douglas-fir and four plots in which the digital terrain model was poorly defined, there were 122 plots available for the modelling. Of these 122 plots, 28, 29, 11 and 54 plots were respectively allocated to the Doon, Pentland, Saddle Peak and Shag River forests. Within each plot diameter (at breast height; standard height of 1.4 m used in New Zealand) and green crown height were recorded for all trees and heights were recorded for at least 10 trees. Stand density and basal area were determined for each plot and mean top height, H m (height of the 100 largest diameter trees per hectare) was derived from plot mean height, H mean , and stand density, S, using the following function, Total volume, V, was determined from H m and basal area, B, using the following function, Analyses Models were developed using SAS software (SAS Institute Inc. 2008), to estimate mean top height, green crown height, total stem volume, mean diameter, basal area and total stand density. Explanatory variables used in the modelling included LiDAR metrics, stand age and stand density (determined from plot measurements). The LiDAR metrics included a range of LiDAR height percentiles from the first (H 01 ) to the ninety-ninth (H 99 ), several metrics describing the LiDAR height distribution through the canopy (skewness, coefficient of variation (CV), standard error, kurtosis, elevation interquartile distance (IQ)) and the percentage of returns reaching within 0.5 m of the ground (PC zero ). Initial exploratory analyses were undertaken to assess the strength of linear relationships between all dependant stand characteristics and the explanatory variables. These models were then extended through the development of multiple regression models to ensure that robust models sensitive to a broad range of LiDAR metrics were created. For each multiple regression model, explanatory variables were sequentially introduced into the model starting with the variable that exhibited the strongest correlation with the dependent variable and subsequently the residuals from the model, until further additions were not significant, or did not improve the model coefficient of determination by more than 0.01. Variable selection was undertaken manually, one variable at a time, and plots of residuals were examined prior to variable addition to ensure that the variable was included in the model using the least biased functional form.
Model precision was determined using the coefficient of determination (R 2 ) and the root mean square error Flying height (m) 750 (RMSE). Model bias was examined by plotting actual values against predictions and predicted values against residual values and dependant variables included in each model. The variance inflation factor was used to ascertain the degree of multicollinearity within the multiple regression models with values over 10 signifying problematic multicollinearity. Predictions from these models were not used to spatially estimate forest characteristics in this study.

Ranges in LiDAR, stand metrics and slope
There was marked variation in LiDAR and stand metrics among plots. Mean top height, green crown height, volume, mean diameter, basal area, and stand density, respectively, ranged, six-fold, 17-fold, 226-fold, six-fold, 73-fold and 12-fold across the plot series (Table 2). Ranges in LiDAR metrics were similarly marked, with the 5 th , 50 th and 95 th LiDAR height percentiles ranging 41-fold, 10-fold, and seven-fold, respectively (Table 2). Variation was also wide for PC zero , which averaged 5.99%, ranging from 0.06 − 66.5% (Table 2). Plot slope averaged 16°and ranged from 2 to 41°.

Bivariate correlations
Green crown height was most strongly related to H 01 with correlation coefficients diminishing as LiDAR height percentile increased (Table 3). The strongest relationships between green crown height and other LiDAR metrics were with the LiDAR height coefficient of variation and skewness (Table 3). Green crown height also exhibited strong positive relationships with both mean top height and volume (Table 3). Mean top height was strongly correlated with the upper LiDAR height percentiles peaking with the seventyfifth percentile, H 75 (Table 3). In contrast, both basal area and total volume were more strongly correlated with the lower LiDAR height percentiles with the highest correlation coefficients occurring at the tenth percentile H 10 for both characteristics (Table 3). For mean diameter, the strength of correlation coefficients increased with LiDAR height percentile reaching a peak value of 0.86 at the ninety-ninth percentile H 95 . Stand density exhibited insignificant to weakly significant correlations with LiDAR height percentiles. The strongest bivariate relationship for stand density was observed with LiDAR height interquartile distance (Table 3).
Simple models describing the strongest linear relationships between LiDAR metrics and each of the stand characteristics were constructed ( Figure 2). The positive linear relationships between LiDAR metrics and mean top height, green crown height, total stem volume, mean diameter and basal area had respective coefficients of determination (R 2 ) of 0.85, 0.71, 0.83, 0.74 and 0.72. The negative linear relationship between total stand density and LiDAR height interquartile distance had an R 2 of 0.33 ( Figure 2).

Multiple regression models Model formulation
Mean top height was best described by a simple linear relationship with H 75 (Table 4). This relationship accounted for 85% of the variance in the data (Table 4).
Green crown height was best described by H 01 and skewness in the LiDAR height distribution (Table 4). Inclusion of H 01 in a quadratic form accounted for 76% of the variation in the data while addition of skewness as a negative linear term increased the R 2 by 0.03 to 0.79 (Table 4).
The best model for volume included H 10 , PC zero and stand age (Table 4). The strongest contributor to the model was H 10 which when included using a quadratic form had a partial R 2 of 0.84. Inclusion of PC zero as a negative linear term added a further 0.01 to the model while the addition of age as a positive linear term contributed a further 0.01 to the overall model. The final model had a R 2 of 0.86 (Table 4). The best model for mean diameter included H 95, stand density, stand age and PC zero ( Table 4). The largest contributor to the model was H 95 and when included in the model as a positive linear term accounted for 74% of the variance in the dataset. Stand density and PC zero were included in the model as negative linear terms and accounted for a further 4 and 5%, respectively, of the variance in the data. Inclusion of age as a positive linear term added a further 0.03 to the model R 2 and the final model R 2 was 0.86. Basal area was best described by a model that included H 10, PC zero , stand density and stand age ( Table 4). The positive linear relationship with H 10 accounted for 72% of the variance, while a negative linear relationship with PC zero added a further 6%. Stand density and age were both included in the model as positive linear terms and contributed a further 2 and 3%, respectively, to the model. The final R 2 was 0.84.
The best model of stand density included LiDAR height interquartile distance, H IQ , (height difference between the 25 th and 75 th height percentiles), PC zero and stand age (Table 4). All variables were included in the model as negative linear terms and H IQ , PC zero

Model precision and bias
With the exception of stand density, the coefficients of determination for the final developed models were high ranging from 0.79 to 0.86. The total amount of variance accounted for by the multiple regression model of stand density was markedly lower at 0.55. The RMSE of the final models was relatively low for all variables apart from stand density. Values of RMSE for the final models of mean top height, green crown height, stem volume, mean diameter, basal area and stand density were, respectively, 1.02 m, 0.427 m, 20.2 m 3 ha -1 , 13.9 mm, 3.81 m 2 and 355 stems ha -1 . All variables included within the six models were highly significant (Table 4). The highest variance inflation factor in for any of the terms in the multiple regression model was 1.97 indicating that multicollinearity was well within acceptable limits. Plots of predicted against actual values for all models, apart from stand density, showed little apparent bias (Figure 3). For these five models there was also little pattern in the residuals when these values were plotted against predicted values or variables included in the final models (data not shown). For stand density, values of stand density were underpredicted at very low and very high stand densities (Figure 3f ).

Discussion
With the exception of stand density, the models developed in this study had a high degree of precision. Despite the scattered nature of the resource from which data were acquired the models appeared to be quite robust. The use of LiDAR provides a viable approach for characterising the resource in young stands where there is little prior information (such as those sampled in this study). The strength of relationships between LiDAR metrics and key stand characteristics were sufficiently high to delineate and prioritise areas for thinning. Thinning within the forests in the study area is typically undertaken when stands reach a mean top height of 14 m and the green crown reaches 3 m in height. Both mean top height and green crown height were predicted with a high degree of precision suggesting that LiDAR information could be useful for scheduling operations.
Of the key measurements used in forest inventory, mean top height is always predicted with the most precision using LiDAR. Although the coefficient of determination for the final mean top height model developed in this study was at the lower end of previously reported values, which range from 0.82 to 0.98 Coops et al. 2007;Means et al. 2000;Naesset 2002;Means et al. 1999;González-Ferreiro et al. 2012;Stone et al. 2011), the RMSE was similar to or lower than reported values for Douglas-fir and other coniferous species (Means et al. 2000;Means et al. 1999). Previous research found LiDAR-derived estimates of height for Douglas-fir to be considerably less precise than those estimated for Pinus ponderosa Douglas ex. C. Lawson (Andersen et al. 2006). LiDAR may predict height less precisely in Douglas-fir as the pronounced conical shape of the species is more likely to be missed by LiDAR than the wider crowns of other coniferous species.
Total stem volume was predicted with a high degree of precision. The coefficient of determination for the final model was within the range cited by previous Shown are the estimated coefficients for each term included in the six models. Model statistics shown include the root mean square error (RMSE), partial R 2 and significance category for each variable included in the model, with asterisks ***, **, * denoting significance at P = 0.001, P = 0.01 and 0.05 respectively. Also shown is the overall model R 2 .
studies for coniferous species that vary from 0.46 to 0.97 (Means et al. 2000;Naesset 1997Naesset , 2002. Compared to previous research, the RMSE of 20.2 m 3 ha -1 for the final model compares favourably to previous values that include 18.3 -31.9 m 3 ha -1 (Holmgren and Jonsson 2004), 26.1 -82.8 m 3 ha -1 (Naesset 1997), 71.6 m 3 ha -1  and 73 m 3 ha -1 (Means et al. 2000). Model predictions of basal area had a comparable precision to similar studies within coniferous forests in boreal and temperate regions. Coefficients of determination ranged from 0.62 to 0.94 for models predicting basal area in coniferous forests in the United States of America (Means et al. 1999), Norway (Naesset 2004b(Naesset , 2005(Naesset , 2002 and Denmark (Nord-Larsen and Schumacher 2012). LiDAR-derived predictions of basal area in young and old growth Douglas-fir stands had a coefficient of determination of 0.95 (Means et al. 2000) that exceeds the value reported in this study. This high R 2 may reflect the relatively large plot size (0.25 ha) and wide range in basal area covered by the dataset as increases in both factors often result in increases in the coefficient of determination (see Zhao et al., 2009 for details on effects of plot size on precision).
Although little previous research has developed plotlevel models of either stem diameter or green crown height from discrete-return LiDAR, both stand metrics were found to be strongly correlated to LiDAR metrics in the current study. Stem diameter was most strongly related to H 95 , which most likely acts as a surrogate for mean top height, reflecting the strong relationship between mean diameter and mean top height in the dataset ( Table 3). The strength of these relationships is broadly consistent with research on coniferous species within Norway where coefficients of determination ranged from 0.39 to 0.78 (Naesset 2002). Green crown height was most strongly related to H 01 which is likely to represent the transition between green needles that intercept light to dead or absent needles at the crown base that allow greater light transmission. The degree of precision of the final model of green crown height exceeds that of models previously developed for stands dominated by Picea abies (L.) Karst. (Norway spruce) and Pinus sylvestris L. (Scots pine) where the coefficient of determination ranged from 0.6 -0.71 (Naesset and Okland 2002).
Stand density was only predicted with a moderate degree of precision from LiDAR and stand level metrics. Previous research has also shown the correlation between LiDAR cloud metrics and stand density to be weaker than that for other stand characteristics reporting correlations that are only of moderate strength (Naesset and Bjerknes 2001;Hall et al. 2005;Naesset 2002).
Inclusion of stand age and stand density as predictive variables significantly improved a number of the developed models. Age featured in the final model for four of the six stand characteristics while stand density was a significant contributor to the models of mean diameter and basal area. These findings are broadly consistent with previous research that showed significant improvements in LiDAR models of stem volume through addition of stand density as an explanatory variable .
The results of this study suggest that an accurate stand density map would provide a useful source of ancillary information to support use of LiDAR in predictions of stand metrics. As direct estimation of stand density from LiDAR cloud metrics is typically only moderately precise, the use of individual tree delineation may provide the most accurate method for deriving a stand density map. Although there have been difficulties in precisely counting trees (Vauhkonen et al. 2012), recent research has shown that plantation stand density can be accurately estimated from LiDAR imagery through use of a calibrated tree detection methodology (D. Pont, Scion pers. comm.).
The very low percentage of ground returns recorded in this dataset provided a robust test of the utility of LiDAR data in predicting stand characteristics. Compared to previous research, the mean percentage of ground returns recorded in this study (3.82%) was markedly lower that the typical range found for other coniferous species (Naesset 2002;Naesset and Bjerknes 2001;Naesset 2004a) but was similar to that recorded for high density unthinned hinoki cypress (Chamaecyparis obtusa Sieb. et Zucc.) and sugi (Cryptomeria japonica D. Don) plantations (Takahashi et al. 2006). The low penetration rates in this study are likely to be attributable to the dense closed canopy of the unthinned plantations and the topographically complex terrain over which the study sites were located. In older thinned stands of Douglas-fir, predictive precision of stand characteristics is likely to be higher due to greater canopy penetration by LiDAR and consequently more accurate definition of the digital terrain model and canopy height model. Further research should be undertaken to investigate how well LiDAR data can be used to predict stand characteristics in older thinned stands of Douglas-fir.

Conclusion
This research demonstrates that strong relationships can be developed for mean top height, green crown height, volume, mean diameter and basal area from LiDAR and stand metrics. Although the final model developed for stand density was only of moderate strength, these results are consistent with previous research linking LiDAR cloud metrics to stand density. Further research should investigate how well stand characteristics are predicted at older ages from LiDAR. It would also be useful to compare the accuracy and cost efficacy of predictions from LiDAR data with standard inventory methods for the scheduling of production thinning operations and later more detailed inventories in stands of Douglas-fir.