Predictive Modeling of Black Spruce ( Picea mariana ( Mill . ) B . S . P . ) Wood Density Using Stand Structure Variables Derived from Airborne LiDAR Data in Boreal Forests of Ontario

Our objective was to model the average wood density in black spruce trees in representative stands across a boreal forest landscape based on relationships with predictor variables extracted from airborne light detection and ranging (LiDAR) point cloud data. Increment core samples were collected from dominant or co-dominant black spruce trees in a network of 400 m2 plots distributed among forest stands representing the full range of species composition and stand development across a 1,231,707 ha forest management unit in northeastern Ontario, Canada. Wood quality data were generated from optical microscopy, image analysis, X-ray densitometry and diffractometry as employed in SilviScanTM. Each increment core was associated with a set of field measurements at the plot level as well as a suite of LiDAR-derived variables calculated on a 20 × 20 m raster from a wall-to-wall coverage at a resolution of ~1 point m−2. We used a multiple linear regression approach to identify important predictor variables and describe relationships between stand structure and wood density for average black spruce trees in the stands we observed. A hierarchical classification model was then fitted using random forests to make spatial predictions of mean wood density for average trees in black spruce stands. The model explained 39 percent of the variance in the response variable, with an estimated root mean square error of 38.8 (kg·m−3). Among the predictor variables, P20 (second decile LiDAR height in m) and quadratic mean diameter were most important. Other predictors describing canopy depth and cover were of secondary importance and differed according to the modeling approach. LiDAR-derived variables appear to capture differences in stand structure that reflect different constraints on growth rates, determining the proportion of thin-walled earlywood cells in black spruce stems, and ultimately influencing the pattern of variation in important wood quality attributes such as wood density. A spatial characterization of variation in a desirable wood quality attribute, such as density, enhances the possibility for value chain optimization, which could allow the forest industry to be more competitive through efficient planning for black spruce management by including an indication of suitability for specific products as a modeled variable derived from standard inventory data.


