Comparison of Canopy Height Metrics from Airborne Laser Scanner and Aerial/Satellite Stereo Imagery to Assess the Growing Stock of Hemiboreal Forests

Abstract

This study compared the canopy height model (CHM) performance obtained from large-format airborne and very high-resolution satellite stereo imagery (VHRSI), with airborne laser scanning (ALS) data, for growing stock (stand volume) estimation in mature, dense Latvian hemiboreal forests. The study used growing stock data obtained by ALS-based individual tree detection as training/reference data for the image-based and ALS CHM height metrics-based growing stock estimators. The study only compared the growing stock species-specific area-based regression models which are based solely on tree/canopy height as a predictor variable applied to regular rectangular 0.25 and 1 ha plots and irregular forest stands. This study showed that ALS and image-based (IB) height metrics demonstrated comparable effectiveness in growing stock prediction in dense closed-canopy forests. The relative RMSEs did not exceed 20% of the reference mean values for all models. The best relative RMSEs achieved were 13.6% (IB) and 15.7% (ALS) for pine 0.25 ha plots; 10.3% (IB) and 12.1% (ALS) for pine 1 ha plots; 16.4% (IB) and 12.2% (ALS) for spruce 0.25 ha plots; 17.9% (IB) and 14.2% (ALS) for birch 0.25 ha plots; 15.9% (IB) and 18.9% (ALS) for black alder 0.25 ha plots. This research suggests that airborne imagery and, accordingly, image-based CHMs collected regularly can be an efficient solution for forest growing stock calculations/updates, in addition to a traditional visual forest inventory routine. However, VHRSI can be the fastest and cheapest solution for monitoring forest growing stock changes in vast and dense forestland under optimal data collection parameters.

1. Introduction

Growing stock or stand volume per hectare (m³ ha⁻¹) is an important forest inventory attribute used to assess wood resources and carbon content, and monitor changes in the forest stands [1]. Furthermore, forest inventories implemented at various spatial scales continuously provide accurate, quick and up-to-date information on forest management and monitoring for national needs [2]. Therefore, besides the traditional labour-intensive and expensive methods of assessing forest resources using classical photointerpretation/ stereoplotting, visual assessment, ‘in situ’ instrumental measurements, modern remote sensing (RS) tools and techniques have come to play a crucial role in forest inventory since the beginning of the 2000s [3,4].

The prediction models based on tree biomass allometry functions (diameter at breast height (DBH), tree height) provide the most accurate forest stand growing stock/biomass estimation using traditional ‘in situ’ methods [5,6]. In their turn, modern RS tools such as Light Detection and Ranging (LiDAR) and image-based structure from motion (SfM) [7] allow deriving a tree/canopy height as a variable with a high predictive value for forest stand/plot growing stock estimation, based on a strong relationship between tree height and stem volume [8]. Generally, individual tree height and canopy height models (CHM) are calculated by subtracting the ground surface height from the digital surface model (DSM) obtained from dense three-dimensional (3D) point clouds using RS techniques [9]. Due to the ability of the laser pulse to penetrate the forest canopy, LiDAR, or ALS, is a primary source of geospatial data for forest inventory. It provides detailed information on the canopy height metrics derived from calculated CHMs [10]. SfM and photogrammetric matching techniques such as semi-global matching (SGM) [11] represent another widely used alternative methodology for deriving dense 3D point clouds and DSM calculations using stereo imagery.

The SfM/photogrammetry approach can be subdivided depending on the RS sensor platform used: unmanned aerial systems (UAS) or drones, large-format digital aerial photogrammetry (DAP) (UltraCam, DMC, ADS family sensors) and very high-resolution (<1 m) satellite-derived stereo imagery (VHRSI) [12]. However, UAS is impractical for forest inventories over large areas (more than 50 km²) and is thus outside the scope of this study. In the last decade, numerous studies have compared the effectiveness of DAP and ALS for assessing forest inventory attributes (height, growing stock, biomass, etc.) through a wide geographical range of forest environments [4]. Despite the conclusion of these studies that forest inventory attributes can be estimated slightly more accurately by ALS than by using DAP point clouds, the DAP results are comparable and acceptable, depending on forest density and structure, species diversity and geographical locations [13,14].

The number of studies using VHRSI data related to forest inventory parameters estimation is still relatively small compared to ALS and DAP (

Table 1. Recent studies used stereo VHRSI data for estimating forest growing stock/biomass. The results refer to canopy height metrics as predictor variables obtained from calculated image-based CHMs. NFI = National Forest Inventory.

Studies	(Pearse et al., 2018) [15]	(Fassnacht et al., 2017) [3]	(Immitzer et al., 2016) [16]	(Kattenborn et al., 2015) [17]	(Persson, 2016) [18]	(Straub et al., 2013) [19]	(Yu et al., 2015) [20]
VHRSI Sensor	Pleiades (FB)	WorldView-2	WorldView-2	WorldView-2	Pléiades (FNB)	WorldView-2	WorldView-2
Study location (terrain range -m)	New Zealand (200–780)	Germany (flat)	Germany (210–490)	Germany (110–140)	Sweden (125–145)	Germany (570–710)	Finland (125–185)
Dominant Tree Species	Pinus radiata	temperate forest: Scots Pine, European beech, oaks	temperate forest: European beech, oaks, Scots Pine	temperate forest: Scots Pine, European beech, oaks	boreal forest: Scots Pine, Norway spruce, birch	temperate forest: spruce, beech	boreal forest: Scots Pine, Norway spruce
Reference field data plots	circular (0.06 ha) plots (n = 195)	Point-Centred Quarter plot-less method (n = 80)	NFI stands (n = 92) probability-based sampling	circular (r = 35 m) plot cluster (n = 101)	circular plots (r = 10 m) (n = 326; 219)	rectangular plots 1 ha (n = 228)	rectangular plots (1024 m2) (n = 91)
Estimated forest parameter	Stem volume (m3 ha−1)	Aboveground biomass (Mg ha−1)	Stem volume (m3 ha−1)	Aboveground Biomass (Mg ha−1)	Stem volume (m3 ha−1)	Stem volume (m3 ha−1)	Stem volume (m3 ha−1)
VHRSI R2	0.7	0.64	0.53	0.57	0.72; 0.70	0.615	0.92
VHRSI RMSE (% of mean)	22.2	20.0	31.7	18.3	27.9; 30.3	44.4	15.9
ALS R2	0.72	n/a	n/a	n/a	0.84; 0.76	0.708	0.92
ALS RMSE (% of mean)	21.1	21	n/a	n/a	21.1; 27.5	38.02	15.91

). Generally, these studies confirm the effectiveness of VHRSI height metrics in growing stock/biomass estimation compared to ALS data. However, the potential application of VHRSI data to measure important forest inventory characteristics needs further investigation [3].

