Combining Multi-Date Airborne Laser Scanning and Digital Aerial Photogrammetric Data for Forest Growth and Yield Modelling

Abstract

The increasing availability of highly detailed three-dimensional remotely-sensed data depicting forests, including airborne laser scanning (ALS) and digital aerial photogrammetric (DAP) approaches, provides a means for improving stand dynamics information. The availability of data from ALS and DAP has stimulated attempts to link these datasets with conventional forestry growth and yield models. In this study, we demonstrated an approach whereby two three-dimensional point cloud datasets (one from ALS and one from DAP), acquired over the same forest stands, at two points in time (circa 2008 and 2015), were used to derive forest inventory information. The area-based approach (ABA) was used to predict top height (H), basal area (BA), total volume (V), and stem density (N) for Time 1 and Time 2 (T1, T2). We assigned individual yield curves to 20 × 20 m grid cells for two scenarios. The first scenario used T1 estimates only (approach 1, single date), while the second scenario combined T1 and T2 estimates (approach 2, multi-date). Yield curves were matched by comparing the predicted cell-level attributes with a yield curve template database generated using an existing growth simulator. Results indicated that the yield curves using the multi-date data of approach 2 were matched with slightly higher accuracy; however, projections derived using approach 1 and 2 were not significantly different. The accuracy of curve matching was dependent on the ABA prediction error. The relative root mean squared error of curve matching in approach 2 for H, BA, V, and N, was 18.4, 11.5, 25.6, and 27.53% for observed (plot) data, and 13.2, 44.6, 50.4 and 112.3% for predicted data, respectively. The approach presented in this study provides additional detail on sub-stand level growth projections that enhances the information available to inform long-term, sustainable forest planning and management.

1. Introduction

Sustainable forest management requires accurate information on both current and projected stand conditions. While the current information on stand composition, structure, and age is collected during forest inventories, long-term planning relies on models that describe forest stand dynamics [1]. We refer to these models as growth simulators throughout this paper.

Growth simulators project forest stand attributes over time and depend on a set of input attributes describing individual trees or a stand [2]. The type and number of input attributes into the models depend on the growth simulator’s level of abstraction with the whole-stand growth models requiring less detailed information than single-tree or size-class models [2]. The choice of growth simulator depends on its availability for the area of interest—e.g., in Canada growth simulators have been developed for each province separately, with whole-stand models being popular in British Columbia and Alberta. In the majority of cases, the core set of input attributes for stand-level growth simulators includes species composition, stand age, and top height or alternatively site index. Additional inputs can improve model accuracy and can include information on stand density, basal area, stocking, canopy cover, information on insect damage, and silvicultural practices such as thinning or fertilization [3,4]. Simulator output typically consists of a list of predicted stand attributes for a specified age sequence (e.g., 10, 20, 30 years in the future), stratified by species. Stand attributes include top height, basal area, merchantable and total volume, and stand density.

Airborne laser scanning (ALS) has demonstrated capacity for characterizing a number of important forest stand attributes [5,6]. Three-dimensional point clouds acquired with ALS precisely characterize the distribution of vegetation and allow for accurate estimation of forest stand attributes, including height, canopy cover, basal area, biomass, and volume [6,7]. Because of the high positional accuracy of the acquired point clouds and the ability of a laser pulse to pass through small openings in forest canopy, some attributes, such as tree height and canopy cover, are estimated directly. Other attributes, such as basal area or volume, require development of models that utilize ALS metrics as predictor variables. A multitude of studies have proven the suitability of ALS-derived forest stand attributes in variety of forest environments e.g., [8,9,10,11].

The area-based approach (ABA) and the individual tree detection (ITD) approach are the two most common approaches for ALS-based predictions of forest-related attributes. ITD requires individual tree tops to be located and tree crowns to be delineated [11,12]. Individual tree attributes are then derived from the properties of the point cloud inside each delineated crown. This method depends on the accuracy of the tree identification and is prone to errors that result from over- or under-segmentation of tree crowns. In the ABA, forest stand attributes are estimated for a grid cell (typically 20–30 m on a side), based on metrics that summarize the distribution of the point cloud within the cell [5,13]. The grid cell is therefore a fundamental unit of the ABA, offering more spatial detail than a traditional polygon-based inventory, while avoiding the bias often introduced by the ITD approach.

Digital aerial photogrammetry (DAP) approaches have been found to be a good alternative data source to ALS for stand characteristics. In some studies, DAP-based point clouds have provided similar accuracies for ABA-derived forest stand attributes like height, basal area or volume [14,15,16,17,18], while the cost of acquiring these data is markedly lower than ALS data. However, the DAP point cloud primarily represents the top of the canopy and the ground surface is often not sufficiently represented, which limits the application of DAP data in moderate to dense forest stands [19,20]. This limitation may be overcome by integrating ALS and DAP data sources; stand heights may be then characterized by normalizing the DAP point cloud elevations with an ALS-based terrain model. Due to lower acquisition costs, DAP data has been proposed as an alternative to ALS data for updating forest inventories, provided an initial ALS acquisition is available to provide the requisite ground elevations under canopy [19].