Introduction
Recent estimates suggest that forests cover approximately 31% of the global land area; however, this area is declining at a rate of approximately 5.2 million ha•year −1 [1].Primary forests account for 36% of total global forest area but decreased by 40 million hectares from 2000 to 2010, largely as the result of reclassification following some form of human intervention [1].Sustainable forest management approaches have initiated a shift to increased harvest of shorter-rotation, second growth forests to satisfy wood supply demands while simultaneously supporting the conservation of old-growth, primary forests.This shift from primary to secondary forests has changed the composition and quality properties of the wood supply currently directed to mills [2].As older forest stocks are depleted, short-rotation, fast-growing crops will dominate the future wood supply [3].Trees harvested before maturity will be smaller and have greater proportions of juvenile wood, which is formed during periods of fast growth, and has inferior quality characteristics (e.g., lower density, shorter fibre length) and greater variability compared to wood formed at maturity [4,5].Thus, the global forest industry will have to operate at much greater efficiency to support increased demands for wood from a less extensive, lower quality resource.
Slow growth rates of boreal tree species provide a competitive advantage in global forest products markets based on the preferred wood quality characteristics [6].The species composition and growing conditions typical of the Canadian boreal forests provide a key and unique benefit (mainly related to fibre morphology and strength) in quality over volumetric production of the fibre resources globally.These advantageous characteristics of the Canadian fibre supply cannot be found anywhere else in the world [7].Canada is therefore ideally placed to embrace the value-chain concept as an alternative strategy to the volume-based paradigm [8], and meet the challenge of increasing demand from a shrinking global wood supply.
The value chain describes the path of activities and processes involved in development of a product; optimization of this value chain in forestry results in products gaining value at each stage by following a wood fibre usage strategy that matches fibre characteristics with their best possible uses at appropriate prices [8].More than half the cost of pulp production is related to acquisition and transport of wood [9], so any value chain optimization (VCO) scheme requires knowledge of the relationship between properties of the raw materials and yield.It is critical that decisions about allocation be made based on optimizing the use of the wood supply while simultaneously meeting criteria for ecological, economic and social sustainability [8].The advantage of taking informed fibre utilization decisions was demonstrated in Sweden by Duchesne et al. [10], by employing a field-based log sorting system derived from basic information on species, canopy position and stem position.This sorting captured key differences in wood supply characteristics (e.g., fibre length and basic density) that influenced paper quality parameters such as tensile and tear indices [10].Zhang [5] suggested that a broader wood quality concept emphasizes the totality of wood characteristics that affect the value chain.For example, wood is a critical component of the global carbon cycle and good estimates of carbon sequestration and storage in forests is another potential value that requires knowledge of key functional traits such as wood density [2,11].
Forest Resource Inventory (FRI) systems are critical tools that support forest management and planning activities across broad spatial scales [2].Advances in remote sensing technologies have improved FRI systems substantially in recent time; however, they still lack information on the characteristics of wood that determine quality, which is critical to achieving VCO goals.Information on wood quality attributes of standing trees prior to harvest would provide significant benefits related to planning harvests and marketing wood products, and could also provide insights into a variety of ecological processes (e.g., wood decay rates could be estimated based on density) that facilitate sustainable management decisions [12].The efficient collection of wood quality information has been made possible by SilviScan technology [13], creating the potential for statistical modeling and prediction of wood and fibre quality on the landscape [2].
Remote sensing data may provide the potential to capture variation in tree growth and wood quality over broad spatial and temporal scales [2].There is extensive literature devoted to description of direct relationships between intrinsic wood quality attributes and extrinsic factors, from both ecological [11] and management [2] perspectives that cover the tree and stand-level.However, very few studies have explored the scaling-up of these direct relationships to understand patterns of variation in wood properties at the landscape-level.Recently, Giroud et al. [14] were able to develop a two-stage model that provided useful estimates of stand-level variation in various wood and fibre attributes, suggesting the possibility of large-scale mapping of black spruce wood quality.However, their model required extensive field data and they also indicated that performance could be improved through the use of LiDAR-derived stand metrics as predictor variables.
Even if some reduction in the accuracy of wood quality predictions occurs when the predictor variables are derived from remote sensing data rather than direct field measurements, the corresponding capacity to scale-up the predictions by using landscape-level source data balances this by creating the potential to characterize wood quality attributes of the resource over a broad spatial scale [12].Van Leeuwen et al. [2] developed a strong case for the use of LiDAR-derived extrinsic indicators of wood quality to support broad scale modeling and mapping of key wood and fibre attributes.Their approach emphasized the use of descriptive indicators such as height, diameter and crown dimensions of trees.They also acknowledged the potential to develop models based on mechanistic indicators (e.g., nutrients, physiography, precipitation and temperature).Recent studies have suggested the potential for using LiDAR-derived predictors in various modeling schemes designed to provide useful estimates of wood quality attributes [12,15].
More precise forest inventory information provides an exciting opportunity for the development of new modeling tools for forest management.Such tools could be used to optimize both growth and value of forest production, and make the Canadian forest industry more competitive in the global market.In addition to traditional volumetric growth and yield modeling, there is a growing interest to develop tools that predict wood quality attributes in standing timber.The main attributes of interest include density, microfibril angle (MFA), modulus of elasticity (MOE), cell wall thickness and fibre length, and estimates of these have been derived based on relationships to easily measurable tree, site and stand level ecological variables [16,17].To address the demand for mapping average wood fibre properties of standing trees, we have explored the potential of using LiDAR derived variables to predict average stem-level fibre attributes in specific stand types of black spruce in the Hearst Forest of Ontario, Canada.The main objectives of the study were to (a) develop relationships between wood quality attributes (e.g., wood density) and metrics derived from remotely sensed LiDAR point cloud data for black spruce to predict average tree-level fibre quality within Enhanced Forest Inventory grid cells (20 m × 20 m); and (b) develop a predictive raster map (20 m × 20 m) of wood density for average black spruce trees in the stands across the forest landscape.

Study Area
Our study was completed in the Hearst Forest (HF), within the boreal forest region of Ontario, Canada (Figure 1).The Hearst forest management unit is an area of 1,231,707 ha and contains elements representative of the boreal forest region [18,19].Most of the area within HF includes part of the northern Claybelt, which is characterized by flat to rolling terrain (elevation 86-435 m) underlain by poorly drained glacial-lacustrine clay deposits.Across the forest landscape, there are also minor deposits (ranging from sand to clay materials) derived from various glacial processes.Climate data from the nearest weather station at Kapuskasing, Ontario, Canada (49.41 • N, 82.50 • W) indicate that the annual mean daily temperature is 0.7 • C, ranging from −18.7 • C in the coldest month (January) to 17.2 • C in the warmest month (July).Total annual precipitation is 831.8 mm, with just over 1/3 falling as snow [19,20].

Sample Plot Networks
The plot network we utilized in the HF consisted of 426 temporary sample plots (TSPs) established in the boreal forest region of northeastern Ontario during 2010.These plots were originally established to provide calibration data for the development of stand-level forest mensuration models from Light Detection and Ranging (LiDAR) data [21].They were circular 400m 2 plots located in a stratified random sampling design intended to cover the full range of variation in species composition, stand age and development, and ecological site conditions across the forest management unit.The scope of this study covered only situations where black spruce forms an important component of the forest, and required very intensive analysis of individual core samples.