In Latvia, corresponding to National Forest Inventory (NFI) data administered by the Latvian State Forest Research Institute (LSFRI), the “Silava” forest land areas occupy 3.39 million hectares (53% of state territory). In Latvia, traditionally, forest management inventories are carried out using visual inspection probability-based sampling, typically supported by instrumental measurements (e.g., relascope) in selected plots, and aggregated to provide stand average conditions [21]. However, a drawback of visual estimation is that it is difficult to standardise the data quality, leading to systematic and random errors [22]. According to Grīnvalds (2014) [21], who compared the growing stock determined by an inventory procedure adopted in Latvia with accurate measurements (callipers and harvester data), the visual-based forest inventory underestimates growing stock by an average of 17.6%.

The main factor constraining the application of ALS data for accurate forest inventory in Latvia is the data acquisition cost and irregular data updating. For example, in Latvia, the LiDAR data were acquired only once between 2013 and 2019. At the same time, since 2007, aerial imagery with a ground sampling distance (GSD) of 0.25 m has been regularly collected (with a three-year cycle) to perform orthophoto cover of the entire territory. However, at the moment, the image-based 3D canopy height information is not used in the forest inventory. Potentially, VHRSI data can also help collect forest inventory data. In addition, there is a lack of research on the effectiveness of using image-based RS data for growing stock estimation, considering Latvia’s dense forest canopy structure and tree species. Thus, the primary aim of this study is to compare the canopy height model (CHM) performance obtained from stereo large-format airborne and very high-resolution satellite imagery (VHRSI), with ALS data, for growing stock estimation in mature, dense Latvian hemiboreal forests.

There are two main approaches to estimating forest parameters from RS data: (1) the area-based approach (ABA) based on the use of statistical canopy height distributions (e.g., mean canopy height and canopy height percentiles) in selected stands or plots; and (2) the single-tree or individual tree detection (ITD) approach, based on tree-level attribute estimation and individual tree crown delineation. For image-based RS data comparison, the study used the ABA approach because the image-based CHMs only represent the non-transparent outer “canopy blanket” cover of dominant trees, making it challenging to use ITD methods [23]. This study compares the RS data performance to estimate growing stock in two ABA scenarios: (1) in regular rectangular sample plot sizes (0.25 and 1 ha) and (2) in forest stands. This study defines forest stands as homogeneous and spatially continuous parts of a forest, uniform in tree species composition, structure, age and distribution. The sample plots are fixed-sized forest parts regularly distributed inside forest stands.

The study only compares the growing stock species-specific ABA models based solely on tree/canopy height metrics as a predictor variable, excluding DBH. This was chosen because the image-based and ALS data do not allow measurement of DBH directly, unlike tree/canopy height. Another reason was that the calculation of growing stock is practically fast and straightforward when using one variable (e.g., canopy height metrics) directly obtained from the RS data.

2. Materials and Methods

2.1. Methodology Overview

Details related to the RS data used (Section 2.3) and a detailed description of the approaches (Section 2.4) used in the given study follow this review of the methodology. The overview of the study workflow is presented in Figure 1.

Three main issues influenced the estimation of the growing stocks in the study area. The first is the lack of accurate reference field data to train/test models for growing stock prediction based on CHM height metrics. The available forest inventory (FI) data based on visual probability-based sampling were rejected due to the complexity of data quality control and systematic errors [21]. Therefore, the FI data were used only to filter forest stand polygons and for the verification of ALS-based individual tree detection results. For this reason, the ABA models’ stands/plots reference values were calculated based on the local canopy maxima ITD method applied on ALS-based CHM, and the allometry between stem volume and tree height [24]. Then, the total growing stock of every established reference stand/plot was calculated as the sum of the stem volumes of all detected trees inside the given stand/plot, Figure 2. Thus, the study used stand/plot growing stock obtained by ALS-based individual tree detection as training/reference data for the image-based and ALS CHM height metric-based growing stock estimators [25,26]. The calculated trained reference stand growing stock values were compared and then verified with the inventory data by linear regression analysis.

For a successful ITD approach, the allometric relationships between stem volume and tree height must exist. Thus, the second problem associated with this study was the lack of allometric equations to estimate stem volume solely by tree height. This challenge has been solved using available harvest data, which automatically calculates the merchantable volume of individual felled tree trunks. Thus, the study established four allometric equations between field-obtained (harvest data) individual tree stem volume and tree height, calculated separately for four dominant tree species.

The main issue is a lack of accurate field/reference data for the growing stock model validation and calibration. Thus, the derived ITD and ABA models have some level of uncertainty that cannot be accurately determined. Mainly, the uncertainty comes from (1) established stem volume–tree height allometry from harvest data due to the stem taper coefficient, (2) systematic errors in ALS-based height predictions for single trees and (3) ‘local maxima’ tree detection performance (omission and commission errors) and the time lag between ALS (2017) and imagery (2020). As a result, the trained reference growing stock values are biased in different directions depending on the species. Therefore, the study’s statistical models were used to infer precision but not bias, only looking at the spread of values relative to the predicted mean [27]. Furthermore, statistical model separation by species allowed the bias to be kept more or less constant (systematic) for species-based regression model comparison during consistency/correlation analysis. As a result, the study does not present the regression equations, as they need to be calibrated against ground truth, concentrating only on RS data comparison analysis.

Since ALS data collection was carried out three years before the satellite imagery acquisition, the ALS-extracted trees’ peak heights were adjusted (increased) using trees’ annual growth rates to compensate for changes in tree heights due to a time shift. The change detection associated with forest clear-cutting was conducted throughout the study area. The clear-cutting mask was generated using GeoEye1 (GE1) by Digital Globe (USA) stereo satellite imagery, orthophoto and DSM data obtained from airborne UltraCam by Vexcel Imaging (Austria) imagery.