The capacity of the ALS and DAP data to characterize stand conditions offers new possibilities for analyzing stand dynamics over large areas when the datasets are acquired at different points in time. ABA-derived layers informing stand attributes at different times should allow forest managers to track stand dynamics between the two acquisition times. Having the cell-level information extended from a single, or two points in time, to the desired length of forest management and planning periods is of great importance and could provide critical information supporting sustainable forest management. However, the development of methodologies to integrate growth simulators with ABA-derived forest inventory layers is in its infancy. First and foremost, there is a discrepancy between the inputs required to parameterize growth simulators and outputs from ABA. Growth simulators typically require information on stand age, species, and top height, at a minimum, with top height being the only variable well demonstrated to be accurately predicted from ALS or DAP and therefore available at cell level. Second, most growth simulators are applied at the stand-level, whereas ABA outputs are provided at the individual grid-cell level; generalizing ABA outputs to the stand-level results in a loss of within-stand information [19,21]. Third, the influence of ABA prediction accuracy on subsequent growth and yield projections is not yet known. Lastly, it is unclear how multi-temporal ABA predictions could be used to increase the accuracy of cell-level growth projections.

Using 3D point clouds to analyze forest dynamics has significant potential for providing detailed and accurate characteristics of forest stand changes through time. Examples rely either on multiple ALS or DAP acquisitions, or combined data from ALS and DAP, and in most cases focus on stand height, as this attribute can be assessed directly with highest accuracy [22,23,24,25]. For example, Stepper et al. [26] derived periodic annual increment (PAI) of forest height from two DAP datasets acquired at two different times. Similarly, Véga & St-Onge [27] analyzed forest stand height growth based on historical aerial photographs and ALS data. Historical analog photographs were used to create digital canopy height models (CHMs) to observe changes in dominant tree height. Cao et al. [28] used multi-temporal ALS data to assess forest biomass dynamics in a mixed subtropical forest in China. Similar to Goodbody et al. [29], the change in forest stand attributes was either modeled directly, or was calculated as a difference between two separate predictions. Recently Socha et al. [30] used repeated laser scanning data to model height growth and site index. They used two point cloud datasets with the second one acquired 5 years after the first, and successfully developed site-specific trajectories of top height.

Studies that demonstrate how multi-temporal ALS or DAP data can be linked to a growth simulator for forest attribute projection are more limited. Falkowski et al. [31] used ALS-based individual tree attributes as input parameters for the Forest Vegetation Simulator (FVS) [32]. FVS is an individual-tree growth and yield simulator that relies on species and diameter at breast height (DBH) for input. The authors compared inventory- and ALS-based projections, which followed the same trends, and showed a correlation coefficient for projected basal area values above 0.91 and root mean square differences (RMSD) of 1.97 m² ha⁻¹ on average. Recently Lamb et al. [33] proposed an approach to linking individual tree growth simulators with ALS data. Instead of using ALS-derived attributes as direct inputs, they used the estimated attributes to assign a tree list to each 20 × 20 m cell using k-nearest neighbor imputation. The tree lists were then used as input to the growth simulator to project stand attributes over time.

Tompalski et al. [21] demonstrated how single-date ALS data and ABA outputs can be integrated into stand-level growth simulators to provide improved cell-level growth projections. Yield curves can be assigned to each ABA cell based on stand age and species information from existing forest inventory data at the stand level, and ABA attributes that correspond to the output of the growth simulator employed. This approach utilizes the existing inventory information to generate a database of yield curves representing all possible stand conditions in the area. Next, candidate yield curves are selected from the database based on each of the ABA attributes and the stand age indicated in the forest inventory. Final yield curves are then assigned to each of the ABA cells based on a weighted mean of the candidate curves, with the weight being the measure of prediction accuracy of the ABA models. Although not statistical in nature, this approach was found to be suitable for large area cell-level growth and yield analysis as it takes advantage of all available attributes estimated with ABA, not only top height. Moreover, the discrepancy between the candidate curves can be used to characterize the uncertainty related to the final assigned yield curve.

