Three-Stage Up-Scaling and Uncertainty Estimation in Forest Aboveground Biomass Based on Multi-Source Remote Sensing Data Considering Spatial Correlation

Abstract

Airborne LiDAR (ALS) data have been extensively utilized for aboveground biomass (AGB) estimation; however, the high acquisition costs make it challenging to attain wall-to-wall estimation across large regions. Some studies have leveraged ALS data as intermediate variables to amplify sample sizes, thereby reducing costs and enhancing sample representativeness and model accuracy, but the cost issue remains in larger-scale estimations. Satellite LiDAR data, offering a broader dataset that can be acquired quickly with lower costs, can serve as an alternative intermediate variable for sample expansion. In this study, we employed a three-stage up-scaling approach to estimate forest AGB and introduced a method for quantifying estimation uncertainty. Based on the established three-stage general-hierarchical-model-based estimation inference (3sGHMB), an RK-3sGHMB inference method is proposed to make use of the regression-kriging (RK) method, and then it is compared with conventional model-based inference (CMB), general hierarchical model-based inference (GHMB), and improved general hierarchical model-based inference (RK-GHMB) to estimate forest AGB and uncertainty at both the pixel and forest farm levels. This study was carried out by integrating plot data, sampled ALS data, wall-to-wall Sentinel-2A data, and airborne P-SAR data. The results show that the accuracy of CMB ($R_{adj}^2$ = 0.37, RMSE = 33.95 t/ha, EA = 63.28%) is lower than that of GHMB ($R_{adj}^2$ = 0.38, RMSE = 33.72 t/ha, EA = 63.53%), while it is higher than that of 3sGHMB ($R_{adj}^2$ = 0.27, RMSE = 36.58 t/ha, EA = 60.43%). Notably, RK-GHMB ($R_{adj}^2$ = 0.60, RMSE= 27.07 t/ha, EA = 70.72%) and RK-3sGHMB ($R_{adj}^2$ = 0.55, RMSE = 28.55 t/ha, EA = 69.13%) demonstrate significant accuracy enhancements compared to GHMB and 3sGHMB. For population AGB estimation, the precision of the proposed RK-3sGHMB (p = 94.44%) is the highest, providing that there are sufficient sample sizes in the third stage, followed by RK-GHMB (p = 93.32%) with sufficient sample sizes in the second stage, GHMB (p = 90.88%), 3sGHMB (p = 88.91%), and CMB (p = 87.96%). Further analysis reveals that the three-stage model, considering spatial correlation at the third stage, can improve estimation accuracy, but the prerequisite is that the sample size in the third stage must be sufficient. For large-scale estimation, the RK-3sGHMB model proposed herein offers certain advantages.

1. Introduction

Forests are a pivotal component of the Earth’s ecosystems, playing an essential role in ecological conservation, water source protection, air purification, and carbon cycling [1,2,3]. Aboveground biomass (AGB) is a critical metric in assessing forest health, and its rapid and precise estimation is imperative in the refined management of forest resources and research on carbon sequestration [4]. Traditionally, the acquisition of forest AGB has relied on design-based inference (DB), which uses plot-based data to estimate AGB at a regional level. However, this approach necessitates that the plot data adhere strictly to the principles of random or systematic sampling to ensure unbiased estimation. Therefore, this method faces significant challenges in terms of high labor and time costs when conducting large-scale surveys [5]. The advancement of remote sensing technology, particularly the utilization of multi-source remote sensing data, has opened up new possibilities in the large-scale inversion of forest AGB. Model-based inference has emerged as a widely adopted method of quantitative estimation in remote sensing, not only providing comprehensive and spatially continuous results but also significantly improving the population estimation precision under the same sample conditions [5,6]. Nonetheless, these methods are still contingent on sample size, and ensuring model accuracy requires a substantial number of samples, which perpetuates the issue of high labor and time costs. Multi-stage approaches have emerged as a research focus within multi-source remote sensing estimation technology, utilizing features strongly correlated with target variables to model and obtain a large number of relative true values, thereby expanding the original sample set, enhancing sample representativeness, and reducing the costs associated with obtaining measured samples [7,8,9].

In multi-stage modeling, intermediate variables play a crucial role in ensuring the precision of results. Light detection and ranging (LiDAR) data can accurately depict forest vertical structure information and have a strong correlation with forest resource parameters, such as average height, diameter at breast height (DBH), volume, and AGB, thus offering a distinct advantage in the estimation of forest resource parameters. Airborne LiDAR data (ALS), in particular, have been extensively applied in the quantitative inversion of forest resource parameters [10,11,12]. Due to the high cost of data acquisition, it is often challenging to achieve the quantitative estimation of forest resources across large areas using ALS data alone [13]. Therefore, ALS data are commonly employed as an intermediate variable in multi-stage modeling to expand the sample size [14,15,16,17,18]. Satellite LiDAR data, such as those from Global Ecosystem Dynamics Investigation (GEDI); Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2); and Terrestrial Ecosystem Carbon Inventory Satellite (TECIS), are widely applied in multi-stage modeling, and they are also easier to obtain and involve relatively lower costs for end-users compared with ALS data [19,20,21]. Given the challenges involved in obtaining measured data for footprint locations due to positioning biases and inaccessibility, ALS data can be used as an intermediate variable to construct models that predict the relative true values of forest AGB at footprint locations [19,22,23,24,25,26]. This approach mitigates the impact of positioning errors and the difficulties associated with tracking footprints, meaning that the combination of sampled airborne LiDAR data and spatially discrete satellite waveform LiDAR data represents a common technical approach to large-area forest parameter estimation.

Currently, the data framework for one-stage inference, such as conventional model-based inference (CMB) [17], typically consists of a limited number of sample data and wall-to-wall data. The data framework for multi-stage inference includes the following three configurations: (1) a limited number of plot data combined with sampling ALS data and wall-to-wall data; (2) a limited number of plot data combined with satellite waveform LiDAR data and wall-to-wall data; (3) utilizing a limited number of plot data, sampling ALS data, satellite waveform LiDAR data, and wall-to-wall data. In the multi-stage approach, except for the first stage, the dependent variables in subsequent stages are derived from the estimated values of the previous stage, which serve as relative true values. Consequently, error propagation can significantly influence the final result. Therefore, quantifying its uncertainty is crucial in practical applications. Scholars have conducted studies on this issue, employing models such as two-stage generalized hierarchical model-based inference (GHMB) [9] and three-stage generalized hierarchical model-based inference (3sGHMB) [8], but neither method accounts for the impact of spatial correlation. Geostatistical methods represent a widely recognized approach to mitigating the effects of spatial correlation and have been extensively employed in the estimation of forest resource parameters [27]. Several studies have integrated interpolation techniques into regression or machine learning algorithms, utilizing the spatial interpolation of residuals and then re-adding them to the original results to minimize the influence of spatial correlation on model outcomes. This integration has led to significant improvements in the accuracy of parameter estimations [28,29,30,31]. However, these methods are less commonly applied to multi-stage modeling, especially for uncertainty estimation. Zhao et al. [32] improved the GHMB model by incorporating spatial autocorrelation and introduced the RK-GHMB estimation method, which applies regression-kriging (RK) at the second stage of the modeling process. Their results demonstrate that integrating geostatistical methods can enhance the estimation accuracy of the model and make the uncertainty estimation more scientific and reasonable. To date, no studies have reported on the effectiveness of integrating geostatistical methods into three-stage inference. Additionally, there is a dearth of research comparing the advantages and disadvantages of three-stage inference relative to one- and two-stage inference [33].

Based on the aforementioned analysis, in this study, we introduce an RK model into the three-stage inference method 3sGHMB for forest AGB estimation, denoted as RK-3sGHMB, addressing the challenges arising from the fact that the current 3sGHMB does not account for spatial autocorrelation in the third stage. We then evaluate this method and compare it with CMB, GHMB, RK-GHMB, and 3sGHMB, analyzing their respective strengths and weaknesses, and assessing the effectiveness of the developed model. We also discuss the applicable scenarios of the developed model, with the goal of providing an efficient and high-precision multi-source remote sensing collaborative estimation technology for forest parameters, to support national forest resource annual monitoring services.

2. Materials and Methods

2.1. Study Area

The study area comprised the whole coverage area of two forest farms named Upper Yangge Forest Farm and Tidal Slag Forest Farm (Figure 1). The study area is located in Genhe City, Hulunbuir City, Inner Mongolia, where the terrain is undulating, with a relative elevation difference of between 100 m and 300 m. The average altitude is above 1000 m, characterizing it as a high-latitude, cold region. It has a cold temperate humid forest climate, with characteristics of a continental monsoon climate, featuring moist and cold air, long winters, short summers, and a rainy season from July to August each year. The vegetation is dominated by forest and grassland, with a forest coverage rate as high as 75%. Typical tree species include Larix gmelinii, Betula platyphylla, Pinus sylvestris var. mongolica Litv, and Populus davidiana [34].

2.2. Data

2.2.1. Plot Data

The forest growth in the study area is relatively slow. To increase the sample size, data from two years, 2021 and 2022, were used for subsequent analysis. A total of 103 plots (25 m × 25 m) were obtained over the two years (Figure 2). Trees with a DBH greater than 5 cm within the plots were measured, and parameters such as DBH, tree height, east–west crown width, north–south crown width, and height to the base of the live crown were measured. The positions of the four corners and the center of each plot were precisely measured using a FindCM GPS-RTK produced by China Qianxun Sl company with a measurement accuracy better than 0.15 m. Subsequently, the tree species-specific allometric growth equations proposed by Zhou et al. [35] were used to calculate the AGB of each tree within the plots. The total biomass and AGB density for each plot were calculated from the AGB of individual tree within the respective plots. The unit of total biomass is ton (t), and the unit of biomass density is tons per hectare (t/ha).

2.2.2. ALS Strip Data Scenario Simulation

Wall-to-wall ALS data were collected in August 2022 (Figure 2) using the CAF-LiCHy data acquisition system, with the LiDAR sensor model being Riegl LMS-Q680i produced by RIEGL that in Horn, Austria [36]. The LiDAR point cloud data underwent preprocessing, including attitude correction, noise point removal, coordinate transformation, flight line stitching, and system error correction. Subsequently, LASTOOLS software [37] was utilized for point cloud classification, denoising, and normalization to obtain normalized vegetation point cloud data. Fusion software [38] was employed to extract features such as the mean height, canopy density, height variance and standard deviation, coefficient of variation, interquartile distance, skewness, maximum height, and height percentiles for subsequent analysis [39,40]. Detailed information is presented in Table 1. The study area was divided into 1,539,092 grids (25 m × 25 m) denoted as N, and the mean of each LiDAR data feature within each grid was calculated for further analysis. In order to simulate real-world application scenarios, 13 strips were extracted to serve as a link to GEDI data (Figure 2b), with both the strip width and spacing set at 1.5 km. There were 80 plots located within the strips and 23 plots outside the strips. Plot data within the strips were used for model training, while plot data outside the strips were used for model validation.

2.2.3. Airborne P-SAR

In this study, we obtained wall-to-wall Airborne P-SAR data for the study area, employing processing methodologies in accordance with those detailed by Fan et al. [41]. The primary steps encompass (1) multi-looking and refined Lee polarization filtering to generate multi-look complex (MLC) PolSAR data, followed by the application of a 5 × 5 moving window filter to the MLC results to further attenuate the impact of speckle noise; (2) the integration of SAR imaging orbit information into digital elevation model (DEM) data and the geocoding of the PolSAR data based on the radar–Doppler geolocation model to establish the correspondence between the SAR image slant range space and the geographic coordinate space while also obtaining relevant angular parameters that describe the local imaging geometry, such as projection angles and local incidence angles; (3) the implementation of a three-stage topographic–radiometric correction method [42] to address the polarization azimuth angle, effective scattering area, and angle effects, in order to compensate for the effects of terrain undulations on polarization characteristics; (4) conducting Yamaguchi decomposition [43] on the topographic–radiometric corrected polarization coherence matrix to extract various polarization decomposition components that, along with backscatter intensity, form the polarization feature set for each band (Table 2). The final processed multi-band polarimetric SAR data and features have a pixel resolution of 25 m × 25 m, ensuring spatial congruence with the sample plots.

