Integration of field sampling and LiDAR data in forest inventories: comparison of area-based approach and (lognormal) universal kriging

Abstract

Key message

We compared (lognormal) universal kriging with the area-based approach for estimation of forest inventory variables using LiDAR data as auxiliary information and showed that universal kriging could be an accurate alternative when there is spatial autocorrelation.

Context

Forest inventories supported by geospatial technologies are essential to achieve a spatially informed assessment of forest structure. LiDAR technology supplies comprehensive and spatially explicit data enabling the estimation of wide-scale forest variables.

Aims

To compare the area-based approach with universal kriging for estimation of the stem density, basal area, and quadratic mean diameter using LiDAR data as auxiliary information.

Methods

We used data from 202 inventory plots, distributed in four Forest Management Units with differences in structure and management, and a 6-points/m² resolution LiDAR dataset from a Pinus sylvestris L. forest in Spain to test the accuracy of the (lognormal) universal kriging and the area-based approaches.

Results

In those Forest Management Units where the analyzed variables showed spatial autocorrelation, kriging showed better results than the area-based approach in terms of RMSE and Pearson coefficient between observed and estimated values, although lognormal universal kriging provided slightly biased estimations (up to 2%).

Conclusion

Universal kriging is an accurate method for estimation of forest inventory variables with LiDAR data as auxiliary information for those variable exhibiting spatial autocorrelation.

1 Introduction

Forest monitoring is essential to design management strategies and conservation policies (Scott and Gove 2002; Köhl et al. 2006). Forest inventories are the main sources of information to assess forest structure and for monitoring forest dynamics. Geographic information systems (GIS) and remote sensing technologies enable forest resources to be mapped through wall-to-wall estimation of variables from plot measurements and have transformed how forest inventory data are employed (Tomppo et al. 2008).

Light Detection and Ranging (LiDAR) provides accurate terrain elevation models and comprehensive and spatially explicit information of vegetation height (Lindberg et al. 2012), enabling the estimation of forest attributes over large areas (Ayrey and Hayes 2018). LiDAR has become an essential tool for forest managers (White et al. 2013), over recent decades, with LiDAR-derived metrics supported by field inventory data proving to be a useful approach for the evaluation of forest structure (Blomdahl et al. 2019), estimation of fire-related variables (Hall et al. 2005; Marino et al. 2018), biomass assessment (Hall et al. 2005), and measurement of individual tree attributes or species characterization (Lin and Hyyppä 2016; Harikumar et al. 2017). LiDAR height metrics are highly correlated with growing stock variables, such as basal area (BA), volume, and aboveground biomass in forests with relatively simple structure (Lefsky et al. 1999; Means et al. 2000). But relationships between height and BA and biomass are not so straightforward for complex forest structures (Palace et al. 2015).

Statistical inference is required to evaluate the field measurements over large areas, usually incorporating LiDAR information through regression models (White et al. 2013). The area-based approach (ABA) (Næsset 2002) for the prediction of forest variables is a successful technique for the integration of LiDAR and field inventory data in forest inventories (White et al. 2013). Under ABA, linear models can be built and applied to obtain estimates of forest attributes over an entire study area using wall-to-wall LiDAR metrics as auxiliary variables. Other modeling approaches have been applied with ABA, such as the k-nearest neighbor method (Tomppo and Halme 2004) or the linear mixed-effects models (Mauro et al. 2017). The ABA has been widely applied for estimating forest stand attributes at fine scale (Means et al. 2000; Magnussen et al. 2012) or for generalized prediction models from forest inventory data (Næsset 2004; Bouvier et al. 2015).

