An Alternative Approach to Using LiDAR Remote Sensing Data to Predict Stem Diameter Distributions across a Temperate Forest Landscape

We apply a spatially-implicit, allometry-based modelling approach to predict stem diameter distributions (SDDs) from low density airborne LiDAR data in a heterogeneous, temperate forest in Ontario, Canada. Using a recently published algorithm that relates the density, size, and species of individual trees to the height distribution of first returns, we estimated parameters that succinctly describe SDDs that are most consistent with each 0.25-ha LiDAR tile across a 30,000 ha forest landscape. Tests with independent validation plots showed that the diameter distribution of stems was predicted with reasonable accuracy in most cases (half of validation plots had R2 ≥ 0.75, and another 23% had 0.5 ≤ R2 < 0.75). The predicted frequency of larger stems was much better than that of small stems (8 ≤ x < 11 cm diameter), particularly small conifers. We used the predicted SDDs to calculate aboveground carbon density (ACD; RMSE = 21.4 Mg C/ha), quadratic mean diameter (RMSE = 3.64 cm), basal area (RMSE = 6.99 m2/ha) and stem number (RMSE = 272 stems/ha). The accuracy of our predictions compared favorably with previous studies that have generally been undertaken in simpler conifer-dominated forest types. We demonstrate the utility of our results to spatial forest management planning by mapping SDDs, the proportion of broadleaves, and ACD at a 0.25 ha resolution.


Introduction
Remote sensing data provide landscape-level information on forest resources, and can be used to inform spatial forest management planning.A number of statistical approaches have been used for estimating stand attributes from airborne LiDAR data in particular [1], and these have subsequently been used to produce detailed maps showing how forest characteristics vary across a landscape.Aboveground carbon density or biomass is commonly estimated [2][3][4] because of its importance to carbon accounting [5].Basal area, stem density and stem volume are often estimated as well for forest management planning [6][7][8][9].Despite the demonstrated utility of airborne LiDAR for mapping several types of stand characteristics, its potential for predicting stem diameter distributions (SDDs-the number of stems per ha across a set of diameter classes) has not yet been fully realized.SDDs represent the fundamental recordings from which other characteristics can be calculated, including basal area, stem density, quadratic mean diameter and (when coupled with wood density information) aboveground carbon density.SDDs also reflect the cumulative effects of growth, disturbance, and regeneration through time, and therefore represent an important indicator of past (and future) stand dynamics and are critical for developing management strategies [10].
Current methods for predicting SDDs from LiDAR include area-based approaches (ABA), which operate at the plot level, and individual tree detection (ITD), which operates at the tree level [11].Parametric ABAs use summary metrics derived from the LiDAR point cloud to predict the parameters of a SDD-usually a Weibull distribution [12]-via a regression model [13][14][15][16].These parameters are often combined with LiDAR-based predictions of basal area, total stem number, or total plot volume so that the number of stems of different sizes can be scaled to match predictions [13,14,16].Alternatively, non-parametric ABAs use nearest neighbor methods to estimate plot attributes as a weighted average of plots with similar LiDAR height return characteristics [17][18][19].Non-parametric ABAs tend not to generalize well to new sites because they are sensitive to within-stand covariance between structural variables and LiDAR acquisition properties [20].
ITD methods predict stem diameters from individual tree crowns that are delineated from the point cloud [21].Stem diameters are estimated from regressions against tree height and crown area [22][23][24], or from other metrics derived from the delineated point cloud [25].Although ITD offers a more mechanistic approach than ABA, its application is limited because ITD requires high density data (>10 points/m 2 versus point densities as low as 0.5 pulses/m 2 for ABA; [26]), it is less effective at capturing broadleaf crowns than coniferous ones [27], it does not detect all trees beneath the canopy [22], and it carries high computational demands.Some of these issues can be addressed through hybrid approaches that combine ITD with ABA [28].For example, the distribution of stem diameters can be extrapolated from an ITD model to encompass smaller stems [29], or diameters can be calibrated to match an ABA-predicted stand metric [30].
A third approach to estimating stand structure uses allometric relationships between tree height and crown size to model the canopy envelope that generates a given LiDAR point cloud [31].This type of model-based approach accounts for how tree size, crown shape, crown overlap, crown permeability, and canopy gaps jointly affect the height distribution of LiDAR return heights [32].It does not infer the locations of individual tree crowns, but does estimate tree size distributions for a localized area [33].The approach benefits from a mechanistic footing, like ITD, whilst maintaining the simpler, spatially-implicit features of an ABA.To date, however, these models have not been widely applied or tested for different forest types.
Our aim was to develop a method for estimating SDDs that has a mechanistic footing, like the ITD approach, whilst maintaining the simpler, spatially-implicit features and lower data demands of an ABA.We previously developed an allometry-based algorithm that describes how the height distribution of LiDAR first returns reflects the size, arrangement, and permeability of individual tree crowns in a mixed temperate forest [32].Here, we couple this algorithm with Bayesian model-fitting techniques to infer the SDDs and proportion of broadleaved species that produced observed LiDAR returns in 0.25 ha tiles across a 30,000 ha landscape.Our objectives were: (1) to evaluate the ability of an allometry-based modelling approach to predict SDDs and stand structural attributes in uneven-aged mixed hardwood stands; and (2) to map variation in stand structure and carbon density across a temperate forest landscape using the predicted stem diameter distributions.

