Modeling the Effect of the Spatial Pattern of Airborne Lidar Returns on the Prediction and the Uncertainty of Timber Merchantable Volume

Lidar data are regularly used to characterize forest structures. In this study, we determine the effects of three lidar attributes (density, spacing, scanning angle) on the accuracy and the uncertainty of timber merchantable volume estimates of balsam fir stands (Abies balsamea (L.) Mill.) in eastern Canada. We used lidar point clouds to compute predictor variables of the merchantable volume in a nonlinear model. The best model included the mean height of first returns, the proportion of first returns below 2 m and the canopy surface roughness index. Our analysis shows a high correlation between lidar and field data of 119 plots (pseudo-R2 = 0.91), however, residuals were heteroscedastic. More precise parameter estimates were obtained by adding to the model a variance function of variables describing the mean height of returns and the skewness of the area distribution of triangulated lidar returns. The residual standard deviation was better estimated (3.7 m3 ha−1 multiplied by the variance function versus 28.0 m3 ha−1). We found no effect of density on the predictions (p-value = 0.74). This suggests that the height and the spatial pattern of returns, rather than the density, should be considered to better assess the uncertainty of merchantable volume estimates.


Introduction
Airborne Laser Scanning (ALS or lidar) has proven to be a useful 3d tool for the measurement of forest attributes over the past decades.It generates point clouds which are a three-dimensional representation of the volumetric interaction between pulse photons and illuminated objects.The spatial distribution of returns within the point cloud depends on the characteristics of the lidar systems and the target object [1,2].The scan mechanisms (e.g., oscillating mirror, rotating mirror, nutating mirror) mounted on lidar scanners produce different scanning patterns of pulses on the ground [1].For example, the bi-directional scan mechanism of oscillating mirrors produces a see-saw scanning pattern where pulses tend to be more homogeneously distributed along the center of the flight line than at the borders.Consequently, an alternation of local clumps of returns and local gaps will be observed within the point cloud [3].Scanning angle is another setting which influences the spatial distribution pattern.For example, high-scanning angles increase the distance traveled by the incident pulse.Objects located further from the lidar sensor are likely to be occluded or undetected [4], thus producing an irregular spatial distribution of returns.
The effects of the surface characteristics of scanned objects on the scanning pattern are not controlled by the user.Gatziolis et al. [1], for example, demonstrated that higher return densities are recorded in forest stands compared to green pastures under similar lidar settings.Korpela et al. [2] noted that leaf size and orientation and foliage density affect the intensity of lidar returns.Balsa-Barreiro et al. [5] noted that topography and land cover influence the density and spacing of returns.
Metrics (e.g., height metrics, canopy cover metrics), derived from point clouds are often used to model forest attributes such as volume [6,7].Hence, highly correlated models between lidar data and field volume have been achieved in many studies [8][9][10][11].Studies on lidar modeling focus on investigating the optimal explanatory variables that best explain the variability of the forest attributes.The relationship between the explanatory and dependent variable is obtained by means of a regression model.Least square procedures are commonly used to estimate the parameters of the regression model [12].These procedures focus on minimizing the sum of squared residuals.Heteroscedasticity occurs in a model when the variance of the residuals is not constant [12,13].Some possible causes of heteroscedasticity can be measurement errors of the data, autocorrelation or misspecification of the model.This can lead to a significant loss of precision where the estimated parameters of the model are inefficient, although unbiased.The standard errors and confidence intervals are also unreliable.In other words, the forest attribute can still be predicted, but the predictions are uncertain.In an operational context, uncertain predictions of a forest attribute such as the merchantable volume have a direct effect on the evaluation of its economic value [14].Nonetheless, heteroscedasticity may be corrected using different methods [12,15,16].Variable transformation for example modifies the measurement scale of variables in a model and heteroscedasticity may then be corrected [12,17].Variance modeling identifies the sources of heteroscedasticity in a model and attempts to better estimate the model uncertainty through a variance function [15].
Few studies can be found on the effects of lidar on the uncertainty of predictive models (e.g., [4,18,19]).Yet, with the increasing use of lidar in forestry, reliable estimates of forest attributes derived from lidar are necessary.There is, therefore, a need to investigate whether factors like the spatial pattern of returns can also influence the predictions or the uncertainty of modeled forest attributes.This study, therefore, aims to analyze the effects of three lidar attributes (return density, return spacing, scanning angle) on the predictions and the uncertainty of merchantable volume estimates.