The study used the Photomod v7.1 (JSC Racurs) software package for photogrammetric chain processing, including image-based DSM generation. All ALS point cloud processing works, including DTM and DSM calculations, local canopy maxima routines and canopy metrics collection, were performed in FUSION v4.30 [28]. All steps related to image-derived CHM grid processing, CHM height metrics collection and Geographic Information System (GIS)-based analysis were fulfilled by SAGA GIS v7.9 [29] and open-source QGIS v3.22 [30] software.

2.2. Study Area

The “Taurkalne large forests” forestland area is located in the central-southern part of Latvia (56°25′N, 24°55′E), Figure 3. About 150 km² of the study area is flat terrain with a relief range of 40 to 70 m (orthometric height) and an average annual rainfall of 690 mm. It represents a typical hemiboreal forest land in Latvia with mostly closed-canopy, dense and mature deciduous and evergreen trees. The focus of this research was on the region-dominant tree species of deciduous Bétula pendula (Silver birch), Álnus glutinósa (Black alder), evergreen Pinus sylvestris (Scots pine) and Pícea ábies (Norway spruce).

2.3. Data Overview

2.3.1. Forest Stand Choice and ALS Data

The JSC “Latvia’s State Forests” provided the study area’s harvest and forest inventory data. This study used visual inspection-based 2019 forest inventory (FI) data to perform the forest stand polygons selection. The FI stand data included information about dominant and co-dominant tree species and proportion, stand area, age, canopy height and estimated growing stock. All FI polygons were classified into four main forest stand groups according to the tree species composition coefficient. Only forest stands with more than 75% of the corresponding dominant tree species were selected (

Table 2. Summary statistics of study-selected forest inventory stands, separated by the four dominant tree species.

Tree Species	No. Stands	Stands Area (ha)			Age (Years)	Canopy Height (m)	Growing Stock	Basal Area (m2 ha−1)
		Min	Max	Mean	Mean	Mean	Mean	Mean
Pine (Pinus sylvestris)	687	0.14	25.6	1.8	88	23	273	24
Spruce (Pícea ábies)	291	0.13	8.1	1.3	55	21	240	24
Birch (Betula pendula)	145	0.08	11.7	1.5	49	20	197	18
European black alder (Álnus glutinósa)	128	0.18	9.4	1.6	66	23	270	23

).

The Latvian Geospatial Information Agency (LGIA) provided open access ALS data acquired by MGGP Aero (Poland) on 20 May 2017 in leaf-on conditions. The Riegl LMS Q680 full-waveform sensor was used for ALS data collection (four returns recorded), with ~400 kHz pulse rate, 800 m average flying height and flying speed of ~240 km/h, scan angle ± 22.5 degrees and sidelap 30%. The obtained ALS data (~6 pts m⁻²) processing and classification of ground points were performed using the Terrasolid software.

2.3.2. Stereo Satellite and Airborne Imagery

The study used two sets of stereo images obtained in the summer of 2020 using an optical satellite (GeoEye-1) and a large-format airborne camera (UltraCam Eagle Mark 1),

Table 3. The main characteristics of stereo pair satellite (GE1) and stereo airborne (UltraCam) imagery used in this study. All angular characteristics are indicated in degrees.

Sensor	Image ID/Overlap	Acquisition Date/Time	Delivered Bands	Product GSD (m)	Elevation Angle	Convergence Angle	Base/Height Ratio	Sun Elevation/Azimuth
GE01	10500500ACDD9F00, 10500500ACDD9E00	7 August 2020 12:27 PM	PAN, R,G,B,NIR	0.5 2	64 73.8	14.7	0.28	48.5/159
UltraCam	80% along-track 40% across-track	13 July 2020 9:30 AM	R,G,B,NIR	0.25	-	6.2	0.11	35/105

. The sensor-corrected in-track GeoEye-1 (GE1) stereo pair with an OrthoReady Stereo (OR2A) processing level was delivered by Digital Globe (USA). In addition to satellite data, twenty 4-band (NIR-R-G-B) large-format aerial images with an average resolution of 0.25 m (GSD) were obtained using the UltraCam camera. The images were acquired by Georeal (Czech Republic) and provided by the LGIA. The photos with a focal length of 100.5 mm were taken at 4600 m AGL. The UltraCam images formed a rectangular two-strip block with an 80% along-track and 40% across-track overlap. While the GE1 stereo satellite images covered the study area entirely, the UltraCam airborne stereo imagery block had only 26% (40 km²) coverage (Figure 3).

2.4. Methods

2.4.1. Stem Volume Allometry Equations from Harvester Data

The harvest data were associated with stands of final-felled forests across the study area. However, the harvester only registered machine positions resulting in stripe-patterned GPS locations on the site, so the exact geolocation positions were unavailable. Thus, the harvest data were only used to establish allometric equations for each species between field-derived individual merchantable tree stem volume and corresponding tree height. Harvesters automate every tree stem’s dimensions and the subsequent merchantable stem volume estimates. However, the unused treetop is not measured [31]. Thus, the study used a stem taper coefficient to predict an entire tree height for every species. The stem taper coefficients were calculated using the over-bark stem taper curve and average felling canopy height obtained by forest inventory data. Finally, the stem volume species-specific allometric models using individual tree height as a solely independent variable were created during regression analysis.

2.4.2. Imagery Orientation and Co-Registration with ALS Data

External GE1 imagery orientation (least-squares bundle adjustment) was performed using provided rapid positioning capability (RPC) data based on an empirical model [32] and measured ground control points (GCP). Altogether, the eighteen well-identified artificial GCPs such as road intersections, concrete slab corners, poles and well-identified natural (e.g., tree stumps) points were carried out from the ALS data to reach the best co-registration of the stereo imagery with ALS. All GCPs were well distributed over the study area and manually measured using stereo mode. However, due to the UltraCam imagery’s partial coverage of the study area, only eight GCPs were used for its final least-squares bundle adjustment. As a result, the optical sensors’ geo-referencing bundle adjustments Appendix A,

Table A1. Optical sensors’ geo-referencing bundle adjustment results based on GCP measurements (root mean square errors (RMSE)).