Study Area
The data used in our study include airborne LiDAR and ground plots from Haliburton Forest and Wildlife Reserve in central Ontario, Canada (Figure 1).This area is a 30,000-ha mixed temperate forest dominated by Acer saccharum (sugar maple), with a significant coniferous component, including Tsuga canadensis (eastern hemlock) and Abies balsamea (balsam fir).There are some lowland regions with approximately 50 lakes and 250 wetlands [34], but the majority of the landscape is characterized by undulating, rough topography typical of Precambian Shield landscapes [35].The forest has been managed with various forms of partial harvesting for the past two centuries.Selection silviculture, which endeavours to maintain a balanced distribution of stems across size classes, has been practiced since the late 1970s.A multiple-use concept has been applied since 1984 to combine sustainable timber management with recreational use such as hunting, fishing, camping, hiking and snowmobiling [34].Stands across the forest are largely uneven-aged, and there are some regions where old-growth remains due to inaccessibility to harvesting machinery.Natural disturbances include frequent small-scale gap formation and infrequent large-scale wind events [36].Stand structure and composition, including SDDs, have important effects on productivity, timber yield, habitat suitability, and recreational opportunities in the forest and its surrounding region.
Remote Sens. 2017, 9, 944 3 of 17 characterized by undulating, rough topography typical of Precambian Shield landscapes [35].The forest has been managed with various forms of partial harvesting for the past two centuries.Selection silviculture, which endeavours to maintain a balanced distribution of stems across size classes, has been practiced since the late 1970s.A multiple-use concept has been applied since 1984 to combine sustainable timber management with recreational use such as hunting, fishing, camping, hiking and snowmobiling [34].Stands across the forest are largely uneven-aged, and there are some regions where old-growth remains due to inaccessibility to harvesting machinery.Natural disturbances include frequent small-scale gap formation and infrequent large-scale wind events [36].Stand structure and composition, including SDDs, have important effects on productivity, timber yield, habitat suitability, and recreational opportunities in the forest and its surrounding region.