Sample Plot Networks
The plot network we utilized in the HF consisted of 426 temporary sample plots (TSPs) established in the boreal forest region of northeastern Ontario during 2010.These plots were originally established to provide calibration data for the development of stand-level forest mensuration models from Light Detection and Ranging (LiDAR) data [21].They were circular 400-m 2 plots located in a stratified random sampling design intended to cover the full range of variation in species composition, stand age and development, and ecological site conditions across the forest management unit.The scope of this study covered only situations where black spruce forms an important component of the forest, and required very intensive analysis of individual core samples.Therefore, we collected samples from a subset of TSPs (n = 75) selected to represent the range of site conditions occupied by black spruce on the HF (Figure 1).The details of the sampling design were discussed by Pokharel et al. [16].

Data on Wood Quality Characteristics
Large (12 mm diameter) increment cores were collected at breast height (1.3 m) from co-dominant black spruce individuals selected at random from the sample population in a given plot.For trees to be eligible for sampling, they had to be free of any visible signs of stress or disturbance relative to the overall population on the plot and composed of solid wood, free of decay in the stem where the core was extracted.We collected and analyzed a total of 111 samples.These samples were initially stored in a freezer, then prepared, processed and analyzed using standard protocols for analysis of wood and fibre attributes using SilviScan™ technology at FPInnovations wood products facility in Vancouver, Canada.Details of the laboratory data collection procedure and derivation of fibre attributes were described by Defo and Uy [22].In brief, the air-dried samples (8 percent moisture content) were cut into strips of 2 mm thickness (tangential) and 7 mm height (longitudinal) using a twin-bladed saw.Strips were sanded with 400 grit and finished with 1200 grit paper.Each sample strip was scanned for radial and tangential fibre dimensions (cell scanner), density (X-ray densitometry), and micro-fibril angle (X-ray difractometry) from pith to bark.Fibre dimension and density were determined at 25 µm resolution, whereas the micro-fibril angle was determined at 5 mm resolution.SilviScan analysis provides a number of wood and fibre attributes including density (kg•m −3 ); microfibril angle; modulus of elasticity (GPa); coarseness (µg•m −1 ); wall thickness (µm); radial diameter (µm); tangential diameter (µm) and specific surface area (m 2 •kg −1 ).A stand, site, tree, ring and sub-ring level hierarchical data-base was created for the SilviScan data in Microsoft Access.Data processing and statistical tests were performed in the R statistical computing environment [23].We sub-set the data to capture representative samples of the wide range of growing conditions and age structure for black spruce in the Hearst Forest.Our chronologies ranged in duration from 25 to 163 years.We estimated the average stem-level basal area-weighted wood density from annual ring-level data for each sample.

LiDAR and FRI Data set
A full suite of LiDAR-derived variables originating from wall-to-wall coverage of the Hearst Forest was made available for the analysis at the plot level for each of 75 black spruce plots.The LiDAR data were collected between 4 July and 4 September 2007.The coverage area of 15,740 km 2 was flown with a Leica ALS50 sensor mounted on a Cessna 310 aircraft.The sensor had a 30 • field of view, with a pulse rate of 119,000 Hz and a scan rate of 32 Hz.The resulting data coverage included 20% overlap, and had an average density of 0.82 points m −2 .The vendor processed the data to ASPRS format using TerraScan [24].The LiDAR data were further processed to trim the variation created by large veteran trees in young stands by removing the top 5% of the point cloud when calculating measures of spread or vertical complexity [25].The LiDAR predictor variables presented here were generated as described in Woods et al. [26].These variables were estimated on a 20 m × 20 m raster for the entire HF, and we extracted the values for the variables of interest associated with the location of each plot in the data set.A description of the variables derived from the LiDAR point cloud that had some degree of correlation with wood fibre properties is provided in Table 1.Some of the variables we used in the analysis were derived from direct measurement in the plots, such as leaf area index (LAI) quadratic mean diameter (QMD) and basal area (BA).LAI estimates were derived based on the method as described by Pope and Treitz [27].

Correlations
To gain some insight into the important variables that might be associated with variation in wood quality, we ran a simple Spearman correlation analysis between the response and predictor variables in the data set.The nonparametric Spearman approach was used to avoid the potential influence of unusual cases (extreme values) or non-normal distributions on the description of associations.All coefficients with a p-value ≤ 0.05 were considered significant.