Optical Sensors	Nr. GCPs	RMSE X (m)	RMSE Y (m)	RMSE Z (m)
GeoEye1	18	0.31	0.35	0.33
UltraCam	8	0.19	0.16	0.18

, were accurate and robust.

2.4.3. Image-Based and ALS-Based CHM Calculations

Two image-based DSM models, GE1 (PAN 0.5 m GSD) and UltraCam (NIR-G-B, 0.25 m GSD), were generated using an SGM algorithm, implemented in Photomod. The derived ALS-based grid DTM (1 m GSD) was used for computing image-based CHMs by subtracting DTM from DSM.

Additionally, two raster-grid CHM models were created from the ALS point cloud. The first ALS-based grid CHM was generated to perform individual tree detection using the local maxima approach (see next chapter). First, the grid CHM with an optimal 0.75 m pixel resolution was chosen, based on the provided ALS point cloud average density, and assigning the highest return of the ALS points within each grid cell. Then, the 3 × 3 median filter was applied to the calculated grid CHM for the optimal result of the local maxima approach. To save the local maxima (localised peak) height values from the filtering, the Fusion software ‘peaks’ switch was used. The second ALS-based grid CHM was calculated without applying any smoothing filters and was used for calculating canopy height metrics during growing stock estimation by the area-based approach.

2.4.4. Individual Tree Detection from ALS Data

The local maxima algorithm, an effective and simple method for estimating the tree height with high accuracy [33], was applied to the ALS-based grid CHM to detect individual trees using the ‘CanopyMaxima’ tool in Fusion [28]. The local maxima algorithm uses a localised circular search window associated with the crown size to identify canopy height maxima (peaks). Therefore, the search radius was modified to be more appropriate for the vegetation density and species of this study. The search window settings were calibrated by comparing the obtained results (number of detected trees) with ~8000 manually stereo-digitized trees from the UltraCam imagery. Thus, the following optimised equations were used for individual tree detection using the ‘CanopyMaxima’ routine in Scots pine and Black alder stands (1) and spruce and birch stands (2), where width (m) = the diameter of the search window in meters, and H_CHM = canopy height:

$$width\left(\mathrm{m}\right)=1.5+0.03\ast H_{CHM}$$

(1)

$$width\left(\mathrm{m}\right)=1.5+0.02\ast H_{CHM}$$

(2)

After determining the most appropriate ‘CanopyMaxima’ settings, the individual tree detection (ITD) with a tree height threshold of 2 m was applied to the selected inventory stands (Table 2). The extracted trees’ peak heights were adjusted (increased) using the trees’ annual growth rates [34] to compensate for changes in tree heights and stem volume, due to a time shift between imagery (2020) and ALS data (2017) acquisitions. The annual growth rates were chosen based on the mean age and height of selected polygons: 0.12 (m year⁻¹) for Scots pine, 0.17 (m year⁻¹) spruce, 0.15 (m year⁻¹) birch and 0.16 (m year⁻¹) black alder. Finally, the tree stem volume of each detected tree was calculated by applying corresponding species-specific allometric equations obtained during harvest data regression analysis. The ITD results were then used to acquire stand/plot growing stock reference data for model comparison based on the area-based approach. For this purpose, the estimated total growing stock by ITD of every stand (Figure 2a) was computed as the sum of all detected tree stem volumes. Then, all growing stock values were calculated to m³ ha⁻¹ using the stand area. Finally, the final forest stand growing stock values were compared/verified with the inventory data by linear regression analysis.

2.4.5. Assessing the Growing Stock with the Area-Based Approach

The regular rectangular 1 ha and 0.25 ha plots were generated inside the corresponding forest stands across the study area covered by UltraCam and GE1 imagery to identify the most suitable canopy height metrics for growing stock estimation based on regular plot patterns (Figure 2b). The 0.25 ha plots were used as primary plots for developing prediction models based on the four main tree species. In addition, the study used the 1 ha plot size to investigate how the increased size of the plot affects the growing stock prediction. The 1 ha plots were used only for Pinus sylvestris since, for other species, the number of such plots was insufficient (less than 60) to perform regression analysis and comprehensive comparison [35].

In addition to the growing stock prediction models based on regular plot patterns, the stand-based models were chosen based on the forest inventory polygons (Table 2). The stand polygons were selected throughout the entire study site. Unfortunately, the UltraCam-based models for stand-based growing stock prediction were excluded from comparison due to the low number of selectable stands (less than 65) across the UltraCam-covered area, leaving only GE1 and ALS-based models.

Altogether, four different CHMs (ALS point cloud, raster-grids: ALS, image-based GE01 and UltraCam) were chosen to calculate corresponding height metrics for further stand/plot-based regression analysis (Figure 4):

Two different ALS point cloud-based height metric sets, including all pulse returns, were calculated: with no cut-off and 0.5 m (ground and low vegetation cut-off) height thresholds. All point cloud and grid-based canopy height metrics included the most commonly used parameters: minimum, maximum, average, median heights, height percentiles and quadratic mean canopy height (QMCH):

$$QMCH=\sqrt{\frac{{\sum }_{i=1}^n{H}_j^2}{n}}$$

(3)

where $H_j$ is the single ALS point CHM height. All growing stock regression results were related to m³ ha⁻¹ using the stand/plot area.

The reference total growing stock of every stand/plot was calculated using the ITD local maxima approach, previously described. Regression analyses were then undertaken using the obtained reference value of the growing stock for every sampled stand/plot as the dependent variable, and image-based and ALS canopy height metrics as the independent variables.

The non-linear stepwise regression was carried out to obtain the relationships between the growing stock and computed CHM height metrics of ALS and image-based models. Only one model’s best height predictor (independent variable) was chosen during stepwise analysis. Then, after testing several non-linear regression functions, the study determined the cubic (3rd degree) polynomial function with zero intercept as the best relationship describer:

$$y=ax+b{x}^2+c{x}^3$$

(4)

In all tested cases, the cubic polynomial function showed slightly better results than the exponential function (y = exp (x)), which was second best. However, the polynomial function was also chosen because it is a monotonic non-decreasing function with one critical/stationary point (0,0), more precisely using only one essential independent variable. The root means square error (RMSE), relative root mean square of the mean (RMSE%) and mean error (ME) were computed as a criterion on which models are to be judged to evaluate the height metrics’ precision and the stand/plot estimates:

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{{\sum }_{i=1}^n{({\widehat{x}}_j-x_j)}^2}{n}}$$