Study Area
The study was conducted at the Montmorency Forest which is a teaching and research forest facility of the Université Laval.The forest is located in central Quebec, Canada (47.3 • N, 71.1 • W), about 70 km north of Quebec City (Figure 1).The study site covered an area of 66 km 2 .Elevation at the site ranged from approximately 460 to 1040 m above sea level.The average annual precipitation was 1589 mm, and the mean annual temperature was 0.3 • C [20].The site lies on the Laurentian Plateau and is part of the boreal forest domain.Balsam fir (Abies balsamea [L.] Miller), black spruce (Picea mariana [Miller] BSP and white spruce (Pinus glauca [Moench] Voss) are the dominant conifers found in the forest.White birch (Betula papyrifera Marshall) and trembling aspen (Populus tremuloides) are also common [21].

Field Plot Data
The Montmorency Forest has a network of circular permanent plots (radius = 11.28 m, area = 400 m 2 ).The plots center positions have been georeferenced using a dual-frequency antenna GPS (Trimble Yuma, accuracy ~30 cm) and were not post-processed.Plots are re-measured on a five-year cycle.In this study, we used field data that were collected between 2007 and 2014.A sample of 119 plots was selected following these criteria: the dominant species was balsam fir, the last cut was before 1996 and the minimum dominant height was 7 m.The merchantable trees with a diameter at breast height (dbh) above 9 cm were measured using a diameter tape.Tree height was measured on eight sample trees per plot using a Vertex hypsometer.A height-diameter relationship was adjusted for each plot (nonlinear mixed effects model) and then used to calculate the height of the other trees.The merchantable volume MV was computed using height and diameter based equations as described in [22].The total MV was computed as the sum of MVs for all trees.Twenty-four field plots had been initially inventoried in 2011, the same year as the lidar survey.Attributes of the other plots were estimated by interpolating field measurements between two dates of inventory.Table 1 presents a summary of the field data for 2011.
(dbh) above 9 cm were measured using a diameter tape.Tree height was measured on eight sample trees per plot using a Vertex hypsometer.A height-diameter relationship was adjusted for each plot (nonlinear mixed effects model) and then used to calculate the height of the other trees.The merchantable volume MV was computed using height and diameter based equations as described in [22].The total MV was computed as the sum of MVs for all trees.Twenty-four field plots had been initially inventoried in 2011, the same year as the lidar survey.Attributes of the other plots were estimated by interpolating field measurements between two dates of inventory.Table 1 presents a summary of the field data for 2011.

Lidar Data
Airborne laser data were acquired in August 2011 with an ALTM 3100 discrete return sensor (Teledyne Optech, Vaughan, Canada).The average flying altitude was 1000 m above ground.The pulse repetition rate was of 100 kHz with a 50% strip overlap.The scanning angle ranged from 0° to 24° relative to the nadir.The sensor could record up to four measurements per pulse.Each field plot was covered by two to four lidar strips.The average above-ground return density was 6.4 m −2 .Information for each return was recorded in las files (version 1.0).The returns classified as ground (by the vendor) were interpolated to construct a digital terrain model using the LAStools software [6].Non-ground returns were normalized by measuring the elevation difference with the underlying terrain model.Table 2 presents a summary of the lidar data of the study site.

Lidar Data
Airborne laser data were acquired in August 2011 with an ALTM 3100 discrete return sensor (Teledyne Optech, Vaughan, ON, Canada).The average flying altitude was 1000 m above ground.The pulse repetition rate was of 100 kHz with a 50% strip overlap.The scanning angle ranged from 0 • to 24 • relative to the nadir.The sensor could record up to four measurements per pulse.Each field plot was covered by two to four lidar strips.The average above-ground return density was 6.4 m −2 .Information for each return was recorded in las files (version 1.0).The returns classified as ground (by the vendor) were interpolated to construct a digital terrain model using the LAStools software [6].Non-ground returns were normalized by measuring the elevation difference with the underlying terrain model.Table 2 presents a summary of the lidar data of the study site.

Data Analysis
The data analysis is organized in six sections: (Section 2.4.1)describes the nonlinear generalized least squares regression; (Section 2.4.2) provides an analysis of the spatial distribution of returns using an example field plot; (Section 2.4.3)calculates the lidar variables from the point clouds; (Section 2.4.4) describes the merchantable volume modelling; (Section 2.4.5)analyses the residual variance; and (Section 2.4.6)provides the model validation.

Generalized Least Squares in Nonlinear Regression
Heteroscedasticity occurs in a nonlinear regression (NLIN-GLS) model when the error variance is not constant over all the observations [15]: where y i is the value of the response variable measured on the ith experimental unit, and is expressed as a known function f (.) of some explanatory variables (x i ) and parameters (β) plus a random error term (ε i ).In this model, the heteroscedasticity is considered by the fact that the variance σ 2 of this latter term is not constant, but depends on a known function g(.) of some covariates z i and parameters δ.Equation ( 1) was used to predict the merchantable volume of our study site.The variance functions g(.) that we have considered had one of the following expressions [15]: The first variance function is a power function of the absolute value of the covariate z i , the second one is an exponential function of this covariate, and the others are various combinations of the previous ones.
The SP group of covariates, which will be described later, have been included in this model either in the fixed part through the f (.) function or in the random part through the g(.) function.

Analysis of the Spatial Distribution of Returns
The return density, return spacing and scanning angle play an important role in the spatial distribution of returns.To provide a visual example of this, we have extracted a sample of the 2d distribution of all the returns of a given field plot that had been covered by three lidar strips (#8, #9 and #10) during the survey (Figure 2).The strips had an average return density of 2.15, 2.19 and 2.04 m −2 respectively, an average return spacing of 6.8, 6.8 and 7.0 dm and an average scanning angle of 16 • , 1.6 • , and 9.4 • .The returns were plotted in a 2d geographic grid and a line segment was drawn between consecutive returns to show the scanning sequence.Strips #8 and #10 had a heterogeneous spatial distribution of returns.Local clumps of returns alternated with local gaps.Conversely, strip #9 had a more homogeneous spatial distribution of returns throughout the plot.When the three strips were pooled together at the plot level, the return spacing drastically decreased (4.0 dm).However, local gaps were still observable within the plot.We triangulated neighboring returns into a 3d network to analyze the spacing distribution.This was done using the Delaunay triangulation algorithm of ArcGIS 10.3 ( ® ESRI).The area of each triangle was calculated.Figure 3 shows the area distribution of triangulated returns for the previous example field plot.Strips #8 and #10 which had a heterogeneous pattern (see Figure 2) produced an asymmetric right-skewed distribution.The areas ranged from 0.2 to 211.44 dm 2 and 0.68 to 217.29 dm 2 respectively.Conversely, strip #9, which had a homogeneous pattern, produced a symmetric leptokurtic distribution.The areas ranged from 0.14 to 67.12 dm 2 .The heterogeneous pattern observed, when pooling all strips together, induced an asymmetric right-skewed mesokurtic distribution.The number of small triangles increased.The areas ranged from 0.06 to 53.14 dm 2 .This example shows that an irregular distribution of returns persisted locally even with three overlapping strips.Strips were always pooled together in the following analyses.We triangulated neighboring returns into a 3d network to analyze the spacing distribution.This was done using the Delaunay triangulation algorithm of ArcGIS 10.3 ( ® ESRI).The area of each triangle was calculated.Figure 3 shows the area distribution of triangulated returns for the previous example field plot.Strips #8 and #10 which had a heterogeneous pattern (see Figure 2) produced an asymmetric right-skewed distribution.The areas ranged from 0.2 to 211.44 dm 2 and 0.68 to 217.29 dm 2 respectively.Conversely, strip #9, which had a homogeneous pattern, produced a symmetric leptokurtic distribution.The areas ranged from 0.14 to 67.12 dm 2 .The heterogeneous pattern observed, when pooling all strips together, induced an asymmetric right-skewed mesokurtic distribution.The number of small triangles increased.The areas ranged from 0.06 to 53.14 dm 2 .This example shows that an irregular distribution of returns persisted locally even with three overlapping strips.Strips were always pooled together in the following analyses.We triangulated neighboring returns into a 3d network to analyze the spacing distribution.This was done using the Delaunay triangulation algorithm of ArcGIS 10.3 ( ® ESRI).The area of each triangle was calculated.Figure 3 shows the area distribution of triangulated returns for the previous example field plot.Strips #8 and #10 which had a heterogeneous pattern (see Figure 2) produced an asymmetric right-skewed distribution.The areas ranged from 0.2 to 211.44 dm 2 and 0.68 to 217.29 dm 2 respectively.Conversely, strip #9, which had a homogeneous pattern, produced a symmetric leptokurtic distribution.The areas ranged from 0.14 to 67.12 dm 2 .The heterogeneous pattern observed, when pooling all strips together, induced an asymmetric right-skewed mesokurtic distribution.The number of small triangles increased.The areas ranged from 0.06 to 53.14 dm 2 .This example shows that an irregular distribution of returns persisted locally even with three overlapping strips.Strips were always pooled together in the following analyses.

Lidar Variable Generation
Lidar returns were clipped to the extent of each plot spatial boundary to extract the corresponding point cloud.Metrics were computed both at the strip and the plot levels using Fusion software [7].They were classified into three standard groups of explanatory variables: canopy height variables (percentiles of height, mean, mode), canopy density variables (proportion of first returns above or below a 2 m, 5 m or 7 m threshold) and canopy structural heterogeneity variables (standard deviation, variance, coefficient of variation, canopy surface roughness index).The canopy surface roughness index, hereafter called rumple index, was computed as the ratio of the canopy outer surface area to the ground surface area as measured by the lidar-derived canopy surface and digital terrain models in a 1 m × 1 m grid [23].
A fourth group of explanatory variables, SP, was computed to describe the spatial distribution pattern of returns.Variables included in this group were: density variable (computed as the number of above ground returns of a given plot divided by its area), area distribution variables (mean, standard deviation, minimum, maximum, skewness and kurtosis of the area distribution of triangulated returns) and scanning angle variables (mean, mode, standard deviation, minimum, and maximum of the scanning angle distribution).

Establishing the Merchantable Volume Model
The basic merchantable volume model was built at the plot level using one variable per standard group iteratively.Nonlinear relationships were observed between the field MV and the lidar canopy density/canopy heterogeneity variables.Hence, each subset of variables was combined into the NLIN-GLS described in Equation ( 1): where MV is the merchantable volume, H the canopy height variable, CD the canopy density variable, CSH the canopy structural heterogeneity variable and ε the residual term.The SP variables were then added to Equation (3): The best subsets of variables were selected based on the corrected Akaike Information Criterion AICc, the residual standard deviation RSD, the level of significance of the model p-value, and the coefficient of determination pseudo-R 2 .Variable multicollinearity in the model was tested by measuring the variance inflator factor VIF with the threshold value set to 10 [24].The statistical analysis was done with R software [25].We used the "nlme" package [26] for the model fitting.

Residual Variance
The absolute residuals were regressed against the predicted MVs to verify the heteroscedasticity of the residuals.A significant p-value of the MV parameter was an indication of heteroscedasticity.The MV uncertainty was then assessed by iteratively adding an SP variable to the error variance function (random part of the model): The regression parameters were estimated using an NLIN-GLS regression.The AICc, RSD, p-value and the 95% confidence intervals (95% CI) were computed to compare models.

Model Validation
The leave-one-out cross-validation was used to assess the accuracy of the merchantable volume models.A new dataset was created by removing one field plot iteratively.The MV model was fitted to the new dataset (training data) and used to predict the MV of the removed field plot (test data).The procedure was repeated until predicted values were obtained for all the field plots.

Results
We found a very good correlation between field and predicted merchantable volumes for balsam fir stands even before the introduction of the SP variables (pseudo-R 2 = 0.91, Table 3, and Figure 4).The best predictors were the mean height of first returns (Mean_h), the proportion of first returns below 2 m (Prop_2m) and the rumple index (Ri).The VIF factor was below 10.The RSD was 28.0 m 3 ha −1 for the regression model and 28.98 m 3 ha −1 after cross-validation.The best equation was:

Results
We found a very good correlation between field and predicted merchantable volumes for balsam fir stands even before the introduction of the SP variables (pseudo-R 2 = 0.91, Table 3, and Figure 4).The best predictors were the mean height of first returns (Mean_h), the proportion of first returns below 2 m (Prop_2m) and the rumple index (Ri).The VIF factor was below 10.The RSD was 28.0 m 3 ha −1 for the regression model and 28.98 m 3 ha −1 after cross-validation.The best equation was: MV = exp(2.04)* Mean_h 1.33 * Prop_2m -0.08 * Ri 0.64 + ε (6)   Figure 5 shows plots of the model residuals versus the predicted MV, Mean_h and Skew_area variables.Residuals had an outward-opening funnel form: they were more variable as the predicted MVs increased.The MV parameter was significant (p = 0.02) when regressing the absolute residuals Figure 5 shows plots of the model residuals versus the predicted MV, Mean_h and Skew_area variables.Residuals had an outward-opening funnel form: they were more variable as the predicted MVs increased.The MV parameter was significant (p = 0.02) when regressing the absolute residuals against the predicted MVs.This indicated the presence of heteroscedasticity and needed to be corrected.
against the predicted MVs.This indicated the presence of heteroscedasticity and needed to be corrected.3)) and the lower scatter plots, after the addition (see Equation ( 5)).
The spatial distribution of returns barely influenced the predicted MV.We found no significant improvement when adding an SP predictor to the fixed part of the MV model.Table 4 shows a model comparison between the basic MV model and 3 MV models built with an additional SP variable: MV = exp (2.01) * Mean_h 1.34 * Prop_2m -0.08 * Ri 0.68 * Skew_area -0.06 + ε MV = exp (2.25) * Mean_h 1.32 * Prop_2m -0.08 * Ri 0.66 * Mean_angle -0.09 + ε MV = exp (2.00) * Mean_h 1.34 * Prop_2m -0.08 * Ri 0.63 * Density 0.02 + ε The additional SP variables did not improve the basic model (AICc = 1137.79,1137.48 and 1139.32 respectively for Equations ( 7) -( 9) versus 1137.23 for Equation ( 6)).The density variable (Equation ( 9)) was the least significant (p = 0.74).RSDs were however similar for all models.Table 5 shows a comparison between MV models before (basic model) and after addition of a variance function (corrected models).Adding SP variables to the random part of the MV model  3)) and the lower scatter plots, after the addition (see Equation ( 5)).
The spatial distribution of returns barely influenced the predicted MV.We found no significant improvement when adding an SP predictor to the fixed part of the MV model.Table 4 shows a model comparison between the basic MV model and 3 MV models built with an additional SP variable: MV = exp (2.25) * Mean_h 1.32 * Prop_2m −0.08 * Ri 0.66 * Mean_angle −0.09 + ε MV = exp (2.00) * Mean_h 1.34 * Prop_2m −0.08 * Ri 0.63 * Density 0.02 + ε The additional SP variables did not improve the basic model (AICc = 1137.79,1137.48 and 1139.32 respectively for Equations ( 7) -( 9) versus 1137.23 for Equation ( 6)).The density variable (Equation ( 9)) was the least significant (p = 0.74).RSDs were however similar for all models.Table 5 shows a comparison between MV models before (basic model) and after addition of a variance function (corrected models).Adding SP variables to the random part of the MV model improved the precision.The corrected MV models had lower AICc value (1110.32,1123.58 and 1106.72 versus 1137.23 for the basic model).Mean_h and Skew_area had the most effect on the unexplained variance as residuals were more variable when these variables increased (Figure 5, Table 6).The MV uncertainty was better assessed by adding the following variance function to Equation (6): var(ε i ) = σ 2 * Mean_h 0.85 * exp (0.31 * Skew_area) 2 (10) The RSD, which was 28.0 m 3 ha −1 initially, was estimated to 3.7 m 3 ha −1 multiplied by the variance function, when considering the heteroscedasticity of residuals.