In this paper, we build upon the approach presented in Tompalski et al. [21] and focus on applying it to enhance forest growth projections when multiple ABA outputs, derived at different points in time are available. The unique contributions of this study are: (i) an analysis of how the accuracy of cell-level yield projections change depending on whether a single or multiple ABA outputs are used during curve assignment, and (ii) determination of the main drivers of the yield curve matching error. We use ALS and DAP point clouds acquired at two points in time (ALS at T1—circa 2008 and DAP at T2—2015) in an ABA. Four key stand attributes were first modelled (top height, basal area, total volume, and number of trees per hectare) for both the T1 and T2 datasets, and were then used to assign a yield curve to each 20 × 20 m grid cell. The differences in the accuracy of the curve assignments were compared by calculating bias and root mean square error (RMSE) between the projected and reference data. The comparison was performed for two different scenarios, depending on whether only the data acquired at T1 were used, or whether a combination of data acquired from T1 and T2 were used. We also considered how the ALS and DAP point clouds could be integrated to analyze forest growth and how the prediction errors associated with the ALS- and DAP-derived models influenced the yield curve assignment.

2. Materials

2.1. Study Area

The study area was approximately 700,000 ha in size, located near the community of Slave Lake, Alberta, Canada (Figure 1). The area is situated in the Central Mixedwood, Lower Foothills, and Upper Foothills natural subregions [34]. The mean summer and winter temperatures are 20 °C and −21 °C, respectively, while the annual precipitation is 600 mm [34]. Forests in the study area are actively managed for lumber and pulpwood and the area is also subject to extensive oil and gas exploration.

2.2. Forest Inventory Data

Alberta Vegetation Inventory (AVI) [35] data were used to provide auxiliary information on stand attributes. AVI data is generated using aerial photo interpretation to delineate homogeneous forest stands (polygons) at least 2 ha in size, supported with field plots and field visits to assign characteristics for the delineated polygons. Polygon attributes include species composition, crown closure class, height, year of origin, structure, condition, and productivity class (Good, Medium, Fair, and Unproductive) [35]. We used the polygon level AVI information on stand age and species composition for yield curve assignment.

The forest inventory on the study area contained 46,913 unique polygons or stands and represented a total area of 587,875 ha. Stands were aggregated into four groups depending on the dominant species: SB—black spruce group, SW—white spruce group, PL—lodgepole pine group, and AW—trembling aspen group (

Table 1. Forest stand characteristics in the study area based on the Alberta Vegetation Inventory.

Species Group and Code	Included Dominant Species	Number of Stands		Total Area		Mean Area	Stand Age
		#	%	ha	%	ha	Mean	σ
Aspen group—AW	balsam poplar (Populus balsamifera) trembling aspen (Populus tremuloides)	17,134	36.52	250886.20	42.68	14.64	81	38
Pine group—PL	jack pine (Pinus banksiana) lodgepole pine (Pinus contorta v. latifolia) tamarack/larch (Larix laricina)	11,830	25.22	162834.58	27.70	13.76	112	55
Black spruce group—SB	black spruce (Picea mariana)	9947	21.20	90498.25	15.39	9.10	117	47
White spruce group—SW	alpine fir (Abies lasiocarpa) balsam fir (Abies balsamea) Douglas fir (Pseudotsuga menziesii) Engelmann spruce (Picea englemannii) white spruce (Picea glauca)	8002	17.06	83656.40	14.23	10.45	120	44
Total		46,913	100.00	587875.43	100.00	12.53	103	49

).

2.3. Plot Data

Permanent sample plot (PSP) data were used to create ABA models that predict forest stand attributes of interest, as well as to evaluate the growth and yield projections based on the remote sensing data. The PSPs were originally established in the study area between 2004 and 2007 (121 plots in total at T1), and then re-visited in 2015 (45 plots at T2, Figure 1). Tree measurements on each plot consisted of species, DBH, height, height to living crown, crown position, tree condition, and stem location. All trees within the boundary of a plot exceeding 7 cm DBH were measured. Plot radius was 11.2 m (400 m²), with centers recorded with a GPS receiver of sub-meter accuracy.

Individual tree measurements were aggregated to plot level summaries that consisted of top height (H), basal area (BA), total volume (V), and number of trees per hectare (N). Plot summaries for all PSPs measured at T1 and T2 were used as reference data. However, due to the time lag between field and remote sensing data collection, which ranged from 1 to 4 years, a number of anthropogenic and non-anthropogenic disturbances occurred. The AVI information was used to flag plots that had recorded disturbances between the two time steps and these plots were then removed from the dataset, reducing the number of plots available to 98 and 35 for T1 and T2, respectively (

Table 2. Characteristics of permanent sample plots (PSPs) used to generate area-based models of forest stand attributes.

Time	Species Group Code	N	H		BA		V		N
			Mean	SD	Mean	SD	Mean	SD	Mean	SD
T1	AW	39	19.6	5.3	32.3	15.7	269.3	179.6	1169.5	629.9
T1	PL	17	12.9	5.9	22.8	19.1	162.5	183.4	1166.4	1193.8
T1	SB	26	15.4	5.2	30.6	17.5	190.4	140.1	1438.5	1058.3
T1	SW	16	21.8	7.1	36.1	18.1	304.4	193.3	873.4	697.0
T1	Total	98	17.7	6.4	30.8	17.4	235.6	178.2	1192.0	887.5
T2	AW	10	21.9	4.2	37.3	10.7	326.1	124.1	1102.0	454.4
T2	PL	9	13.8	5.7	25.6	24.0	202.1	257.4	1122.9	1061.3
T2	SB	11	16.2	4.0	37.5	16.8	230.2	123.7	1846.8	1171.7
T2	SW	5	18.7	5.7	33.7	18.5	232.2	132.7	955.7	714.5
T2	Total	35	17.6	5.5	33.9	17.8	250.6	169.5	1320.6	957.6

).

