Sample strategies for bias correction of regional LiDAR-assisted forest inventory Estimates on small woodlots

Abstract

Key message

This study presents an easy-to-apply variable probability sample design that is an efficient and cost-effective method to correct for local bias in regional LiDAR-assisted forest inventory estimates. This design is especially useful for small woodlot owners.

Context

Light detection and ranging (LiDAR)-derived forest inventory estimates are generally unbiased at landscape levels but may be biased locally. One solution to correct local bias is to use ground-based double sampling with ratio estimation where the LiDAR estimates form the large sample covariate and the ground plots are used to estimate a correction or calibration ratio.

Aims

Our objectives were to test the performance of different sample strategies, to correct for local bias, and to determine the most efficient and cost-effective sampling design.

Methods

We compared five sample selection methods and four plot types using simulation. Sample sizes and inventory costs required to achieve 5% standard error were calculated to assess sampling efficiency.

Results

The results showed that bias can be corrected successfully using a doubling sampling approach with ratio estimation, and that variable probability selection methods were more efficient than equal probability selection methods. A big basal area factor (BAF) plot was the most cost-effective on-the-ground plot type.

Conclusion

The most efficient and cost-effective sampling design was list sampling with big BAF plots. This combination can be used to calibrate LiDAR-derived forest inventory estimates for a variety of forest attributes.

1 Introduction

Airborne light detection and ranging (LiDAR) scanning (ALS) is a remote sensing technique that can provide three-dimensional spatial information about forest canopies and the underlying terrain. Based on this information, ALS is used to produce estimates of forest attributes, such as tree height, diameter, crown width, and volume (Næsset 2002; Corona and Fattorini 2008; Asner et al. 2011; Hayashi et al. 2015). These estimates are commonly referred to as Enhanced Forest Inventory (EFI) in Canada (Canada 2013).

In New Brunswick, the provincial government is in the process of producing freely available EFI estimates for the entire province (Dick 2019). The advantage of EFI is that it provides high-resolution spatial forest inventory estimates across the landscape (Andersen et al. 2011). However, these estimates are derived from relative sparse ground data sets that may have questionable sample selection properties; as a result, although landscape-level estimates may be unbiased, there is no evidence that these estimates are accurate and unbiased for smaller forest areas within the forest zone covered by EFI (Hayashi et al. 2015). The need to evaluate and potentially calibrate these estimates is required to assure appropriate forest management decisions are made on a local basis.

Localized bias can arise due to differences in past management treatments, microsites, and disturbance regimes, which can all produce substantial differences in stand structure and area-based forest attributes (Clutter et al. 1983; Oliver and Larson 1996). The high spatial resolution provided by EFI may connotate an accuracy that is not realized on the ground. Evaluation and calibration of these estimates is required to make appropriate forest management decisions on a local scale. One solution is to locally calibrate EFI estimates using a sampling with covariates approach (Iles 2003, p. 443; Kershaw et al., 2016, pp. 343–353).

Sampling with covariates (Kershaw et al., 2016, pp. 343–353), double sampling, or “sampling to correct” (cf. Iles 2003, p. 443) is a widely use sampling approach. It is different than small area estimation where models are built for target areas with limited sample observations (Goerndt et al. 2011). With sampling with covariates, EFI estimates for the target area are used as the “large” sample and a “small” subsample of ground plots is obtained to correct EFI estimates. The sampling process is separated into two components. First is the sample selection method. In practice, systematic sampling is most common, since sample points can easily cover the target area (Freese 1962; Payandeh and Ek 1971; Iles 2003, pp. 171–182). However, applying variable probability selection methods, where selection probabilities are dependent on parameters of interest (e.g., EFI volume or biomass estimates), has potential to be a more efficient sampling strategy (Basu 1969). The second component is the sampling unit (plot type). Fixed area plots and horizontal point samples are common in North America (Kershaw et al., 2016 Chap.9). Controlling the number of trees per plot results in trade-offs between increasing variability between plots (i.e., increasing required sample sizes) and decreasing time spent on each plot (Freese 1961; Gambill et al. 1985; Lynch 2017; Yang et al. 2017). Probabilistic subsampling for detailed tree measurements can improve sample efficiencies in many situations (Grosenbaugh 1979; Bell et al. 1983; Iles 2003, pp. 564–565; Marshall et al. 2004; Kershaw et al., 2016; Yang et al. 2017).

This paper explores the efficiency of different sample selection methods and different plot types on developing cost-effective and locally calibrated LiDAR-assisted inventories for small woodlots. The specific objectives include the following: (1) assessment of bias calibration; and (2) assessment of sampling efficiency, required sample sizes, and associated inventory costs of different sample selection method × plot type combinations.