Parametric Approach
A multiple regression model was fitted using all potential LiDAR-derived variables and wood density as a response variable.We used backward selection in the base package of the R statistical computing environment.We estimated the variance inflation factor and the majority of these variables had values greater than 10, indicating that the model was overfit.We further eliminated one variable at a time, then evaluated the significance of the remaining variables in the model and assessed their variance inflation factor simultaneously.A model with the most parsimonious set of variables was produced and its fit was evaluated.

Nonparametric Approach
Our preliminary correlation analysis suggested that there could be multicollinearity among the LiDAR derived predictor variables; therefore, we further explored the potential of using nonparametric approaches (i.e., regression tree and random forests) for predicting wood quality characteristics such as density for black spruce.Nonparametric techniques such as classification and regression trees or random forests have gained popularity in ecological modeling because these approaches can generalize model predictability to a scale appropriate to management applications, and they can make use of a large number of predictor variables without over-fitting [29,30].These models require conceptually meaningful explanatory variables that are measured across the entire range of conditions over which model projection is desired, otherwise the resulting model will not represent any condition not identified in the fitting data.In the case of forest growth and yield modeling, many predictor variables can be generated easily at the stand-level from various spatial data sources including the FRI and LiDAR coverage.
Classification and regression trees provide an output tree that is easy to interpret and apply within a GIS environment, allowing great flexibility in the choice of explanatory variables (even with complex interactions or multicollinearity) and generally produce high classification accuracy [29,30].A regression tree is built based on a binary partitioning algorithm which successively splits a response variable (e.g., wood density) into mutually exclusive homogeneous groups.For each split, the algorithm considers each candidate explanatory variable and selects the one that reduces the greatest residual sum of squares.Eventually, a tree is built by successively repeating the splitting process until the user-defined parameters, (typically the critical cross-validation error value and number of observations at each terminal node) are satisfied.Regression trees provide a set of classification rules derived from a dendrogram, estimates of means for the response variable at each node, and a coefficient of determination to assess the explanatory value of the model.We fit a regression tree using our sample density data and suite of predictor variables; however, the analysis was limited to a data set of 75 plots and the general applicability of the RT model was unknown.
We employed the random forests (RF) [31] approach to overcome some of the limitations of regression tree analysis, including high sensitivity of the output tree to small changes in the fitting data set.One of the objectives of this project was to identify important LiDAR-derived variables while predicting wood quality attributes such as density for black spruce in the boreal forest region of Ontario; therefore, we identified the importance of variables based on our RF simulation.RF is a machine learning approach that we employed to generate five thousand regression trees from the fitting data set through a boot-strapped sampling procedure (selecting cases and predictor variables for each tree at random).The technique uses ensemble strategies to generate predictions for the out-of-bag fitting data (i.e., cases not selected for a given tree) and the targeted population (i.e., the forest area the model will be applied to).RF was chosen as the approach because it offers numerous advantages including computational efficiency, the requirement for only a few tuning parameters (e.g., number of randomly selected variables, number trees, and tree size), clear estimates of the reliability of model prediction based on out-of bag observations, the potential to predict missing values, the ability to handle thousands of predictor variables, the lack of concern related to over-fitting of the model and a quantitative listing of predictor variable importance in the classification/regression process [30][31][32][33].
The bootstrapping, bagging and ensembling make RF predictions robust, unbiased, and generalizable across the population.

Correlation Analysis
Many of the wood quality variables were associated with wood density with relatively strong correlation coefficients (e.g., ≥ ±0.48).Furthermore, our preliminary analysis indicated that variables such as cell wall thickness (WT) and specific surface area (SSA) were highly correlated (r > 0.87) with wood density, and had similar model fit statistics as wood density.Alternatively, other wood quality variables such as microfibril angle (MFA), modulus of elasticity (MOE), radial diameter (RD) and tangential diameter (TD) exhibited weak relationships with predictor variables.Thus, we focused on wood density as the single response variable that captured the most variance in wood quality that could be modelled using LiDAR-derived predictors.Several of the predictor variables such as QMD, D1, D2, CC2, CC6, P20, P60, VCI, LAI, DA, top height and basal area all had Spearman correlation coefficients to stem-level average wood density greater than ±0.30(Table 1).Most of the associations were positive, with the exception of D1, D2 and first returns (Figure 2).fitting of the model and a quantitative listing of predictor variable importance in the classification/regression process [30][31][32][33].The bootstrapping, bagging and ensembling make RF predictions robust, unbiased, and generalizable across the population.