2.4. ALS Data

ALS data were acquired for the study area between 2006 and 2008, with the majority acquired in 2008. Optech ALTM 3100 was used in all three data acquisitions, with very similar flight and acquisition parameters (

Table 3. Airborne laser scanning (ALS) data acquisition parameters.

Year of Acquisition	2006	2007 and 2008
Sensor	Optech ALTM 3100	Optech ALTM 3100
Flying height	1250 m AGL	1400 m AGL
Flight speed	160 kts	160 kts
Pulse repeatition rate	50 kHz	70 kHz
Scanning frequency	30 Hz	33 Hz
Scan angle	50 deg	50 deg
Beam divergence	0.3 mrad	0.3 mrad
Average point density	1.53	1.52 (2007)/1.68 (2008)

). The average point density for all point clouds was 1.63 points/m².

2.5. DAP Data

The digital image data were acquired on 26 April, 9 May and 13 May 2015 with a Z/I DMC II 230 camera, and consisted of 1527 images. Images included blue, green, red, and near-infrared spectral bands. The ground sample distance (GSD) was 0.3 m. The along track and across track overlap were 60% and 30%, respectively. Images were processed following the standard routines for DAP that included block alignment, and dense cloud building, using Agisoft Photoscan. A total of 134 ground control points (GCPs) were used to ensure accurate registration of the data. The dense point clouds were built using the “high” setting, which resulted in an average density of 0.82 points/m².

2.6. Point Cloud Data Processing

The ALS point clouds were processed following standard processing routines, which included tiling, ground classification, and height normalization. Processing was performed using the LAStools software package [36], with a tile size of 1000 × 1000 m and a 20 m buffer. Ground classification was based on adaptive TIN models [37]. Points classified as ground were then used to normalize the ALS point cloud elevations relative to the ground surface.

DAP point clouds were processed using routines similar to the routines used to process the ALS data. Data were tiled using a similar tile size and tile overlap. However, because the DAP data primarily characterize the top of the canopy, height normalization for each tile was performed by incorporating the ALS-based points classified as ground. After normalization the ALS ground points were removed.

Point cloud metrics were calculated using both the ALS and the DAP data. Metrics included measures of central tendency (mean, median, mode), measures of dispersion (variance, standard deviation, interquartile range), percentiles, proportions, and densities of point heights above ground. All metrics were calculated using all returns above 2 m and were calculated at both the plot- and grid cell-level using the FUSION software package. A complete list of point cloud metrics generated by FUSION can be found in McGaughey [38].

3. Methods

The overall design of the study is shown on Figure 2. Using ground plot data and point cloud metrics, we first developed predictive models for four selected stand attributes, for the T1 and T2 datasets. We then generated plot-level and wall-to-wall estimates of the selected attributes and used them to match a yield curve at both the plot- and wall-to-wall level, using two different approaches. In approach 1 only T1 data was used, while approach 2 utilized both T1 and T2 data. At the plot level, both curve matching approaches were then validated by comparing the projected values with the values observed in the field. We demonstrated a result of the curve matching technique on a 2 × 2 km subset of the study area.

3.1. ABA Modeling

Multiple linear regression was used to estimate H, BA, V, and N from ALS data at T1 and DAP data at T2. Eight separate regression models were created using a stepwise variable selection approach. Models for H were based on a single ALS- or DAP-based metric describing canopy height. Modeling of BA, V, and N used a maximum of four variables following the approach of Bouvier et al. [39]. These variables should characterize stand height, heterogeneity of stand height, canopy cover, and stand vertical complexity. We tested different combinations of variables from each of these groups, as well as different variable transformations, selecting final models based on the AIC coefficient. A log-log transformation was required to achieve linearity for models predicting BA, V, and N. The modeling results were assessed using leave-one-out cross validation. We calculated the R² coefficient, bias and RMSE (absolute and relative) and used a paired Wilcoxon signed rank test to evaluate the null hypothesis that the median difference between the compared observed and predicted values was zero at α = 0.05. In case of models that required data transformation, model statistics were calculated using the back-transformed variables with bias correction factors [40].

3.2. The Growth and Yield Projection System (GYPSY)