2 Material and methods

2.1 Study area

The study area was the 80-ha Femelschlag Research Site on the Noonan Research Forest (NRF; N 45°59′12″, W 66°25′15″), managed by the University of New Brunswick’s Faculty of Forestry and Environmental Management in New Brunswick, Canada. The NRF has a permanent grid of big BAF sample plots (BBAF) located on a North-South/East-West grid at 100-m spacing. The forest inventory is conducted every 10 years, with the inventory used in this study collected in 2012. The count BAF was 2 M (i.e., each tallied tree represented 2 m² ha⁻¹) and the measure BAF was 27 M. For each count tree, diameter at breast height (DBH, nearest 0.1 cm), status (live/dead), and species were measured and recorded. On the measure-trees, total height (HT, nearest 0.1 m) was measured. There were 1735 count trees and 100 measure-trees recorded on the 83 sample points contained within the Femelschlag Research Site.

The Femelschlag Research Site was established in 2014 with 83 fixed area plots (FAP) established on the same plot centers as the BBAF. The FAPs had a radius of 11.28 m (0.04 ha), and trees with >= 6 cm DBH were sampled. Tree species, status, DBH, and HT were recorded for 6252 sample trees. Table 1 summarizes the plots and tree measurements for both the BBAF and FAP data.

Table 1 Minimum, average, and maximum of number of trees per plot, number of species per plot, number of trees per ha (density), basal area (BA, m² ha⁻¹), diameter at breast height (DBH, cm), and total height (HT, m) for fixed area plots (FAP) and big BAF sample plots (BBAF) data

EFI estimates were provided by the Department of Energy and Resource Development, New Brunswick. EFI estimates were derived from ALS collected in 2015. The cell resolution was 20 m × 20 m, and a total of 2181 cells covered the Femelschlag Research Site. While EFI provides estimates of several forest attributes for each EFI cell, only gross total volume (GTV, m³ ha⁻¹) was used in this study. The time differences between data collections and LiDAR estimates were ignored due to relatively slow tree growth (10–15 rings per cm) and minor natural or human disturbances in 2012–2015.

Operationally, aligning field samples with EFI cells will produce more accurate results. However, because our sampling grid and sample points were established prior to EFI, we did not have perfect spatial alignment with EFI cells. We assumed that our grid was a representative sample of the area around the sample locations. Thus, in the simulation, each EFI cell was assigned one of the 83 sample points based on distance between cell center and sample point location. This will add variability to our results but also will demonstrate the power of sampling with covariates.

2.2 Field sample point selection

In this study, the large sample was EFI cells and five sample selection methods were considered to select the small samples: (1) simple random sampling (SRS); (2) spatial systematic sampling (SYS); (3) systematic ordered-list sampling (SOL); (4) sampling with probability proportional to prediction (3P); and (5) list sampling (LIST).

SRS and SYS are two equal probability selection methods that are frequently used to select both large and/or small samples in double sampling. SYS is typically applied spatially, but can be applied to ordered lists of factors (Iles 2003, p. 176). The idea of SOL is to sort sample units from small to large by their covariate (i.e., GTV), then select sample units at a regular interval across the covariate range. Since SOL ensures equal distribution over the covariate range, it may provide a more efficient double sampling scheme than either SRS or SYS (Iles 2003, pp. 171–174).

3P and LIST are two variable probability selection methods. 3P was introduced by Grosenbaugh (1964, 1965). Following Iles’ (2003, pp. 450–451) procedure, a cruiser predicts (or estimates) an attribute related to the variable of interest for each individual in the population (tree, plot, ALS cell, etc.), and compares this prediction to a uniform random number drawn from the expected range of predictions. If the prediction ≥ the random number, then the individual is measured. Although the sample size in 3P cannot be precisely controlled, Iles (2003) suggested controlling the range of random numbers to obtain the desired sample size. The maximum random number, KZ, is calculated using:

$$ \mathrm{KZ}=\frac{K}{p}=K\cdot \left(\frac{N}{n}\right) $$

(1)

where K = the expected value (mean) of the predictions; p = proportion of population to be sampled; N = population size; and n = desired sample size. In our application, since the EFI values are known and sample selection can be performed in advance, the number of units selected will be known in advance of field work despite being a random variable.