Correlation Analysis
Many of the wood quality variables were associated with wood density with relatively strong correlation coefficients (e.g., ≥ ±0.48).Furthermore, our preliminary analysis indicated that variables such as cell wall thickness (WT) and specific surface area (SSA) were highly correlated (r > 0.87) with wood density, and had similar model fit statistics as wood density.Alternatively, other wood quality variables such as microfibril angle (MFA), modulus of elasticity (MOE), radial diameter (RD) and tangential diameter (TD) exhibited weak relationships with predictor variables.Thus, we focused on wood density as the single response variable that captured the most variance in wood quality that could be modelled using LiDAR-derived predictors.Several of the predictor variables such as QMD, D1, D2, CC2, CC6, P20, P60, VCI, LAI, DA, top height and basal area all had Spearman correlation coefficients to stem-level average wood density greater than ±0.30(Table 1).Most of the associations were positive, with the exception of D1, D2 and first returns (Figure 2).

Parametric Model
Although many predictor variables were correlated with wood density, they were also correlated with one another, leading us to suspect there would be multicollinearity among the predictor variables (variance inflation factor-VIF > 10) when we fit the multiple regression model.For instance, top height was one of the variables that entered the model of wood density, but it was also highly correlated with QMD and VCI (r > 0.8; VIF > 24).In such cases, we chose the variable to enter the model from the group of correlated predictors that explained the most variance, which in this case was QMD, as it explained the most variance compared to the top height or VCI.
The coefficient of determination (r 2 ) for the multiple regression model of wood density was moderate (0.48).The most parsimonious model is presented in Equation ( 1), and includes QMD, LAI and P20 as predictors (Table 2).This parsimonious model had variance inflation factor values of less than 2 for all predictor variables.A residuals plot and Q-Q normal plot showed that the model residuals were distributed normally and homogenously with constant variance, validating the assumptions related to parametric model.
where DENSITY = wood density in kg•m −3 , QMD is quadratic mean diameter (cm), LAI is leaf area index and P20 is second decile LiDAR height (m).

Parametric Model
Although many predictor variables were correlated with wood density, they were also correlated with one another, leading us to suspect there would be multicollinearity among the predictor variables (variance inflation factor-VIF > 10) when we fit the multiple regression model.For instance, top height was one of the variables that entered the model of wood density, but it was also highly correlated with QMD and VCI (r > 0.8; VIF > 24).In such cases, we chose the variable to enter the model from the group of correlated predictors that explained the most variance, which in this case was QMD, as it explained the most variance compared to the top height or VCI.The coefficient of determination (r 2 ) for the multiple regression model of wood density was moderate (0.48).The most parsimonious model is presented in Equation ( 1), and includes QMD, LAI and P20 as predictors (Table 2).This parsimonious model had variance inflation factor values of less than 2 for all predictor variables.A residuals plot and Q-Q normal plot showed that the model residuals were distributed normally and homogenously with constant variance, validating the assumptions related to parametric model.
where DENSITY = wood density in kg•m −3 , QMD is quadratic mean diameter (cm), LAI is leaf area index and P20 is second decile LiDAR height (m).

Regression Tree and Random Forests Modeling
The regression tree analysis identified P20 (second decile LiDAR height in m) as the most important predictor variable while classifying wood density for black spruce.Regression trees split the data set into two classes that minimize internal heterogeneity and provide a predicted mean for each class, as well as a set of rules to map the class on the landscape.Notably, the first split in the regression tree, which explained the most variation in wood density, separated out the plots with low wood density from high wood density plots (Figure 3), based on whether the value of P20 was ≥1.08.The regression tree further split plots into two sub-classes with a cut-off quadratic mean diameter of ≥12.6 cm.The entire regression tree produced 11 terminal nodes that separated along a gradient of wood density, and explained 75% of the total variance of wood density in the sample population.Such a regression tree could be easily implemented in a GIS to provide average estimates of wood fibre attributes for a given polygon or even at the raster level.For example, an average co-dominant tree in a given polygon or raster with less than 1.08 m of second decile LiDAR height and in a stand where the QMD is ≤12.63 cm would have a predicted average wood density of 614.7 kg•m −3 (Figure 3).