Discussion
This study aimed to develop an efficient lidar model to predict merchantable volume in balsam fir stands and to better assess the model uncertainty.As uncertainty increases, the model parameters can still be estimated but the confidence intervals are unreliable [12,27].Uncertain predictions induce measurement errors in the volume of stands and consequently on the evaluation of their economic value [14].The uncertainty of merchantable volume should, therefore, be correctly estimated.
As expected, highly significant model for predicting the MV could be developed using standard lidar metrics (pseudo-R 2 = 0.91, RSD = 28.0 m 3 ha −1 ).This is consistent with other studies indicating that lidar variables can accurately predict forest structure attributes in the eastern Canadian boreal forest [11,28,29].The optimal explanatory variables of the model were the mean height of first returns, the proportion of first returns below 2 m and the rumple index.Mean height is a good predictor of volume as shown in many studies (e.g.: [8,11]).The variable was positively correlated to the MV in the model as shown in Table 3.
The proportion of first returns below 2 m variable could be correlated with the fraction of canopy gaps.Canopy gaps are small openings of stands where there are smaller or no trees [30].They are often found in mature stands, and reduce the overall stand volume.Several plots of our study site were within mature forest stands.The variable was negatively correlated to the MV in our model as shown in Table 3.
The rumple index measures the canopy surface roughness or structural complexity.It has also been defined as a 3d measure of canopy heterogeneity [23].It can be used to characterize multi-cohorts stands such as those present in our study site.Mature trees, which have bigger volume, are more likely to be present in multi-aged stands [31,32].The variable was positively correlated to the MV in the model as shown in Table 3.
Figure 2 shows that the spatial distribution of returns can remain irregular despite a high return density or the presence of overlapping strips.The consequences of an irregular spacing pattern are directly evident when building, for example, canopy surface models or digital elevation models [3,33,34].In our case study, the irregular pattern induced an alternation of local clumps (overscanned areas) and local gaps (areas without data) within the point cloud.This pattern did not affect the accuracy of the model as adding an SP variable to the fixed part did not improve the predictions.However, the residuals were affected by the pattern.This shows that an irregular spatial distribution of returns can increase the uncertainty of predictive models.
The variance function helped to identify the sources of heteroscedasticity of residuals.Residuals were more heteroscedastic with increases of the mean height of returns (Mean_h) and the skewness of the area distribution of triangulated returns (Skew_area); see Figure 5 and Table 6.This suggests a higher MV uncertainty in high stands (e.g., mature forest stands) and with irregular lidar scanning patterns.Natural high stands tend to have a more heterogeneous canopy height structure compared to low stands [35].An irregular scanning pattern over these stands could over or under characterize the height distribution substantially.Yet these stands are preferentially cut during forest operations.Their MV, therefore, need to be precisely predicted.Conversely, low stands, tend to have a more homogeneous canopy height structure.The over or under characterization of the height distribution due to the irregular scanning pattern would be less substantial.
We further tested the robustness of the corrected MV model for a low-density scenario like in some operational context (pulse density ≤ 2 m −2 ).We obtained similar results to the previous analyses.The RSD was 27.4 m 3 ha −1 initially and was estimated to 3.7 m 3 ha −1 multiplied by the variance function when considering the heteroscedasticity of residuals.Mean_h and Skew_area remained the best combination of covariates for the variance function.The spacing distribution of returns should also be considered when planning lidar surveys or, at least, when analyzing acquired lidar data.
Pooling lidar strips together at the plot level increased the return density (6.4 m −2 on average).However, our analysis shows that the return density did not influence the prediction of the merchantable volume estimates for the study site (p = 0.74).This result is in accordance with other studies (e.g., [10,11,36]) done even in low return density contexts.
We tested variable transformation to correct for heteroscedasticity.The response and explanatory variables of the MV model were strictly positive (see Tables 1 and 2), and three variables (MV, Prop_2m, and Ri) had a right-skewed distribution.We thus tested a square root, a cube root, and a logarithmic transformation.The heteroscedasticity was best reduced with the logarithmic transformation.The p-value for the predicted MV parameter was improved, however, it remained significant (0.01).The model had the following form: log(MV) = 1.99 + 1.18 * log (Mean_h) − 0.08 * log (Prop_2m) + 0.95 * log (Ri) + ε Acquiring lidar data with a small scanning angle would also be an ideal solution to obtain more homogeneous data but the acquisition costs would also be greatly increased.Other authors have chosen to exclude from the point clouds either pulses transmitted at high scanning angles [18,37] or irregularly spaced sectors within strips [3].However, these changes would substantially alter the point cloud structure.A statistical approach, as the one we have proposed, has the advantage of alleviating these problems while entailing no additional costs.We recommend the use of an additional spatial distribution variable when modeling forest attributes models from lidar point clouds having an irregular spatial distribution of returns.
A practical application of our study is the establishment of reliable uncertainty maps of predicted merchantable volumes.Users can then confidently assess the reliability of the estimates and consequently better plan their harvesting operations.