LIST is based on the probability proportional to contribution to the population total. In this procedure, as illustrated by Kershaw et al. (2016, pp. 353–356), the listed covariate (GTV, in this study) of sample units is firstly cumulated across the whole population. A sample unit is selected if a random number drawn between 0 and the total accumulation lies within the range of the cumulated sums. Thus, a sample unit that has a larger covariate will have higher probability of being selected and measured in detail. With LIST, the exact number of units selected is controlled and known in advance of field work.

2.3 Plot compilation

For each of the 83 sample points, FAP and BBAF data were available. To further explore the performance of tree subsampling within plots, four sample tree selection systems (plot types) were considered in this study: (1) fixed area plots with all trees measured (FAP); (2) fixed area plots with 3P subsampling of measured trees (3PN); (3) horizontal point samples with all trees measured (HPS); and (4) big BAF sample plot (BBAF).

The 3PN was simulated using the FAP data. Height predictions, based on a height-diameter model (MacPhee et al. 2018), were compared with uniform random numbers between 0 and KZ (Eq. 1). If the predicted height ≥ the random number, the measured tree height would be used. The desired tree sample size (n) was set at the same level of effort as observed in the big BAF samples (n = 100 trees across the 83 field plots). KZ (Eq. 1) was calculated using average predicted height (K = 11.51 m), and N = 6252 sample trees. Volume per ha (m³ ha⁻¹) for a full FAP plot was calculated by estimating total volume for each tree using the measured DBH, HT, and Honer et al.’s (1983) metric volume equations, and then summing over all trees on the plot and multiplying by the expansion factor (25 for 11.28-m radius FAP). For the 3PN plots, the mean ratio of measured HT to predicted HT (${\overline{R}}_{\mathrm{HT}}$) for all measure trees was calculated and multiplied by each predicted HT, and the adjusted HTs and measured DBHs were used to estimate volume using the same process used on the full FAP.

To expand the BBAF data to represent HPS data, the NRF height-diameter model (MacPhee et al. 2018) was used to predict HT on all count trees. Volume was then estimated using Honer et al.’s (1983) metric volume equations and multiplied by each tree’s associated expansion factor (varies with tree DBH) and summed to obtain volume per ha (m³ ha⁻¹) for each sample point. For the BBAF, volume and volume to basal area ratio (VBAR) were calculated for each measured tree across the entire sample, and the mean VBAR was then multiplied by each plot’s basal area per ha (m² ha⁻¹) to estimate volume per ha (m³ ha⁻¹) for each sample point (Iles 2003, pp. 314–322, 564–565; Marshall et al. 2004). The EFI GTV values and the four different plot-level summaries are archived in the University of New Brunswick’s DataVerse data archive (Hsu and Kershaw 2020).

2.4 Sample simulation

As described above, EFI cells formulate our large sample of covariates and we will select a smaller subsample of these cells to obtain simulated field measurements. Sample sizes between 10 and 50, in steps of 5, were simulated for each of the five sample selection methods and each plot type was used at each selected cell. Each sample selection method × sample size was simulated 100 times. Sampling was conducted without replacement except for LIST which, as typically implemented, was conducted with replacement (Freese 1962; Kershaw et al., 2016, pp. 353–356). Preliminary analyses showed that the GTV for our study area did not have appreciable bias (approximately 2–9%, depending on plot type) relative to the field plots (bias = cell GTV – Field Volume). To test our hypothesis that sampling with covariates and ratio estimation can be used to correct bias in EFI cells and to provide meaningful comparisons between different sample design combinations, we added multiplicative biases into each EFI cell by generating random numbers between 1.05 and 1.15 for 10% bias and 1.15 and 1.25 for 20% bias (the adjustment was in one direction because we were testing effects of bias, not error). All sample simulations were carried out in the R statistical program (R Development Core Team 2019).

2.5 Sample compilation

In this study, ratio estimation (Freese 1962; Iles 2003, pp. 314–322; Kershaw et al., 2016, pp. 351–353) was used to calibrate (or correct) GTV using the small sample field volume estimates (VOL). Ratio estimation can be considered within either a design-based or a model-based framework (e.g., Gregoire 1998; Magnussen 2015). Here, we choose the model-based (or model-dependent) approach. Because GTV is expected to be very similar to VOL, the relationship between GTV and VOL should be linear and pass through 0. For the ith selected EFI cell in the jth replicate sample in the kth sample design combination (sample selection method × plot type × sample size), the VOL to GTV ratio (R_ijk) was calculated using:

$$ {R}_{\mathrm{ijk}}=\frac{{\mathrm{VOL}}_{\mathrm{ijk}}}{{\mathrm{GTV}}_{\mathrm{ijk}}} $$

(2)