2.2.4. GEDI Data

The GEDI LiDAR system, launched aboard the International Space Station (ISS) in December 2018, comprises three lasers, with one laser divided into two coverage beams and the other two operating at full power. The coverage beams and full-power beams are able to detect the ground through 95% and 98% forest canopy cover, respectively [44]. GEDI creates spot trajectories on the Earth’s surface through beam jitter, with each trajectory consisting of circular footprints that are 60 m along-track and 25 m in diameter, and the distance between adjacent trajectories is approximately 600 m, with a scanning swath width of about 4.2 km [45,46]. These data have been widely applied to the estimation of AGB in forests [47,48]. The data used in this study are cover, canopy height, and height quantiles, i.e., rh10, rh20, rh30, rh40, rh50, rh60, rh70, rh75, rh80, rh85, rh90, rh95, and rh100 from the GEDI Level 2A and Level 2B data products. The data acquisition period was from 1 July 2022 to 30 September 2022, and some invalid spots were filtered based on the parameters inherent in the data [49]. The data were downloaded and preprocessed using the rGEDI package in R [50], with outliers removed, resulting in data features from 13,010 location spots, as shown in Figure 2a, of which 5715 GEDI spots are located within the strips (Figure 2b).

2.2.5. Sentinel-2A Data

Two scenes of sentinel-2A multispectral data were downloaded from the Copernicus Open Access Hub website (https://dataspace.copernicus.eu/ (accessed on 10 May 2023)). Due to the influence of cloudy and wet weather, it was not possible to obtain imagery that is perfectly synchronized with the ground plot data acquisition time; therefore, only available imagery from the growing season with a close time match was acquired, dated 8 July 2022. The European Space Agency’s open-source SNAP software (https://step.esa.int/main/download/snap-download/ (Version 8.0.0.0, accessed on 20 May 2023)) was utilized to perform radiometric calibration, atmospheric correction, and orthorectification on each scene. After these preprocessing steps had been performed, the pixel size of the Sentinel-2A data were set to 25 m × 25 m, consistent with the size of the sample plots; then, the two preprocessed scenes were mosaicked to obtain multispectral remote sensing data covering the entire study area. Common vegetation indices were selected as candidate features (Table 3), and an additional novel vegetation index, the kernel-normalized difference vegetation index (KNDVI), was included. This index is derived from the concept of kernel functions, and previous studies have shown that it can significantly improve the estimation accuracy of forest carbon stocks [51]. The corresponding kernel function vegetation indices, KNDVIre1, KNDVIre2, and KNDVIre3, were calculated with the methods shown in [51].

2.3. Methods

2.3.1. Overview

We propose the application of RK at the third stage of a three-stage model, and in this study, we aim to confirm whether it can promote accuracy and reduce uncertainty in estimation results. We compare it with CMB, GHMB, RK-GHMB, and 3sGHMB to analyze their strengths and weaknesses in practical application scenarios. The flowchart for this study is depicted in Figure 3, and the data frameworks for the five methods are as follows:

Case A: CMB using plot data and wall-to-wall airborne P-SAR and Sentinel-2A data.

Case B: GHMB using plot data, sampled ALS strip data, and wall-to-wall airborne P-SAR and Sentinel-2A data.

Case C: RK-GHMB using the same data as in Case B.

Case D: 3sGHMB using plot data, sampled ALS strip data, GEDI data, and wall-to-wall airborne P-SAR and Sentinel-2A data.

Case E: RK-3sGHMB using the same data as in Case D.

In this study, the N grid of the study area is denoted as population U. The grids for Tidal Slag Forest Farm constitute the subpopulation U₁, with a total number of 821,909, denoted as N₁. The grids for Upper Yangge Forest Farm constitute the subpopulation U₂, with a total number of 717,183, denoted as N₂. For the sake of clarity in subsequent discussions, the model within the CMB is designated as $F_{ms}$; the first-stage models for the GHMB, RK-GHMB, 3sGHMB, and RK-3sGHMBs are represented by $F_{ml}$; the second-stage models for GHMB and RK-GHMB are denoted by $G_{ls}$ and $G_R$ (in Figure 3, $G_R$ is denoted as “model G_ls + RK”), respectively; the second-stage models for 3sGHMB and RK-3sGHMB are denoted by $Q_{lg}$; and their third-stage models are labeled as $Q_{gs}$ and $Q_R$ (in Figure 3, $Q_R$ is denoted as “model $Q_{gs}$ + RK”). The dataset constituted by the sample plots is represented by $S_I$, which is located in the airborne LiDAR strips and includes the measured forest AGB, ALS features, and wall-to-wall features. The sample size n for $S_I$ is 80. The sample plots outside the ALS strips are used for validation, with a count of 23. From the GEDI spots within the ALS strips, a random selection is made to construct the $S_{II}$ dataset, which includes ALS features, GEDI features, and wall-to-wall features, and the sample size k is 5000. Since there are no measured data at the GEDI spot locations, data for an additional GEDI spot for which LiDAR-based biomass is available are randomly selected to form a validation dataset to verify the $G_{lg}$ model; the sample size of the dataset is 1000. From all GEDI spots, random selection is performed to create the $S_{III}$ dataset for model training, which includes GEDI features and wall-to-wall features. The sample size of $S_{III}$ is h = 10,000.

2.3.2. Feature Selection

We employed a genetic algorithm for feature selection, which is a heuristic method that simulates natural selection to solve optimization problems and has a wide range of uses in the field of remote sensing [60]. The R² and RMSE of internal and external validation samples were used as evaluation metrics. We implemented genetic algorithm feature selection for each model, utilizing the gafs function within the caret package in R [61].

2.3.3. Case A: The CMB

This case follows the well-established model-based inference [9,17]. The general form of the model is as follows:

(1)$F_{ms}:\mathit{y}=f({\mathit{X}}_{S_I},{\mathit{\alpha }}_{S_I})+\mathit{\epsilon },\mathit{\epsilon }~{\rm N}(0,\mathit{\Omega })$

where $\mathit{y}$ is the vector of measured forest AGB data, ${\mathit{X}}_{S_I}$ is an n × q dimension matrix of wall-to-wall features with respect to $S_I$, ${\mathit{\alpha }}_{S_I}$ is the corresponding parameter vector, and $\mathit{\epsilon }$ is the random error following ${\rm N}(0,\mathit{\Omega })$. f represents the model form for fitting. In this study, we used logistic regression for model fitting and estimation due to its ability to ensure non-negative values and provide a certain physical significance [15]. The same model form was also used in Case B, Case C, Case D, and Case E. The pixel uncertainty is calculated based on the methods proposed by Saarela et al. [9] and McRoberts et al. [62], using the following formula:

(2)$RMSE({\hat{y}}_{F_{m{s}_i}})=\sqrt{{\tilde{\mathit{X}}}_i^{T}Cov({\hat{\mathit{\alpha }}}_{S_I}){\tilde{\mathit{X}}}_i+V\left({\epsilon }_i\right)}$

(3)${\tilde{\mathit{X}}}_i={\displaystyle \frac{\partial F_{ms}({\hat{\mathit{\alpha }}}_{S_I},{\mathit{X}}_i)}{\partial {\hat{\mathit{\alpha }}}_{S_I}}}$

where $RMSE({\hat{y}}_{F_{m{s}_i}})$ represents the uncertainty of the model $F_{ms}$ at map unit i; ${\hat{y}}_{F_{m{s}_i}}$ is the predicted AGB using model $F_{ms}$ for map unit i; ${\tilde{\mathit{X}}}_i$ is a (q + 1)-length vector composed of the partial derivatives of the model $F_{ms}$ at unit i with respect to ${\hat{\mathit{\alpha }}}_{S_I}$, where q is the number of parameters. $Cov({\mathit{\alpha }}_{S_I})={({\tilde{\mathit{X}}}_{S_I}^T\mathit{\Omega }{\tilde{\mathit{X}}}_{S_I})}^{-1}$, and ${\tilde{\mathit{X}}}_{S_I}$ is an n × q matrix of partial derivatives of the model $F_{ms}$ based on the dataset $S_I$ with respect to ${\hat{\mathit{\alpha }}}_{S_I}$, and $\mathit{\Omega }$ is the variance–covariance matrix of the residuals. Without considering spatial autocorrelation, the diagonal elements are constituted by $V\left({\epsilon }_i\right)$, which is the variance of ${\epsilon }_i$. It can be calculated by the method proposed by McRoberts [62], with the formula $V\left({\epsilon }_i\right)=a{{\overline{\hat{y}}}_{S_I}}_i^b+\delta$, where a and b are the model parameters, $\delta$ is the residuals term, ${\overline{\hat{y}}}_{S_I}$ and is the mean of the predicted forest AGB for each group after stratification.

The population mean is as follows:

(4)${\mu }_{F_{ms}}={\displaystyle \frac{1}{N}}{{\displaystyle {\sum }_{i=1}^{N}y}}_{F_{m{s}_i}}$

with the variance of ${\mu }_{F_{ms}}$ being

(5)$Var({\mu }_{F_{ms}})={\mathit{l}}_U^T{\tilde{\mathit{X}}}_{U}Cov({\hat{\mathit{\alpha }}}_{S_I}){\tilde{\mathit{X}}}_U^T{\mathit{l}}_U$

where ${\mathit{l}}_U$ is an N-length vector with each element being 1/N, and ${\tilde{\mathit{X}}}_U$ is a partial derivative matrix of $F_{ms}$ with respect to the population U. When estimating the population mean and variance for the two sub-regions, U and N are replaced by U₁, N₁ and U₁, N₂.

2.3.4. Case B: The GHMB

In the case of GHMB, the forest AGB values estimated by sampling ALS data serve as extended samples for the estimation of the full-coverage forest AGB [9]. Similarly to Equation (1), the general form of the model for the first stage is as follows:

(6)$F_{ml}:\mathit{y}=f({\mathit{P}}_{S_I},{\mathsf{\beta }}_{S_I})+\mathsf{o},\mathsf{o}~{\rm N}(0,{\mathit{\Omega }}^{*})$

where ${\mathit{P}}_{S_I}$ is an n × p dimensional matrix of ALS features with respect to $S_I$; ${\mathit{\beta }}_{S_I}$ is the corresponding model parameter vector; and $\mathsf{o}$ is the random error following ${\rm N}(0,\mathit{\Omega })$.

The second model general form is as follows:

(7)$G_{ls}:{\mathit{y}}_{F_{ml}}=f({\mathit{X}}_{S_{II}}^{*},{\mathit{\alpha }}_{S_{II}})+\mathit{v},\mathit{v}~{\rm N}(0,\mathit{C})$

where ${\mathit{y}}_{F_{ml}}$ is the predicted AGB estimated by the model $F_{ml}$, ${\mathit{X}}_{S_{II}}^{*}$ is a k × g matrix of wall-to-wall data features with respect to $S_{II}$; ${\mathit{\alpha }}_{S_{II}}$ is a vector of the corresponding model parameters; and $\mathit{v}$ is a vector of random error following ${\rm N}(0,\mathit{C})$.

The uncertainty is calculated using the same method as Equation (2):

(8)$RMSE\left({\hat{y}}_{G_{l{s}_i}}\right)=\sqrt{{\tilde{\mathit{X}}}_i^{*T}Cov({\hat{\mathit{\alpha }}}_{S_{II}}){\tilde{\mathit{X}}}_i^{*}+V\left(v_i\right)}$

(9)${\tilde{\mathit{X}}}_i^{*}={\displaystyle \frac{\partial G_{ls}\left({\hat{\mathit{\alpha }}}_{S_{II}},{\mathit{X}}_i^{*}\right)}{\partial {\hat{\mathit{\alpha }}}_{S_{II}}}}$

where ${\hat{y}}_{G_{l{s}_i}}$ is the predicted AGB using $G_{ls}$ at the unit i; ${\tilde{\mathit{X}}}_i^{*}$ is a (g + 1)-length vector of partial derivatives of the $G_{ls}$ model with respect to ${\hat{\mathit{\alpha }}}_{S_{II}}$; $Cov({\hat{\mathit{\alpha }}}_{S_{II}})$ is the estimated covariance matrix for ${\hat{\mathit{\alpha }}}_{S_{II}}$; $V(v_i)$ is the variance of random error $v_i$ at unit i, with the calculation approach being identical to that of $V\left({\epsilon }_i\right)$ and the model form being $V(v_i)=a{\overline{\hat{y}}}_{S_{I{I}_i}}^2+b{\overline{\hat{y}}}_{S_{I{I}_i}}+c+{\delta }_2$ with a, b, and c as corresponding parameters; and ${\delta }_2$ is the random error. According to Saarela et al. [9], $Cov({\hat{\mathit{\alpha }}}_{S_{II}})$ in the two-stage method can be expressed as follows:

(10)$Cov({\hat{\mathit{\alpha }}}_{S_{II}})={({\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{-1}{\tilde{\mathit{X}}}_{S_{II}}^{*})}^{-1}+{\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{-1}Cov({\mathit{y}}_{F_{m{l}_{S_{II}}}}){\mathit{C}}^{-1}{\tilde{\mathit{X}}}_{S_{II}}^{*}{({\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{-1}{\tilde{\mathit{X}}}_{S_{II}}^{*})}^{-1}$

where ${\tilde{\mathit{X}}}_{S_{II}}^{*}$ is the k × g matrix of partial derivatives of the model $F_{ms}$ based on dataset $S_{II}$ with respect to ${\hat{\mathit{\alpha }}}_{S_{II}}$. $Cov({\mathit{y}}_{F_{m{l}_{S_{II}}}})={\tilde{\mathit{P}}}_{S_{II}}Cov({\hat{\mathit{\beta }}}_{S_I}){\tilde{\mathit{P}}}_{S_{II}}^T$, ${\mathit{y}}_{F_{m{l}_{S_{II}}}}$ is the predicted AGB using the $F_{ml}$ model with respect to $S_{II}$; ${\tilde{\mathit{P}}}_{S_{II}}$ is the k × p matrix of partial derivatives of the model $F_{ml}$ based on dataset $S_{II}$ with respect to ${\hat{\mathit{\beta }}}_{S_I}$; $Cov({\hat{\mathit{\beta }}}_{S_I})={({\tilde{\mathit{P}}}_{S_I}^T{\mathit{\Omega }}^{*}{\tilde{\mathit{P}}}_{S_I})}^{-1}$, with ${\tilde{\mathit{P}}}_{S_I}$ being the n × p matrix of partial derivatives of the model $F_{ml}$ based on dataset $S_I$ with respect to ${\hat{\mathit{\beta }}}_{S_I}$.

The population mean estimated by GHMB can be written as follows:

(11)${\mu }_{G_{ls}}={\displaystyle \frac{1}{N}}{\displaystyle {\sum }_{i=1}^N{\hat{y}}_{G_{l{s}_i}}}$

The corresponding variance of ${\mu }_{G_{ls}}$ is

(12)$Var({\mu }_{G_{ls}})={\mathit{l}}_U^T{\tilde{\mathit{X}}}_U^{*}[{({\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{-1}{\tilde{\mathit{X}}}_{S_{II}}^{*})}^{-1}+{\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{-1}Cov({\mathit{y}}_{F_{m{l}_{S_{II}}}}){\mathit{C}}^{-1}{\tilde{\mathit{X}}}_{S_{II}}^{*}{({\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{-1}{\tilde{\mathit{X}}}_{S_{II}}^{*})}^{-1}]{\tilde{\mathit{X}}}_U^{*T}{\mathit{l}}_U$

where ${\tilde{\mathit{X}}}_U^{*}$ is a partial derivative matrix of $G_{ls}$ with respect to the population U. When estimating the population mean and variance for the two sub-regions, U and N are replaced with U₁, N₁ and U₂, N₂.

2.3.5. Case C: The RK-GHMB

Zhao et al. [27] proposed an improved GHMB model named RK-GHMB, which uses an RK model in the second stage of GHMB, considering the spatial correlation of residual errors. The general form of the first stage is identical to Equation (6). Assuming that the residual, which cannot be predicted using the trend model of (1), is second-order stationarity, the model for the second stage is as follows:

(13)$G_R:{\hat{y}}_{G_{ls}-R{K}_i}={\hat{y}}_{G_{l{s}_i}}+V_{kriging}$

where ${\hat{y}}_{G_{ls}-R{K}_i}$ is the predicted AGB of the model $G_R$, and $V_{kriging}$ is the residual result interpolated by ordinary kriging, which can be calculated as follows:

(14)$V_{kriging}={\displaystyle \sum _{j=1}^k{d}_j}{v}_j,{\displaystyle \sum _{j=1}^k{d}_j}=1$

where $d_j$ is the kriging weight of the residual $v_j$ at the jth index; the total weights must add up to one. The $d_j$ can be calculated as follows:

(15)$\mathit{d}={\mathit{C}}^{*-1}\cdot {\mathit{c}}_{G_{l{s}_i}}$

(16)${\mathit{c}}_{G_{l{s}_i}}=\left[\begin{array}{c}Cov(v_i,v_1)\\ Cov(v_i,v_2)\\ ⋮\\ Cov(v_i,v_k)\end{array}\right]$

where $\mathit{d}$ is an m-length vector of $d_j$ kriging weights; ${\mathit{C}}^{*}$ is a k × k covariance matrix of residuals for $S_{II}$, with its diagonal elements constituted by $V(v_i)$ and the off-diagonal elements calculated using $V(v_i)$, $V(v_j)$, and the spatial correlation ${\rho }_{ij}$; i.e., $Cov(v_i,v_j)=\sqrt{V(v_i)\cdot V(v_j)}\cdot {\rho }_{ij}$. ${\rho }_{ij}$ is expressed as ${\rho }_{ij}=1-{\displaystyle \frac{\gamma (d_{ij})}{C_0+C_1}}$, where $\gamma (d_{ij})$ is a semi-variogram, and its form is either exponential, gaussian, spherical, or linear. $C_0$ and $C_1$ are the nugget and partial sill, respectively. ${\mathit{c}}_{G_{l{s}_i}}$ is a k-length vector of covariance for the unit map i.

Thus, the uncertainty of RK-GHMB can be written as follows:

(17)$RMSE(y_{G_{ls}-R{K}_i})=\sqrt{\begin{array}{c}{({\tilde{\mathit{X}}}_i^{*}-{\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{*-1}{\mathit{C}}_{G_{l{s}_i}})}^T\\ \{{({\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{*-1}{\tilde{\mathit{X}}}_{S_{II}}^{*})}^{-1}+{\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{*-1}Cov({\mathit{y}}_{F_{m{l}_{S_{II}}}}){\mathit{C}}^{*-1}{\tilde{\mathit{X}}}_{S_{II}}^{*}{({\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{*-1}{\tilde{\mathit{X}}}_{S_{II}}^{*})}^{-1}\}\\ ({\tilde{\mathit{X}}}_i^{*}-{\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{*-1}{\mathit{C}}_{G_{l{s}_i}})+V(v_i)-{\mathit{c}}_{G_{l{s}_i}}^T{\mathit{C}}^{*-1}{\mathit{c}}_{G_{l{s}_i}}\end{array}}$

The population mean estimated by RK-GHMB is

(18)${\mu }_{G_{ls}-R{K}_i}={\displaystyle \frac{1}{N}}{{\displaystyle {\sum }_{i=1}^N\hat{y}}}_{G_{ls}-R{K}_i}$

with the variance estimator as

(19)$Var({\mu }_{G_{ls}-R{K}_i})=\begin{array}{c}{l}_U^T\{({\tilde{\mathit{X}}}_U^{*}-{\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{*-1}{\mathit{c}}_{G_{l{s}_U}})\\ \{{({\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{*-1}{\tilde{\mathit{X}}}_{S_{II}}^{*})}^{-1}+{\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{*-1}Cov({\mathit{y}}_{F_{m{l}_{S_{II}}}}){\mathit{C}}^{*-1}{\tilde{\mathit{X}}}_{S_{II}}^{*}{({\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{*-1}{\tilde{\mathit{X}}}_{S_{II}}^{*})}^{-1}\}\\ {({\tilde{\mathit{X}}}_U^{*}-{\tilde{\mathit{X}}}_{S_{II}}^{*T}{\mathit{C}}^{*-1}{\mathit{c}}_{G_{l{s}_U}})}^T\}{\mathit{l}}_U\end{array}$

where the parameters are the same as in Equation (12).

2.3.6. Case D: The 3sGHMB

The 3sGHMB model includes three stages [8]. The general form of the first stage is identical to the GHMB, and the form of the second stage is as follows:

(20)$Q_{lg}:{\mathit{y}}_{F_{ml}}=f({\mathit{Z}}_{S_{II}},{\mathit{\eta }}_{S_{II}})+\mathit{e},\mathit{e}~{\rm N}(0,\mathit{\Sigma })$

where ${\mathit{Z}}_{S_{II}}$ is a k × l dimension matrix of GEDI features with respect to $S_{II}$; ${\mathit{\eta }}_{S_{II}}$ is the corresponding model parameter vector; and $\mathit{e}$ is the random error following ${\rm N}(0,\mathit{\Sigma })$. The general form of the third stage can be written as follows:

(21)$Q_{gs}:{\mathit{y}}_{G_{ls}}=f({\mathit{X}}_{S_{III}}^{**},{\mathit{\alpha }}_{S_{III}})+\mathit{\theta },\mathit{\theta }~{\rm N}(0,\mathit{\Lambda })$

where ${\mathit{X}}_{S_{III}}^{**}$ is an h × r dimension matrix of wall-to-wall features with respect to $S_{III}$; ${\hat{\mathit{\alpha }}}_{S_{III}}$ is the corresponding model parameter vector; and $\mathit{\theta }$ is the random error following ${\rm N}(0,\mathit{\Lambda })$.

The estimation of uncertainty is analogous to Equation (8):

(22)$RMSE\left(y_{Q_{g{s}_i}}\right)=\sqrt{{\tilde{\mathit{X}}}_i^{**T}Cov\left({\hat{\mathit{\alpha }}}_{S_{III}}\right){\tilde{\mathit{X}}}_i^{**}+V\left({\theta }_i\right)}$

(23)${\tilde{\mathit{X}}}_i^{**}={\displaystyle \frac{\partial Q\left({\hat{\mathit{\alpha }}}_{S_{III}},{\mathit{X}}_i^{**}\right)}{\partial {\hat{\mathit{\alpha }}}_{S_{III}}}}$

where $y_{Q_{g{s}_i}}$ is the predicted AGB using $Q_{gs}$ at the map unit i; ${\tilde{\mathit{X}}}_i^{**}$ is a (r + 1)-length vector of partial derivatives of the $Q_{gs}$ model with respect to ${\hat{\mathit{\alpha }}}_{S_{III}}$; $Cov({\hat{\mathit{\alpha }}}_{S_{III}})$ is the estimated covariance matrix for ${\hat{\mathit{\alpha }}}_{S_{III}}$; $V\left({\theta }_i\right)$ is the variance of random error ${\theta }_i$ at the map unit i, with a calculation approach identical to that of $V\left({\epsilon }_i\right)$ and the model form $V\left({\theta }_i\right)=a{\overline{\hat{\mathit{y}}}{S}_{III}}_i^b+{\delta }_4$ with a and b as corresponding parameters; and ${\delta }_4$ is the random error. Analogous to $Cov({\hat{\mathit{\alpha }}}_{S_{II}})$, $Cov({\hat{\mathit{\alpha }}}_{S_{III}})$ can be expressed as follows:

(24)$Cov({\hat{\mathit{\alpha }}}_{S_{III}})={({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}+{\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}Cov({\mathit{y}}_{Q_{l{g}_{S_{III}}}}){\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}^{**}{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}$

where ${\tilde{\mathit{X}}}_{S_{III}}^{**}$ is an h × r matrix of partial derivatives of the model $Q_{gs}$ based on dataset $S_{III}$ with respect to ${\hat{\mathit{\alpha }}}_{S_{III}}$. $\mathit{\Lambda }$ is a diagonal matrix; each element is estimated using the model form presented as $V\left({\theta }_i\right)$. The covariance of ${\mathit{y}}_{Q_{l{g}_{S_{III}}}}$ can be expressed using ${\tilde{\mathit{Z}}}_{S_{III}}$ and ${\hat{\mathit{\eta }}}_{S_{II}}$; i.e., $Cov({\mathit{y}}_{Q_{l{g}_{S_{III}}}})={\tilde{\mathit{Z}}}_{S_{III}}Cov({\hat{\mathit{\eta }}}_{S_{II}}){\tilde{\mathit{Z}}}_{S_{III}}^T$. ${\tilde{\mathit{Z}}}_{S_{III}}$ is an h × r matrix of partial derivatives of the model $Q_{lg}$ based on dataset $S_{III}$ with respect to ${\hat{\mathit{\eta }}}_{S_{II}}$. The covariance matrix of ${\hat{\mathit{\eta }}}_{S_{II}}$ was estimated as follows:

(25)$Cov({\hat{\mathit{\eta }}}_{S_{II}})={({\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}})}^{-1}+{\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}Cov({\mathit{y}}_{F_{m{l}_{S_{II}}}}){\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}}{({\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}})}^{-1}$

where ${\tilde{\mathit{Z}}}_{S_{II}}$ is a k × r matrix of partial derivatives of the model $Q_{lg}$ based on dataset $S_{II}$ with respect to ${\hat{\mathit{\eta }}}_{S_{II}}$; the calculation method employed in $Cov({\mathit{y}}_{F_{m{l}_{S_{II}}}})$ is consistent with the approach utilized in Equation (10). By replacing $Cov({\mathit{y}}_{F_{m{l}_{S_{II}}}})$ and $Cov({\hat{\mathit{\eta }}}_{S_{II}})$ with its estimators, the covariance matrix estimator of ${\hat{\mathit{\alpha }}}_{S_{III}}$ can be expressed as

(26)$Cov({\hat{\mathit{\alpha }}}_{S_{III}})=\begin{array}{c}{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}+\\ {\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}{\tilde{\mathit{Z}}}_{S_{III}}{({\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}})}^{-1}{\tilde{\mathit{Z}}}_{S_{III}}^T{\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}+\\ {\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}{\tilde{\mathit{Z}}}_{S_{III}}{\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\hat{\mathit{Z}}}_{S_{II}}{({\tilde{\mathit{P}}}_{S_I}^T\mathit{\Omega }{\tilde{\mathit{P}}}_{S_I})}^{-1}{\hat{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}}\\ {({\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}})}^{-1}{\tilde{\mathit{Z}}}_{S_{III}}^T{\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}^{**}{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}\end{array}.$

The population mean estimated by 3sGHMB is

(27)${\mu }_{Q_{gs}}={\displaystyle \frac{1}{N}}{{\displaystyle {\sum }_{i=1}^N\hat{y}}}_{Q_{gs}}$

with the variance estimator being

(28)$Var({\mu }_{Q_{gs}})=\begin{array}{c}{\mathit{l}}_U^T{\tilde{\mathit{X}}}_U^{**}\{{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}+\\ {\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}{\tilde{\mathit{Z}}}_{S_{III}}{({\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}})}^{-1}{\tilde{\mathit{Z}}}_{S_{III}}^T{\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}^{**}{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}+\\ {\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}{\tilde{\mathit{Z}}}_{S_{III}}{\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\hat{\mathit{Z}}}_{S_{II}}{({\tilde{\mathit{P}}}_{S_I}^T\mathit{\Omega }{\tilde{\mathit{P}}}_{S_I})}^{-1}{\hat{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}}\\ {({\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}})}^{-1}{\tilde{\mathit{Z}}}_{S_{III}}^T{\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}^{**}{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}\}{\tilde{\mathit{X}}}_U^{**T}{\mathit{l}}_U\end{array}$

where ${\tilde{\mathit{X}}}_U^{**}$ is a partial derivative matrix of $Q_{gs}$ with respect for the population U.

2.3.7. Case E: The Proposed RK-3sGHMB

Similarly to Equation (13), we replaced the third general form $Q_{gs}$ as follows:

(29)$Q_R:y_{Q_{gs}-R{K}_i}=y_{Q_{g{s}_i}}+V_{kriging}^{*}$

where $y_{Q_{g{s}_i}}$ is the predicted AGB of the model $Q_{gs}$, and $V_{kriging}^{*}$ is the residual result interpolated by ordinary kriging, which can be calculated as follows:

(30)$V_{kriging}^{*}={\displaystyle \sum _{j=1}^h{w}_j}{\theta }_j,{\displaystyle \sum _{j=1}^h{w}_j}=1$

where $w_j$ is the kriging weight of the residual ${\theta }_j$ at the jth index; the total weights must add up to one. The $w_j$ can be calculated as follows:

(31)$\mathit{w}={\mathit{\Lambda }}^{*-1}\cdot {\mathit{c}}_{Q_{g{s}_i}}$

(32)${\mathit{c}}_{Q_{g{s}_i}}=\left[\begin{array}{c}Cov({\theta }_i,{\theta }_1)\\ Cov({\theta }_i,{\theta }_2)\\ ⋮\\ Cov({\theta }_i,{\theta }_h)\end{array}\right]$

where $\mathit{w}$ is an h-length vector of the $w_j$ kriging weights; ${\mathit{\Lambda }}^{*}$ is an h × h covariance matrix of residuals for $S_{III}$, with its diagonal elements constituted by $V({\theta }_i)$, its off-diagonal elements calculated using $V({\theta }_i)$, $V({\theta }_j)$, and the spatial correlation ${\rho }_{ij}^{*}$; i.e., $Cov({\theta }_i,{\theta }_j)=\sqrt{V({\theta }_i)\cdot V({\theta }_j)}\cdot {\rho }_{ij}^{*}$. ${\rho }_{ij}^{*}$ is expressed as ${\rho }_{ij}^{*}=1-{\displaystyle \frac{\gamma (w_{ij})}{C_0^{*}+C_1^{*}}}$, where $\gamma (w_{ij})$ is a semi-variogram in which the form is determined in the same manner as in $\gamma (d_{ij})$. $C_0^{*}$ and $C_1^{*}$ are the nugget and partial sill, respectively. ${\mathit{c}}_{Q_{g{s}_i}}$ is an h-length vector of covariance for the unit map i.

Thus, the uncertainty of RK-3sGHMB for the map unit i can be written as follows:

(33)$RMSE(y_{Q_{gs}-R{K}_i})=\sqrt{\begin{array}{c}{({\tilde{\mathit{X}}}_i^{**}-{\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\mathit{c}}_{Q_{g{s}_i}})}^T\\ \left\{\begin{array}{c}{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}+\\ {\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{Z}}}_{S_{III}}{({\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}})}^{-1}{\tilde{\mathit{Z}}}_{S_{III}}^T{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{X}}}_{S_{III}}{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}+\\ {\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{Z}}}_{S_{III}}{\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\hat{\mathit{Z}}}_{S_{II}}{({\tilde{\mathit{P}}}_{S_I}^T\mathit{\Omega }{\tilde{\mathit{P}}}_{S_I})}^{-1}{\hat{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}}\\ {({\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}})}^{-1}{\tilde{\mathit{Z}}}_{S_{III}}^T{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{X}}}_{S_{III}}^{**}{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}\end{array}\right\}\\ ({\tilde{\mathit{X}}}_i^{**}-{\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\mathit{c}}_{Q_{g{s}_i}})+V({\theta }_i)-{\mathit{c}}_{Q_{g{s}_i}}^T{\mathit{\Lambda }}^{*-1}{\mathit{c}}_{Q_{g{s}_i}}\end{array}}$

The population mean estimated by RK-3sGHMB can be expressed as

(34)${\mu }_{Q_{gs}-RK}={\displaystyle \frac{1}{N}}{{\displaystyle {\sum }_{i=1}^N\hat{y}}}_{Q_{gs}-R{K}_i}$

with the variance estimator being

(35)$Var({\mu }_{Q_{gs}-R{K}_i})=\begin{array}{c}{\mathit{l}}_U^T\{({\tilde{\mathit{X}}}_U^{**}-{\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\mathit{c}}_{Q_{g{s}_U}})\\ \left\{\begin{array}{c}{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}+\\ {\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{Z}}}_{S_{III}}{({\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}})}^{-1}{\tilde{\mathit{Z}}}_{S_{III}}^T{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{X}}}_{S_{III}}{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}+\\ {\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{Z}}}_{S_{III}}{\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\hat{\mathit{Z}}}_{S_{II}}{({\tilde{\mathit{P}}}_{S_I}^T\mathit{\Omega }{\tilde{\mathit{P}}}_{S_I})}^{-1}{\hat{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}}\\ {({\tilde{\mathit{Z}}}_{S_{II}}^T{\mathit{\Sigma }}^{-1}{\tilde{\mathit{Z}}}_{S_{II}})}^{-1}{\tilde{\mathit{Z}}}_{S_{III}}^T{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{X}}}_{S_{III}}^{**}{({\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\tilde{\mathit{X}}}_{S_{III}}^{**})}^{-1}\end{array}\right\}\\ {({\tilde{\mathit{X}}}_U^{**}-{\tilde{\mathit{X}}}_{S_{III}}^{**T}{\mathit{\Lambda }}^{*-1}{\mathit{c}}_{Q_{g{s}_U}})}^T\}{\mathit{l}}_U\end{array}$