Aerial Input Data
Discrete-return airborne LiDAR data were collected in August 2009 using a four-pass system (Optech ALTM 3100: fight altitude = 1500 m, pulse repetition frequency = 70 kHz, pass overlap = 30%), which can record up to four returns for each pulse.The resultant point cloud is made up of returns that each include an x, y and z coordinate and a return number.The z coordinates were converted to the height above ground by subtracting the ground height estimated by a digital elevation model (carried out using Optimal Geomatics' proprietary methods).First returns are least affected by acquisition properties [37], and therefore methods using only first returns are considered more robust to applications in new areas [38].Consequently, our analysis considered only first returns, which had an average point density of 2 points/m 2 .
Across the area of LiDAR coverage, ADS40 LEICA imagery was captured for the purpose of photo-interpretation of the forest resources of Ontario.The stereoscopic imagery was interpreted by

Aerial Input Data
Discrete-return airborne LiDAR data were collected in August 2009 using a four-pass system (Optech ALTM 3100: fight altitude = 1500 m, pulse repetition frequency = 70 kHz, pass overlap = 30%), which can record up to four returns for each pulse.The resultant point cloud is made up of returns that each include an x, y and z coordinate and a return number.The z coordinates were converted to the height above ground by subtracting the ground height estimated by a digital elevation model (carried out using Optimal Geomatics' proprietary methods).First returns are least affected by acquisition properties [37], and therefore methods using only first returns are considered more robust to applications in new areas [38].Consequently, our analysis considered only first returns, which had an average point density of 2 points/m 2 .
Across the area of LiDAR coverage, ADS40 LEICA imagery was captured for the purpose of photo-interpretation of the forest resources of Ontario.The stereoscopic imagery was interpreted by trained forest managers to classify the forest according to Ontario's Forest Resources Inventory program guidelines and implemented by certified forest inventory interpreters in conjunction with an extensive ground-based plot network [39].We used information on the proportion of the canopy occupied by individual species (i.e., aerial estimates of species composition: SPECIES AERIAL ) in combination with the ground-based plot data to develop a binary ecosite classification: sugar maple dominated stands and mixed stands.
We first performed a Principal Component Analysis (PCA) on the relative abundance of the eight most prevalent species recorded in 114 calibration plots (described below), with all other species grouped as broadleaves or conifers (SPECIES GROUND ).We then developed rule-based conditions for ecosites (forest types) that best separated the calibration plots along the first and second PCA axes.We defined the sugar maple ecosite to be where SPECIES AERIAL contained low proportions of balsam fir, white spruce and eastern hemlock (<20%) and a high proportion of broadleaves (>60%), which were predominantly sugar maple.Areas that did not meet the conditions for the sugar maple ecosite were defined to belong to the mixture ecosite.This ecosite definition separated the majority of plots on the first two PCA axes based on their species composition as a proportion of both stem numbers and basal area (Figure 2).
HD LIDAR for each 0.25-ha LiDAR tile was defined by the proportion of first returns falling into a set of 1-m height bins.We also calculated gap fraction and top canopy height as LiDAR summary metrics for each tile.Gap fraction (GF) was calculated as the proportion of first returns recorded below 2 m [20].Top canopy height (TCH) was calculated by averaging the maximum height recorded from each complete 1 m 2 section of the tile, excluding gaps where all returns were below 2 m.The tile-level data that we used to predict stem diameter distributions (SDD PREDICT ) therefore consisted of HD LIDAR , TCH, GF and ecosite.[39].We used information on the proportion of the canopy occupied by individual species (i.e., aerial estimates of species composition: SPECIESAERIAL) in combination with the ground-based plot data to develop a binary ecosite classification: sugar maple dominated stands and mixed stands.
We first performed a Principal Component Analysis (PCA) on the relative abundance of the eight most prevalent species recorded in 114 calibration plots (described below), with all other species grouped as broadleaves or conifers (SPECIESGROUND).We then developed rule-based conditions for ecosites (forest types) that best separated the calibration plots along the first and second PCA axes.We defined the sugar maple ecosite to be where SPECIESAERIAL contained low proportions of balsam fir, white spruce and eastern hemlock (<20%) and a high proportion of broadleaves (>60%), which were predominantly sugar maple.Areas that did not meet the conditions for the sugar maple ecosite were defined to belong to the mixture ecosite.This ecosite definition separated the majority of plots on the first two PCA axes based on their species composition as a proportion of both stem numbers and basal area (Figure 2).
HDLIDAR for each 0.25-ha LiDAR tile was defined by the proportion of first returns falling into a set of 1-m height bins.We also calculated gap fraction and top canopy height as LiDAR summary metrics for each tile.Gap fraction (GF) was calculated as the proportion of first returns recorded below 2 m [20].Top canopy height (TCH) was calculated by averaging the maximum height recorded from each complete 1 m 2 section of the tile, excluding gaps where all returns were below 2 m.The tile-level data that we used to predict stem diameter distributions (SDDPREDICT) therefore consisted of HDLIDAR, TCH, GF and ecosite.The biplot shows the first two axes of a Principal Components Analysis (PCA), and the percentage of variation explained by each axis (in parentheses).The PCA was performed on the main species in each calibration plot as a proportion of (a) the total stem number and (b) the total basal area; the shape and colour of the points denote the ecosite classification for these plots.

Ground Plot Data
In 2010-11, 1540.25-hacircular ground plots were inventoried across the study area (Figure 1).The plots were stratified to span a gradient in both species composition and canopy openness based on the aerial imagery described above (see Table 1 for a summary of stand characteristics).We randomly selected 114 of these plots (the same calibration plots from [32]) and used them to calculate the mean, variance, and covariance of SDD parameters and structural attributes across the landscape.The biplot shows the first two axes of a Principal Components Analysis (PCA), and the percentage of variation explained by each axis (in parentheses).The PCA was performed on the main species in each calibration plot as a proportion of (a) the total stem number and (b) the total basal area; the shape and colour of the points denote the ecosite classification for these plots.

Ground Plot Data
In 2010-11, 1540.25-hacircular ground plots were inventoried across the study area (Figure 1).The plots were stratified to span a gradient in both species composition and canopy openness based on the aerial imagery described above (see Table 1 for a summary of stand characteristics).We randomly selected 114 of these plots (the same calibration plots from [32]) and used them to calculate the mean, variance, and covariance of SDD parameters and structural attributes across the landscape.The remaining 40 plots (validation plots) were used to evaluate model performance.Data for each plot consisted of the species and diameter at breast height (dbh) of all stems with a dbh ≥ 8 cm.For each calibration plot, we fit a truncated Weibull distribution to the stem diameters to estimate shape (k) and scale (λ) parameters of its SDD.We calculated the proportion of broadleaves (p B ) based on total basal area.We also calculated quadratic mean diameter (QMD, in cm; QMD C for conifers and QMD B for broadleaves), aboveground carbon density (ACD, in Mg C/ha), basal area (BA, in m 2 /ha; BA C for conifers and BA B for broadleaves), and the density of stems with dbh ≥ 8 cm (N, in ha −1 ; N C for conifers and N B for broadleaves) for each plot.For ACD, individual stem biomass was calculated from dbh using Canadian biomass equations [40,41], then summed to obtain the total plot biomass.Biomass was halved to give the carbon content [42].The biomass equations used to estimate ACD from field data were species-specific, but generic allometries for conifers and broadleaves were used in calculating ACD from SDD PREDICT .These attributes were variously used for model priors and to measure model performance, as detailed below.

Overview of the Modelling Approach
We sought to predict SDDs, along with the relative abundance of broadleaves and conifers, from the height distribution of first returns in 0.25 ha LiDAR tiles (HD LIDAR ).Briefly, our approach estimates four parameters (described in the following section) for each tile that specify a SDD that, when converted into a predicted distribution of LiDAR first return heights (HD PREDICT ) using the allometry-based algorithm described in [32], most closely approximates HD LIDAR for a given tile.The allometry-based algorithm uses relationships between stem diameter, tree height, crown diameter, and crown taper to estimate the canopy height profile for a given SDD, then uses correction factors that account for crown overlap, crown permeability, and gaps to compute a corresponding HD PREDICT .The resulting SDD is used to calculate a suite of structural attributes (BA, N, QMD, and ACD) for each LiDAR tile across the full landscape (an overview of the approach used is depicted in Figure 3).
We used Bayesian parameter estimation methods, for two reasons.Firstly, a Bayesian approach allows us to specify prior distributions for model parameters that discourage the model from predicting SDDs that are inconsistent with those observed across the set of calibration plots, or that do not fit simple empirical relationships with LiDAR summary statistics.Secondly, a Bayesian approach allows us to formally propagate parameter uncertainty into model predictions.Parameter and model uncertainty is considered inherently within a single framework, which can help with identifying and challenging the source of model uncertainties.We evaluate the performance of our model by comparing SDD PREDICT and related stand attributes to observations from the validation plots.The prior and likelihood together inform the selection of the next set of parameters using the Metropolis-Hastings criterion.This process is iterated many times to obtain a chain of parameter samples from the posterior distribution.The bottom panel describes the five parameters being estimated: k, λ, N, p , σ .

Description of the Predictive Model
Our model (Figure 3) describes the SDD of a given LiDAR tile with four parameters: the total number of stems with a diameter at breast height (dbh) ≥ 8 cm (N), the proportion of these stems that are broadleaved species (p ), and the shape (k) and scale (λ) parameters of a truncated Weibull distribution for stem diameters.We only considered stems with a dbh ≥ 8 cm to reflect the ground plot measurements.To generate a list of broadleaf and conifer stem diameters (D and D ), N is first partitioned into the number of broadleaves ( N = N × p ) and the number of conifers (N = N -N ).D and D are then assumed to correspond to evenly-spaced quantiles of a Weibull distribution.For sets of i = 1 … N broadleaf and i = 1 … N conifer stems, we define quantiles q = (i − 0.5)/N , where N is either N or N .The set of broadleaf and conifer stem diameters are then obtained from the Weibull's inverse cumulative distribution function:

Description of the Predictive Model
Our model (Figure 3) describes the SDD of a given LiDAR tile with four parameters: the total number of stems with a diameter at breast height (dbh) ≥ 8 cm (N), the proportion of these stems that are broadleaved species (p B ), and the shape (k) and scale (λ) parameters of a truncated Weibull distribution for stem diameters.We only considered stems with a dbh ≥ 8 cm to reflect the ground plot measurements.To generate a list of broadleaf and conifer stem diameters (D B and D C ), N is first partitioned into the number of broadleaves (N B = N × p B ) and the number of conifers (N C = N -N B ). D B and D C are then assumed to correspond to evenly-spaced quantiles of a Weibull distribution.For sets of i = 1 . . .N B broadleaf and i = 1 . . .N C conifer stems, we define quantiles q i = (i − 0.5)/N X , where N X is either N B or N C .The set of broadleaf and conifer stem diameters are then obtained from the Weibull's inverse cumulative distribution function: where D X represents either D B or D C , q is taken from the set of quantiles for either broadleaf or conifer stems (as appropriate), and d min is the minimum stem diameter measured in the plots (8 cm).This procedure yields a list of stems that matches the number, composition, and diameter distribution specified by the given parameters.
We next apply Spriggs et al.'s (2015) allometry-based algorithm (we use the two-allometry version for simplicity) to the list of stems in order to generate the predicted height distribution of LiDAR first returns (HD PREDICT ) [32].Given an observed distribution of LiDAR first returns for a particular tile (HD LIDAR ), our goal is to infer the Bayesian posterior distribution of SDD parameters, θ = {k, λ, N, p B , σ L }, that could have generated the observed LiDAR returns.In addition to the four parameters described above, θ also includes a noise term (σ L ) for the residual error between HD LIDAR and HD PREDICT .Using Bayes' theorem: The term on the left represents the joint posterior distribution of SDD parameters for a given LiDAR tile (P( θ|HD LIDAR )), and the terms on the right represent the model likelihood (P( HD LIDAR |θ)), prior parameter distribution (P(θ)), and a normalisation constant (η) that ensures the posterior is a valid probability distribution.We used Markov chain Monte Carlo (MCMC) sampling to generate samples of parameter values from the posterior distribution using: (1) a likelihood function measuring the goodness-of-fit between HD LIDAR and HD PREDICT for given values of θ, and (2) informative priors, which specify the probability of a given set of θ from related (but statistically independent) analyses of the LiDAR point cloud data, an aerial imagery-derived ecosite classification, and the set of calibration plots.(The normalization constant, which ensures the posterior distribution integrates to one, is ignored when using MCMC sampling.) Parameter estimation was carried out in C# using the Filzbach library for Bayesian inference [43].For each tile, three independent Markov chains were run where each chain consisted of a burn-in phase of 1500 iterations and a sampling phase of 5000 iterations.We thinned the samples from each chain by retaining every 50th value.All chains converged within the prescribed number of burn-in iterations.We describe the formulations of our model likelihood and prior distributions in detail in the sections that follow.

Model Likelihood
For a given tile, the likelihood of HD LIDAR given the HD PREDICT generated from θ is calculated from differences in the proportion of returns falling within each 1-m height interval (h where h = H at the highest interval).Specifically, the likelihood for a given set of parameters (θ) is the product of the probability densities of each of these differences under a normal distribution, where f represents the normal probability density function with mean zero and standard deviation σ L .

Prior Distributions
We incorporated informative priors within our modelling approach to help address the issue that a single LiDAR distribution could be generated by a number of different stand structures.The priors were chosen to ensure that the sampled parameters represented realistic stands and to include the additional information available from the ecosite classification data.As a result, the joint prior distribution for the model parameters incorporated tile-level information from the ecosite classification and LiDAR summary metrics (TCH and GF), as well as landscape-level information on stand attributes from the 114 calibration plots.We used these different sources of information to define two separate components of the joint prior distribution: where P S is a prior for the SDD parameters, and P M is a prior for stand attributes calculated from the SDD parameters.For both P S and P M , we derived a prior probability distribution from the mean, variance, and covariance of SDD parameters and stand attributes across the set of calibration plots.The prior for P S also used tile-level estimates of TCH and GF to help predict stem density.Both components of the prior were specified separately for sugar maple and mixture ecosites.P S was based on information from LiDAR and the ground-plot data.Firstly, we used the summary LiDAR metrics TCH and GF to help estimate stem density as part of P S .Using linear regression, we found that the total number of stems in a given tile was related to TCH and GF for each ecosite: We used this relationship to define residual variation in stem density as ε N = N − N. Next, we reasoned that estimated values for the SDD parameters should fall within the landscape-level distribution observed across the set of calibration plots.We therefore defined P S based on the mean, variance, and covariance of the parameters k, λ, p B , as well as of ε N , across the calibration plots: where MVN represents a multivariate normal distribution with a mean vector µ and covariance matrix Σ. Parameter values for P S are listed in Table S1 (see also Figure S2).Finally, we expected stand attributes derived from the SDD parameters to also follow the landscape-level distribution observed in the calibration plots.We chose five structural attributes (BA C , BA B , QMD C , QMD B and TCH) that could be calculated for a given SDD, and specified their contribution to the joint prior distribution as: P M helps ensure that a given set of SDD parameters yields stands that have a similar overall structure to those observed across the landscape.Rather than being a probability distribution for the parameters themselves, P M is defined as a probability distribution for quantities that are calculated from the list of stems associated with given SDD parameters (somewhat analogous to defining priors on the predicted response from a regression-type model [44]).The parameter values for P M are listed in Table S1.It is important to note that since we estimate a separate set of SDD parameters for each individual tile, the fitted relationships that inform the priors only serve to constrain the SDD parameters and are not used to predict the SDD parameters directly, as in standard area-based approaches.

Model Output and Performance
Model performance was evaluated with the validation plots by comparing the predicted SDDs and associated stand attributes with the corresponding plot measurements.The fit of the predicted stem diameter distributions (SDD PREDICT ) to the true observed distributions (SDD TRUE ) was quantified by binning the SDDs into 3 cm diameter intervals (i), where T is the total number of bins, then calculating the coefficient of determination (R 2 ) as: where SDD TRUE is the mean number of stems in each bin of the observed SDD.A positive R 2 indicates that model predictions are an improvement over the simple mean number of stems, while R 2 = 1 indicates perfect agreement between the predicted and observed number of stems in each diameter bin.We also used the Reynolds error index (e) [45] and root mean square error (RMSE) to compare our predictions for SDDs and stand attributes with previous studies: .100, ( 8) where P i is the predicted value and O i is the measured value from a given plot i and V is the number of validation plots (V = 40).

SDD Predictions
The model showed moderate overall success in predicting SDDs (Figure 4).SDDs were predicted well (R 2 ≥ 0.75) in half of the 40 validation plots and moderately well (0.5 ≤ R 2 < 0.75) in another nine plots.The remaining 11 plots did not demonstrate particularly good SDD predictions (R 2 < 0.50), with six of these being worse than a simple mean.The mean Reynolds error index was 53.3 ± 24.7 across all plots (Table 2).Overall, sugar maple plots were predicted more accurately than mixture plots (19 of 24 vs. 10 of 16 plots having R 2 > 0.5), a finding also supported by the Reynolds error index.
There was no significant relationship between R 2 and the number of stems or basal area (Pearson correlation coefficient = −0.05 and 0.12, respectively).Small-diameter stems were the most poorly predicted (see Figure S3), being generally over-predicted in the lower basal area plots and under-predicted in higher basal area plots.As an indication of the impact on predictive performance, if the smallest bin is excluded from the R 2 calculation (i.e., considering only stems with dbh ≥ 11 cm), the number of moderately well-fitting plots (R 2 ≥ 0.5) increases from 29 to 34 (85%).correlation coefficient = −0.05 and 0.12, respectively).Small-diameter stems were the most poorly predicted (see Figure S3), being generally over-predicted in the lower basal area plots and underpredicted in higher basal area plots.As an indication of the impact on predictive performance, if the smallest bin is excluded from the R 2 calculation (i.e., considering only stems with dbh ≥ 11 cm), the number of moderately well-fitting plots (R 2 ≥ 0.5) increases from 29 to 34 (85%).

Stand Attributes
Among the stand attributes considered (Figure 5), our model was most effective in predicting quadratic mean diameter (RMSE = 3.64 cm), aboveground carbon density (RMSE = 21.4Mg C/ha), and basal area (RMSE = 6.99 m 2 /ha).Weibull shape and scale parameters were generally over-predicted (RMSE: shape = 0.504 and scale = 7.87).Total number of stems in sugar maple plots was often over-predicted, whilst in mixture plots stem numbers tended to be under-predicted (RMSE = 68.0).Overall, sugar maple plots were better predicted than the mixture plots (BA RMSE: 4.84 vs. 9.32 m 2 /ha; N RMSE: 62.4 vs. 75.6).Partitioning basal area and the number of stems into broadleaf and conifer components revealed that broadleaves were estimated more accurately than conifers (BA RMSE: broadleaves = 5.93 m 2 /ha, conifers = 7.39 m 2 /ha; N RMSE: broadleaves = 38.6,conifers = 67.3).The errors were spatially uniform in the validation plots across the landscape.

Applying the Model Across the Landscape
We used the predicted SDDs and associated stand attributes to map size distributions, forest composition, and aboveground carbon at a landscape scale (Figure 6).Predicted ACD had a

Applying the Model Across the Landscape
We used the predicted SDDs and associated stand attributes to map size distributions, forest composition, and aboveground carbon at a landscape scale (Figure 6).Predicted ACD had a maximum of 176 Mg C/ha, a mean of 74 Mg C/ha, and a standard deviation of 27 Mg C/ha (Figure 6a).Our best estimate for the total ACD of the landscape was 1.611 ± 0.002 Tg C (±standard deviation, representing predictive uncertainty).Over 75% of tiles are predicted to be at least half broadleaved, with just 10% predicted to have less than 30% broadleaf stems (Figure 6b).The Weibull size distribution's parameters were strongly related to one another (Pearson correlation coefficient = 0.82; see Figure S2b), and therefore the form of the SDD could largely be inferred from the shape parameter alone.Figure 6c shows a slightly higher prevalence of the "reverse J" distribution in the southwest, northwest, and northeast, but that there are large areas exhibiting a unimodal distribution in the central part of the forest.Though a full geospatial analysis is beyond the scope of this work, our maps reveal strong spatial relationships between the three attributes (ACD, proportion of broadleaves and Weibull distribution of stem diameters), with conifer-dominated areas supporting a lower carbon density and being most often dominated by small-diameter stems.
with just 10% predicted to have less than 30% broadleaf stems (Figure 6b).The Weibull size distribution's parameters were strongly related to one another (Pearson correlation coefficient = 0.82; see Figure S2b), and therefore the form of the SDD could largely be inferred from the shape parameter alone.Figure 6c shows a slightly higher prevalence of the "reverse J" distribution in the southwest, northwest, and northeast, but that there are large areas exhibiting a unimodal distribution in the central part of the forest.Though a full geospatial analysis is beyond the scope of this work, our maps reveal strong spatial relationships between the three attributes (ACD, proportion of broadleaves and Weibull distribution of stem diameters), with conifer-dominated areas supporting a lower carbon density and being most often dominated by small-diameter stems.In (c) red regions correspond with a "reverse J" distribution and yellow regions with a unimodal distribution when predicted SDDs are reduced to the single bestfit axis of variation depicted in the inset figure (also see Figure S2).White regions represent nonforested areas, such as lakes, or land not under the same ownership as the forest.In (c) red regions correspond with a "reverse J" distribution and yellow regions with a unimodal distribution when predicted SDDs are reduced to the single best-fit axis of variation depicted in the inset figure (also see Figure S2).White regions represent non-forested areas, such as lakes, or land not under the same ownership as the forest.
The model is founded upon a mechanistic interpretation of the interactions between LiDAR first returns and the canopy, yet still operates at the plot level rather than on individual tree crowns.While it shares some of the allometric characteristics of ITD models, it is not spatially explicit and does not carry such high point density requirements.All components of the underlying algorithm relating stem diameters to LiDAR return heights are interpretable and can be re-parameterized or modified as required for different forests.
Testing with validation plots showed that model performance was good in most cases, despite being applied in a relatively challenging study area compared with other temperate or boreal forests.Most previous studies that have estimated SDDs from airborne LiDAR have been conducted in forests dominated by conifers, particularly Norway spruce and Scots pine [13][14][15]17,18,46,47].We have tested our model on an uneven-aged mixed temperate forest, where fewer studies have attempted to predict SDDs from airborne LiDAR [16].We aimed to mechanistically incorporate general features of a forest crowns to enable the model to be modified for use in new regions.However, the accuracy of the model is likely to be reduced in areas where the understory layer contributes a significant proportion of the ACD, such as in tropical forests.