where VOL_ijk = field estimated volume using one of the four plot types; and GTV_ijk = the GTV of selected EFI cell. The mean ratio (${\overline{R}}_{\mathrm{jk}}$) was then calculated for the jth replicate sample in the kth sample design combination:

$$ {\overline{R}}_{\mathrm{jk}}=\frac{\sum \limits_{i=1}^m{R}_{\mathrm{ijk}}}{m_{\mathrm{jk}}} $$

(3)

where m_jk = the number of samples. We note that under the design-based approach, different estimators of the population mean would be necessary, depending on the selection probabilities of the individual cells under the design (Thompson 2012 Chap. 7). From a model-based perspective, the sampling probabilities are irrelevant so long as selection is “noninformative” in that it is not influenced by the residual value itself (e.g., Chambers and Clark 2012); a single estimating equation can be used, depending on the characteristics of the error term for individual cells. For further discussion on noninformative samples in forest inventory, see Magnussen (2015). Equation 3, a mean-of-ratios estimator, is the minimum-variance model-unbiased estimator when the residual variance is proportional to GTV², or the residual standard deviation is proportional to GTV (Brewer 1963). Adjusted GTV (${\mathrm{VOL}}_{\mathrm{ijk}}$) was then obtained using:

$$ {\mathrm{VOL}}_{\mathrm{ijk}}={\overline{R}}_{\mathrm{jk}}\cdot {\mathrm{GTV}}_i $$

(4)

The mean adjusted GTV (${\overline{\mathrm{VOL}}}_{\mathrm{jk}}$) was then calculated as:

$$ {\overline{\mathrm{VOL}}}_{\mathrm{jk}}=\frac{\sum \limits_{i=1}^N{\mathrm{VOL}}_{\mathrm{ijk}}}{N} $$

(5)

The overall error of the sample can be considered from two different perspectives. If we view the LiDAR-derived GTV as a random variable associated with its own standard error, then overall sample error can be calculated using Bruce’s formula (Iles 2003, p. 588; Marshall et al. 2004; Kershaw et al., 2016, p. 380) as derived from Goodman (1960):

$$ \mathrm{se}\%\left({\overline{\mathrm{VOL}}}_{\mathrm{jk}}\right)=\sqrt{\mathrm{se}\%{\left({\overline{R}}_{\mathrm{jk}}\right)}^2+\mathrm{se}\%{\left(\overline{\mathrm{GTV}}\right)}^2} $$

(6)

where $\mathrm{se}\%\left({\overline{\mathrm{VOL}}}_{\mathrm{jk}}\right)$ = the percent standard error of the mean adjusted GTV; $\mathrm{se}\%\left({\overline{R}}_{\mathrm{jk}}\right)$ = the percent standard error of the mean ratio; and $\mathrm{se}\%\left(\overline{\mathrm{GTV}}\right)$ = the percent standard error of the mean original GTV. The percent standard error was calculated using:

$$ \mathrm{se}\%\left(\overline{X}\right)=100\cdot \left(\frac{s(X)/\sqrt{m}}{\overline{X}}\right)=\left(\frac{\mathrm{CV}(X)}{\sqrt{m}}\right) $$

(7)

where X = R or GTV; $\overline{X}$= $\overline{R}$ or $\overline{\mathrm{GTV}}$ (averages); s(X) = standard deviation of X, m = sample size; and CV(X) = coefficient of variation of X. The use of Eq. 7 for R conforms to the usual treatment of the plots as independent estimates. However, for GTV, one might argue that ignoring spatial autocorrelation leads the resulting $\mathrm{se}\%\left(\overline{\mathrm{GTV}}\right)$ to be a slight underestimate. On the other hand, a strong case can be made that the LiDAR-derived GTV, for a given property and a given LiDAR collection, is a fixed value and that randomness only enters through variation in the VOL_ijk of the individual sampled cells. In that case, the term $\mathrm{se}\%\left(\overline{\mathrm{GTV}}\right)$ in Eq. 6. equals zero, and spatial autocorrelation becomes a non-issue. Ignoring spatial autocorrelation, $\mathrm{se}\%\left(\overline{\mathrm{GTV}}\right)$ calculated using Eq. 7 is about two orders of magnitude smaller than $\mathrm{se}\%\left({\overline{R}}_{\mathrm{jk}}\right)$in all cases studied, such that $\mathrm{se}\%\left({\overline{\mathrm{VOL}}}_{\mathrm{jk}}\right)$is identical to within 5 or 6 decimal places. Even with a generous allowance for spatial autocorrelation (e.g., an increase of an order of magnitude), the influence of $\mathrm{se}\%\left(\overline{\mathrm{GTV}}\right)$ on $\mathrm{se}\%\left({\overline{\mathrm{VOL}}}_{\mathrm{jk}}\right)$ is negligible and can hardly be represented meaningfully or discerned in graphical or tabular presentation of the results; it is certainly too small to be of management importance, so we consider it no further and use Bruce’s formula (Eq. 6) in all analyses presented here.