where the parameters are the same as in Equation (28). If there is no spatial correlation in the dataset $S_{III}$, then ${\mathit{c}}_{Q_{g{s}_i}}$ and ${\mathit{c}}_{Q_{g{s}_U}}$ are approximately equal to the zero vector, and then Equations (34) and (35) are approximately equal to Equations (22), (27) and (28).

2.3.8. Evaluation Criteria for Modeling and Up-Scaling

We selected the adjusted coefficient of determination ($R_{adj}^2$), root mean square error (RMSE), and estimation accuracy (EA) as the evaluation indicators for the model fitting. The calculation formulas are given as (36), (37), (38), and (39).

(36)$R^2=1-{\displaystyle \frac{{{\displaystyle {\sum }_{i=1}^n\left(y_i-{\hat{y}}_i\right)}}^2}{{{\displaystyle {\sum }_{i=1}^n\left(y_i-\overline{y}\right)}}^2}}$

(37)$R_{adj}^2=1-\left(1-R^2\right){\displaystyle \frac{n-1}{n-1-P}}$

(38)$RMSE=\sqrt{{\displaystyle \frac{{{\displaystyle {\sum }_{i=1}^n\left(y_i-{\hat{y}}_i\right)}}^2}{n-1}}}$

(39)${\rm E}A=\left(1-{\displaystyle \frac{RMSE}{\overline{y}}}\right)\times 100\%$