GYPSY (Growth and Yield Projection System) is a forest stand-based growth simulator developed by the Alberta Agriculture and Forestry Ministry of the Government of Alberta, Canada [4] for projecting the growth of the four main species groups found in Alberta. GYPSY consists of several sub-models that predict top height, density, basal area increment, and gross total and merchantable volume, with an optional model for percent stocking used for young, post-harvest stands. Detailed information on the models and accuracies can be found in [41]. The sub-models utilize stand-level input variables to project forest stand attributes. The required input data includes species group, site index (or alternatively top height and total age), stem density, and basal area. Outputs consist of estimated top height, basal area, volume, stem density, site index, and increments in these attributes, at specified times intervals over a specified period. We assumed that errors introduced by GYPSY’s internal sub-models are equal for each projection and consequently we did not account for them in our study.

Following the approach of Tompalski et al. [21], we used the simulator to create a database of yield curve templates based on all possible input combinations (e.g., every combination of species groups, top height, total age, density, and basal area). The database represented all possible stand conditions in the study area and was based on the range of stand attributes in the existing forest inventory. The yield curves were generated for the four specified species groups, from 1 to 200 years, by 1 year increments. An individual yield curve template consisted of four sequences of values representing top height, basal area, volume, and stem density, between 1 and 200 years, and being estimated by the simulator based on species group, top height, total age, and basal area.

3.3. Yield Curve Matching

A yield curve from the database was assigned to each of the ground plots and to each of the 20 m cells corresponding to the wall-to-wall ABA raster data outputs for H, BA, V, and N (Figure 3). The selection process to assign the best possible yield curve that matched stand attributes was based either on observations recorded at T1 (approach 1), or based on the combination of observations at both T1 and T2 (approach 2). First, candidate curves were selected from the database, based on stand age, species group and the minimal difference between the stand attribute and a value of a yield curve. This resulted in four candidate curves for approach 1, and eight candidate curves for approach 2 (four from T1 plus four from T2). The final yield curve was derived by calculating a weighted mean of the candidate curves, with the percent of explained variance in the ABA model used as a weight. This allowed assigning more weight to the attributes that were modeled with greater accuracy and increased the reliability of the final curve. When T1 and T2 data were used (approach 2), the weighted means were calculated separately for each data set, and then averaged. The range of attribute values projected to 80 years on each of the candidate curves was used as a measure of uncertainty for the curve assignment. The uncertainty of each attribute x was calculated as the absolute value of the relative difference between a value of attribute x on a matched curve with a value of an attribute on the curve with the largest difference (x_max):

$$\mathsf{\Delta }x=\left|\frac{x_{max}-x_C}{x_C}\right|\ast 100$$

(1)

The curve assignment procedure resulted in a yield curve and uncertainty value assigned to each of the plots or each of the grid cells. This allowed us to project the stand attributes to a chosen age of interest.

3.4. Validation of the Yield Curve Assignment

To compare the curve matching accuracy based on the different sets of input data, the curve matching procedure was applied to four scenarios. Each scenario was based either on observed or predicted plot attributes, following either approach 1 or approach 2. In each of the four cases, the matched yield curve was used to project the stand attributes from T1 to T2, and then compared to the values observed at T2. We calculated absolute and relative bias, and absolute and relative RMSE, and determined whether the differences between the two compared plot attributes were significant with the paired Wilcoxon test (α = 0.05) under each scenario.

We performed an exploratory data analysis of the ABA modeling error, the error in curve assignment, and the uncertainty of the yield curve assignment to better understand how these errors were related to one another. Spearman correlation coefficients were calculated for each variable pair (e.g., curve matching error of H and prediction error of H, or curve matching error of V and uncertainty of curve assignment of V), and relationships between variable pairs assessed using scatterplots. Multiple linear regression was used to determine the drivers of the curve matching error, which involved initially identifying a set of factors that could potentially influence the accuracy of curve matching. These factors were the prediction error, uncertainty of the yield curve assignment, stand age, and species. The absolute value of the curve matching error was regressed against the uncertainty of curve assignment, absolute prediction error of the ABA models, age, and species group, for each stand attribute and for both approaches.

3.5. Wall-to-Wall Growth and Yield Projections

The wall-to-wall raster data predicted by ABA, together with auxiliary inventory information on stand age and species group, were used to assign a yield curve to each 20 × 20 m ABA cell. The auxiliary forest inventory data were first converted to raster layers representing stand age and species code, with a cell size and extent that matched the ABA raster layers. Yield curve matching was done using the ABA_T1 predictions only (approach 1) and a combination of ABA_T1 and ABA_T2 (approach 2). Yield curves for each of the attributes (H, BA, V, N) were stored as a separate raster stack, with layers representing attributes at stand ages between 1 and 200 years.

4. Results

4.1. ABA Modeling Results