2.6 Data analyses

2.6.1 Sample design effects

Effects of different sample design combinations on the distributions of ${\overline{\mathrm{VOL}}}_{\mathrm{jk}}$ and $\mathrm{se}\%\left({\overline{\mathrm{VOL}}}_{\mathrm{jk}}\right)$were examined using horizontal error bar plots and bean plots (Kampstra 2008), respectively. To assess overall sampling efficiency and determine the variance contributions of selection methods and plot types, a nonlinear mixed effects model was fitted:

$$ \mathrm{se}\%\left({\overline{\mathrm{VOL}}}_k\right)=\frac{b_0}{m^{b_1}} $$

(8)

where m = sample size; and b_i = regression coefficients. Random effects were fitted using sample selection method and plot type nested within sample selection method to both coefficients. With the use of Eq. 8, minimum sample size requirements were determined for the 5% error level, and the design effect (Särndal et al. 1992) was calculated using:

$$ {\mathrm{DE}}_{\mathrm{st}}=\frac{m_{5\%,s}}{m_{5\%,t}} $$

(9)

where DE_st = the design effect between sample design s and t; and m_{5%, i} = the sample size requirement to achieve 5% error from sample design s or t.

2.6.2 Influence of cost constraints on optimal sample design

A field survey design usually needs to consider both statistical efficiency and time or cost (Smith 1938). Both Lynch (2017) and Yang et al. (2017) showed that optimal sample design varies in relationship to cost constraints. Building upon the work of Lynch (2017), Yang et al. (2017) expressed total inventory cost as the sum of overhead costs (e.g., the planning, preparation and compilation), total measurement costs (C_Meas), and plot-to-plot travel costs (C_Travel). Overhead costs are fixed. C_Meas is composed of plot establishment costs, tree identification costs, and tree measurement costs:

$$ {C}_{\mathrm{Meas}}={C}_P+k\cdot E\left[c\right]\cdot m+r\cdot M $$

(10)

where C_P = sample point establishment costs; k = per tree costs to determine an “in” tree, identify species, and measure DBH; E[c] = expected number of trees per sample point; m = required sample size; r = cost to measure other tree attributes (HT, in this study); and M = the total number of measure trees. Zeide (1980) estimated travel distance between two sampling points using a square grid: $d=\sqrt{A/m}$ (A = total area of forest in square of units of distance). The equation for C_Travel can then be expressed as:

$$ {C}_{\mathrm{Travel}}={C}_d\cdot m\cdot \sqrt{A/m}={C}_d\cdot \sqrt{A\cdot m} $$

(11)

where C_d = travel cost per meter.

Given Eqs. 10 and 11, costs for a given sample size and sample design were estimated. Finally, the design effect was expressed as a ratio of costs:

$$ {\mathrm{DE}}_{st}=\frac{C_{5\%,s}}{C_{5\%,t}} $$

(12)

The inventory costs used in this study were based on ancillary data collected on the NRF (Yang et al. 2017): C₀ = 1300, C_p = 2.70, k = 0.675, r = 4.41, and C_d = 0.0945 (2015 Canadian dollars). Crew size was assumed constant (two people) for all designs.

3 Results

Differences in mean adjusted volumes ${\overline{\mathrm{VOL}}}_{\mathrm{jk}}$ and associated standard errors across the three bias levels (no bias, 10%, and 20% added bias) were negligible, except for SRS. It is important to note the negligible differences demonstrate the efficiency of sampling with covariates to correct for local bias. However, because of the small differences in results, we only present results based on 10% bias to illustrate the ability of the ratio estimation and the efficiency of sampling with covariates.

3.1 Sample design effects