Assessing SDD Predictions
The model predicted SDD distributions reasonably well (R 2 > 0.5) for 85% of the validation plots, with sugar maple plots predicted slightly more accurately than mixture plots (Figure 4).The smallest stems (8 ≤ dbh < 11 cm) were the least well predicted, as has also been reported in a comparison of parametric and non-parametric ABAs [20].The smallest trees are most likely problematic because they are not captured well by LiDAR pulse data [48], particularly when considering first returns from low density data [49].Issues also arise from the stems falling just below the measurement threshold of 8 cm, but which occupy the same space as those falling just above the threshold.
Our approach yielded a mean Reynolds error index value of 53.3 in predicting SDDs, which is at the low end of estimates from coniferous forests: e.g., 55.1 and 59.2 [20], 69.5-184.3[50], 59.7-111.3[51], 70-152 [52] and 619.9-943.8[47].Our validation plots were larger (0.25 ha) than those used in many other studies, however increased plot size would act to reduce model error at a larger spatial scale (all else being equal), thereby biasing comparisons in our favor.

Assessing Predictions of Stand Attributes
Few studies have predicted SDDs from airborne LiDAR and then estimated stand attributes from the resultant list of diameters.Mean dbh has been presented in some cases [15,46], but more commonly stand attributes are estimated directly using a standard ABA, which is then used to calibrate stem diameter frequencies [14,17].Vauhkonen and Mehtätalo (2015) [47] used a combination of an ABA and ITD to provide estimates for stem number per hectare, basal area and quadratic mean diameter, which allows for the most comprehensive comparison with our methods (their 0.04-ha plots were considerably smaller than ours, however).Our RMSE of 6.99 m 2 /ha for basal area and 3.64 cm for quadratic mean diameter fell within their recorded ranges (2.9-18.8m 2 /ha and 2.1-4.6 cm, respectively), whilst our RMSE of 272 (scaled to per hectare) for the total number of stems was slightly better than their results (299-1030 stems/ha).
To our knowledge, no studies have converted broadleaf and conifer stem diameters predicted from LiDAR to aboveground carbon density.Previous estimates of temperate forest biomass using purpose-built methods have produced a RMSE of 58.0 Mg/ha [53], 28.9 Mg/ha [54] and 31.6 Mg/ha [55].For comparison, doubling our RMSE for ACD (since we take ACD to be half the biomass) gives a value of 42.8 Mg/ha, which suggests that our model can estimate ACD or biomass reasonably well from SDD predictions.
Broadleaves were predicted more accurately than conifers both in terms of basal area and total stem number (e.g., the difference in BA RMSE for conifers vs. broadleaves was 1.46 m 2 /ha).This is contrary to the findings of Peuhkurinen et al. (2008) [18] who, in their study aimed at predicting species-specific diameter distributions, predicted the volume of deciduous species with only half the accuracy of conifers in boreal forest plots.Broadleaves were minor species in their study and dominant in ours, and estimation methods are expected to best fit the dominant species.Sugar maple plots were also generally better predicted than mixture plots, which follows from their higher proportion of broadleaf stems.Conifers may be captured less accurately by our model because the reflective properties of coniferous crowns are difficult to capture mechanistically [56].We partly accounted for this by incorporating a crown permeation term in the algorithm that predicted return heights, but crown permeability would likely increase further if this algorithm were re-tuned to coniferous stand data.