(5)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}\left(\%\right)=\frac{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}{\stackrel{-}{x}}$$

(6)

$$\mathrm{M}\mathrm{E}=\frac{{\sum }_{i=1}^n({\widehat{x}}_j-x_j)}{n}$$

(7)

where ${\widehat{x}}_j$ is the predicted value for stand/plot j, $x_j$ is the observed value and $\stackrel{-}{x}$ is the mean of observed values in n stand/plots.

3. Results

3.1. Species-Specific Allometry between Stem Volume and Tree Height

Species-specific allometry between merchantable stem volume and observed tree height was established using the polynomial function Equation (4),

Table 4. Regression results of species-specific allometry between merchantable stem volume and single independent variable (tree height) obtained from harvester data. R² is the coefficient of determination and RMSE% is the RMSE percentage of the species-specific stem volume mean of all observations.

	Stem Volume
	No. Obs	R2	RMSE (m3)	RMSE %
Pine	47084	0.82	0.25	30%
Spruce	42588	0.89	0.13	31%
Birch	33059	0.77	0.21	35%
Black alder	22983	0.80	0.16	33%

. The results show a strong correlation (R² > 0.77) to be sufficient for individual tree stem volume detection solely by tree height. The parameters of corresponding regression equations and correlation trendline graphs are presented in Appendix A,

Table A2. List of estimated coefficients (a,b,c) of a single independent variable (tree height) on tree merchantable stem volume prediction based on non-linear cubic polynomial regression Equation (4) with zero intercept.

	Tree Height
	a	b	c
Scots Pine	2.334 × 10−2	−2.41 × 10−3	1.10 × 10−4
Spruce	1.450 × 10−2	−1.47 × 10−3	8.00 × 10−5
Birch	3.400 × 10−4	4.00 × 10−4	3.45 × 10−5
Black alder	1.413 × 10−2	−1.29 × 10−3	7.20 × 10−5

and Figure A1.

3.2. Growing Stock Estimation Comparison by Individual Tree Detection Approach

Using allometric equations obtained during the previous step, every detected tree stem volume was calculated after applying the ‘CanopyMaxima’ approach. The comparison of the ‘CanopyMaxima’ ITD-predicted species-specific stand-based growing stock with the observed values of the forest inventory data (Table 2) is presented in

Table 5. The verification results of the linear regression of the predicted stand-based growing stock (m³ ha⁻¹) obtained using ALS-based ITD, with the observed values from the forest inventory data. Positive ME means growing stock overestimation by the ALS-based ITD model.

	No. Stands	R2	RMSE (m3 ha−1)	RMSE % of Mean	ME (m3 ha−1)	ME % of Mean
Pine	687	0.79	68	25%	25	9%
Spruce	291	0.74	73	30%	−24	−10%
Birch	145	0.84	55	28%	3.6	2%
Black alder	128	0.71	91	34%	3.2	1%

. Since the given forest inventory data was based on visual estimation (relascope based), the growing stock estimation by ALS individual tree detection approach was found acceptable for further use in the area-based method.

3.3. Growing Stock Estimation Comparison by Area-Based Approach

The ALS-based ITD results were used to acquire reference growing stock of every stand/plot used for the ABA. The ABA growing stock estimates, based on the correlation with canopy height metrics obtained in species-specific rectangular (0.25 ha and 1 ha) plots, are presented in

Table 6. Comparison of growing stock regression results in species-specific regular rectangular plots, based on a single height variable obtained from the corresponding five CHMs height metrics sets (Equation (4)). The second-best variable is in brackets. The two best models per species are marked by grey shading. R²: coefficient of determination, RMSE%: RMSE percentage of the growing stock mean across the study area, p50 … p80: height percentiles, mean: mean of canopy height, QMCH: quadratic mean of canopy height.

Height Metrics Models	Independent Variable	R2	RMSE (m3 ha−1)	RMSE %
Scots pine −0.25 ha 1087 plots			s
ALS cloud metrics–0.5 m cut-off	p50 (p60)	0.93	47.4	16.9%
ALS cloud metrics–no cut	QMCH (p80)	0.92	50.8	18.1%
ALS-based GRID 0.75 m	p60 (Mean)	0.95	44.0	15.7%
GeoEye-1 GRID 0.5 m	Mean (p50)	0.95	44.3	15.8%
UltraCam GRID 0.25 m	Mean (p50)	0.96	38.2	13.6%
Scots pine −1 ha 217 plots
ALS cloud metrics–0.5 m cut-off	p50 (p60)	0.95	38.4	13.7%
ALS cloud metrics–no cut	QMCH (p80)	0.94	39.8	14.2%
ALS-based GRID 0.75 m	Mean (p60)	0.96	33.9	12.1%
GeoEye-1 GRID 0.5 m	Mean (p50)	0.97	33.2	11.9%
UltraCam GRID 0.25 m	Mean (p50)	0.97	28.8	10.3%
Norway spruce −0.25 ha 379 plots
ALS cloud metrics–0.5 m cut-off	Mean (QMCH)	0.91	30.7	15.3%
ALS cloud metrics–no cut	p75 (p70)	0.92	27.2	13.6%
ALS-based GRID 0.75 m	p60 (p50)	0.94	24.4	12.2%
GeoEye-1 GRID 0.5 m	Mean (p50)	0.91	32.8	16.4%
UltraCam GRID 0.25 m	Mean (p50)	0.90	34.2	17.1%
Birch −0.25 ha 223 plots
ALS cloud metrics–0.5 m cut-off	Mean (QMCH)	0.89	33.3	17.5%
ALS cloud metrics–no cut	QMCH (p75)	0.90	32.4	17.1%
ALS-based GRID 0.75 m	Mean (p50)	0.93	26.9	14.2%
GeoEye-1 GRID 0.5 m	Mean (p50)	0.89	34.0	17.9%
UltraCam GRID 0.25 m	Mean (p40)	0.88	35.8	18.8%
Black Alder −0.25 ha 338 plots
ALS cloud metrics–0.5 m cut-off	p60 (p70)	0.90	53.8	19.2%
ALS cloud metrics–no cut	p80 (p90)	0.90	54.5	19.5%
ALS-based GRID 0.75 m	p75 (p80)	0.91	52.8	18.9%
GeoEye-1 GRID 0.5 m	Mean (p50)	0.94	44.4	15.9%
UltraCam GRID 0.25 m	Mean (p50)	0.93	46.4	16.6%