The proportion of explained variance (R²) for the ABA models derived from the ALS (T1) and DAP (T2) data was similar, with the lowest values for ABA models predicting N (0.275 and 0.246, for T1 and T2, respectively) (

Table 4. Predictive models for top height (H), basal area (BA), total plot volume (V), and stem density (N). Dependent variables denoted with T1 were modeled using ALS data acquired at T1 (2008), while T2 variables were modeled using DAP data from T2 (2015).

Dependent Variable	Predictive Model	R2
HT1	3.344 + 0.814 ∗ P99	0.76
BAT1	${10}^{0.96+0.81\ast SD+0.15\ast STRATA\_10\_20}$	0.71
VT1	${10}^{1.47+1.24\ast SD+0.13\ast STRATA\_10\_20}$	0.78
NT1	${10}^{3.37-0.91\ast P99+0.34\ast perc\_above\_2+0.22\ast STRATA\_10\_20}$	0.26
HT2	$5.22+0.81\ast P99$	0.83
BAT2	${10}^{0.26+1.30\ast P90-0.13\ast perc\_above\_2}$	0.60
VT2	${10}^{0.46+1.88\ast P95-0.17\ast perc\_above\_2}$	0.72
NT2	${10}^{2.23+0.86\ast P90-0.08\ast perc\_above\_2}$	0.25

R²—adjusted coefficient of determination; P90, P95, P99—90th, 95th, and 99th percentile of normalized point heights, respectively; SD—standard deviation of point heights; STRATA_10_20—proportion of points between 10 and 20 m above ground; perc_above_2—percentage of points above 2 m above ground.

, Figure 4). A power regression approach was used for modeling BA, V, and N, and the predicted values were back-transformed incorporating a bias correction factor to correct for the systematic bias introduced by the log transformation. Model statistics (bias, RMSE, R²) were calculated using the back-transformed values in case of models that required variable transformation Models for H, BA, and V, for T1 and T2, had R² values over 0.7, with the highest for H_T2 (0.827). The lowest R² value among these three variables was for BA_T2 (0.596). The correlation coefficients between independent variables were always lower than 0.7 for models with more than one independent variable. All independent variables remaining in each of the models were significant (α = 0.05).

There were similar prediction accuracies for the T1 and T2 models (Figure 4). In both cases, the predictions of H had the lowest absolute and relative RMSE and no bias, while predictions of N had the poorest accuracy (RMSE% > 65%). Predictions of BA and V showed small bias and had a relative RMSE of around 40%. Results of the paired Wilcoxon test (p-values reported in Figure 3) indicated that the null hypothesis that the median value of the difference between field-measured and predicted plot volume was 0 could not be rejected in all but one case (for V_T1 p-value = 0.044).

4.2. Yield Curve Matching

The four approaches to yield curve matching resulted in four different curves being assigned to each plot. Yield curves are shown, together with the T1 and T2 values used during the weighting process for three representative plots, in Figure 5. The three selected plots demonstrate cases when the error in curve matching was low or medium (plots 80 and 88), and when there was a larger discrepancy between the curves matched based on observed and predicted data (plot 58).

The greatest agreement between plot-level stand attributes and matched yield curves was achieved when the observed ground plot data were used as the input information for the matching (Figure 6A,C). The accuracy decreased when predicted attributes were used instead, with the relative RMSE increasing in some cases by more than 25% (Figure 6B,D). The accuracy was higher when both T1 and T2 data were used; however, the differences in bias or RMSE values were not as large as that between scenarios that differ by data type. The relative RMSE values for the scenario that used predicted data at T1 and T2 were lower by 7.28%, 0.11%, 1.55%, for H, BA, and V, respectively, and higher by 8.21% for N, when compared to relative RMSE for the scenario that used predicted data at T1 only. The largest discrepancies were observed for N in scenarios that used predicted values during curve matching, while the smallest discrepancies were found for H when T1 and T2 data were used.

4.3. Uncertainty in the Yield Curve Matching

There was a strong correlation between the model prediction errors and the error in curve matching, with Spearman correlation coefficient values exceeding 0.6 for the majority of the stand attributes and data types (Figure 7). Conversely, there were no strong relationships between prediction error and uncertainty, and between curve matching error and uncertainty (R < 0.3 for the majority of stand attributes).