Implementing the Model across the Landscape
We have presented a model that translates LiDAR data on canopy structure into spatially explicit estimates of stand-level attributes that are widely-used by forest managers.Mapping these attributes reveals interesting landscape patterns; for example, Figure 6 shows that coniferous stands tend to correspond with high small stem frequencies and low carbon densities.Further spatial analyses could reveal how these patterns relate to factors such as elevation, slope, aspect, and natural or anthropogenic disturbances.
There was considerable variation in SDDs across the landscape, with a higher prevalence of unimodal distributions towards the central area of the forest.This type of information is particularly valuable to forest managers, who rely on maps of stand structure for planning harvest schedules.Under ecosystem-based forest management, it is important that harvests balance timber yields against regeneration potential and non-timber values, such as the provision of wildlife habitat and aesthetics.Carbon accounting further demands quantitative estimates of current timber stocks and their uncertainty [57].Spatial estimates of stem numbers and sizes from LiDAR can enable predictions for expected timber yield across the landscape, while also ensuring that forest age and compositional structure meets a range of targets for non-timber values.
Applying the model to a repeat LiDAR dataset would provide further insight into how the landscape is changing over time.The broadleaf and conifer stem diameters outputted from the model may also be used to both test and initialize simulation models of forest dynamics, which predict how the structure of individual stands will evolve through tree-level growth and mortality.The combined application of our model to future LiDAR surveys and forest simulator projections could enable a new understanding of the landscape-scale response of forest ecosystems to changes in environmental conditions [58].