Spatial autocorrelation is a common characteristic of many forest inventory variables linked to site and climate factors as well as ecological processes driving the stand development (Legendre 1993). In this sense, kriging exploits the spatial autocorrelation of forest inventory variables to increase the accuracy of estimations (Montes et al. 2005). When there is an underlying variable distribution explained by auxiliary variables, universal kriging (UK) incorporates auxiliary information as a mean function (Montes and Ledo 2010), and the residuals between the observed values and the mean function exhibit spatial autocorrelation. Nowadays, many forest inventories combine information from two phases: information from remotely sensed data, such as satellite images, aerial photogrammetry, or LiDAR, is collected in the first phase, whereas the second phase comprises terrestrial plots (Mandallaz and RongHua 1999). Universal kriging employing variables derived from aerial photogrammetry as auxiliary information has proved to be effective in reducing estimation error in forest inventories (Mandallaz 2000). However, kriging estimator optimality is based on normality assumptions (Matheron 1971) that the target variables in a forest inventory often do not meet. Even though logarithmic transformations are often required for these variables, the difficulties involved in back-transforming lognormal kriging interpolations (Dowd 1982) have prevented the use of lognormal kriging in operational forest inventories. However, robust lognormal kriging estimators as proposed by Cressie (2006) or Nussbaum et al. (2014) have demonstrated their viability for interpolation of soil variables.

Our aim in this work is to evaluate and compare the performance of two techniques for estimation of stem density per ha (N), basal area (BA), and quadratic mean diameter (QMD) in dense stands of pines: (i) ABA, where field-based response variables are linked to LiDAR metrics through linear models and (ii) the UK approach for interpolation of field-based variables with LiDAR metrics as auxiliary variables. We also aim to assess the effect of spatial autocorrelation on the UK estimation error and to provide scientific rationale for the choice of ABA or UK in LiDAR-supported forest inventories.

2 Materials and methods

2.1 Study area

The study took place in Valsaín pinewood, a 7622-ha Pinus sylvestris L. forest located in Sierra de Guadarrama, Spain (Fig. 1) at an altitude between 1200 and 2130 m. Valsaín yields good-quality timber and plays an important role as a carbon sink and in the protection of a high-mountain area (Donés and Garrido 2001; Moreno-Fernández et al. 2015).

Fig. 1

Location of Valsaín pinewood in Europe and study area delimitation (Valsaín pinewood). Forest Management Units (FMUs) are defined, and field inventory plots are represented with triangles

The study area has an average N of 604 stems/ha with a minimum of 75 stems/ha and a maximum of 2147 stems/ha. The average BA is 42 m²/ha with a minimum of 14 m²/ha and a maximum of 91 m²/ha. QMD shows an average of 33 cm/ha, a minimum of 14 cm/ha, and a maximum of 63 cm/ha. As described in the forest management plan, the area is divided into three Forest Management Units (FMUs). The group shelterwood system is implemented during a 40-year regeneration period in those forest compartments that have reached the rotation age (120 years), while in the other compartments, a moderate thinning regime is applied. Regeneration areas and polewoods predominate in FMU 1, whereas dense mature stands and polewoods predominate in FMUs 2 and 3, respectively. We considered a fourth FMU comprising the pine stands located at altitude > 1800 m since they show specific high-mountain stand structure and are managed for protection. The distribution and characteristics of each FMU are shown in Fig. 1 and Table 1.

Table 1 Characteristics of Forest Management Units (FMUs) in the study area

2.2 Data

2.2.1 Field inventory

Two hundred two plots of 13-m radius were positioned using a Leica SR530 GPS receiver with sub-metric accuracy. Within the plots, the diameter at breast height (DBH) and tree height (H) of all trees with DBH ≥ 7.5 cm were manually measured with a caliper and a Vertex III dendrometer. When modeling with LiDAR, circular plots are recommended (White et al., 2013) since they are easier to establish (Adams et al. 2011) and the edge effects are reduced (Wulder et al. 2012). We evaluated N, BA, and QMD stand variables. The inventory was carried out at the same time as the LiDAR acquisition survey (2009).

2.2.2 LiDAR data

The LiDAR flight covered the entire study area with an average density of 6 points/m² and a minimum of 1.6 points/m² which has proven to be appropriate to derive a digital terrain model (DTM), stand volume, or canopy metrics (Jakubowski et al. 2013). The flight details are summarized in Table 2.

Table 2 LiDAR flight plan parameters

Matching and adjustment of LiDAR flights were carried out using TerraMatch software (Terrasolid 2021a), and point classification was performed using TerraScan software (Terrasolid 2021b). A DTM and a digital surface model (DSM) were built with a 1-m cell grid following two steps: (i) automatic and manual detection of atypical points and (ii) automatic ground/not-ground point classification. For the second step, the point cloud was filtered out, bare ground points were automatically extracted with TerraScan software and visually checked for verification, and the point cloud was classified into two single classes: ground and vegetation points.