Multiple regression analysis of the curve matching error confirmed that the prediction error of the ABA models was the main predictor of the matching accuracy (

Table 5. Multiple regression analysis of the curve matching error.

	H		BA		V		N
Approach	1	2	1	2	1	2	1	2
R2	0.29	0.12	0.66	0.72	0.64	0.81	0.83	0.81
model term:
prediction error	0.21	0.47 *	0.79 ***	0.72 ***	0.88 ***	0.81 ***	0.85 ***	0.85 ***
uncertainty	0	0	−0.01	−0.03 *	−0.01	−0.02	−0.06	0.19
age	−0.02	−0.01	0.03	0.04 *	0.26	0.11	0.3	1.7
species-PL	1.27	0.09	1.62	7.96 **	32.56	25.38	19.14	153.8
species-SB	−0.37	−0.19	2.95	5.66 *	17.33	35.91 *	−168.11	−42.93
species-SW	−0.51	0.24	1.61	5.69 *	12.35	16.84	−320.85*	−196.21

H—top height; BA—basal area; V—total volume; N—stem density; R²—adjusted coefficient of determination, * p < 0.05; ** p < 0.01; *** p < 0.001.

). The prediction error was significant for all but one model, with regression coefficients between 0.47 and 0.88. The coefficient value for these models indicated the relationship between curve matching error and prediction error. For example, a value of 0.5 indicated that an increment of 1 unit in prediction error results in 0.5 unit increment of the curve matching error. Other explanatory variables were not significant in most cases, with uncertainty and age being significant only in the case of BA (T1 and T2 approach). The species variable showed significance in only two cases (V in approach 2, and N in approach 1). The coefficient of determination was high for BA, V, and N, with highest value of 0.83 for N (approach 1).

The relationship between each variable pair (i.e., each variable projected to T2 with approach 1 and with approach 2) was strong, with the correlation coefficient between 0.89 for N and 0.95 for V (

Table 6. Comparison of stand attributes projected to T2 based on approach 1 and approach 2.

Variable	R	MD	RMSD	RMSD%	p-Value
H [m]	0.93	0.02	2.45	14.49	0.93
BA [m2/ha]	0.91	−1.01	3.98	13.99	0.98
V [m3/ha]	0.95	−11.50	44.36	20.75	0.57
N [stems/ha]	0.89	−30.72	151.96	17.67	0.31

H—top height; BA—basal area; V—total volume; N—stem density; R—Spearman’s correlation coefficient; MD—mean difference; MD%—relative MD, RMSD—root mean squared difference; RMSD%—relative RMSD; paired Wilcoxon test used to calculate the p-value; n = 35 for each variable.

). The mean differences indicated that the values of BA, V, and N projected with approach 2 were lower. The relative RMSD was highest for V (20.75%). For all comparisons, the results of the Wilcoxon test indicated that the difference between the medians of the compared variables (e.g., H projected to T2 based on T1 data only, and based on T1 and T2 data) was not significant.

4.4. Analysis of the Wall-to-Wall Projections

The results of the cell-based, wall-to-wall attribute projections derived with both approaches were saved as raster stacks. A representative 2 × 2 km subset for approach 2 is shown in Figure 8. To further illustrate the growth information assigned to each cell, four exemplar cells were selected and their corresponding yield curves shown.

5. Discussion

ALS and DAP point cloud data were used to derive forest stand attributes with ABA, and then used to assign yield curves at the plot- and cell-level. It was found that the yield curve matching approach based on the combined T1 and T2 data resulted in more accurate projections of forest attributes and that the main driver of the curve matching error was the prediction error of the ABA models. The differences between the two approaches to curve matching were more pronounced when the observed ground plot data was used and less when the matching relied on predicted data. However, the relative RMSE was virtually the same in both approaches for BA and V. We did not find significant relationships between stand age, dominant species, or the uncertainty in curve assignment with the curve matching accuracy. These results indicate that the accuracy of ABA models was the most important factor in determining the accuracy of curve matching. This result confirms that there were only modest gains in accuracy from the use of ABA outputs from two time periods in this environment.

The characteristics of the GYPSY simulator, used to create the yield curve template database, dictated the choice of the ABA variables to be modelled at T1 and T2. The accuracy of the models differed, with more accurate estimates for H, and less accurate for N, which is typical for ABA. Importantly, when used as input data during curve matching, model predictions were not treated equally, but rather weighted by the amount of explained variance. This allowed the influence of potential outliers to be limited and more weight to be assigned to predictions for which there was greater confidence.

There were similar levels of accuracy for both ALS- and DAP-based ABA predictions. In both cases, the three-dimensional point data was of low density (1.5 and 0.82 points m⁻² for ALS and DAP point clouds, respectively). Perhaps as a consequence of the data density, combined with the complexity of forest environment, the accuracy of the attribute predictions based on both sources of point cloud data were comparable or slightly lower than reported in other studies e.g., [15,16,39]. Thus, it appears that future enhanced forest inventories can rely on DAP data, with the digital terrain model required for point cloud normalization provided by ALS [42,43,44,45].