. All results are based on a non-linear polynomial regression (Equation (4)) with one best independent variable from the extracted canopy height metrics. All models showed a high correlation between the forest growing stock and variables based on canopy height metrics (R² > 0.88, RMSE% < 20%). The image-based plot height metrics were the best for the Scots pine and Black alder forest species, while the ALS metrics were the best for spruce and birch, respectively. The Scots pine case shows noticeably improved estimation accuracy of growing stock with increased plot size, from 0.25 ha to 1 ha.

In addition to the regular rectangular (0.25; 1 ha) plot-based evaluation, the forest stand-based growing stock estimation was chosen for comparison (Figure 2a, Table 2). The stand-based growing stock estimation comparison results are shown in

Table 7. Growing stock regression results in species-specific stands (Table 2) based on a single height variable obtained from four CHM height metrics (Equation (4)). The second-best variable is in brackets. The two best models per species are marked by grey shading. R²: coefficient of determination, RMSE%: RMSE percentage of the growing stock mean across the study area, p50 … p80: canopy height percentiles, mean: mean of canopy height, QMCH: quadratic mean of canopy heights, numbers in parentheses: the total number of stands.

Height Metrics Models	Independent Variable	R2	RMSE (m3 ha−1)	RMSE %
Scots pine—stands (687)
ALS cloud metrics—0.5 m cut-off	p60 (p50)	0.86	52.2	18.6%
ALS cloud metrics—no cut-off	p70 (QMCH)	0.89	46.8	16.7%
ALS-based GRID 0.75 m	p60 (p70)	0.90	44.8	16.0%
GeoEye-1 GRID 0.5 m (677 stands)	p50 (p60)	0.89	46.3	16.6%
incl.: Scots pine—stands > 1 ha (400)
ALS cloud metrics—no cut-off	QMCH (p75)	0.90	46.4	16.6%
ALS-based GRID 0.75 m	p60 (p70)	0.92	42.3	15.1%
GeoEye-1 GRID 0.5 m (400 stands)	p50 (p60)	0.91	45.1	16.1%
incl.: Scots pine—stands < 1 ha (287)
ALS cloud metrics—no cut-off	p70 (p75)	0.87	51.6	18.4%
ALS-based GRID 0.75 m	p60 (p70)	0.87	50.8	18.1%
GeoEye-1 GRID 0.5 m (277 stands)	p50 (p60)	0.86	56.5	20.2%
Norway spruce—stands (291)
ALS cloud metrics—0.5 m cut-off	p50 (p60)	0.92	33.7	16.9%
ALS cloud metrics—no cut-off	p80 (p75)	0.93	31.4	15.7%
ALS-based GRID 0.75 m	p70 (p60)	0.94	29.2	14.6%
GeoEye-1 GRID 0.5 m (288 stands)	p60 (p50)	0.92	32.7	16.4%
Birch—stands (145)
ALS cloud metrics—0.5 m cut-off	Mean (p60)	0.92	33.9	17.8%
ALS cloud metrics—no cut-off	QMCH (p80)	0.94	28.4	14.9%
ALS-based GRID 0.75 m	p70 (p75)	0.94	29.0	15.3%
GeoEye-1 GRID 0.5 m (143 stands)	p50 (p60)	0.92	33.7	17.7%
Black Alder—stands (128)
ALS cloud metrics—0.5 m cut-off	p60 (p70)	0.88	52.3	19.4%
ALS cloud metrics—no cut-off	p80 (p75)	0.90	48.8	18.1%
ALS-based GRID 0.75 m	p70 (p75)	0.92	46.5	17.2%
GeoEye-1 GRID 0.5 m (127 stands)	p50 (p60)	0.87	54.6	19.5%

. In addition, in the pine case, the stands were separated into two parts based on area values (less than 1 ha and more than 1 ha) to study what effect the stand area has on the growing stock estimation efficiency. Finally, the UltraCam-based models were excluded due to the low number of selected stands covered by UltraCam imagery.

In general, the forest stand-based ABA showed slightly worse or similar (birch case) results for forest growing stock estimation compared to the 0.25 ha plots-based ABA. Due to the incompleteness in canopy detection rates using GeoEye-1 stereo imagery [36], the number of stands used in GeoEye-1 regression analysis is slightly decreased compared to ALS-based models.

4. Discussion

4.1. Growing Stock Estimation Based on Individual Tree Detection

Besides errors in stem volume–tree height allometry from harvest data by using the stem taper coefficient, tree detection omission (missing trees) and commission (false, non-existent trees) errors, using the ITD local maxima, significantly affect growing stock estimation [10]. This is primarily due to incomplete detection of the understorey and overlapping trees of the second forest floor. Therefore, it is necessary to consider that canopy maxima applied to gridded CHMs can adequately identify only dominant and co-dominant trees in the upper canopy. To evaluate the impact of ITD omission-related errors, the ~8000 manual stereo-based delineated trees were compared with the number of harvested trees in corresponding polygons after clearcut logging. As manually extracted trees mainly belong to the first forest floor (dominant and co-dominant trees), the harvested trees with heights less than 14 m (1st percentile of manually delineated trees) were excluded from the comparison as belonging to the second forest floor (intermediate and suppressed trees). As a result, the number of manually detected trees was similar to the remaining harvested trees in the ±10% range for all species. Therefore, it suggests that manually extracted trees used for ITD calibration explain ±10% of all trees on the forest’s first floor. In this turn, the number of excluded trees (omission errors) varied from 30% for coniferous (pine, spruce) to 40% for deciduous (birch, black alder), from which the overwhelming majority belong to the second forest floor.