2.2.3 LiDAR metrics

A set of elevation metrics were derived from the LiDAR point cloud using FUSION software (McGaughey 2018). LiDAR point cloud metrics corresponding to field inventory plots were extracted based on the circular areas of 13-m radius.

We initially defined 43 LiDAR metrics related to height, height variability, vegetation cover (Lefsky et al., 2005), and number of returns (Hawbaker et al. 2010; Bater et al. 2011; Woods et al. 2011) (Table 3). A 2-m threshold was applied to remove ground and understory (Bouvier et al. 2015), considering an appropriate approach when studying mature, coniferous forests (White et al. 2013). A post-mission quality check assessment to evaluate the final product was performed for LiDAR point-cloud metrics as recommended in White et al. (2013).

Table 3 LiDAR metrics considered for modeling

2.3 Methods

2.3.1 Area-based approach

For the prediction of forest variables in the ABA, we applied linear regression models with N, BA, and QMD as Z variables, while 43 preselected LiDAR metrics were tested as auxiliary variables (Table 3). Linear relationships were built with no variable transformation or by using (i) the exponential model (Eq. 1), which assumes a linear relationship between the log transformation of the Z variable and the auxiliary variables (Næsset 2002; Frazer et al. 2011; Bouvier et al. 2015):

$$ \ln \left(\hat{Z}\right)={\beta}_0+\sum \limits_{j=1}^p{\beta}_j{X}_j $$

(1)

and (ii) the potential model, which assumes a linear relationship between the log transformation of both Z and auxiliary variables (Eq. 2).

$$ \ln \left(\hat{Z}\right)={\beta}_0+\sum \limits_{j=1}^p{\beta}_j\ln \left({X}_j\right) $$

(2)

In each FMU, we carried out a preselection of those LiDAR metrics that present a priori relationship with each response variable based on the literature and the expert knowledge of stand structure and dynamics. We then analyzed the Pearson correlation coefficient between LiDAR metrics and the Z variables to identify those metrics with the highest correlation with each response variable. This preselection step provided a set of ten LiDAR variables for each variable and FMU. To obtain the definitive set of predictors for each response variable (N, BA, and QMD), different combinations of no more than four variables were tested by comparing the Akaike information criterion (AIC) calculated by the StepAIC function implemented in MASS R package (Venables and Ripley 2002), discarding redundant variables (Ayrey and Hayes 2018). Seven outliers were detected and filtered out in six models.

In order to determine the final model to be used, we studied the performance of the three models, selected for each response variable and FMU: without variable transformation, with log transformation of Z variable, and with log transformation of Z and LiDAR variables, by comparing the coefficient of determination (R²), Pearson coefficient (ρ_Z,X) and bias. These coefficients were calculated as follows:

$$ {R}^2=1-\frac{\sum \limits_{i=1}^n{\left({Z}_i-{\hat{Z}}_i\right)}^2}{\sum \limits_{i=1}^n{\left({Z}_i-\overline{Z_i}\right)}^2} $$

(3)

$$ {\rho}_{\mathrm{Z},\mathrm{X}}=\frac{\sigma_{\mathrm{Z}\mathrm{X}}}{\sigma_{\mathrm{Z}}{\sigma}_{\mathrm{X}}} $$

(4)

$$ \mathrm{bias}\left(\%\right)=\frac{\frac{1}{n}\times \sum \limits_{i=1}^n\left({\hat{Z}}_i-{Z}_i\right)}{\overline{Z_i}}\cdotp 100 $$

(5)

where n is the sampling size; Z_i is the observed value of variable Z; ${\hat{Z}}_{\mathrm{i}}$is the predicted value of variable Z; σ_Z,X is the covariance of Z and auxiliary variable X; σ_Z is the standard deviation of Z variable; and σ_X is the standard deviation of X.

Back-transformation was necessary previous to validation. As a systematic bias is introduced with log-transformation, the application of a model correction factor (CF) was necessary (Bouvier et al. 2015):