Projections based on the two approaches were not statistically different when validated at the plot level. Accurate estimates of the ABA-predicted stand attributes were more important for curve matching than the use of two separate datasets representing stand conditions at T1 and T2. This finding is based on the much larger discrepancies between attribute projections based on observed and predicted data, than between projections based on observed or predicted data only. As demonstrated, there was a high correlation between the ABA prediction error and the error of curve matching, and this error was the main component driving the accuracy of the yield projections. This demonstrates that accurate wall-to-wall growth and yield analysis at the cell level can be performed based on ABA predictions representing only a single point in time, as long as the accuracy of the ABA predictions is high. Therefore, the initial investment in field data collection and ALS acquisition (e.g., higher density) is beneficial not only for the developed ABA predictions at T1 but also for growth and yield analysis based on those predictions. Including additional predictions of stand attributes at T2 during curve matching may provide more confidence in the projections as any discrepancies between stand attributes at T1 and T2 can be easily identified; however, the accuracy of projected stand attributes will increase only slightly. We suspect that the small improvement in the yield curve matching accuracy for scenarios based on combining T1 and T2 data was related, in part, to the relatively short time period (for this forest environment) between the ALS and DAP data acquisitions. Improvement in the accuracy of projected attributes relative to what is reported herein may be greater in forests with more rapid growth rates. In the mixedwood boreal, improvements in the accuracy of projected stand attributes should increase when either the time between two acquisitions is longer, or when additional datasets that extend the time series of point cloud data are incorporated.

The flexibility of the yield curve matching approach allows for the use of multiple ABA-predicted forest inventory layers representing several points in time of stand development. The results presented herein suggest that a longer chronosequence of cell-level stand attributes estimated with the ABA would likely result in a higher accuracy of matched yield curves. Initial ALS acquisition at T1 and multiple DAP datasets collected every 3–5 years would not just update the inventory; it would also provide more information on stand dynamics and therefore improve growth and yield projections by increasing the accuracy of yield curve assignment. In addition, a data assimilation approach could be used to increase the accuracy of forest attributes at each step by combing the previous and new field observations [46,47]. However, this would require a growth model to be developed that in addition to projecting forest stand attribute would also inform on the accuracy of that projection. When longer time series of point cloud-based inventories are available, cells located in undisturbed and healthy stands may be used to improve existing growth and yield projections or to describe response to thinning, fertilization, disease, or insect infestation, by characterizing discrepancies from the matched curves. Most importantly, the yield curve matching approach allows the same level of spatial detail provided in the ABA to be maintained. The cell-level information on the growth patterns of the key stand characteristics provides opportunities to analyze dynamics at a sub-stand level and can be valuable information for sustainable forest management and planning.

Finally, it is important to recognize that there are more factors to be taken into consideration when using forest growth simulators and utilizing the derived projections in forest management and planning. This includes information related to forestry (e.g., fertilization, diseases, thinning) but also related to socio-economic values (e.g., value of wood products, job growth, harvest rates). Characterization of these additional pieces of information is beyond the capacity of ALS or remotely sensed data in general. However, the accurate, timely, and detailed forest inventory information that can be generated from ALS data, when used in growth simulators, can provide improved projections of forests in the future, which has value for forest planning and management. The approach demonstrated herein is a step towards an improved integration of ALS-derived forest inventory attributes with growth simulators. The cell-level enhanced growth and yield projections can be used to generate dynamic treatment units using spatial optimization methods to aggregate cells into larger cutting areas more suitable for forest planning than individual cells [48]. The increased level of spatial detail in the generated projections allows for the optimization of the size and shape of the treatment units depending on the desired forest planning outcome and management objectives.

6. Conclusions

Detailed information on stand dynamics is of high importance for forest resource management. Given the increasing availability of pre-existing ALS data to complement lower cost acquisitions of DAP data, integrating the two datasets to enhance forest growth and yield modeling becomes possible. This study explored the use of point cloud-derived forest inventory information, representing two points in time, to find the best matching yield curve at a 20 by 20 m cell level to support growth simulation modelling. ALS- and DAP-based point clouds were used to represent stand conditions at T1 and T2, respectively. ABA was then used to predict H, BA, V, and N, for both datasets. The yield curve assignment relied on GYPSY; however, the yield-curve matching was based on all predicted forest stand attributes, not just top height, age, and species. We analyzed how the error in curve matching changed, depending on whether T1 only, or both T1 and T2 ABA-derived stand attributes were used. The results showed that when the two datasets were used, the yield curves were matched with higher accuracy; however, because the accuracy of curve matching depended on ABA prediction error, accurate estimates of cell-level stand attributes were critical. Projections based on yield curves matched with the multi-date approach 2 were lower than those using a single date of ALS data exclusively (approach 1). Enhanced growth and yield analyses based on combined sources of point cloud data allow forest managers to benefit from additional detail and greater confidence in yield projections. Sub-stand analysis of forest dynamics provides opportunities to improve forest operations, specify more accurately areas that require management interventions, characterize uncertainty, and improve long-term planning.