Regression Tree and Random Forests Modeling
The regression tree analysis identified P20 (second decile LiDAR height in m) as the most important predictor variable while classifying wood density for black spruce.Regression trees split the data set into two classes that minimize internal heterogeneity and provide a predicted mean for each class, as well as a set of rules to map the class on the landscape.Notably, the first split in the regression tree, which explained the most variation in wood density, separated out the plots with low wood density from high wood density plots (Figure 3), based on whether the value of P20 was ≥1.08.The regression tree further split plots into two sub-classes with a cut-off quadratic mean diameter of ≥12.6 cm.The entire regression tree produced 11 terminal nodes that separated along a gradient of wood density, and explained 75% of the total variance of wood density in the sample population.Such a regression tree could be easily implemented in a GIS to provide average estimates of wood fibre attributes for a given polygon or even at the raster level.For example, an average co-dominant tree in a given polygon or raster with less than 1.08 m of second decile LiDAR height and in a stand where the QMD is ≤12.63 cm would have a predicted average wood density of 614.7 kg•m −3 (Figure 3).Even though a regression tree would be easily implemented in a GIS to provide average estimates of wood fibre attributes for a given polygon, the classification is sensitive to small changes in the fitting data set.To overcome this limitation, we applied RF to the same set of fitting data.The variable importance matrix generated by random forests indicated that P20, QMD, CC4, D1, CC2, D2, TPH, P60, and TOPHT were the most important variables (reduction in MSE at least by 20 percent for each variable) for classifying wood density at the landscape level (Figure 4).The random forests model for average stem-level wood density explained 39 percent of variability in the sample population with an estimated root mean square error of 38.8 kg•m −3 .The model under-predicted with respect to high density wood and over-predicted with respect to low density wood (Figure 5).
Even though a regression tree would be easily implemented in a GIS to provide average estimates of wood fibre attributes for a given polygon, the classification is sensitive to small changes in the fitting data set.To overcome this limitation, we applied RF to the same set of fitting data.The variable importance matrix generated by random forests indicated that P20, QMD, CC4, D1, CC2, D2, TPH, P60, and TOPHT were the most important variables (reduction in MSE at least by 20 percent for each variable) for classifying wood density at the landscape level (Figure 4).The random forests model for average stem-level wood density explained 39 percent of variability in the sample population with an estimated root mean square error of 38.8 kg•m −3 .The model under-predicted with respect to high density wood and over-predicted with respect to low density wood (Figure 5).

Discussion
Despite differences in approach, predictor variables and geographic region, our model performance (r 2 = 0.39; RMSE = 38.8kg•m −3 ) compared well to other models of wood fibre attributes that made use of LiDAR-derived predictor variables.In one of the first studies of airborne LiDAR metrics and their relation to wood quality, Hilker et al. [12] extracted two components of wood