Conclusions
In this study, we have developed an area-based approach to translating low density airborne LiDAR data into stem diameter distributions and the proportion of broadleaves vs. conifers using an allometry-based algorithm and Bayesian model-fitting techniques.We have demonstrated that the method can provide comparative estimates for stand structure to previous studies and has been implemented across a structurally heterogonous temperate forest.We made our model as general as possible by incorporating canopy features mechanistically and by using first returns to minimize the effect of LiDAR acquisition properties.A natural progression would be to test the model in boreal forests as well as other temperate forests to assess its broader applicability.The model could be tested without re-calibration, since it is based on general conifer and broadleaf allometries, using 0.25-ha tile data for HD LIDAR , TCH, GF and an ecosite classification (alternatively, a forest could be classified entirely as sugar maple or a mixture).Further refinements might use new plot data to re-parameterize the model priors and stem number relationship, or to include other known predictive relationships, to better constrain parameter estimates for a new area.If a suitable network of plots already exists (as is common where LiDAR is used for forest inventories [11]), then a final step would be to re-parameterize the allometry-based algorithm for translating plot inventory data to

Figure 1 .
Figure 1.Locations of 154 calibration and validation ground plots in Haliburton Forest and Wildlife Reserve in Ontario, Canada.Black dots represent calibration plots and light grey dots represent validation plots.