where $y_i$ represents the reference data; ${\hat{y}}_i$ represents the predicted data; $\overline{y}$ refers to the mean value of the reference data; and n is the number of reference data.

The precision of the population prediction for the five methods was evaluated using standard error (SE) and estimation precision (P) indicators, as shown in Formulas (40) and (41) [63].

(40)$SE=\sqrt{V({\overline{\mu }}_{model})}$

(41)$P=(1-{\displaystyle \frac{t\cdot SE}{\overline{y}}})\times 100\%$

where $V({\overline{\mu }}_{model})$ is the variance estimation of mean value by model; t is a reliability indicator, with a 95% confidence interval being t = 1.96.

3. Results

3.1. Feature Selection Results

Genetic algorithms were employed to perform feature selection for each model. To prevent overfitting, each model’s feature selection was subjected to 5-fold and 10-fold cross-validation, conducted 50 times. The optimal number of iterations and features for the final genetic algorithm was determined by calculating the RMSE and R² of the external samples. For the model $F_{ml}$, the optimal number of iterations was 46, with the optimal features being Hvar, H98, and HCD (Figure 4a). For the model $F_{ms}$, the optimal number of iterations was 49, with the optimal features being PHV, DVIre1, EVI, and EVIre1 (Figure 4b). For the model $Q_{lg}$, the optimal number of iterations was 39, with the optimal features being rh90, rh10, rh50, and cover (Figure 4c). For the model $G_{ls}$, the optimal number of iterations was 11, with the optimal features being KNDVIre1, PHV, EVI, EVIre1, and Vol (Figure 4d). For the model $Q_{gs}$, the optimal number of iterations was 20, with the optimal features being PHV, Hlx, EVIre1, Vol, EVIre3, KNDVIre1, and PHH (Figure 4e).