Figure 1 shows the mean and range of estimated ${\overline{\mathrm{VOL}}}_{\mathrm{jk}}$ by different sample designs at 10% added bias. Field measured volume varied by plot type (Fig. 1): FAP averaged 272.8 (± 6.69) m³ ha⁻¹; 3PN averaged 270.5(± 5.78) m³ ha⁻¹; HPS averaged 247.0 (± 7.97) m³ ha⁻¹; and BBAF averaged 256.2 (± 7.21) m³ ha⁻¹. Although the range of ${\overline{\mathrm{VOL}}}_{\mathrm{jk}}$ for equal probability selections was greater than the 95% confidence interval (CI) of 83 field plot measured volumes and greater than variable probability selections, mean ${\overline{\mathrm{VOL}}}_{\mathrm{jk}}$ for every combination was close to the mean field measured volume (Fig. 1), except for SOL. As sample size increased, 3P and LIST generally produced narrower ranges in estimated ${\overline{\mathrm{VOL}}}_{\mathrm{jk}}$ than equal probability selection methods (Fig. 1). SOL consistently had the most extreme ${\overline{\mathrm{VOL}}}_{\mathrm{jk}}$ when sample size was 10 (Fig. 1a); however, this difference was minimal as sample size increased. 3P generally had the smallest range of estimated ${\overline{\mathrm{VOL}}}_{\mathrm{jk}}$ for most plot types and sample sizes, except LIST had smaller ranges than 3P for HPS and BBAF plots for sample size 20 (Fig. 1b).

Fig. 1