$$ \mathrm{CF}=\exp \left(\frac{{\mathrm{SEE}}^2}{2}\right) $$

(6)

where SEE is the standard error of regression estimates obtained with lm function in stats R package (Wilkinson and Rogers 1973; Chambers 1992; R Core Team 2019). Therefore, the back-transformation would be (Eq. 9) for the exponential model, and (Eq. 10) for the potential model.

$$ \hat{Z}=\left({e}^{\beta_0}\prod \limits_{j=1}^p{e}^{\beta_j\cdotp {X}_j}\right)\cdotp CF $$

(7)

$$ \hat{Z}=\left({e}^{\beta_0}\prod \limits_{j=1}^p{X_j}^{\beta_j}\right)\cdotp CF $$

(8)

A leave-one-out cross-validation was carried out to assess model accuracy, estimating the Pearson coefficient, the bias, and the RMSE (%) through Eqs. 4, 5, and 9.

$$ \mathrm{RMSE}\left(\%\right)=\frac{\sqrt{\frac{1}{n}\sum \limits_{i=1}^n{\left({\hat{Z}}_i-{Z}_i\right)}^2}}{\overline{Z_i}}\cdotp 100 $$

(9)

2.3.2 Universal kriging approach

In the UK model, the variable value Z(s₀) at the location s₀ is expressed as the sum of a mean function of the auxiliary variables X_j (s₀) and a spatially autocorrelated residual δ(s₀):

$$ Z\left({s}_0\right)=\sum \limits_{\mathrm{j}=0}^p{\beta}_j{X}_j\left({s}_0\right)+\delta \left({s}_0\right) $$

(10)

We used as Z variable the field-based N, BA, or QMD (or the log transformation of these variables) and as auxiliary variables, X_j (s₀), the intercept and the LiDAR-derived variables. As the mean function establishes linear relationships between the Z variable and the auxiliary variables, analogous to those of ABA models, the same variable selection and transformations applied for ABA were applied as well for UK (Table 4).

Table 4 Final models for stand structure variables (Z variables) applying area-based approach (ABA) for stem density (N), basal area (BA), and quadratic mean diameter (QMD)

The UK prediction of variable Z at location s₀ is a linear combination of values Z in sampling locations s_i corresponding to the field plots:

$$ P\left(Z;{\boldsymbol{s}}_0\right)=\sum \limits_{\mathrm{i}=1}^n{\lambda}_{\mathrm{i}}Z\left({\boldsymbol{s}}_{\mathrm{i}}\right) $$

(11)

The coefficients λ_i that weight each observed value Z(s_i) are calculated minimizing the prediction error under the unbiasedness condition (Cressie 1993). The kriging variance can be calculated as (Cressie 1993):

$$ {\sigma}_{\mathrm{k}}\left({\boldsymbol{s}}_0\right)=\sum \limits_{\mathrm{i}=1}^n{\lambda}_{\mathrm{i}}\gamma \left({\boldsymbol{s}}_0-{\boldsymbol{s}}_{\mathrm{i}}\right)+\sum \limits_{\mathrm{j}=0}^p{m}_{\mathrm{j}}{X}_{\mathrm{j}}\left({\boldsymbol{s}}_0\right) $$

(12)

To calculate λ_i, the variogram (d) is needed (where d is the distance between locations). Under UK (d) can be estimated from the residual semivariance:

$$ \gamma (d)=\frac{1}{2n(d)}\sum \limits_i^{n(d)}{\left(\delta \left({s}_i\right)-\delta \left({s}_i+d\right)\right)}^2 $$

(13)

Restricted maximum likelihood (REML) was applied for the simultaneous estimation of the coefficients of the mean function and the parameters of the variogram (Harville 1974; Mardia and Marshall 1984). Visual inspection of the residual variogram was used for selection of the variogram model. Finally, the spherical variogram model was fitted for all the variables analyzed.

In those models where log-transformations were applied, lognormal kriging (Dowd 1982) was used. The unbiased estimator of Z(s_i) when the UK prediction is carried out using the log-transformed variable is given in Eq. 14 (Cressie 2006):