Figure 1 .
Figure 1.Locations of 154 calibration and validation ground plots in Haliburton Forest and Wildlife Reserve in Ontario, Canada.Black dots represent calibration plots and light grey dots represent validation plots.

Figure 2 .
Figure 2. Separation of calibration plots according to an ecosite classification defined from aerial imagery.The biplot shows the first two axes of a Principal Components Analysis (PCA), and the percentage of variation explained by each axis (in parentheses).The PCA was performed on the main species in each calibration plot as a proportion of (a) the total stem number and (b) the total basal area; the shape and colour of the points denote the ecosite classification for these plots.

Figure 2 .
Figure 2. Separation of calibration plots according to an ecosite classification defined from aerial imagery.The biplot shows the first two axes of a Principal Components Analysis (PCA), and the percentage of variation explained by each axis (in parentheses).The PCA was performed on the main species in each calibration plot as a proportion of (a) the total stem number and (b) the total basal area; the shape and colour of the points denote the ecosite classification for these plots.

Figure 3 .
Figure 3. Overview of the Markov chain Monte Carlo approach used for estimating the parameters required to predict stem diameter distributions (SDDs).In the top panel, the HDLIDAR-SDD model includes five parameters: four to generate SDDPREDICT (from which the allometry-based algorithm produces the accompanying HDPREDICT), and one for the residual error.The likelihood quantifies the goodness-of-fit between HDPREDICT and the LiDAR tile inputted into the model (HDLIDAR).The prior reflects the probability of the current set of parameters based on conditions found across a set of calibration plots within the study area combined with a LiDAR regression and ecosite classification.The prior and likelihood together inform the selection of the next set of parameters using the Metropolis-Hastings criterion.This process is iterated many times to obtain a chain of parameter samples from the posterior distribution.The bottom panel describes the five parameters being estimated: k, λ, N, p , σ .