Comparison of field measured volumes with 95% confidence intervals based on the full 83 samples and the mean and range of $\overline{\mathit{\mathsf{VOL}}}$ (with 10% bias added) for the 100 replicates by five sample selection methods and four plot types. (Selection method and plot type acronyms are as defined in Table 2; ave.$\overline{\mathit{\mathsf{VOL}}}$ = mean of the 100$\overline{\mathit{\mathsf{VOL}}}\mathit{\mathsf{s}}$, $\overline{\mathit{\mathsf{GTV}}}$ = mean of uncorrected EFI volume with 10% bias, $\overline{\mathit{\mathsf{Field}}.\mathit{\mathsf{VOL}}}$ = mean of field measured volume

The $\mathrm{se}\%\left({\overline{\mathrm{VOL}}}_{\mathrm{jk}}\right)$ decreased with increasing sample size for all sample design combinations (Fig. 2). While all selection methods produced many sample instances with low $\mathrm{se}\%\left({\overline{\mathrm{VOL}}}_{\mathrm{jk}}\right)$, the equal probability methods often produce extreme values which increased their mean $\mathrm{se}\%\left({\overline{\mathrm{VOL}}}_{\mathrm{jk}}\right)$. 3P and LIST had lower mean $\mathrm{se}\%\left({\overline{\mathrm{VOL}}}_{\mathrm{jk}}\right)$ than all equal probability selection methods. Plots using subsampling methods (3PN and BBAF) had smaller mean $\mathrm{se}\%\left({\overline{\mathrm{VOL}}}_{\mathrm{jk}}\right)$ than plots where all trees were measured (FAP and HPS; Fig. 2).

Fig. 2

The averages and distributions of $\mathsf{se}\%\left(\overline{\mathit{\mathsf{VOL}}}\right)$ with 10% bias by five sample selection methods and four plot types across the range of sample sizes (10 to 50). (Selection method and plot type acronyms are as defined in Table 2)

The fixed effects component of Eq. 8 only explained 13% of the variation in $\mathrm{se}\%\left({\overline{\mathrm{VOL}}}_k\right)$ and the full model including the random effects explained 19% with an overall root mean square error of 3.5%. The fixed parameter estimates were as follows: b₀ = 25.0915 (1.98981) and b₁ = 0.4223 (0.02181). The standard errors for the random effects associated with sample selection method were as follows: s(b₀) = 4.2605 and s(b₁) = 0.0431; and those associated with plot type nested within sample selection method were as follows: s(b₀) = 3.0578 × 10⁻⁸ and s(b₁) = 3.1263 × 10⁻². While the model accounts for a small percentage of the variation due to high variability in some sample designs (especially in SOL), the results are still informative and useful for sample design comparisons.

Table 2 shows the combined fixed + random effects coefficient estimates by sample selection methods and plot types. The b₀ is the theoretical population standard deviation, and the b₁ represents the rate of standard error reduction with increasing sample size. Larger b₁ represent more efficient designs while larger b₀ represent more variable designs. Theoretically, b₁ should equal 0.5. The b₀ of SRS and SYS were similar; SOL had the largest b₀, and 3P had the lowest. The b₀ of LIST was higher than SRS and SYS. The b₁ for SRS and SYS were usually lower than the ones obtained for variable probability selection methods. The b₁ were slightly lower than 0.5 in general, except for LIST with FAP and 3PN. Based on the same sample selection method, the use of subsampling plot types (3PN and BBAF) also had higher b₁ than all trees measured within plots (FAP and HPS).

Table 2 Coefficient estimates (fixed + random effects) of the nonlinear mixed effect model (Eq. 8) for estimating minimum sample size requirements for the scenario with 10% added bias by selection method and plot type. (Sample selection methods: SRS = simple random sampling, SYS = systematic sampling, SOL = systematic ordered-list sampling, 3P = sampling with probability proportional to prediction, LIST = list sampling. Plot types: FAP = fixed area plot, 3PN = fixed area plot with 3P subsampling of measured trees, HPS = horizontal point sample, BBAF = big BAF sampling)

Table 3 shows the required sample sizes for 5% standard error derived from Eq. 8. The required sample sizes did not increase substantially across the different bias levels, except for SRS. For SRS, the required sample sizes increased proportionally to the bias level increment (Table 3). LIST required the lowest sample sizes for all plot types and bias levels. Figure 3 shows the design effect based on minimum required sample sizes. FAP under SRS and SYS were more efficient than HPS and BBAF, except when these plot types were selected using 3P or LIST. 3P and LIST for all plot types were generally more efficient than other selection methods, especially for the two fixed area plot types (FAP and 3PN). LIST with 3PN was more efficient than all other sample selection method × plot type combinations.

Table 3 The required sample sizes under 5% sampling error for each combination of sample selection methods and plot types by three bias levels (No bias; 10% bias, and 20% bias). (Selection method and plot type acronyms are defined in Table 2)

Fig. 3

Sample size–based design effect comparisons by five sample selection methods and four plot types. Sample size calculations were based on a desired standard error of 5%. Upward pointing (positive) bars indicate that the named panel sample design is more effective than the x-axis sample design (DE = x-axis sample size/panel sample size). Downward pointing (negative) bars indicate the x-axis design is more efficient than the named panel design (DE = panel design sample size /x-axis sample size). (Selection method and plot type acronyms are as defined in Table 2)

3.2 Costs

Unlike the higher sample size design effects (i.e., more sample efficiency) associated with FAP (Fig. 3), when costs were considered, FAP no longer had any positive design effects except when used with variable probability selection methods (Fig. 4). As with sample size (Fig. 3), variable probability selection was more efficient than equal probability designs (Fig. 4). Subsampling plot types had larger cost-based design effects than the plot types with all trees measured (Fig. 4). As a result, the total costs associated with 3PN and BBAF were much less than FAP and HPS (Fig. 5). Because FAP had more expected “in” trees per plot than HPS (Table 3), even though 3PN had lower required sample sizes than BBAF (Table 3, Fig. 3), the costs for 3PN were still higher than BBAF (Figs. 4 and 5).

Fig. 4

Cost-based design effect comparisons by five sample selection methods and four plot types. Sample size calculations were based on a desired standard error of 5%. Upward pointing (positive) bars indicate that the named panel sampling design is more cost-effective than the x-axis sample design (DE = x-axis cost/panel cost). Downward pointing (negative) bars indicate the x-axis design is more cost-efficient than the named panel design (DE = panel design cost/x-axis cost). (Selection method and plot type acronyms are as defined in Table 2)

Fig. 5

Total costs by five sample selection methods and four plot types across three bias levels (No bias; 10% bias, and 20% bias). (Selection method and plot type acronyms are as defined in Table 2)

4 Discussion

In forestry, inventory errors in the 10–20% range are common and acceptable, while the estimates should be unbiased (Schreuder et al. 1993; Gregoire and Valentine 2008). Under the approach of sampling with covariates and ratio estimation, auxiliary variables which support a sampling process should not be a source of bias (Iles 2003) as long as the primary variable (small sample) is unbiased. Typically, bias is less than 5%, and biases in our corrected results ranged from 0.1 to 5.0% (except SOL which bias was around 8% in small sample size). As a result, bias can generally be corrected using ratio estimation without substantially increasing sampling effort. Our results support this opinion, except for SRS (Table 3). Because SRS does not consider any information about the population and the added bias increased population variance, SRS required more samples as bias increased (Table 3). SYS and SOL, which also are equal probability selection methods, ensure spatial or distributional coverage across the population range, and, as a result, the added bias did not have a substantial impact on sample sizes (Table 3).

In theory, SOL should be more efficient than SRS and SYS because SOL samples across the population range (Iles 2003, pp. 175–183). However, in our case, SOL produced the most extreme volume estimates (Fig. 1) and required more samples to achieve a given standard error (Table 3, Fig. 2). It is because the ratios of field volume to EFI volume were calculated using the nearest field plot and the selected EFI cells were not spatially aligned. A large-valued field plot, coupled with an EFI cell centered on a small gap or at an edge, can result in a very large ratio, which could bias the overall average ratio and subsequent corrected volume estimate. Given that SOL ensures sampling across the range of EFI values, these extreme estimates were more likely due to our data limitations, especially at the lowest sample sizes (Fig. 1a). At the larger sample sizes (Fig. 1c, d), SOL was on par with SYS and slightly better than SRS.

Most ALS inventory studies align field plots and ALS cells (White et al. 2013). This was not done here because the field data sets were established independently of EFI cells. We used all EFI cells and matched closest field plots, which was purposely done to increase population size (large sample) and variability for our sample simulations. Despite this, the sampling with covariates and ratio estimation approach was still robust, and produced unbiased estimates with relatively small, cost-effective sample sizes. This study demonstrates that it is possible to use existing ground data that does not exactly align directly with EFI cells for the purpose of locally calibrating LiDAR-assisted predictions. The fundamental outcome of this study is the requirement for an underlying, probabilistic sample at the second stage (Iles 2003, pp. 314–322). The small sample enables us to estimate a total, and the large sample enables us to spread that total over the area of interest (Freese 1962; Iles 2003, pp. 314–322).

Variable probability selection, which allows sampling with probability proportional to the variable of interest (Basu 1969; Iles 2003 Chap.11), minimized extreme value estimates (Fig. 1), reduced overall sample size requirements (Table 3, Figs. 2 and 3), and reduced inventory costs (Figs. 4 and 5). While West (2017) found that variable probability selection was not superior to SRS, this study shows both 3P and LIST had better performance than any of the equal probability selection methods for all criteria evaluated, similar to the results in Yang et al.’s (2019) study. It should be noted that West (2017) only implemented 3P and a variant of 3P that did not require a census of covariates. In this study, LIST typically outperformed 3P.

In traditional forest inventory, LIST and 3P methods were not widely implemented because of the census demands (either a priori for LIST or concurrent for 3P). As shown here, estimates derived from ALS provide suitable lists of covariates for LIST and 3P. Besides the EFI estimated volume used in this study, covariates extracted directly from ALS point cloud data (e.g., LiDAR metrics) also can be used (Yang et al. 2019). Similar to Yang et al.’s (2019) work, Parker and Evans (2004, 2009) showed that double sampling could be an effective approach to derive LiDAR estimates for forest and individual tree attributes, but they only used systematic selection. Other authors have explored additional sample selection criteria such as space filling designs (Junttila et al. 2013) and local pivotal methods (LPM; Grafström et al. 2012; Grafström and Ringvall 2013). These designs are similar to SOL in that they attempt to balance sample allocation across multiple spatial and statistical dimensions. Yang et al.’s (2019) results found that LPM performed as well as 3P and LIST when developing forest-level LiDAR-assisted inventories. As spatial scale increases, stratified sampling or cluster sampling may be more effective (Schreuder et al. 1993; Gregoire and Valentine 2008). However, even with these designs, the principles of variable probability selection can be considered to select sampling units.

The inventory cost of a sample procedure is influenced by two factors: (1) required sample size, which is mainly influenced by sample selection methods (Table 3, Figs. 2 and 3); and (2) the number of “in” trees and measure-trees, which are dependent on plot types. Sample size generally is the main factor used to assess sample efficiency. However, because the spatial extent in this study was only 80 ha, which can be considered as a stand-level or small woodlot inventory, the cost for tree measurement is much higher than for traveling between plots (Yang et al. 2017). Thus, plot type was the major factor impacting total costs rather than selection methods. HPS, as a variable probability selection method, had lower numbers of “in” trees than FAP (Table 3). Although in our study area, HPS and BBAF required larger sample sizes than FAP and 3PN, they still had lower costs. Plot types with subsampling saved much effort in measuring trees (Table 3). Ratios used in 3PN and BBAF were effective for reducing variation between individual measurements (Iles 2003, pp. 564–565; Marshall et al. 2004).

5 Conclusion

Our study shows how LiDAR-assisted forest inventory estimates can be corrected for localized conditions, and the potential of using sampling with covariates and ratio estimation. This approach can be applied to not only LiDAR-assisted estimates but to other remote sensing techniques, such as satellite data and even ground-based optical sensors (Hsu 2019). And it can be used to estimate any parameter of interest, including biomass or carbon content (Chen et al. 2020). As ALS point cloud data and EFI estimates are becoming publicly available (especially in Atlantic Canada), the applications presented in this study become feasible for both large industrial forest landholders and smaller private woodlot owners, giving them high-resolution forest information that is locally calibrated to support management decisions.