3.2. Model Fitting Results and Accuracy Evaluation

Independent samples were employed to validate each model. The results indicate that models $F_{ml}$ and $Q_{lg}$ exhibit higher accuracy; both $R_{adj}^2$ values exceed 0.80, with estimation accuracies close to 80%. This suggests that using ALS-based AGB and GEDI-based AGB for multi-stage methods is feasible (Figure 5a,b). For the single-stage model using the CMB method, the model $R_{adj}^2$ is 0.37, with an RMSE of 33.95 t/ha and an accuracy of 63.28% (Figure 5c). In the case of multi-stage models, as the sample size increases, the precision of the GHMB model is improved, with an $R_{adj}^2$ reaching 0.38, an RMSE of 33.72 t/ha, and an accuracy reaching 63.53% (Figure 5d). However, the model accuracy of the 3sGHMB method decreases, with an $R_{adj}^2$ dropping to 0.27, an RMSE of 36.58 t/ha, and an accuracy of 60.43% (Figure 5f). The ratio of the nugget to the sill (C₀/(C₀ + C₁)), which represents the spatial dependence structure, is estimated for RK-GHMB and RK-3sGHMB. In general, a ratio below 25% indicates a strong spatial dependence structure. A ratio between 25% and 75% suggests a moderate spatial dependence structure, while a ratio exceeding 75% implies a weak spatial dependence structure [32]. In this study, all the ratios are less than 75% (Table 4, Figure 6). The global Moran’s I of the two models is 0.20 and 0.40, with both p-values being 2.2 e16, also indicating the presence of spatial autocorrelation effects. The optimal residual variogram model forms for both models are exponential, with a range of 2625.28 m and 16,069.84 m, respectively, which means that the RK model could have a limited effect beyond these distances (Table 4). After considering spatial correlation, the model accuracy improved. The $R_{adj}^2$ of the RK-GHMB method increased to 0.60, with an RMSE reduction to 27.07 t/ha and an improvement in accuracy to 70.72% (Figure 5e). The $R_{adj}^2$ of the RK-3sGHMB method increased to 0.55, with an RMSE reduction to 28.55 t/ha and an accuracy improvement to 69.13% (Figure 5g). The fitting formulas for each model are presented in Table 5.

3.3. Model Estimation and Uncertainty Analysis

We estimated the uncertainty at the pixel level using the five inferences; the spatial distribution is shown in Figure 7a–e. It can be seen that the uncertainty range of the CMB is 2–48.3 t/ha. The uncertainty of the forest AGB low-value area from CMB is lower than that of other models, while the uncertainty in high-value areas is larger. The uncertainty range of the GHMB is 10–32.22 t/ha, and that of the 3sGHMB model is 26.09–41.40 t/ha. The uncertainty of the RK-GHMB and RK-3sGHM, which, when considering spatial correlation, ranges from 0.61 to 32.69 t/ha and from 18.06 to 38.63 t/ha, respectively. Their uncertainty is significantly lower than that of GHMB and 3sGHMB, especially in areas near the samples; the uncertainty of RK-3sGHMB is the smallest, but the uncertainty in the low-value area is larger (Figure 7f).

From the population mean estimation results and uncertainty, it can be seen that the precision of the multi-stage method is superior to that of single-stage inference, considering that the spatial correlation of the multi-stage method shows better accuracy than that of the original models. Among the models, the precision of the population mean by RK-3sGHMB is the highest, with the precision of the Upper Yangge and Tidal Slag Forest Farms reaching 93.48% and 95.39%, respectively, with the combined population accuracy of the two farms being 94.44% (Table 6).

4. Discussion

The results of this study indicate that although multi-stage inferences are trained with a larger number of samples in the second or third modeling stage, only the accuracy of GHMB is improved, while the accuracy of 3sGHMB is lower than that of CMB (Figure 5). The primary reason for this may be related to the propagation of errors. Although two-stage inference also experiences error propagation, the features of ALS are strongly correlated with plot AGB. By utilizing ALS features to build a model to plot AGB, to predict the corresponding AGB at other locations in the sampled ALS strips, we thereby expanded the sample size for the second stage of modeling, while keeping the expanded samples more uniformly distributed in space and more representative of the population. The errors introduced by error propagation can be compensated for by the increased number of units in the expanded samples; hence, the model accuracy is still improved, which is consistent with the conclusions of existing research [14]. In the case of three-stage inference that does not consider spatial relationships, the final model is affected by the error propagation from the first two models. Whether this is the cause of the reduced model accuracy requires further investigation.

For population mean estimation, as the AGB sample size used in modeling increases, the representativeness of the sample improves, leading to higher precision and more stable estimation results for the population (Figure 8), which is consistent with existing research [14,17,64,65]. It is important to highlight that all the samples are obtained through random sampling, as shown in Figure 8. For two-stage and three-stage inferences, the change in sample quantity pertains exclusively to the final stage [33]. Concurrently, when estimating models in the other stages, all available samples are utilized to ensure optimal model accuracy [14]. Models that account for spatial correlation exhibit higher accuracy (Table 6, Figure 8), which is closely related to the sample size employed in modeling. Specifically, when the sample size is less than approximately 489, methods that consider spatial relationships demonstrate a comparable accuracy in population mean estimation to those that do not account for spatial correlation. In contrast, when the sample size exceeds approximately 489, the accuracy of models incorporating spatial relationships gradually improves with the increasing sample size. This trend is attributed to the inherent characteristics of geostatistical methods [66]. Under identical sample size conditions, the CMB has a lower population estimation uncertainty than other inferences; this is due to the error propagation in multi-stage models [17]. In this study, the sample size used for the CMB model was only 103, and the curve for CMB model variation and sample size only reached 103. However, Saarela et al. [67] also showed that under the same sample size conditions, the uncertainty of the CMB is less than that of multi-stage inferences. When the sample size is less than approximately 2998, the precision of population estimation results of two-stage inference (GHMB and RK-GHMB) is always better than that of three-stage inferences (3sGHMB and RK-3sGHMB). When the sample size exceeds approximately 2998, as the sample size increases, the estimation results for RK-3sGHMB begin to surpass those of GHMB but remain lower than those of RK-GHMB. When the sample size exceeds a certain value, the RK-3sGHMB model estimation results are optimal among all methods, which is also the reason why the precision of RK-3sGHMB for population estimation is the highest in this study, indicating that geostatistical methods can effectively improve model estimation accuracy. Although the accuracy of the 3sGHMB model is lower than that of CMB, when the sample size exceeds approximately 1000, the population estimation precision is better than that of the CMB model with only 103 samples, but it can never compare with GHMB, RK-GHMB, and RK-3sGHMB.

Based on the aforementioned analysis, we can categorize the application scenarios for inferences in this study as follows. (1) The integration of limited sample plots into comprehensive wall-to-wall data, including, but not limited to, ALS data, optical imagery, and SAR data, as exemplified in Case A, represents a widely adopted framework for AGB estimation. This framework is particularly well suited to CMB. However, the costs associated with field plot surveys are substantial, which can represent a significant barrier to implementation. In an ideal scenario without budgetary constraints, wall-to-wall ALS data would be the optimal choice due to their high accuracy and detailed information. Nonetheless, obtaining full-coverage ALS data over extensive areas is often impractical due to logistical and financial limitations. (2) For smaller study areas (counties, townships, forest farms, etc.), the strategic sampling of ALS data (collected by flying over some strips or blocks in the study area) can reduce the sample size of field plot data required. Utilizing sampled ALS data to estimate target variables can expand the sample size, potentially reducing the cost associated with field plot acquisition and enhancing the representativeness of the samples, particularly for regions that are inaccessible to humans. This framework is amenable to both GHMB and RK-GHMB. It is imperative to note that the RK-GHMB method can only manifest its advantages when there is an adequate number of ALS data to expand the modeling samples in the second stage. Furthermore, if actual AGB data at the satellite waveform LiDAR footprint locations can be procured, then the ALS data used for sample modeling could be substituted with more cost-effective satellite LiDAR data, such as GEDI, ICESat-2, TECIS, and GF-7. However, obtaining actual values for footprint locations in practice remains challenging. (3) In the context of larger study areas (such as at the national, provincial, or municipal level) and with budgetary constraints on ALS data acquisition, employing a limited number of field plots in conjunction with partial sampling ALS data linking the field plots to satellite footprints to augment sample size can improve the uniformity and representativeness of sample distribution. In such scenarios, the RK-3sGHMB exhibits certain advantages.

The results of this study demonstrate that RK-3sGHMB can improve the estimation accuracy of model fitting and the precision of population prediction and can enhance the rationality of the uncertainty estimation results; however, geostatistical methods were only applied in the final stage. Given that the second-stage model also exhibits spatial correlation, future research will explore a three-stage inference estimation approach that accounts for spatial correlation in both the second and third stages. However, the derivation of the uncertainty estimation formula for the second stage using the RK model is too complex and challenging to implement. The spatial filtering estimation method for feature vectors could serve as a primary direction for future research [68,69]. This method has the potential to reduce the model formula complexity and improve the operational efficiency of the model.

5. Conclusions

In the present study, we successfully implemented a three-stage forest AGB estimation method, RK-3sGHMB, by integrating RK into the third stage of 3sGHMB, and conducted a comparative analysis with the CMB, GHMB, RK-GHMB, and 3sGHMB methods, leading to the following conclusions. (1) When the sample data show no spatial correlation or the model estimation does not take spatial correlation into account, multi-stage inference does not demonstrate a distinct advantage, yielding a lower estimation accuracy under equivalent sample size conditions compared with single-stage inference. (2) The integration of geostatistical methods can significantly enhance the model accuracy, which is correlated with the number of samples. As the sample size in the last stage increases, the precision of the estimation results consistently improves. With an adequate number of GEDI samples, the RK-3sGHMB model proposed in this study exhibits superior population estimation precision compared to other models, offering certain advantages in the estimation of forest resource parameters across broader areas.