Discussion
Despite differences in approach, predictor variables and geographic region, our model performance (r 2 = 0.39; RMSE = 38.8kg•m −3 ) compared well to other models of wood fibre attributes that made use of LiDAR-derived predictor variables.In one of the first studies of airborne LiDAR metrics and their relation to wood quality, Hilker et al. [12] extracted two components of wood quality from a principal components analysis of several lodgepole pine wood fibre attributes, and produced a model with LiDAR-derived predictors that explained just over half of the variance (r 2 = 0.55) in the component they described as wood strength (which was related to density).In a more direct comparison to our study, Luther et al. [15] modelled black spruce wood density using multiple linear regression with airborne LiDAR-derived predictors in Newfoundland, and produced results that were similar (r 2 0.40, RMSE = 30.9kg•m −3 ) to ours.Blanchette et al. [34] collected terrestrial LiDAR data and the extracted canopy metrics from a higher resolution point cloud, which resulted in an improvement in performance of models of black spruce density in Newfoundland (r 2 = 0.63; RMSE = 23.87 kg•m −3 ).Our results from the nonparametric random forests modeling approach in Ontario black spruce make a different type of prediction, essentially describing the density of a typical or average tree that is representative of specific conditions in the landscape (Figure 6).Nevertheless, it is interesting that the model emerging from the ensemble strategy we employed produced fit statistics that compare favourably to the models derived from more traditional parametric approaches.
The consistent performance of models that attempt to predict wood fibre attributes from LiDAR-derived predictors using different systems and approaches likely arises from the consistency in the influence of the specific drivers of wood quality that LiDAR metrics are able to capture.Van Leeuwen et al. [2] identified tree height and growth, crown dimensions, vertical foliage distribution, stocking and branchiness as the characteristics of forest canopies that could be extracted from LiDAR data and related to various wood fibre characteristics.In our study, the main predictors of importance among the various modeling approaches were consistently P20 (second decile LiDAR height in m) and quadratic mean diameter.Other predictors such as D1, D2, P10, P30 and LAI were of secondary importance and differed according to the modeling approach.One group of these predictors (P20, D1, D2, P10, P30) provides different quantities that represent the proportion of returns that occur in the lower bins of the canopy profile (a higher proportion representing a canopy of greater depth), which had an inverse relationship to wood density.The other group of predictors (QMD, BA) represents the level of competition on the site, which also had a negative relationship to density.Thus, our models reflect the negative influence of canopy depth and high site occupancy on black spruce wood density.For lodgepole pine forests in Alberta, Hilker et al. [12] found that the relative thickness of the upper canopy layer was the only significant LiDAR-derived predictor of the wood strength principal component they extracted; however, they also noted a strong relationship between average diameter increment and the strength component, suggesting an important competition effect.Luther et al. [15] identified the percentage of all returns above ground (a forest cover metric), the ratio of the canopy height model area to ground area (a canopy complexity measure) and the maximum value of the bare earth terrain surface model (a terrain complexity metric) as the primary predictors of black spruce wood density; however, these predictors were selected for model inclusion based on parsimony, and it could be that other predictors would also provide similar performance.Nevertheless, the canopy and forest cover metrics suggest similar relationships between complexity, occupancy and wood density as revealed in our analyses.Furthermore, Blanchette et al. [34] identified the nearest-neighbour ratio (an index of occupancy/competition) as a strong predictor of density in their models based on terrestrial LiDAR data.Therefore, it appears that the models we produced identified predictor variables that are consistent with other approaches and systems that have the common thread of using LiDAR data to quantify relationships between canopy characteristics and wood quality attributes.A large amount of unexplained variability still exists in our model.This is due to the fact that there is a large variation in wood properties between individual trees.These variations could be related to differences in age, site, genotypes, and management [35,36].There is substantial variation A large amount of unexplained variability still exists in our model.This is due to the fact that there is a large variation in wood properties between individual trees.These variations could be related to differences in age, site, genotypes, and management [35,36].There is substantial variation in fibre characteristics along the chronology from pith to bark.In our previous study of ecosite-based variation in wood fibre properties, we identified earlywood:latewood proportions as an effective response variable for capturing the pith to bark signal caused by variable growing conditions.Therefore, in black spruce, variation in latewood percentage could be largely influencing age-related changes in wood properties such as density, based on the transition from the formation of juvenile to mature wood.In addition to age, there are other sources of unexplained variability that are largely dependent on site, genotypes and social position of the tree in the stand.We have already explored site effects based on ecosite classification; however, some structural variables may also be useful in identifying site quality, such as top height and LAI.The age differences in our data have obscured the ability to detect the differences between sites, as site quality and age in the Hearst Forest are confounded (better growing sites tend to be younger).In our previous ecosite-based analyses, we eliminated this confounding effect by truncating our response variables at a common 50-year index value; however, age is primary driver of wood quality attributes and varies considerably over the landscape.The LiDAR-based approach could provide a key element of predictive wood quality models if it can capture age effects through structural variables such as VCI and top height.
Very few forest growth and yield models have integrated wood quality information into their modeling framework.Furthermore, models incorporating wood quality are typically temporal, with the objective of predicting the future characteristics of trees or stands derived from a set of initial measurements, such as MELA (Metsälaskelma), an individual tree model developed by the Finnish Forest Research Institute [37], and the "TreeRing" model, which generates predictions on a daily time step [38].The incorporation of value characteristics into forest inventory systems, however, requires large-scale spatial predictive modeling.Since the spatial distribution of wood properties is shaped by vertical structure and environmental factors, it is necessary to explore the relationship among these variables and to select the modeling approach that could best capture the complex relationships among wood quality attributes and broad-scale drivers of tree growth.
The majority of research on wood quality modeling has focused on plantation grown species such as radiata pine in New Zealand [39,40] or loblolly pine in the United States [41,42].These modeling efforts have been largely concentrated at the ring or stem level due the fact that there is a large variation in wood quality characteristics within and between individual trees [35,36].Age is one of the variables that explains the most ring-to-ring variation in wood quality characteristics and has been widely used.Models projected at the ring-level that use age as one of the covariates are useful tools to test research hypotheses but have little or no utility for management of natural stands or developing a predictive landscape-level map to enhance value-based management planning and product optimization, unless age can be accurately predicted from remote-sensing data [43].
In this study, we limited our analysis to the use of only LiDAR-derived variables to predict wood density.In addition to LiDAR-derived variables, there is a potential opportunity to characterize fibre properties based on site (i.e., ecosite); individual tree level characteristics such as crown height, crown width and crown volume; FRI-based variables such as species composition, stand history and management; and spatial variables such as moisture availability slope, aspect, elevation and latitude.In our earlier study [16], ecosite was one of the most important predictor variables in regression tree and RF analyses, which reduced a substantial amount of the residual sum of squares for predicting wood density as a response variable.We could possibly improve the ultimate version of a predictive wood quality model by combining LiDAR, ecosite, FRI and spatial data in order to predict key wood quality characteristics.
Enhanced forest resource inventories (eFRI) have demonstrated the potential of LiDAR data [26,44] and individual tree classification (ITC) approaches to provide extremely accurate spatial inventory information.Individual tree classification along with characterization of crowns would increase the accuracy and reliability of predicting wood quality attributes at the landscape level.A multivariate approach was used elsewhere to employ LiDAR-derived variables to quantify wood properties and scale them up to the stand level [2,12,45].Our approach here is unique and has demonstrated some of the correlation between wood quality attributes and some of the easily obtainable stand level attributes that could be derived from LiDAR data.