$$ \hat{Z}\left({\boldsymbol{s}}_0\right)=\exp \left(\mathrm{P}\left(\ln (Z);{\boldsymbol{s}}_0\right)+{\sigma}_{\mathrm{k}}\left({\boldsymbol{s}}_0\right)/2-\sum \limits_{\mathrm{j}=0}^p{m}_{\mathrm{j}}{X}_{\mathrm{j}}\left({\boldsymbol{s}}_0\right)\right) $$

(14)

where $\hat{Z}\left({\boldsymbol{s}}_0\right)$ is the unbiased back transformation of P(ln(Z); s₀), which is the UK leave-one-out prediction at validation plot s₀, and m_j are the Lagrange multipliers to ensure the condition $\sum \limits_{\mathrm{i}=1}^n{\lambda}_{\mathrm{i}}{X}_{\mathrm{j}}\left({\boldsymbol{s}}_{\mathrm{i}}\right)={X}_{\mathrm{j}}\left({\boldsymbol{s}}_0\right)$.

Leave-one-out cross-validation of the UK model was carried out in order to compare the results of the ABA and the UK approaches. The geostatistical analysis was implemented using the software package Geostat developed in Matlab® by the authors.

2.3.3 Comparison of ABA and UK

The bias (Eq. 5), Pearson coefficient between estimated and observed values (Eq. 4), and RMSE (Eq. 9) for the leave-one-out cross-validation residuals of the ABA and UK models were used to compare the results from both approaches. Box plots were used to compare the prediction errors of both ABA and UK models for each FMU and variable, estimated through SEE and σ_k, respectively.

3 Results

3.1 Area-based approach

Exponential models were fitted for N in FMUs 1, 2, and 3, and BA in FMUs 1 and 3. QMD models were linear in all the FMUs, whereas in FMU 4, the models were linear for all the variables analyzed (Table 4). In every N model, regardless of the FMU, the LiDAR variables included were related with cover fraction (PFR_A2 or PAR_A2) and mid to high height elevation percentiles (missing in the FMU 2 where we found height mode instead). In BA models, the most frequent variables included were h₂₀ and PFR_A2 (FMUs 1 and 3). In the case of QMD, the variables included were h₄₀ in forest units I and II and h₇₀ in FMUs 3 and 4 (Table 4).

3.2 Universal kriging approach

REML estimations of variogram parameters for N, BA, and QMD for each FMU are shown in Table 5. The ratio from the residual variogram partial sill to the empirical variogram sill represents the variance linked to spatial autocorrelation (Montes and Ledo 2010), which is maximum in FMU 4 for all the analyzed variables and minimum in FMU 3 (Fig. 2 and Table 5). The spatial correlation ranged from 141 m in FMU 1 to 2101 m in FMU 4 for N, from 978 m in FMU 2 to 11556 m in FMU 4 for BA and from 1299 in FMU 1 to 4576 m in FMU 2 for QMD. The variance explained by the LiDAR variables (estimated as a ratio from the linear combination of variograms and cross-variograms of the auxiliary variables weighted by their β coefficients to the empirical variogram sill) is lower than the variance linked to spatial autocorrelation, for QMD in FMUs 1 and 3, and BA in FMU 3, whereas the LiDAR metrics stand for more than 50% of the variance of N in FMUs 2 and 3, QMD in FMU 2, and BA in FMU 4 (Fig. 2).

Table 5 Spherical variogram model parameters and trend function coefficients with their respective p-values for stem density (N), basal area (BA), and quadratic mean diameter (QMD) in universal kriging models

Fig. 2

Universal kriging variogram model γ (discontinuous line) fitted to the residual variogram (hollow circles) and linear combination of the fitted variogram model and $\sum \limits_{\mathrm{i}}^{n(d)}\sum \limits_{\mathrm{m}=1}^p\sum \limits_{\mathrm{j}=1}^p{\beta}_{\mathrm{j}}{\beta}_{\mathrm{m}}{\left({X}_{\mathrm{j}}\left({\boldsymbol{s}}_{\mathrm{i}}\right)-{X}_{\mathrm{m}}\left({\boldsymbol{s}}_{\mathrm{i}}+h\right)\right)}^2$ (continuous line) vs empirical semivariogram (filled circles) of the stem density (N), basal area (BA), and quadratic mean diameter (QMD) in a Forest Management Unit 1 (FMU 1), b FMU 2, c FMU 3, and d FMU 4