Based on the most recent research on Latvian forests, the overall growing stock of the second forest floor accounts for an average of 13.1% of the total forest growing stock [37]. Thus, only considering that omission and commission errors in tree detection of the forest’s first floor usually do not fully compensate for each other, this study’s ‘canopy maxima’ ITD approach underestimates the growing stock by at least 13% due to the intermediate and suppressed trees. After a comparison of the observed growing stock of inventory stands with ITD-based estimates (Table 5), the results showed that overestimation (except spruce case) of ITD-based prediction varied from 1–2% (deciduous trees) to 9% (pine). Inaccuracies can explain this in the growing stock estimation of the Latvian forest inventory data collected by the visual-based probability-based sampling method. According to Grīnvalds (2014) [21], who compared the growing stock determination from the inventory procedure adopted in Latvia with accurate measurements (callipers and harvester data), the forest inventory underestimates growing stock by an average of 17.6% (pine 18.2%, spruce 9.5%, birch 19.2% and black alder 26%). Thus, the underestimated inventory data partially compensates for underestimating the forest growing stock calculated using the ALS-based ITD in this study. This fact can also partly clarify the growing stock underestimation by ALS-based ITD in the spruce case (Table 5), as the forest inventory most accurately explains the spruce stock. Nevertheless, in this study, conducting a comprehensive analysis and propagating the growing stock uncertainty from stem volume–tree height allometry and ITD to a stand/plot level without accurate reference data (at least a few hectares) was impossible.

4.2. ALS and GeoEye1 Image-Based CHM Performance in Growing Stock Estimation

In general, growing stock estimations using the ABA confirmed the conclusions of previous studies (Table 1), including Balenović et al. (2017) [38] and Goodbody et al. (2019) [4], such that DAP and VHRSI image-based CHMs show a similar performance to ALS-based height metric estimators. However, one unanticipated finding was that for all stand/plot-based models, the best growing stock predictors were based on grid-based (ALS 0.75 m or image-based) canopy height models, except for only one case for birch stands, where ALS-cloud height metrics were higher. Moreover, the image-based UltraCam canopy height metrics were ~2–3% (RMSE of mean) better than ALS metrics for pine and black alder growing stock estimations in the plot-based (0.25; 1 ha) approach. From previous studies in Table 1, only Yu et al. (2015) [20] reached a similar accuracy level for WorldView-2 VHRSI compared with ALS point cloud estimators in growing stock and AGB estimation in Finland boreal forest. In the stands case, for all models, the ALS grid-based metrics explained the growing stock estimation 2% (RMSE % of mean) better than the GE1 image-based height metrics.

The high performance of image-based models could be explained by the closed-canopy, tree height homogeneity of pine and black alder forest species without canopy cover gaps in the study area. In turn, the occlusions and shadows of conic spruce and multi-treetop “loose” crowns of birch reduced the accuracy of image-based models compared to ALS-based models. The previous research by Mielcarek et al. (2020) [14] and Straub et al. (2013) [19] showed that canopy inhomogeneity, gaps and significant height variability in canopy cover strongly influence the efficiency of image-based matching techniques in accurate outer canopy reconstruction. The relatively poor performance of ALS point cloud metrics, especially in comparison with the same ALS data grid-based metrics, could be that models based just on one variable of tree height are more suitable for forest parameter estimation linked to canopy height. As in our case, the growing stock data were calculated/trained based on individual tree detection related solely to canopy maxima/height. In turn, the ALS point-cloud metrics, especially when using more than one predictor variable (e.g., height + basal area, stand age, crown size, etc.), are more robust in models related to vertical vegetation distribution and crown/canopy density [39].

A crucial aspect of image-based CHMs is that they cannot accurately estimate canopy heights compared with ALS point clouds. A recent study that used the same imagery data sets showed that all image-based models underestimated canopy height [36]. For example, the canopy height estimation accuracy of GE1-based models differed from −1.49 m (pine) to −0.31 m median errors (black alder) compared with ALS data.

The occurrence of poor image matching, especially GE1 VHRSI imagery, results in incomplete canopy detection (CHM completeness), as described in Table 7. Regarding the recent research using the same GE1 imagery [36], the CHM completeness rates ranged from 97% for pine to 99% for other tree species in the study area. This specificity should always be taken into account using stereo imagery. In practice, this problem does not occur when collecting ALS data, except for deciduous trees during the leaf-off season.

Factors such as a base-to-height (B/H_VHRSI) ratio, sun-to-sensor geometry, sensor-to-target geometry, vertical vegetation structure and species composition directly impact the completeness and vertical accuracy of image-based CHMs. From a practical point of view, it means that regression equations based on any accurate data sets (e.g., ALS of field data), including tree/canopy height as a variable, cannot be automatically used for image-based models (e.g., growing stock estimation). This is because canopy height underestimation compared with ALS or field data will lead to a biased, incorrect output. This is also true in any case of applying available regression equations obtained from one imagery set to another. For example, two imagery sets obtained by the same sensor but with different acquisition parameters (e.g., B/H_VHRSI ratio) will lead to two different results using the same regression equation. In other words, any image-based regression model (any satellite stereo or large format airborne imagery > 0.2 m GSD) is site-species-specific, sensor-specific and data set acquisition parameters specific. Therefore, this study does not present the regression equations, concentrating only on RS data comparison analysis. From all the above, model calibration/verification is needed before applying imagery data to the existing regression model.

In this study, the GE1 imagery was acquired with an ideal sun-to-sensor (summer noonday) and optimal base-to-height ratio geometry, which impacted the outcome positively. Considering all the factors affecting the image-based model performance, this study confirms that an increase in image resolution (from 0.5 m GE1 to 0.25 m UltraCam) is not a significant advantage for the ABA, which agrees with White et al. (2013) [23].

4.3. Plot Size Choice and Stand/Plot Specificity

The study results confirm that plot size strongly influences the growing stock estimation based on regular plot use in ABA. In the case of pine, the increment of the rectangular plot size from 0.25 ha to 1 ha led to an improvement of 25%, or ~10 m³ ha⁻¹ for the pine model. Thus, by decreasing the perimeter-to-area ratio, the plot size increment reduces the “edge effect” when a treetop detected outside the plot boundary has part of a crown described by height metrics within the plot boundary [40]. However, the choice of plot size will require a balance between minimising edge effects and keeping a sufficiently large scale of growing stock variability related to species-specific and forest heterogeneity in an area of interest. This study is consistent with the conclusion by Zolkos et al. (2013) [41] that a one-hectare plot is optimal for forest inventory parameters mapping, simultaneously capturing forest heterogeneity and reducing “edge effect” uncertainty in mature, closed-canopy and homogeneous forests.