Author Contributions

Conceptualization, E.C. and X.D.; methodology, X.D. and E.C.; validation, X.D.; formal analysis, X.D.; investigation, X.D., E.C., L.Z. and Y.M.; data curation, X.D., Y.F. and J.W.; writing—original draft preparation, X.D.; writing—review and editing, E.C., L.Z. and X.D.; visualization, X.D.; funding acquisition, E.C. All authors contributed to interpreting the results and to the improvement of the article. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The field data, ALS data, and P-SAR data are not publicly available to protect the privacy of private landowners.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Location and coverage of the study area.

Figure 2. (a) Spatial distribution of sample plots, data from Sentinel-2A, P-SAR, and GEDI in the study area; and (b) the spatial distribution of ALS strips and GEDI data within the strips.

Figure 3. Overall workflow of the study; the green square blocks represent field sample plots; the purple square blocks represent ALS sample plots; the yellow circles represent the GEDI sample plots; the black bars represent the ALS strips; the elements enclosed within the blue boxes collectively constitute a unified whole.

Figure 4. Variation in the feature selection iteration curves using genetic algorithms according to [Forumla omitted. See PDF.] (a), [Forumla omitted. See PDF.] (b), [Forumla omitted. See PDF.] (c), [Forumla omitted. See PDF.] (d), and [Forumla omitted. See PDF.] (e).

Figure 5. Model accuracy evaluation of [Forumla omitted. See PDF.] (a), [Forumla omitted. See PDF.] (b), MB (c), GHMB (d), RK-GMB (e), 3sGHMB (f), and RK-3sGHMB (g) using independent samples.

Figure 6. Residual variogram and fitted model with exponential, gaussian, spherical, and linear fitting using RK-GHMB (a) and RK-3sGHMB (b).

Figure 7. The results for forest AGB and uncertainty: the spatial distribution of forest AGB and uncertainty estimated by 3sGHMB (a), GHMB (b), RK-3sGHMB (c), RK-GHMB (d), and CMB (e); (f) pixel-level uncertainty (RMSE) with respect to predicted AGB.

Figure 8. The influence of sample size used for [Forumla omitted. See PDF.] estimation using MB, GHMB, RH-GHMBa, 3sGHMB, and RK-3sGHMB.

Table 1

Features of LiDAR.

Feature Name	Feature Symbol	Description
Mean height	Hmean	The mean height of the point cloud within a 25 m × 25 m statistical unit
Forest canopy density	CD	The ratio of canopy backscatter points within a 25 m × 25 m statistical unit to the total number of backscatter points
Variance and standard deviation	Hvar, Hstd	The variance and standard deviation of the point cloud within a 25 m × 25 m statistical unit
Coefficient of variation	Hcv	Coefficient of variation in the point cloud within a 25 m × 25 m statistical unit
Interquartile distance	Hint	Height interquartile distance within a 25 m × 25 m statistical unit for point cloud data
Skewness	Hskw	The skewness of the points cloud within a 25 m × 25 m statistical unit
Maximum and minimum	Hmax, Hmin	The maximum and minimum values of the point cloud within a 25 m × 25 m statistical unit
Percentiles	H10, H20, H30, H40, H50, H60, …, H95	Percentiles of the point cloud at different heights within a 25 m × 25 m statistical unit

Table 2

Polarimetric features extracted from P-SAR data in this study.

Feature Class	Feature Name	Feature Symbol
Backscatter intensity	HH polarization	PHH
	HV polarization	PHV
Yamaguchi decomposition	Surface scattering	Odd
	Dihedral scattering	Dbl
	Volume scattering	Vol
	Helix scattering	Hlx

Table 3

Features of Sentinel-2A.

Feature	Symbol	Formula	Reference
Spectral reflectance	B2, B3, B4, B5, B6, B7, B8a, B11, B12	/	/
Normalized difference vegetation index	NDVI	$\left(B_8-B_4\right)/\left(B_8+B_4\right)$	[52]
Red-edge vegetation index	NDVIre1	$\left(B_8-B_5\right)/\left(B_8+B_5\right)$	[53]
	NDVIre2	$\left(B_8-B_6\right)/\left(B_8+B_6\right)$
	NDVIre3	$\left(B_8-B_7\right)/\left(B_8+B_7\right)$
Kernel-normalized difference vegetation index	KNDVI	$tanh\left({\left({\displaystyle \frac{B_8-B_4}{2{\sigma }_1}}\right)}^2\right)$	[51]
Kernel red-edge vegetation index	KNDVIre1	$tanh\left({\left({\displaystyle \frac{B_8-B_5}{2{\sigma }_2}}\right)}^2\right)$
	KNDVIre2	$tanh\left({\left({\displaystyle \frac{B_8-B_6}{2{\sigma }_3}}\right)}^2\right)$
	KNDVIre3	$tanh\left({\left({\displaystyle \frac{B_8-B_7}{2{\sigma }_4}}\right)}^2\right)$
Difference vegetation index	DVI	$B_8-B_4$	[54]
Red-edge difference vegetation index	DVIre1	$B_5-B_4$	[55]
	DVIre2	$B_6-B_4$
	DVIre3	$B_7-B_4$
Enhanced vegetation index	EVI	$2.5\times \left(B_8-B_4\right)/\left(B_8+6\times B_4-7.5\times B_2+1\right)$	[56]
Red-edge enhanced vegetation index	EVIre1	$2.5\times \left(B_5-B_4\right)/\left(B_5+6\times B_4-7.5\times B_2+1\right)$	[55]
	EVIre2	$2.5\times \left(B_6-B_4\right)/\left(B_6+6\times B_4-7.5\times B_2+1\right)$
	EVIre3	$2.5\times \left(B_7-B_4\right)/\left(B_7+6\times B_4-7.5\times B_2+1\right)$
Simple ratio	RVI	$B_8/B_4$	[57]
Soil-adjusted vegetation index	SAVI	$1.5\times \left(B_8-B_4\right)/\left(B_8+B_4+0.5\right)$	[58]
Green normalized difference vegetation index	GNDVI	$\left(B_8-B_3\right)/\left(B_8+B_3\right)$	[59]
Narrow enhanced vegetation index	EVInirn	$2.5\times \left(B_{8A}-B_4\right)/\left(B_{8A}+6\times B_4-7.5\times B_2+1\right)$	[55]

$B_2$, $B_3$, $B_4$, $B_5$, $B_6$, $B_7$, $B_8$, $B_{8a}$, $B_{11}$, and $B_{12}$ represent the bands corresponding to sentinel-2A; the value of ${\sigma }_1$, ${\sigma }_2$, ${\sigma }_3$ and ${\sigma }_4$ is $\left(B_8+B_4\right)/2$, $\left(B_8+B_5\right)/2$, $\left(B_8+B_6\right)/2$, and $\left(B_8+B_7\right)/2$, respectively.

Table 4

Fitting accuracy of the semivariogram for RK-3sGHMB and RK-GHMB.

Residual Source	Fitting Method	$\mathbf{Nugget} ({\mathit{C}}_0/{\mathit{C}}_0^{*})$	$\mathbf{Partial} \mathbf{Sill} ({\mathit{C}}_1/{\mathit{C}}_1^{*})$	Range (m)	Ratio (%)	${\mathit{R}}_{\mathit{a}\mathit{d}\mathit{j}}^2$	RMSE (t/ha)	EA (%)
RK-3sGHMB	Exponential	328.35	525.55	16,069.84	38.45	0.55	28.55	69.13
	Spherical	414.00	420.82	7109.70	49.59	0.49	30.67	66.83
	Gaussian	516.99	302.11	10,615.19	63.12	0.48	30.70	66.79
	Linear	747.71	747.71	/	/	0.38	38.22	58.7
RK-GHMB	Exponential	120.02	539.06	2625.28	18.21	0.60	27.07	70.72
	Spherical	338.20	533.07	2019.45	38.82	0.53	29.21	68.41
	Gaussian	380.09	533.08	2288.70	41.62	0.41	32.84	64.48
	Linear	610.96	610.96	/	/	0.41	32.85	64.47

Table 5

Fitted forms for models $F_{ms}$, $F_{ml}$, $G_{ls}$, $Q_{lg}$, and $Q_{gs}$.

Model	Model Formula	${\mathit{R}}_{\mathit{a}\mathit{d}\mathit{j}}^2$	RMSE (t/ha)	EA(%)
Fms	$\begin{array}{l}201.09464/(1+exp(-0.70266\times EVI-0.09913\times PHV+\\ 21.18109\times DVIre1-0.14215\times EVIre1-9.69895))\end{array}$	0.37	33.95	63.28
Fml	$\begin{array}{l}299.30254/(1+exp(-0.074878\times Hvar+0.026074\times H98\\ -0.006368\times \mathit{C}D+1.947972))\end{array}$	0.80	19.20	79.24
Gls	$\begin{array}{l}189.4068/(1+exp(-29.7195\times EVI-0.1698 \times PHV+30.0420\times \\ KNDVIre1+41.6752\times Vol+18.6999\times EVIre1+17.9654))\end{array}$	0.38	33.72	63.53
Qlg	$\begin{array}{l}288.317726/(1+exp(-0.094833\times rh90-0.046170\times rh50+\\ 0.059449\times rh10+0.003274\times \mathit{C}over+3.446146))\end{array}$	0.81	15.84	83.78
Qgs	$\begin{array}{l}199.54404/(1+\exp (20.57206\times Hlx-9.80876\times EVIre3+\\ 7.25995\times KNDVIre1+4.68176\times EVIre1-0.18229\times PHV\\ -0.08263\times PHH+0.02713\times Vol+5.37781))\end{array}$	0.27	36.58	60.43

Table 6

The population results and uncertainties for regions.

Region	Methods	Population Mean (t/ha)	Variance in the Population Mean (t/ha)	SE (t/ha)	p (%)
Upper Yangge Qi Forest Farm	CMB	80.43	30.27	5.50	86.60
	GHMB	83.63	19.43	4.41	89.67
	RK-GHMB	81.50	9.99	3.16	92.40
	3sGHMB	84.20	29.33	5.42	87.39
	RK-3sGHMB	79.16	6.94	2.63	93.48
Tidal Slag Forest Farm	CMB	98.66	28.98	5.38	89.31
	GHMB	102.75	17.20	4.15	92.09
	RK-GHMB	100.79	8.77	2.96	94.24
	3sGHMB	101.07	24.35	4.93	90.43
	RK-3sGHMB	103.72	5.96	2.44	95.39
Total	CMB	89.55	29.63	5.44	87.96
	GHMB	93.19	18.32	4.28	90.88
	RK-GHMB	91.15	9.38	3.06	93.32
	3sGHMB	92.64	26.84	5.18	88.91
	RK-3sGHMB	91.44	6.45	2.54	94.44

References