Figure 3 .
Figure 3. Overview of the Markov chain Monte Carlo approach used for estimating the parameters required to predict stem diameter distributions (SDDs).In the top panel, the HD LIDAR-SDD model includes five parameters: four to generate SDD PREDICT (from which the allometry-based algorithm produces the accompanying HD PREDICT ), and one for the residual error.The likelihood quantifies the goodness-of-fit between HD PREDICT and the LiDAR tile inputted into the model (HD LIDAR ).The prior reflects the probability of the current set of parameters based on conditions found across a set of calibration plots within the study area combined with a LiDAR regression and ecosite classification.The prior and likelihood together inform the selection of the next set of parameters using the Metropolis-Hastings criterion.This process is iterated many times to obtain a chain of parameter samples from the posterior distribution.The bottom panel describes the five parameters being estimated: k, λ, N, p B , σ L .

Figure 4 .
Figure 4. Stem diameter distributions for the 40 validation plots.Histograms show the measured stems, points show the predicted number of stems within each bin, and the solid line shows the Weibull distribution fitted to the data.The plots are ordered from low to high total stem numbers; note the changing y-axis.

Figure 4 .
Figure 4. Stem diameter distributions for the 40 validation plots.Histograms show the measured stems, points show the predicted number of stems within each bin, and the solid line shows the Weibull distribution fitted to the data.The plots are ordered from low to high total stem numbers; note the changing y-axis.

Figure 5 .
Figure 5. Predictions of stand attributes versus the measured values in validation plots for (a) Weibull shape parameter for the SDD, (b) Weibull scale parameter for the SDD, (c) quadratic mean diameter, (d) aboveground carbon density, (e) total basal area and (f) total number of stems.Solid lines denote 1:1 relationships.The shape and colour of the points denote the ecosite classification of the plot.

Figure 5 .
Figure 5. Predictions of stand attributes versus the measured values in validation plots for (a) Weibull shape parameter for the SDD, (b) Weibull scale parameter for the SDD, (c) quadratic mean diameter, (d) aboveground carbon density, (e) total basal area and (f) total number of stems.Solid lines denote 1:1 relationships.The shape and colour of the points denote the ecosite classification of the plot.

Figure 6 .
Figure 6.Mapped predictions of stand characteristics across the forest landscape.Purple regions correspond with (a) low ACD and (b) low broadleaf density, and green regions correspond with (a) high ACD and (b) high broadleaf density.In (c) red regions correspond with a "reverse J" distribution and yellow regions with a unimodal distribution when predicted SDDs are reduced to the single bestfit axis of variation depicted in the inset figure (also see FigureS2).White regions represent nonforested areas, such as lakes, or land not under the same ownership as the forest.

Figure 6 .
Figure 6.Mapped predictions of stand characteristics across the forest landscape.Purple regions correspond with (a) low ACD and (b) low broadleaf density, and green regions correspond with (a) high ACD and (b) high broadleaf density.In (c) red regions correspond with a "reverse J" distribution and yellow regions with a unimodal distribution when predicted SDDs are reduced to the single best-fit axis of variation depicted in the inset figure (also see FigureS2).White regions represent non-forested areas, such as lakes, or land not under the same ownership as the forest.

Table 1 .
Mean (±s.d.)values for stand structural attributes across the set of calibration plots.

Table 2 .
Reynolds error index values for all validation plots and for sugar maple and mixture ecosites separately.Lower values indicate a closer match between the predicted and observed stem distributions.