3.3 Comparison of ABA and UK approach

FMU 4 shows the lowest correlation between estimates and field observations in both approaches (Fig. 3, Table 6) whereas FMU 2 showed the highest correlation, reaching R² values up to 73% in QMD. UK obtains better R² and RMSE results (and similar bias) than ABA for those variables exhibiting strong spatial autocorrelation in those FMUs where linear models were fitted (QMD in FMU 2 and all variables—N, BA, and QMD—in FMU 4). ABA and UK showed very similar results for the remaining combination of FMU and variables where linear relationships were fitted. Lognormal kriging for BA in FMU 1, where this variable also exhibited spatial autocorrelation, improved R² and RMSE over ABA (from R² = 0.66 and RMSE = 0.25 with ABA to 0.72 and 0.23 respectively with lognormal universal kriging) at the expense of an increase of bias (from 0.31 with ABA to 1.89 with lognormal universal kriging). Lognormal kriging point prediction bias was less than 2% for all the variables.

Fig. 3

Scatterplot of predicted vs observed values in every Forest Management Unit (FMU) for stem density (N), basal area (BA), and quadratic mean diameter (QMD) under area-based approach (ABA) (a, c, e) and universal kriging (UK) approach (b, d, f). The continuous line represents the 1:1 line

Table 6 Comparison of area-based approach and universal kriging approach for stem density (N), basal area (BA), and quadratic mean diameter (QMD), in terms of Pearson coefficient between estimated and observed values, RMSE, R², and bias

Regarding the prediction errors (Fig. 4), we observe that in every studied model, UK presents a lower error average than ABA, as well as a wider interquartile range of prediction variances due to the effect that the varying distance between plots produces in the final UK prediction.

Fig. 4

Boxplot of the prediction errors for stem density (N), basal area (BA), and quadratic mean diameter (QMD) under area-based approach (ABA) and universal kriging (UK) approach, estimated through the standard error of regression estimates (SEE) and kriging error (σ_k), respectively. a N exponential models. b N linear models. c BA exponential models. d BA linear models. e QMD linear models

4 Discussion

Both the ABA and UK approaches have proved to be useful techniques for integrating field sampling and LiDAR data in forest inventories showing reasonable bias, RMSE, and Pearson coefficient between observed and estimated values for modeling N, BA, and QMD. Our results show that the UK approach may increase the accuracy for those variables exhibiting spatial autocorrelation beyond the distance between sampling plots. The spatial correlation of the analyzed variables depends on the spatial variations of site conditions, disturbance regime, stand regeneration, and development processes and management. FMU 4 presents the greater amount of variance explained by autocorrelation for the analyzed variables, possibly due to the lack of silvicultural treatments in this FMU covering high altitude protection areas, where variability in forest structure is due to changes in the ecological conditions and the natural dynamics, in contrast to the man-induced forest structure in the other managed FMUs, where group shelterwood is applied at forest compartment level, leading to a shorter spatial autocorrelation range.