1. Dixon, R.K.; Brown, S.; Houghton, R.A.; Solomon, A.M.; Trexler, M.C.; Wisniewski, J. Carbon pools and flux of global forest Ecosystems. Science; 1994; 263, pp. 185-190. [DOI: https://dx.doi.org/10.1126/science.263.5144.185] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17839174]

2. Hese, S.; Lucht, W.; Schmullius, C.; Barnsley, M.; Dubayah, R.; Knorr, D.; Neumann, K.; Riedel, T.; Schröter, K. Global biomass mapping for an improved understanding of the CO₂ balance-the earth observation mission carbon-3D. Remote Sens. Environ.; 2005; 94, pp. 94-104. [DOI: https://dx.doi.org/10.1016/j.rse.2004.09.006]

3. Chen, Q. Modeling aboveground tree woody biomass using national-scale allometric methods and airborne lidar. ISPRS J. Photogramm. Remote Sens.; 2015; 106, pp. 95-106. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2015.05.007]

4. Askne, J.I.H.; Soja, M.J.; Ulander, L.M.H. Biomass estimation in a boreal forest from TanDEM-X data, lidar DTM, and the interferometric water cloud model. Remote Sens. Environ.; 2017; 196, pp. 265-278. [DOI: https://dx.doi.org/10.1016/j.rse.2017.05.010]

5. Ståhl, G.; Saarela, S.; Schnell, S.; Holm, S.; Breidenbach, J.; Healey, S.P.; Patterson, P.L.; Magnussen, S.; Næsset, E.; McRoberts, R.E. et al. Use of models in large-area forest surveys: Comparing model-assisted, model-based and hybrid estimation. For. Ecosyst.; 2016; 3, 5. [DOI: https://dx.doi.org/10.1186/s40663-016-0064-9]

6. McRoberts, R.E.; Nasset, E.; Gobakken, T.; Chirici, G.; Condés, S.; Hou, Z.; Saarela, S.; Chen, Q.; Ståhl, G.; Walters, B.F. Assessing components of the model-based mean square error estimator for remote sensing assisted forest applications. Can. J. For. Res.; 2018; 48, pp. 642-649. [DOI: https://dx.doi.org/10.1139/cjfr-2017-0396]

7. Marshak, C.; Simard, M.; Duncanson, L.; Silva, C.A.; Denbina, M.; Liao, T.-H.; Fatoyinbo, L.; Moussavou, G.; Armston, J. Regional tropical aboveground biomass mapping with L-Band Repeat-Pass Interferometric Radar, Sparse lidar, and multiscale superpixels. Remote Sens.; 2020; 12, 2048. [DOI: https://dx.doi.org/10.3390/rs12122048]

8. Saarela, S.; Varvia, P.; Korhonen, L.; Yang, Z.; Patterson, P.L.; Gobakken, T.; Næsset, E.; Healey, S.P.; Ståhl, G. Tree-phase hierarchical model-based and hybrid inference. MethodsX; 2023; 11, 102321. [DOI: https://dx.doi.org/10.1016/j.mex.2023.102321]

9. Saarela, S.; Wästlund, A.; Holmström, E.; Mensah, A.A.; Holm, S.; Nilsson, M.; Fridman, J.; Ståhl, G. Mapping aboveground biomass and its prediction uncertainty using LiDAR and field data, accounting for tree-level allometric and LiDAR model errors. For. Ecosyst.; 2020; 7, 43. [DOI: https://dx.doi.org/10.1186/s40663-020-00245-0]

10. Magnussen, S.; Nord-Larsen, T.; Nielsen, T.R. Lidar supported estimators of wood volume and aboveground biomass from the Danish national forest inventory (2012–2016). Remote Sens. Environ.; 2018; 211, pp. 146-153. [DOI: https://dx.doi.org/10.1016/j.rse.2018.04.015]

11. Silva, C.A.; Hudak, A.T.; Vierling, L.A.; Loudermilk, E.L.; O’Brien, J.J.; Hiers, J.K.; Jack, S.B.; Gonzalez-Benecke, C.; Lee, H.; Falkowski, M.J. et al. Imputation of individual longleaf pine (pinus palustris mill.) tree attributes from field and LiDAR data. Can. J. Remote Sens; 2016; 42, pp. 554-573. [DOI: https://dx.doi.org/10.1080/07038992.2016.1196582]

12. Xu, C.; Manley, B.; Morgenroth, J. Evaluation of modelling approaches in predicting forest volume and stand age for small-scale plantation forests in New Zealand with RapidEye and LiDAR. Int. J. Appl. Earth Obs. Geoinf.; 2018; 73, pp. 386-396. [DOI: https://dx.doi.org/10.1016/j.jag.2018.06.021]

13. Lim, K.; Treitz, P.; Wulder, M.; St-Onge, B.; Flood, M. LiDAR remote sensing of forest structure. Prog. Phys. Geogr.; 2003; 27, pp. 88-106. [DOI: https://dx.doi.org/10.1191/0309133303pp360ra]

14. Xie, B.; Cao, C.; Xu, M.; Bashir, B.; Singh, R.P.; Huang, Z.; Lin, X. Regional forest volume estimation by expanding LiDAR samples using multi-sensor satellite data. Remote Sens.; 2020; 12, 360. [DOI: https://dx.doi.org/10.3390/rs12030360]

15. McRoberts, R.E.; Næsset, E.; Gobakken, T. Inference for lidar-assisted estimation of forest growing stock volume. Remote Sens. Environ.; 2013; 128, pp. 268-275. [DOI: https://dx.doi.org/10.1016/j.rse.2012.10.007]

16. McRoberts, R.E.; Næsset, E.; Saatchi, S.; Quegan, S. Statistically rigorous, model-based inferences from map. Remote Sens. Environ.; 2022; 279, 113028. [DOI: https://dx.doi.org/10.1016/j.rse.2022.113028]

17. Saarela, S.; Holm, S.; Grafström, A.; Schnell, S.; Næsset, E.; Gregoire, T.G.; Nelson, R.F.; Ståhl, G. Hierarchical model-based inference for forest inventory utilizing three sources of information. Ann. For. Sci.; 2016; 73, pp. 895-910. [DOI: https://dx.doi.org/10.1007/s13595-016-0590-1]

18. Saarela, S.; Holm, S.; Healey, S.P.; Patterson, P.L.; Yang, Z.; Andersen, H.E.; Dubayah, R.O.; Qi, W.; Duncanson, L.I.; Armston, J.D. et al. Comparing frameworks for biomass prediction for the Global Ecosystem Dynamics Investigation. Remote Sens. Environ.; 2022; 278, 113074. [DOI: https://dx.doi.org/10.1016/j.rse.2022.113074]

19. Holm, S.; Nelson, R.; Ståhl, G. Hybrid three-phase estimators for large-area forest inventory using ground plots, airborne lidar, and space lidar. Remote Sens. Environ.; 2017; 197, pp. 85-97. [DOI: https://dx.doi.org/10.1016/j.rse.2017.04.004]

20. Pang, Y.; Yu, T.; Jia, W.; Liang, X.; Li, Z.; Fu, A.; Wu, F.; Liu, X.; Zhang, X.; Huang, J. et al. TECIS: The first mission towards forest carbon mapping by combination of lidar and multi-angle optical observations. IGARSS 2023-2023 IEEE International Geoscience and Remote Sensing Symposium; IEEE: Pasadena, CA, USA, 2023; pp. 1890-1893. [DOI: https://dx.doi.org/10.1109/IGARSS52108.2023.10282849]

21. Varvia, P.; Saarela, S.; Maltamo, M.; Packalen, P.; Gobakken, T.; Næsset, E.; Ståhl, G.; Korhonen, L. Estimation of boreal forest biomass from ICESat-2 data using hierarchical hybrid inference. Remote Sens. Environ.; 2024; 311, 114249. [DOI: https://dx.doi.org/10.1016/j.rse.2024.114249]

22. Wulder, M.A.; White, J.C.; Bater, C.W.; Coops, N.C.; Hopkinson, C.; Chen, G. Lidar plots-a new large-area data collection option: Context, concepts, and case study. Can. J. Remote Sens.; 2012; 38, pp. 600-618. [DOI: https://dx.doi.org/10.5589/m12-049]

23. Margolis, H.A.; Nelson, R.F.; Montesano, P.M.; Beaudoin, A.; Sun, G.; Andersen, H.E.; Wulder, M.A. Combining satellite lidar, airborne lidar, and ground plots to estimate the amount and distribution of aboveground biomass in the boreal forest of North America. Can. J. For. Res.; 2015; 45, pp. 838-855. [DOI: https://dx.doi.org/10.1139/cjfr-2015-0006]

24. Narine, L.L.; Popescu, S.C.; Malambo, L. Using ICESat-2 to estimate and map forest aboveground biomass: A first example. Remote Sens.; 2020; 12, 1824. [DOI: https://dx.doi.org/10.3390/rs12111824]

25. Varvia, P.; Korhonen, L.; Bruguière, A.; Toivonen, J.; Packalen, P.; Maltamo, M.; Saarela, S.; Popescu, S.C. How to consider the effects of time of day, beam strength, and snow cover in ICESat-2 based estimation of boreal forest biomass?. Remote Sens. Environ.; 2022; 280, 113174. [DOI: https://dx.doi.org/10.1016/j.rse.2022.113174]

26. Guerra-Hernández, J.; Narine, L.L.; Pascual, A.; Gonzalez-Ferreiro, E.; Botequim, B.; Malambo, L.; Neuenschwander, A.; Popescu, S.C.; Godinho, S. Aboveground biomass mapping by integrating ICESat-2, SENTINEL-1, SENTINEL-2, ALOS2/PALSAR2, and topographic information in mediterranean forests. GISci. Remote Sens.; 2022; 59, pp. 1509-1533. [DOI: https://dx.doi.org/10.1080/15481603.2022.2115599]

27. Zakeri, F.; Mariethoz, G. A review of geostatistical simulation models applied to satellite remote sensing: Methods and applications. Remote Sens. Environ.; 2021; 259, 112381. [DOI: https://dx.doi.org/10.1016/j.rse.2021.112381]

28. Fayad, I.; Baghdadi, N.; Bailly, J.S.; Barbier, N.; Gond, V.; Herault, B.; El Hajj, M.; Fabre, F.; Perrin, J. Regional scale rain-forest height mapping using regression-kriging of spaceborne and airborne LiDAR data: Application on French Guiana. Remote Sens.; 2016; 8, 240. [DOI: https://dx.doi.org/10.3390/rs8030240]

29. Pouladi, N.; Møller, A.B.; Tabatabai, S.; Greve, M.H. Mapping soil organic matter contents at field level with cubist, random forest and kriging. Geoderma; 2019; 342, pp. 85-92. [DOI: https://dx.doi.org/10.1016/j.geoderma.2019.02.019]

30. Silveira, E.M.; Santo, F.D.E.; Wulder, M.A.; Júnior, F.W.A.; Carvalho, M.C.; Mello, C.R.; Mello, J.M.; Shimabukuro, Y.E.; Terra, M.C.N.S.; Carvalho, L.M.T. et al. Pre-stratified modelling plus residuals kriging reduces the uncertainty of aboveground biomass estimation and spatial distribution in heterogeneous savannas and forest environments. For. Ecol. Manag.; 2019; 445, pp. 96-109. [DOI: https://dx.doi.org/10.1016/j.foreco.2019.05.016]

31. Chen, L.; Ren, C.; Zhang, B.; Wang, Z. Multi-sensor prediction of stand volume by a hybrid model of support vector machine for regression kriging. Forests; 2020; 11, 296. [DOI: https://dx.doi.org/10.3390/f11030296]

32. Zhao, J.; Zhao, L.; Chen, E.; Li, Z.; Xu, K.; Ding, X. An improved generalized hierarchical estimation framework with geostatistics for mapping forest parameters and its uncertainty: A case study of forest canopy height. Remote Sens.; 2022; 14, 568. [DOI: https://dx.doi.org/10.3390/rs14030568]

33. Chen, F.; Hou, Z.; Saarela, S.; McRoberts, R.E.; Ståhl, G.; Kangas, A.; Packalen, P.; Li, B.; Xu, Q. Leveraging remotely sensed non-wall-to-wall data for wall-to-wall upscaling in forest inventory. Int J Appl Earth Obs Geoinf.; 2023; 119, 103314. [DOI: https://dx.doi.org/10.1016/j.jag.2023.103314]

34. Li, C.; Zhang, W.; Li, Z.; Chen, E.; Tian, X. Retrieval of forest above-ground biomass using multi-source data in Genhe, Inner Mongolia. J. Beijing For. Univ. (Chin. Ed.); 2016; 38, pp. 64-72. [DOI: https://dx.doi.org/10.13332/j.1000-1522.20150209]

35. Zhou, G.; Yin, G.; Tang, X. Biomass Equation for Forest Ecosystems in China-Carbon Storage; China Science Publishing: Beijing, China, 2018; pp. 41-80.

36. Pang, Y.; Li, Z.; Ju, H.; Lu, H.; Jia, W.; Si, L.; Guo, Y.; Liu, Q.; Li, S.; Liu, L. et al. LiCHy: The CAF’s LiDAR, CCD and hyperspectral integrated airborne observation system. Remote Sens.; 2016; 8, 398. [DOI: https://dx.doi.org/10.3390/rs8050398]

37. Isenburg, M. LAStools—Efficient Tools for LiDAR Processing. Version 2.0.1 2022; Available online: https://rapidlasso.de/lastools-220107 (accessed on 25 March 2022).

38. McGaughey, R.J. FUSION/LDV: Software for LiDAR Data Analysis and Visualization; January 2021-FUSION Version 4.20 United stated of Department of Agriculture: Washington, DC, USA, 2021; Available online: http://forsys.cfr.washington.edu/software/fusion/FUSION_manual.pdf (accessed on 16 June 2021).

39. Cao, L.; Xu, T.; Shen, X.; She, G. Mapping biomass by integrating Landsat OLI and airborne LiDAR transect data in subtropical forests. J. Remote Sens.; 2016; 20, pp. 665-678. [DOI: https://dx.doi.org/10.11834/jrs.20165210]

40. Donoghue, D.N.M.; Watt, P.J.; Cox, N.J.; Wilson, J. Remote sensing of species mixtures in conifer plantations using LiDAR height and intensity data. Remote Sens. Environ.; 2007; 110, pp. 509-522. [DOI: https://dx.doi.org/10.1016/j.rse.2007.02.032]

41. Fan, Y.; Zhao, L.; Chen, E.; Xu, K.; Zhang, W.; Ma, Y. Evaluation of forest stock estimation ability of high resolution airborne multi-band PolSAR in cold temperate coniferous forests. Natl. Remote Sens. Bull.; 2024; 28, pp. 2525-2539. [DOI: https://dx.doi.org/10.11834/jrs.20243200]

42. Zhao, L.; Chen,; Li, Z.; Zhang, W.; Gu, X. Three-step semi-empirical radiometric terrain correction approach for PolSAR data applied to forested areas. Remote Sens.; 2017; 9, 269. [DOI: https://dx.doi.org/10.3390/rs9030269]

43. Yamaguchi, Y.; Moriyama, T.; Ishido, M.; Yamada, H. Four-component scattering model for polarimetric SAR image decomposition. IEEE T. Geosci. Remote.; 2005; 43, pp. 1699-1706. [DOI: https://dx.doi.org/10.1109/TGRS.2005.852084]

44. Spracklen, B.; Spracklen, D.V. Detenllination of structuml chacteristics of old-growth forest in Ukmine using spaceborne LiDAR. Remote Sens.; 2021; 13, 1233. [DOI: https://dx.doi.org/10.3390/rs13071233]

45. Lin, X.; Xu, M.; Cao, C.; Dang, Y.; Bashir, B.; Xie, B.; Huang, Z. Estimates of forest canopy height using a combination of ICESat-2/ATLAS data and stereo-photogmmmetry. Remote Sens.; 2020; 12, 3649. [DOI: https://dx.doi.org/10.3390/rs12213649]

46. Liu, L.; Wang, C.; Nie, S.; Zhu, X.; Xi, X.; Wang, J. Analysis of the influence of different algorithms of GEDI L2A on the accuracy of ground elevation and forest canopy height. J. Univ. Chin. Acad. Sci.; 2022; 39, pp. 502-511. [DOI: https://dx.doi.org/10.7523/j.ucas.2021.0076]

47. Dorado-Ronda, I.; Pascual, A.; Godinho, S.; Silva, C.A.; Botequim, B.; Ro-dríguez-Gonzálvez, P.; González-Ferreiro, E.; Guerra-Hernández, J. Assessing the accuracy of GEDI data for canopy height and aboveground biomass estimates in Mediterranean forests. Remote Sens.; 2021; 13, 2279. [DOI: https://dx.doi.org/10.3390/rs13122279]

48. Qi, W.; Saarela, S.; Armston, J.; Ståhl, G.; Dubayah, R. Forest biomass estimation over three distinct forest types using TanDEMX InSAR data and simulated GEDI lidar data. Remote Sens. Environ.; 2019; 232, 111283. [DOI: https://dx.doi.org/10.1016/j.rse.2019.111283]

49. Dubayah, R.; Blair, J.B.; Goetz, S.; Fatoyinbo, L.; Hansen, M.; Healey, S.; Hofton, M.; Hurtt, G.; Kellner, J.; Luthcke, S. et al. The Global Ecosystem Dynamics Investigation: High-resolution laser ranging of the Earth’s forests and topography. Sci. Remote Sens.; 2020; 1, 100002. [DOI: https://dx.doi.org/10.1016/j.srs.2020.100002]

50. Silva, C.A.; Hamamura, C.; Valbuena, R.; Hancock, S.; Cardil, A.; Broadbent, E.N.; Almeida, D.R.A.; Silva, J.; Klauberg, C. rGEDI: NASA’s Global Ecosystem Dynamics Investigation (GEDI) Data Visualization and Processing. Version 0.1.9. Available online: https://CRAN.R-project.org/package=rGEDI (accessed on 22 October 2020).

51. Camps-Valls, G.; Campos-Taberner, M.; Moreno-Martínez, Á.; Walther, S.; Duveiller, G.; Cescatti, A.; Mahecha, M.D.; Muñoz-Marí, J.; García-Haro, F.J.; Guanter, L. A unified vegetation index for quantifying the terrestrial biosphere. Sci. Adv.; 2021; 7, eabc7447. [DOI: https://dx.doi.org/10.1126/sciadv.abc7447]

52. Rouse, J.; Haas, R.H.; Deering, D.; Schell, J.; Harlan, J. Monitoring the Vernal Advancement and Retrogradation (Green Wave Effect) of Natural Vegetation; Technical Report NASA/GSFC Type III Final Report NASA/GSFC: Greenbelt, MD, USA, 1974; Corpus ID: 129198382

53. Shoko, C.; Mutanga, O. Examining the strength of the newly-launched Sentinel 2 MSI sensor in detecting and discriminating subtle differences between C3 and C4 grass species. ISPRS J. Photogramm. Remote Sens.; 2017; 129, pp. 32-40. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2017.04.016]

54. Richardson, A.J.; Wiegand, C.L. Distinguishing vegetation from soil background information. Photogramm. Eng. Remote Sens.; 1977; 43, pp. 1541-1552. Corpus ID: 126604551

55. Zhang, X.; Shen, H.; Huang, T.; Wu, Y.; Guo, B.; Liu, Z.; Luo, H.; Tang, J.; Zhou, H.; Wang, L. et al. Improved random forest algorithms for increasing the accuracy of forest aboveground biomass estimation using Sentinel-2 imagery. Ecol. Indic.; 2024; 159, 111752. [DOI: https://dx.doi.org/10.1016/j.ecolind.2024.111752]

56. Liu, H.Q.; Huete, A. A feedback based modification of the NDVI to minimize canopy background and atmospheric noise. IEEE T. Geosci. Remote.; 1995; 33, pp. 457-465. [DOI: https://dx.doi.org/10.1109/TGRS.1995.8746027]

57. Jordan, C.F. Derivation of leaf-area index from quality of light on the forest floor. Ecology; 1969; 5, pp. 663-666. [DOI: https://dx.doi.org/10.2307/1936256]

58. Huete, A.R. A soil-adjusted vegetation index (SAVI). Remote Sens. Environ.; 1988; 25, pp. 295-309. [DOI: https://dx.doi.org/10.1016/0034-4257(88)90106-X]

59. Gitelson, A.A.; Kaufman, Y.J.; Merzlyak, M.N. Use of a green channel in remote sensing of global vegetation from EOS-MODIS. Remote Sens. Environ.; 1996; 58, pp. 289-298. [DOI: https://dx.doi.org/10.1016/S0034-4257(96)00072-7]

60. McRoberts, R.E.; Domke, G.M.; Chen, Q.; Næsset, E.; Gobakken, T. Using genetic algorithms to optimize k-Nearest Neighbors configurations for use with airborne laser scanning data. Remote Sens. Environ.; 2016; 184, pp. 387-395. [DOI: https://dx.doi.org/10.1016/j.rse.2016.07.007]

61. Kuhn, M. Building predictive models in R using the caret package. J. Stat. Softw.; 2008; 28, pp. 1-26. [DOI: https://dx.doi.org/10.18637/jss.v028.i05]

62. McRoberts, R.E.; Næsset, E.; Hou, Z.; Ståhl, G.; Saarela, S.; Esteban, J.; Travaglini, D.; Mohammadi, J.; Chirici, G. How many bootstrap replications are necessary for estimating remote sensing-assisted, model-based standard errors?. Remote Sens. Environ.; 2023; 288, 113455. [DOI: https://dx.doi.org/10.1016/j.rse.2023.113455]

63. Tao, S.; Labrière, N.; Calders, K.; Fischer, F.J.; Rau, E.; Plaisance, L.; Chave, J. Mapping tropical forest trees across large areas with lightweight cost-efective terrestrial laser scanning. Ann. For. Sci.; 2021; 78, 103. [DOI: https://dx.doi.org/10.1007/s13595-021-01113-9]

64. Saarela, S.; Schnell, S.; Grafström, A.; Tuominen, S.; Nordkvist, K.; Hyyppä, J.; Kangas, A.; Ståhl, G. Effects of sample size and model form on the accuracy of model-based estimators of growing stock volume. Can. J. For. Res.; 2015; 45, pp. 1524-1534. [DOI: https://dx.doi.org/10.1139/cjfr-2015-0077]

65. Li, C.; Yu, Z.; Dai, H.; Zhou, X.; Zhou, M. Effect of sample size on the estimation of forest inventory attributes using airborne LiDAR data in large-scale subtropical areas. Ann. For. Sci.; 2023; 80, 40. [DOI: https://dx.doi.org/10.1186/s13595-023-01209-4]

66. Voss, S.; Zimmermann, B.; Zimmermann, A. Detecting spatial structures in throughfall data: The effect of extent, sample size, sampling design, and variogram estimation method. J. Hydrol.; 2016; 540, pp. 527-537. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2016.06.042]

67. Saarela, S.; Holm, S.; Healey, S.P.; Andersen, H.E.; Petersson, H.; Prentius, W.; Patterson, P.L.; Næsset, E.; Gregoire, T.G.; Ståhl, G. Generalized hierarchical model-based estimation for aboveground biomass assessment using GEDI and landsat data. Remote Sens.; 2018; 10, 1832. [DOI: https://dx.doi.org/10.3390/rs10111832]

68. Griffith, D.A. Spatial autucorrelation and eigenfunctions of the geographic weights matrix accompanying geo-referenced data. Can. Geogr.; 1996; 40, pp. 351-367. [DOI: https://dx.doi.org/10.1111/j.1541-0064.1996.tb00462.x]

69. Getis, A.; Griffith, D.A. Comparative spatial filtering in regression analysis. Geogr. Anal.; 2002; 34, pp. 130-140. [DOI: https://dx.doi.org/10.1111/j.1538-4632.2002.tb01080.x]