Conclusions
While lidar return density and scanning angle have little impact on the timber merchantable volume parameter estimates, an irregular spacing among returns influences a model uncertainty.Residuals of a standard merchantable volume predictive model were found to be correlated with the mean height of returns and a variable accounting for the irregularity of spatial distribution of returns.We included these variables in a variance function to correct for the heteroscedasticity.The residual standard deviation was better estimated (3.7 m 3 ha −1 multiplied by a variance function versus 28.0 m 3 ha −1 without the variance function.Therefore, the mean height of returns and the return spatial distributions should be considered for reliable estimates of forest attributes.

Figure 1 .
Figure 1.Locations of the Montmorency Research Forest study site and the field plots over a lidar DTM.

Figure 1 .
Figure 1.Locations of the Montmorency Research Forest study site and the field plots over a lidar DTM.
Remote Sens. 2017, 9, 808 5 of 13 spatial distribution of returns.Local clumps of returns alternated with local gaps.Conversely, strip #9 had a more homogeneous spatial distribution of returns throughout the plot.When the three strips were pooled together at the plot level, the return spacing drastically decreased (4.0 dm).However, local gaps were still observable within the plot.

Figure 2 .
Figure 2. Spatial distribution pattern of returns for the field plot #28 surveyed by three lidar strips.Dark areas depict zones of clumps of returns and white areas depict gap zones where there are no returns.

Figure 3 .
Figure 3. Area distribution of triangulated returns for the field plot #28 surveyed by three lidar strips.

Figure 2 .
Figure 2. Spatial distribution pattern of returns for the field plot #28 surveyed by three lidar strips.Dark areas depict zones of clumps of returns and white areas depict gap zones where there are no returns.
Remote Sens. 2017, 9, 808 5 of 13 spatial distribution of returns.Local clumps of returns alternated with local gaps.Conversely, strip #9 had a more homogeneous spatial distribution of returns throughout the plot.When the three strips were pooled together at the plot level, the return spacing drastically decreased (4.0 dm).However, local gaps were still observable within the plot.

Figure 2 .
Figure 2. Spatial distribution pattern of returns for the field plot #28 surveyed by three lidar strips.Dark areas depict zones of clumps of returns and white areas depict gap zones where there are no returns.

Figure 3 .
Figure 3. Area distribution of triangulated returns for the field plot #28 surveyed by three lidar strips.

Figure 3 .
Figure 3. Area distribution of triangulated returns for the field plot #28 surveyed by three lidar strips.

Figure 4 .
Figure 4. Scatter plot of predicted merchantable volume versus field merchantable volume (observed) of plots.

Figure 4 .
Figure 4. Scatter plot of predicted merchantable volume versus field merchantable volume (observed) of plots.

Figure 5 .
Figure 5. Scatter plots of standardized model residuals versus predicted merchantable volume (MV), observed mean heights of first returns (Mean_h) and skewness of the area distribution of triangulated returns (Skew_area) of plots.The upper scatter plots show the relationship before addition of the variance function to the model (see Equation (3)) and the lower scatter plots, after the addition (see Equation (5)).
(2.04) * Mean_h 1.33 * Prop_2m -0.08 * Ri 0.64 + ε Mean_h = mean height of first returns, Prop_2m = proportion of first returns below 2 m, Ri = rumple index, Density = return density, Mean_angle = mean scanning angle of returns, Skew_area = skewness of the area distribution of triangulated returns, RSD = Residual standard deviation.* The shown value represents the level of significance of the additional variable.

Figure 5 .
Figure 5. Scatter plots of standardized model residuals versus predicted merchantable volume (MV), observed mean heights of first returns (Mean_h) and skewness of the area distribution of triangulated returns (Skew_area) of plots.The upper scatter plots show the relationship before addition of the variance function to the model (see Equation (3)) and the lower scatter plots, after the addition (see Equation (5)).
(2.04) * Mean_h 1.33 * Prop_2m −0.08 * Ri 0.64 + ε Mean_h = mean height of first returns, Prop_2m = proportion of first returns below 2 m, Ri = rumple index, Density = return density, Mean_angle = mean scanning angle of returns, Skew_area = skewness of the area distribution of triangulated returns, RSD = Residual standard deviation.* The shown value represents the level of significance of the additional variable.

Table 2 .
Summary of the lidar data for the Montmorency Forest study site (119 plots) a .Std. = standard deviation, Mean_h = mean height of first returns, Skew_area = skewness of the area distribution of triangulated returns, Prop_2m = proportion of first returns below 2 m, Ri = rumple index. a

Table 3 .
Parameters of the merchantable volume models before and after the addition of the variance function a .

Table 3 .
Parameters of the merchantable volume models before and after the addition of the variance function a .

Table 4 .
Comparison of four predictive merchantable volume MV models before (basic model) and after addition of a spatial distribution variable a .
a Basic model = exp

Table 4 .
Comparison of four predictive merchantable volume MV models before (basic model) and after addition of a spatial distribution variable a .
a Basic model = exp

Table 5 .
Comparison of the merchantable volume MV models before (basic model) and after addition of a variance covariate a .Basic model = exp (2.04) * Mean_h 1.33 * Prop_2m −0.08 * Ri 0.64 + ε Mean_h = mean height of first returns, Prop_2m = proportion of first returns below 2 m, Ri = rumple index, Skew_area = skewness of the area distribution of triangulated returns, RSD = Residual standard deviation.* The shown value represents the level of significance of the additional variance function. a

Table 6 .
Effects of Mean_h and Skew_area variables on the variance of predicted merchantable volume a .
a RSD = Residual standard deviation.