The Pearson coefficients from the cross-validation for N, BA, and QMD attained through ABA are in general within the range reported in other studies (from 0.38 to 0.67 for N (Silva et al., 2018; Hall et al. 2005), from 0.78 to 0.93 for BA (Hall et al. 2005; Reutebuch et al., 2005; Silva et al., 2017), or from 0.39 to 0.78 for QMD (Næsset 2002). However, in FMU 4, which covers the high elevation areas with a more complex forest structure (due to the natural disturbance regime), ABA showed poorer correlations. In this FMU, UK performs better than ABA in terms of Pearson correlation and RMSE (Table 6). LiDAR cover and height percentiles are often the metrics that show the greatest correlation with those variables related to stand growing stock (Means et al. 2000). In our case, h₂₀ had the greatest explanatory power in the BA models and medium-to-high percentiles (h₅₀, h₇₀, and h₉₅) in N models. Height percentiles are usually correlated, and the inclusion of one or another may depend on the stand structure and the LiDAR characteristics. Furthermore, metrics associated with crown cover, such as PFR_A2, PAR_A2, PFR_Ahmean, and PAR_Ahmean, were also relevant for N and BA. The inclusion of spatial pattern measures of pulse returns by height strata may improve the estimations of N and QMD.

In most studies, N estimates show the weakest correlation with LiDAR metrics (Næsset 2002; Goerndt et al. 2010; Hayashi et al. 2014). This may explain the lower Pearson coefficient and RMSE both for ABA and UK (Table 6). The lack of spatial autocorrelation may explain the poor performance of UK for in FMUs 1 and 3, where managed polewoods and regenerating stands predominate (Table 6). On the contrary, in FMU 4, a high mountain forest under natural disturbance regime, the variograms show a noticeable spatial autocorrelation of the analyzed variables beyond 2000 m, and the UK performs much better than the ABA. The bias problems with lognormal kriging have been widely treated in the literature (Journel 1980; Dowd 1982; Cressie and Pavlicová 2005; Cressie 2006; Yamamoto 2007). The RMSE based on cross-validation residuals reported includes the effect of bias, showing that this effect is compensated by an estimation variance reduction when there is spatial autocorrelation (BA in FMU 2). It should be noticed that bias is less than 2% in all the lognormal UK models indicating the suitability of Eq. 14 to bound bias for lognormal UK predictions in our dataset. Other approaches as proposed in Tolosana-Delgado and Pawlowsky-Glahn (2003) may overcome the problems linked to error sensitivity in the variance estimation of the adjustment used to preserve unbiasedness when the kriging prediction is transformed back to the original scale (Clark and Harper 2000).

In the UK models (as seen in Fig. 2), the variance increases with the distance between plots; the plots nearest to the prediction points have significant weight in the prediction (Oliver and Webster 2014). The effect of spatial autocorrelation of both field sampling data and LiDAR metrics over the UK estimators depends on the plot size, the cell size, and the distance between plots (Mauro et al. 2017). Therefore, increasing the number of sampling plots may improve the estimation of the coefficients of both ABA models and UK mean function, but for UK, there should be an additional uncertainty reduction because the kriging variance decrease as the number of nearby sampling plots within the range of spatial correlation increases (Mauro et al. 2017). Increasing the number of plots in those FMUs where the target variables show strong spatial autocorrelation, with a range proximal to the minimum distance between plots, would produce a greater reduction of the UK error. This increase in the number of plots shall be accompanied by a plot size reduction to avoid an inventory cost rise. Another option to reduce the distance between plots would be a two-stage sampling design (Mandallaz and Ronghua 1999) combining a denser network of terrain data points assessed through an inexpensive methodology, such as 2D Terrestrial Laser scanner (Ringdahl et al. 2013) or ForeStereo (Montes et al. 2019) and a second stage of sparser accurately measured sampling plots. It must be taken into account that point estimates are analyzed in our study, whereas operational forest inventories usually estimate variables at forest compartment or forest unit level. Block predictions (referred to UK predictions at forest compartment or forest unit level) should reduce substantially the estimation error concerning the point predictions, especially in small area units, due to the spatial correlation with nearby sampling plots (Lappi 2001). Similarly, using linear mixed-effects models with ABA may reduce the uncertainty of estimations for small-scale (5 to 50 ha) management unit (Mauro et al. 2017).

5 Conclusions

Universal kriging (UK) integrating field sampling and LiDAR data increased modeling accuracy with respect to the area-based approach (ABA) in terms of RMSE and Pearson coefficient for stem density (N), basal area (BA), and quadratic mean diameter (QMD) in Forest Management Units (FMUs) where these variables exhibited spatial autocorrelation (e.g., protection areas subjected to natural disturbance regime). In those cases, the UK approach reduced the estimation error up to 16%.

Lognormal universal kriging was superior to ABA in terms of RMSE and Pearson correlation between cross-validation estimations and observed values when modeling BA in FMU 1, where this variable shows spatial auto-correlation, but bias was, in general, higher with lognormal UK than with ABA.

In those FMU where regeneration areas and managed polewoods predominate, we generally found low autocorrelation and UK did not improve the ABA results.

Under the UK approach, an increment of the sampling intensity along with plot size reduction would increase cost-efficiency in those FMUs where the spatial autocorrelation range gets closer to the plot distance.