Conclusions
Our analysis using SilviScan data evaluated the relationships between an important fibre attribute (wood density) and LiDAR-derived variables using both parametric and nonparametric approaches.There is a strong correlation between wood density and LiDAR-derived variables; however, there also exists a multicollinearity among the predictor variables derived from LiDAR.Random forests analysis identified that P20, QMD, D1, TOPHT, LAI and some other attributes that characterize the crown cover were the important variables for predicting wood density.Given the nature of our data set, in which we have quantified fibre properties at the tree-level and are relating them to LiDAR-derived plot-level metrics, there may be some noise in the predictions, which cannot be removed without using additional predictors.We foresee an opportunity to use LiDAR-derived variables in conjunction with other variables such as ecosite, moisture/wetness index, slope, aspect, elevation, growing degree days, mean annual temperature, mean annual precipitation and latitude to produce a more broadly applicable model.The RF approach provides the potential utility of mapping the average fibre attributes of individual co-dominant black spruce across the forest landscape.In order to increase the reliability and build confidence in the model for the end users, a continuous grid map of wood quality attributes validated using independent set of data is needed.

Figure 1 .
Figure 1.Sample plot location within the Hearst Forest in Ontario, Canada.

Figure 1 .
Figure 1.Sample plot location within the Hearst Forest in Ontario, Canada.

Figure 3 .
Figure 3. Regression tree analysis to predict the wood density of an average black spruce tree in a given stand using LiDAR-derived predictor variables for a sample tree population representing black spruce in the boreal forest of northeastern Ontario, Canada.Note: for each decision node, if true, move to the left.Variables used in regression trees are: HWD.P = hardwood proportion in the stand; MEDHT = median height in m; QMD = quadratic mean diameter in cm; P20 = second decile LiDAR height in m; DA = first

Figure 3 .
Figure 3. Regression tree analysis to predict the wood density of an average black spruce tree in a given stand using LiDAR-derived predictor variables for a sample tree population representing black spruce in the boreal forest of northeastern Ontario, Canada.Note: for each decision node, if true, move to the left.Variables used in regression trees are: HWD.P = hardwood proportion in the stand; MEDHT = median height in m; QMD = quadratic mean diameter in cm; P20 = second decile LiDAR height in m; DA = first returns/all returns (%); CC12 = crown closure > 12 m (%); D2 = cumulative % of the number of returns found in bin 2 of 10 (%); D6 = cumulative % of the number of returns found in bin 6 of 10 (%).

Figure 4 .
Figure 4. Variable importance plot of LiDAR-derived variables generated from the random forests model.%IncMSE = Percent increase in mean-squared error when the variable is removed from the model.

Figure 4 .
Figure 4. Variable importance plot of LiDAR-derived variables generated from the random forests model.%IncMSE = Percent increase in mean-squared error when the variable is removed from the model.

Forests 2016, 7 , 311 12 of 18 Figure 5 .
Figure 5. Actual versus predicted wood density in the fitting data set.Dashed line is the line of best fit; solid line represents a 1:1 correspondence.

Figure 5 .
Figure 5. Actual versus predicted wood density in the fitting data set.Dashed line is the line of best fit; solid line represents a 1:1 correspondence.

Figure 6 .
Figure 6.Predictive raster map (20 m × 20 m) of average wood density (kg⋅m −3 ) of black spruce trees in representative stands across the Hearst Forest, Ontario, Canada.The upper panel represents the total 1.3 million ha area within the Hearst Forest, the lower panel is an enlarged portion of the predicted density raster.White-colored polygons are non-forested area.

Figure 6 .
Figure 6.Predictive raster map (20 m × 20 m) of average wood density (kg•m −3 ) of black spruce trees in representative stands across the Hearst Forest, Ontario, Canada.The upper panel represents the total 1.3 million ha area within the Hearst Forest, the lower panel is an enlarged portion of the predicted density raster.White-colored polygons are non-forested area.

Table 1 .
Description of selected wood quality attributes and LiDAR derived variables.Correlation was estimated between the variable of interest and wood density.

Table 2 .
Parameter estimates and variance inflation factor of the parametric multiple regression model (model 1) fit using some selected predictor variables estimated from airbore LiDAR data.

Table 2 .
Parameter estimates and variance inflation factor of the parametric multiple regression model (model 1) fit using some selected predictor variables estimated from airbore LiDAR data.