In all plot-based models, the mean canopy height was the essential image-based height predictor for growing stock estimation (Table 6). Furthermore, the same results regarding the mean canopy height as an optimal predictor of plot-level growing stock were achieved by Mielcarek et al. (2020) [14] and Stepper et al. (2015) [42]. In turn, the GE1 canopy height estimators varied from p50 (median height) to p60 percentiles for stand-based growing stock estimation. It means that growing stock estimation based on regular plots correlates best with mean canopy height, which does not change depending on the tree species in this study, vertical vegetation structure and plot size. Thus, a plot-based approach is better than stand-based growing stock estimation, where canopy height predictors vary (are not stable) and depend on stand sizes. Study results show that stand-based growing stock estimation may be better for homogeneous forests with larger species-specific areas (more than 1 ha), as shown in Table 7. In comparison, the plot-based growing stock predictions are better for heterogeneous forests with smaller species-specific areas (less than 1 ha).

In terms of using the ALS point cloud estimators, variables based on all points (no cut-off height threshold) showed the best results for all stand-based models. In contrast, the ALS point cloud models with ground cut-off (0.5 m height cut) showed better results only in the pine plot-based growing stock estimation.

Interestingly, in most cases of the area-based approach, the quadratic mean canopy height (QMCH) (3) was, among plot-based models, the best ALS point cloud predictor variable of the forest growing stock estimation. Unlike the percentile heights, which in many cases showed better results, the QMCH includes all tree heights, putting more weight on the ALS points in the upper forest canopy [43]. Since the share of dominant and co-dominant trees in growing stock is the largest, the QMCH could be a more stable predictor variable than percentile heights, as it is less dependent on the ALS penetration rate.

4.4. Study Limitations and Practical Implementations

Besides the uncertainty in forest growing stock estimation by canopy maxima-based individual tree detection, the framework lacks a preliminary step related to tree species composition determination and classification. By default, the study used the forest inventory data for polygon selection, which is connected to the corresponding forest dominant species. The problems are that FI data are not always accurate and could not be available for some areas, thereby emphasising the need to add a tree species-specific classification routine. The reader is referred to Fassnacht et al. (2016) [44], who reviewed the most frequently applied remote sensing-assisted tree species classification algorithms, comparing their advantages and disadvantages.

Because the study performed the regression analysis only for areas with the dominant four tree species (>75%), further research should investigate the species-specific uncertainty in mixed (multispecies) forest areas without the dominance of one of the tree species. In addition, this study showed relatively large differences between growing stock estimation species-specific models (e.g., RMSE for GE1 stand ABA, from 32 in spruce to 48 m³ ha⁻¹ for black alder). Thus, it might conclude that there are no significant reasons to use or develop generic (non-species-specific) equations for growing stock estimation in dense Latvian forests, especially if the tree species composition is known.

Concerning the choice of remote sensing data for frequent growing stock estimation in dense forestland, large-format aerial imagery could be the optimal solution. However, all benefits of using high-density ALS data compared to imagery are offset by limited temporal resolution due to the high cost compared to image-based data and data acquisition time over large areas. Thus, airborne imagery and image-based DSMs collected during the three-year cycle (e.g., in Latvia) can be alternatively used for regular forest growing stock calculations and updates. Despite the high effectiveness in the growing stock estimation of stereo satellite imagery demonstrated in this research, VHRSI use in Latvia is limited due to the inability to accurately predict weather conditions (cloud cover) during data acquisition. Furthermore, the short time for optimal sun-to-sensor data collection (midday during the summer season) and the lack of an existing stereo satellite image archive (database) also limit the use of VHRSI in Latvia. However, the VHRSI could be a cheap alternative for operative data collection and monitoring forest inventory changes in dense forestland to support the decision-making process on forest management at the national or regional level. Additionally, it should be noted that image-based sensors (especially VHRSI) have an advantage in tree species composition detection over ALS data as they provide additional spectral information (e.g., infra-red, red-edge, etc.). It is also necessary to remember the need for ALS data as a DTM source for CHM creation from image-based DSMs.

The study showed that where reference field data are missing or limited, growing stock estimations obtained from ALS data by the ITD approach can be used as training/reference data to develop corresponding image-based regressions by ABA. However, the study results should be interpreted with caution when transferring them to other sites. Different species-specific vegetation structure/compositions (e.g., forest canopy density, canopy height variance), study methodology and accuracy of reference data can lead to different results (see Table 1).

5. Conclusions

This study compared the effectiveness of canopy height metrics obtained from ALS and stereo imagery in estimating the growing stock of hemiboreal forests. It showed that ALS and image-based height metrics are comparable in growing stock prediction in dense closed-canopy forestlands. The relative RMSEs did not exceed 20% of the reference mean values obtained by the ALS-based canopy maxima individual tree detection approach for all models. The best relative RMSEs achieved were 13.6% image-based (IB) and 15.7% (ALS) for pine 0.25 ha plots; 10.3% (IB) and 12.1% (ALS) for pine 1 ha plots; 16.4% (IB) and 12.2% (ALS) for spruce 0.25 ha plots; 17.9% (IB) and 14.2% (ALS) for birch 0.25 ha plots; and 15.9% (IB) and 18.9% (ALS) for black alder 0.25 ha plots. The high performance of image-based plot models in pine and black alder could be explained by their closed-canopy, tree height homogeneity with minimum canopy cover gaps and occlusions.

Finally, this research suggests that airborne imagery, VHRSI and, accordingly, image-based DSMs collected regularly can be an efficient solution for typical forest growing stock calculations/updates in addition to a traditional visual forest inventory routine. Furthermore, despite the high effectiveness of stereo satellite imagery in forest growing stock estimation demonstrated in this research, the VHRSI use in Latvia is limited due to a large number of cloudy days and the short summer season of optimal sun-to-sensor geometry. However, subject to the optimal parameters during data acquisition, the VHRSI could be the cheapest alternative for operative data collection and monitoring forest inventory changes in vast and dense forestland.