Improving Plot-Level Model of Forest Biomass: A Combined Approach Using Machine Learning with Spatial Statistics

: Estimating the aboveground biomass (AGB) at the plot level plays a major role in connecting accurate single-tree AGB measurements to relatively difﬁcult regional AGB estimates. However, AGB estimates at the plot level suffer from many uncertainties. The goal of this study is to determine whether combining machine learning with spatial statistics reduces the uncertainty of plot-level AGB estimates. To illustrate this issue, this study evaluates and compares the performance of different models for estimating plot-level forest AGB. These models include three different machine learning models [support vector machine (SVM), random forest (RF), and a radial basis function artiﬁcial neural network (RBF-ANN)], one spatial statistic model (P-BSHADE), and three combinations thereof (SVM & P-BSHADE, RF & P-BSHADE, and RBF-ANN & P-BSHADE). The results show that the root mean square error, mean absolute error, and mean relative error of all combined models are substantially smaller than those of any individual model, with the RF & P-BSHADE combined method generating the smallest values. These results indicate that a combined approach using machine learning with spatial statistics, especially the RF & P-BSHADE model, improves the accuracy of plot-level AGB models. These research results contribute to the development of accurate large-forested-landscape AGB maps.


Introduction
Estimates of forest plot-level aboveground biomass (AGB) serve to connect single-tree AGB measurements to regional-scale AGB maps. Uncertainties in plot-level estimates can ultimately propagate to regional AGB maps, degrading the quality and credibility of decision making in sustainable forest management [1,2]. Uncertainties in plot-level estimates stem from many sources, such as the non-local allometric equations used to calculate tree-level AGB [3] or insufficiently robust plot-level AGB-estimation models. Thus, improving the plot-level model of forest AGB is a key issue for producing accurate AGB maps.
More recently, a variety of prediction models have been applied to make accurate AGB estimates, including linear models [4], nonlinear machine learning models [5], spatial statistical models [6][7][8], and hierarchical spatial Bayesian models [9,10], regardless of the data source. Investigators have compared the performance of different models for estimating forest biomass and volume. Some report that the random forest and regression models predict similar forest volumes [11], and others report that machine learning outperforms linear regression models for predicting forest attributes [12]. In addition, a spatial model (geographically weighted regression) produced more accurate estimates than a nonspatial multiple regression model [8]. The performance of different models depends on the forest type and data source.
These different models have their own advantages and disadvantages. Traditional parametric models have difficulty characterizing nonlinear relationships between AGB and multiple environmental covariates. In comparison, nonparametric models are advantageous because they are more elastic and are good at fitting different functional forms without prior knowledge [13]. Moreover, nonlinear nonparametric models are advantageous for dealing with nonlinear fitting and have significant potential for applications involving nonlinear systems, such as forest systems. Some nonparametric machine learning algorithms (e.g., K-nearest neighbor, nonlinear support vector machine, and random forest) offer high prediction accuracy for estimating forest AGB [5,14]. Nevertheless, these nonparametric models also have shortcomings, such as overfitting and poor performance beyond the range of training data.
Spatial statistical models form another group of models frequently used to estimate forest AGB [8,15,16]. For example, kriging was used to construct a pantropical forest carbon stock map [17], and geographically weighted regression was used to predict forest AGB integrated with dominant variables that explained the spatial variation in AGB [8]. Spatial statistical approaches consider the possible spatial autocorrelation and heterogeneity of forest structure. Compared with traditional statistical methods, spatial methods integrate spatial information that affects the model response, thus overcoming the constraints of traditional statistical methods that assume sample independence [7,18] and improving our understanding of spatial autocorrelation [19]. However, kriging or geographically weighted regression methods have well-known disadvantages, such as not considering the uncertainty of spatial covariance parameters derived from variograms [9].
The accuracy of current AGB estimates still suffers from significant uncertainty, regardless of which models are applied [20][21][22]. From the model perspective, other reasons in addition to the shortcomings mentioned above may cause this uncertainty: First, the relationship between forest AGB and covariates may not be fully described [9]. Second, the data distribution of forest AGB is not consistent with the assumption of an independent and identical distribution, which is required by the traditional spatial statistical model when it has spatial autocorrelation and spatial heterogeneity. In addition, nonparametric models often do not adequately address nested or hierarchical observations. An interesting question that thus arises is whether integrating different models (i.e., machine learning and spatial statistics) improves the accuracy of AGB estimates.
Different approaches complement the advantages of different models and may yield more accurate AGB estimates than would otherwise be obtained by using a single method. The objective of the present study is to develop and evaluate a combined approach that combines machine learning with spatial statistics to improve the accuracy of plot-level AGB estimates. The proposed method integrates the nonlinear mapping capabilities of machine learning algorithms with the spatial autocorrelation and stratified heterogeneous advantages of a spatial statistical model. Our aim is to answer two specific questions: (1) What are the differences in the accuracy of forest AGB estimates based on the different methods? (2) Can the integration of spatial statistics and machine learning methods improve the accuracy of AGB estimation models at the plot level? We explore these two questions by studying an empirical case for predicting plot-level AGB in a Eucalyptus plantation in Nanjing County, China.

Study Area
The study area was in Fujian Province, Nanjing County, China (117 • 00 -117 • 36 E, 24 • 26 -25 • 00 N, Figure 1). The region has a South Asian tropical monsoon climate. The average annual temperature in Nanjing County is 21 • C, with an annual precipitation of 1700 mm and 340 frost-free days per year. The major soil type is red soil. The study area has a complex topography with significantly varying elevation (0-1566 m). Seventy-four percent (145,009 ha) of the county comprises forests, and 79,346 ha are plantations. The main tree species are Eucalyptus grandis x urophylla, Pinus massoniana, and Cunninghamia lanceolata (Lamb.) Hook. Forest composition, structure, and biomass are spatiotemporally heterogeneous. Over the past decade, the area of Eucalyptus plantations increased by 10,862 ha, reaching 13,338 ha.

Study Area
The study area was in Fujian Province, Nanjing County, China (117°00′-117°36′ E, 24°26′-25°00′ N, Figure 1). The region has a South Asian tropical monsoon climate. The average annual temperature in Nanjing County is 21 °C, with an annual precipitation of 1700 mm and 340 frost-free days per year. The major soil type is red soil. The study area has a complex topography with significantly varying elevation (0-1566 m). Seventy-four percent (145,009 ha) of the county comprises forests, and 79,346 ha are plantations. The main tree species are Eucalyptus grandis x urophylla, Pinus massoniana, and Cunninghamia lanceolata (Lamb.) Hook. Forest composition, structure, and biomass are spatiotemporally heterogeneous. Over the past decade, the area of Eucalyptus plantations increased by 10,862 ha, reaching 13,338 ha.

Reference AGB Data
This dataset contains the reference AGB of all 30 inventory plots in the study area. A total of 30 inventory plots were selected in 2012. The 30 inventory plots included ten Eucalyptus plantation age groups, and each group included three inventory plots. The plots

Reference AGB Data
This dataset contains the reference AGB of all 30 inventory plots in the study area. A total of 30 inventory plots were selected in 2012. The 30 inventory plots included ten Eucalyptus plantation age groups, and each group included three inventory plots. The plots were in the eastern section of the study area ( Figure 1). In each plot (0.04 ha, 20 m × 20 m), we measured the diameter at breast height (DBH) and tree height (H) of all living stems. In addition, we investigated the environmental data and forest attributes of inventory plots, including stand age, stand density, longitude, latitude, and altitude. Table 1 lists the statistics of the inventory plots. The AGB of each plot was obtained by summing the AGB of all trees within the plot. The AGB of each tree was estimated using H, DBH, and three allometric models built in this study.
Obtaining reference AGB data involved the following steps: (1) destructive sampling in inventory plots (i.e., tree harvest), (2) construction of tree-level allometric models, and (3) calculation of reference AGB of inventory plots. Details on obtaining reference AGB data are provided in Section A of the Supplementary Material.

Destructive Sampling in Inventory Plots: Tree Harvest
Trees were harvested from the 30 inventory plots. Three trees with a DBH close to the mean DBH of trees in each plot were cut down, for a total of 90 trees harvested from the 30 plots. We then measured the H and DBH of each harvested tree and weighed the biomass of each organ (foliage, stems, and branches) to obtain the AGB of each harvested tree. Details on the selection of the standard wood and cutting process are provided in Section A of the Supplementary Material. Due to environmental impacts, trees of the same forest age may also have different chemical, physical, and energy characteristics that may affect their biomass [23]. To harvest as many representative trees as possible, we harvested three trees per plot.

Construction of Tree-Level Allometric Models
We divided the 90 harvested trees into three age groups (1-2 years, 3-5 years, and 6-10 years) to construct tree-level allometric models: where a and b are constant coefficients. Table 2 lists the parameters used in these models. The tree-level allometric models (Table 2) were then applied to each tree in each inventory plot according to tree age, DBH, and H. The resulting tree AGB was aggregated on the plot level, thereby producing a reference AGB for each of the 30 inventory plots.
The reference AGB for the 30 inventory plots ranged from 1.02 to 135.79 Mg·ha −1 , with an average value of 47.34 Mg·ha −1 and a standard deviation of 34.46 Mg·ha −1 . The coefficients of variation of the AGB for all inventory plots and for the ten age categories were 0.73 and 0.07-0.37, respectively.

Spatial Characteristic Test of Forest Reference AGB
The spatial-characteristic test of forest reference AGB is the premise for using the P-BSHADE spatial statistical model. We used the R software to calculate Moran's I [24] to evaluate the spatial autocorrelation of the reference AGB between inventory plots. Spatially stratified heterogeneity refers to the within-strata variance being less than the betweenstrata variance. The spatially stratified heterogeneity of the reference AGB was evaluated by using a q-statistic generated by the GeogDetector model. GeogDetector is a software tool that analyzes the spatial variation in the geographic strata (or hierarchical) of variables [25]. First, we used the K-means algorithm to obtain the strata of the reference AGB. Next, we regarded the reference AGB as Y and the strata of the reference AGB as X and inserted them into the GeogDetector model to obtain the q statistics [25,26].

Selection of Variables
To create the plot-level model, we first identified predictor variables. Based on our previous work [27], we selected environmental covariates, including longitude, latitude, and altitude, and forest attribute variables, including stand density, plot mean DBH, plot mean H, and forest age. Pearson's correlation coefficient was used to investigate the correlation between these variables and the reference AGB of the inventory plots.

Model Development
Seven models, including three machine learning models (Figure 2a The spatial statistical model P-BSHADE requires AGB-related variables. In this case study, we used the estimated AGB data of the plots as the AGB-related variables (see "estimated AGB data of plots" in Figure 2). For the combined machine learning and spatial statistical models, the estimated AGB data of the plots were obtained from the results of SVM (Figure 2a The models were trained by using the R Development Core Team. The spatial statistical model P-BSHADE requires AGB-related variables. In this case study, we used the estimated AGB data of the plots as the AGB-related variables (see "estimated AGB data of plots" in Figure 2). For the combined machine learning and spatial statistical models, the estimated AGB data of the plots were obtained from the results of SVM [ Figure 2a], RBF-ANN [ Figure 2b], or RF [ Figure 2c]. For the P-BSHADE model only [ Figure 2d], we established a plot-level allometric model to obtain the estimated AGB data of the plots.
The models were trained by using the R Development Core Team.

Machine Learning
SVM, RF, and ANN are robust machine learning models for estimating forest biomass and volume, mapping soil organic carbon stocks, and geological mapping [11,[28][29][30]. Therefore, we selected SVM, RBF-ANN, and RF machine learning models for this research.
SVM is a method of supervised machine learning that is often used to solve classification problems or regression problems. The basic principle of SVM for classification is to find a hyperplane in the feature space and separate the positive and negative samples with the minimum misclassification rate [31]. The principle of SVM for regression is very

Machine Learning
SVM, RF, and ANN are robust machine learning models for estimating forest biomass and volume, mapping soil organic carbon stocks, and geological mapping [11,[28][29][30]. Therefore, we selected SVM, RBF-ANN, and RF machine learning models for this research.
SVM is a method of supervised machine learning that is often used to solve classification problems or regression problems. The basic principle of SVM for classification is to find a hyperplane in the feature space and separate the positive and negative samples with the minimum misclassification rate [31]. The principle of SVM for regression is very similar to that of SVM for classification [32]. SVM for regression retains the main ideas of minimizing error and individualizing the hyperplane that maximizes the margin, and part of the error is acceptable [32].
RBF-ANN is a three-layer neural network model that includes an input layer, a hidden layer (a Gaussian RBF is used here), and an output layer. The transformation from input space to hidden space may be linear or nonlinear, whereas the transformation from hidden space to output space is linear. The function of the hidden layer is to map the vector from the indivisible low-dimensional linear state to the separable high-dimensional linear state to greatly accelerate the learning and convergence speed and avoid becoming stuck in a local optimum [33].
RF is a machine learning technique and data mining method that combines selflearning technologies and was developed by Breiman in 2001 [34]. RF combines tree predictors such that each tree depends on the values of a random vector that is sampled independently and with the same distribution for all trees in the forest. RF provides accurate predictions [34].
The schematic function for machine learning is where y j is the AGB of the jth inventory plot predicted by a machine learning model (j = 1, . . . , n, n = 30), f (. . .) is a machine learning model represented by a function of x j,k (k = 1, . . . , 4), and x j,1 , x j,2 , x j, 3 , and x j,4 are the central longitude, mean DBH, mean H, and forest age of the jth inventory plot, respectively.

Spatial Statistical Model: P-BSHADE
P-BSHADE is an optimal linear unbiased estimate-interpolation method based on the assumption of the simultaneous existence of the spatial autocorrelation and heterogeneity of the target object. This method is empirically superior for solving the problem of an unrepresentative sample when estimating [35,36].
The core of the model is to minimize the variances between predicted error and unbiased estimation. The prediction process of the P-BSHADE model requires strong spatio-temporal coordination between the predictive variable and the related variable to spatially interpolate the predictive variable. P-BSHADE differs markedly from the kriging and inverse distance weighting (IDW) algorithms. Compared with kriging and IDW, the application of P-BSHADE to forest AGB interpolation has the following advantages: (1) The spatial distribution of forest AGB is characterized by spatial autocorrelation and heterogeneity, which are taken into account in the P-BSHADE model. Considering spatial heterogeneity eliminates the difference in forest AGB distribution caused by varying terrains or different geographic locations. However, kriging and IDW only consider the spatial correlation between plots. (2) In addition, P-BSHADE considers strongly correlated inventory plots as neighboring plots, whereas the kriging and IDW algorithms consider sites that are close in proximity.
In brief, the P-BSHADE model includes three steps: First, the model obtains reference AGB for all inventory plots by using the allometric model. Second, it uses the reference AGB of the target inventory plot and the reference AGB of other inventory plots to obtain the weight relationship between the target inventory plot and the other inventory plots. Third, it uses the reference AGB of the other inventory plots and weights in Equation (3) to predict the AGB of the inventory plots. The specific mathematical formula for the P-BSHADE model is (see Supplementary Material [35,36]) whereŷ j is the estimated AGB of the jth inventory plot using the P-BSHADE model (j = 1, 2, . . ., n, n = 30), y i is the reference AGB of the ith inventory plot (i = 1, 2, . . ., n, n = 30), w ij is calculated by dividing the reference AGB of the ith inventory plot by the (allometric model) estimated AGB of the jth inventory plot (when j = 1, i = 2, 3, . . ., 30; when j = 2, i = 1, 3, 4, . . ., 30).

Combination of Machine Learning and Spatial Statistical Models
Considering the inherent advantages and disadvantages of P-BSHADE and machine learning, this study investigates whether they can improve the accuracy of forest AGB estimates. Therefore, P-BSHADE was separately integrated with the three machine learning methods (SVM, RBF-ANN, and RF) to form three combined models (SVM & P-BSHADE, RBF-ANN & P-BSHADE, and RF & P-BSHADE). The specific and most important combination step was to use the estimated results from the machine learning models as reference data. Therefore, the equation of each combined model is the same as Equation (3), with the description of w ij replaced by "w ij is calculated by dividing the reference AGB of the ith inventory plot by the machine-learning-estimated AGB of the jth inventory plot (when j = 1, i = 2, 3, . . . , 30; when j = 2, i = 1, 3, 4, . . . , 30)".

Split Datasets
Considering that the 30 inventory plots were relatively centralized and that spatial autocorrelation of the data was possible, we used a spatial block cross-validation strategy to split datasets to avoid overfitting of the model when using machine learning. The major difference between standard random k-fold cross-validation and spatial block crossvalidation is how the datasets are split into folds. Numerous studies have shown that spatial block cross-validation is more suitable than random k-fold cross-validation in the field of ecology when using environmental data because it reduces model overfitting [37][38][39][40]. More information on this topic and an R package of the spatial block k-fold cross-validation are available in [40].
After the spatial block, the 30 inventory plots were split into 12 clusters (see the 12 squares of Figure S1 of the Supplementary Material), and the 12 clusters were then split tenfold to apply tenfold cross-validation ( Figure S2 shows one instance of tenfold cross-validation).

Calculate Indicators
To evaluate the accuracy of the AGB estimates of the seven models (SVM, RBF-ANN, RF, P-BSHADE, SVM & P-BSHADE, RBF-ANN & P-BSHADE, and RF & P-BSHADE), the estimated AGB results were compared with the reference AGB of the inventory plots. As performance indicators, we used the mean absolute error (MAE), mean relative error (MRE), root mean square error (RMSE), and normalized root mean square error (nRMSE), which are given by where y p i is the predictive value of the different models, y i is the reference AGB of the ith inventory plot, and n is the number of training datasets. We then used the calculated MAE, MRE, and RMSE to identify the optimal model.

Robustness of Combined Models
To evaluate the robustness of the combined machine learning and spatial statistical models for different time series, we selected 22 independent sample plots and made nondestructive measurements of each tree in July 2019. We repeated the process used for constructing the plot-level model and evaluated the models. We then used the accuracyassessment indexes (MAE, MRE, RMSE, and nRMSE) to determine whether the combined models were more accurate than the single models.

Spatial Characteristic Test of Forest Reference AGB
The spatial autocorrelation test produced a Moran's I value of 0.36, a z score of 4.78, and a p value less than 0.01 (p < 0.01). The spatial distribution of the reference AGB revealed a pattern of aggregation (see red regions in Figure S3 in the Supplementary Material). Less than 1% of the AGB data were randomly distributed (see blue regions in Figure S3 in the Supplementary Material), and the possibility of an aggregated distribution was greater than that of a random distribution. These results suggest that the spatial distribution of the AGB data reveals aggregation and a pattern of strong spatial autocorrelation.
The spatially stratified heterogeneity test produced a q value of 0.87 and a p value less than 0.01, which indicate that the within-layer variances were far less than the sum of variances between different strata. The results indicate that the reference AGB of the 30 inventory plots can be described by spatially stratified heterogeneity. Figure 3 shows the correlation-coefficient matrix of variables. The following variables were strongly correlated with AGB: longitude (r = −0.56), DBH (r = 0.79), H (r = 0.84), and forest age (r = 0.82). Thus, we selected four variables (longitude, DBH, H, and forest age) as covariates for the AGB plot-level models. Table 1 lists the statistical descriptions of these covariates and the AGB statistics for the 30 inventory plots. models for different time series, we selected 22 independent sample plots and made nondestructive measurements of each tree in July 2019. We repeated the process used for constructing the plot-level model and evaluated the models. We then used the accuracy-assessment indexes (MAE, MRE, RMSE, and nRMSE) to determine whether the combined models were more accurate than the single models.

Spatial Characteristic Test of Forest Reference AGB
The spatial autocorrelation test produced a Moran's I value of 0.36, a z score of 4.78, and a p value less than 0.01 (p < 0.01). The spatial distribution of the reference AGB revealed a pattern of aggregation (see red regions in Figure S3 in the Supplementary Material). Less than 1% of the AGB data were randomly distributed (see blue regions in Figure  S3 in the Supplementary Material), and the possibility of an aggregated distribution was greater than that of a random distribution. These results suggest that the spatial distribution of the AGB data reveals aggregation and a pattern of strong spatial autocorrelation.
The spatially stratified heterogeneity test produced a q value of 0.87 and a p value less than 0.01, which indicate that the within-layer variances were far less than the sum of variances between different strata. The results indicate that the reference AGB of the 30 inventory plots can be described by spatially stratified heterogeneity. Figure 3 shows the correlation-coefficient matrix of variables. The following variables were strongly correlated with AGB: longitude (r = −0.56), DBH (r = 0.79), H (r = 0.84), and forest age (r = 0.82). Thus, we selected four variables (longitude, DBH, H, and forest age) as covariates for the AGB plot-level models. Table 1 lists the statistical descriptions of these covariates and the AGB statistics for the 30 inventory plots.

Performance of Plot-Level Models
An allometric model was applied to compare with our seven models. The indicators of the allometric model are greater than those of the seven models developed in this study (Figure 4). The forest AGB estimates obtained by the P-BSHADE method are similar to those obtained by the three machine learning methods.
( Figure 4). The forest AGB estimates obtained by the P-BSHADE method are similar to those obtained by the three machine learning methods.
Of the three machine learning methods, the SVM was the most accurate.    Figure 4).
We now compare the machine learning methods with combined methods over two periods, 2012 and 2019 ( Figure 5). Whether considering the data from 2012 or the new data collected in 2019, the combined model is more accurate than the single machine learning model, although p-values change with the methods and with predictive indicators. (Table 3). These results suggest that the combined models are more accurate than single machine learning models, and these improvements are robust. 3). These results suggest that the combined models are more accurate than single machine learning models, and these improvements are robust.  Table 3. The prediction indicator differs between the three machine learning methods and the combined methods. A paired t-test was used to examine if the predictive indicators (i.e., MAE and RMSE) of the combined methods are greater than those of the machine learning method alone based on the results of the 10-fold cross-validation.

Significance and Challenges of Accurate Plot-Level AGB Estimates
Plot-level AGB models link tree-level AGB measurements to regional-scale AGB estimates. Ignoring the uncertainty of plot-level models underestimates the total uncertainty of pixel-level estimates by 6% [2]. In the field of forest biomass estimation, the accurate Figure 5. Improved accuracy assessment indexes of three combined machine learning and spatial statistical methods are revealed by comparison with three corresponding machine learning methods. Panels (a-d) show the MAE, MRE, RMSE, and nRMSE, respectively; S1 vs. S5 compares S5 with S1, S2 vs. S6 compares S6 with S2, and S3 vs. S7 compares S7 with S3 (S1 = SVM, S2 = RBF-ANN, S3 = RF, S5 = SVM & P-BSHADE, S6 = RBF-ANN & P-BSHADE, S7 = RF & P-BSHADE). Table 3. The prediction indicator differs between the three machine learning methods and the combined methods. A paired t-test was used to examine if the predictive indicators (i.e., MAE and RMSE) of the combined methods are greater than those of the machine learning method alone based on the results of the 10-fold cross-validation.

Significance and Challenges of Accurate Plot-Level AGB Estimates
Plot-level AGB models link tree-level AGB measurements to regional-scale AGB estimates. Ignoring the uncertainty of plot-level models underestimates the total uncertainty of pixel-level estimates by 6% [2]. In the field of forest biomass estimation, the accurate estimation of forest AGB or structure at the plot level is very important for calibrating and validating large-scale forest biomass. However, the distribution of most AGB is either non-Gaussian, skewed, or multimodal, especially in tropical and subtropical regions [41]. Different intensities and factors are coupled together, resulting in high heterogeneity and clear nonlinearity in the spatial distribution of forest AGB. These spatial characteristics of forest AGB make it more difficult to accurately estimate AGB at the plot level.
Here, we integrate the advantages of machine learning and spatial statistics to construct a plot-level AGB model for a subtropical region. Combining the advantages of machine-learning-based quantification of complex nonlinear relationships between AGB and the multiple covariates, in conjunction with the P-BSHADE model, allows the spatial autocorrelation and heterogeneity of forest AGB to be incorporated into the model.

Why a Combined Model Outperforms a Single Machine Learning or Spatial Statistical Model
As expected, the combined methods were more accurate than any single method (either machine learning or spatial statistics). This may be due to the advantages of machine learning, which compensates for the inherent defects of the P-BSHADE model, and vice versa.
On the one hand, the P-BSHADE model has its own merits: (1) It considers the spatial autocorrelation and spatial heterogeneity of the distribution of the target objects, not only to solve the difference between target objects caused by the different terrain or geographic location but also to solve the problem of strong correlation between target objects with remote geographic locations due to similar terrain conditions. (2) The P-BSHADE model calculates the covariance of the reference data of research objects, that is, the reference AGB of inventory plots in our study. This method is more reliable because it avoids the secondorder stationary hypothesis (i.e., when using the kriging algorithm, semivariograms require this hypothesis), which does not correspond with the actual situation. (3) P-BSHADE regards strongly correlated plots as neighboring plots. However, the P-BSHADE model is also handicapped by the fact that the founding assumption does not conform to reality. The assumption is that the estimated AGB is accurate in all inventory plots except in the target inventory plot. In other words, the premise behind using the P-BSHADE model is that the reference AGB data are accurate or strongly correlated with AGB. In reality, the AGB of each inventory plot has a varying degree of uncertainty when using the single P-BSHADE model because the reference AGB data were obtained from the allometric model. As the P-BSHADE model combined with machine learning uses the results optimized by machine learning as the reference AGB data, it further improves the accuracy of AGB mapping.
On the other hand, machine learning also has advantages and disadvantages. Machine learning has the advantage of being able to handle complex, potentially nonlinear relationships between forest AGB and other variables. However, the initial samples of machine learning are randomly selected, which may lead to differences in the results of each operation of the model. In addition, one of the important shortcomings is overfitting, which may perform poorly beyond the training data range. In contrast with machine learning, the P-BSHADE model considers the spatial autocorrelation and spatial heterogeneity of forest AGB and of environmental covariates and the bias of the observed values of the inventory plots, which better corresponds to actual situations. A combined model uses the results of machine learning as the reference AGB data of P-BSHADE, so that the fitting process of the combined model better accounts for spatial relationships than does the single machine learning model. The result is improved accuracy.
Machine learning models or the P-BSHADE model have been used to model the uncertainty of temperature measurements obtained by weather stations [20,29,36]. However, the methods used in these studies were adopted independently. Conversely, the combination of machine learning and spatial statistics can improve the prediction accuracy of AGB maps, which in turn can be used as criteria for improving the accuracy of LiDAR remote-sensing technology and the results of ecological process models. Eventually, these improvements can promote process-oriented projects that require dynamic AGB predictions for large-scale forests in different forest management scenarios.
Our combined methods produce a very small RMSE for the prediction accuracy of AGB, which we explain as follows: (1) The reference AGB of the 30 inventory plots were obtained by summing the AGB of each tree, which was calculated using allometric models constructed from accurate measurements of harvested trees. (2) Machine learning methods were used to quantify the complex nonlinear relationship between AGB and multiple environmental covariates. (3) We applied a spatial statistical method based on the hypothesis of spatial heterogeneity. Although the nRMSE index is calculated in different studies using different datasets and prediction methods in different locations, most studies use nRMSE as an indicator for quantifying the AGB prediction errors of models [5,6,9]. In contrast with other studies, the present work not only focuses on subtropical forests but also uses methodological differences to mitigate uncertainty, especially by comprehensively addressing the sources of uncertainty caused by multiple spatial and environmental covariates.

Why the RF & P-BSHADE Method Outperforms Other Combined Methods
The three combined machine learning and spatial statistical methods produced more accurate AGB predictions than any individual method. The RF & P-BSHADE and SVM & P-BSHADE methods were significantly more accurate than the individual methods, whereas the RBF-ANN & P-BSHADE method was only slightly more accurate. The accuracies of the combined methods depend on the accuracy of the reference AGB data (machinelearning-predicted results) [36]. The combined method adds spatial information based on machine learning. In other words, more accurate predicted machine learning results lead to a greater accuracy of the combined method. This is an important scientific contribution of this work; namely, that combining machine learning and spatial statistics improves the accuracy of forest AGB estimation. Therefore, the different improvements offered by the three combined methods may be attributed to the following two mechanisms: (1) The RF has superior properties, and SVM models are easier to use and optimize than RBF-ANN [42]. The number of training samples determines the number of nodes in the hidden layer of the RBF-ANN model, and the number of nodes significantly affects the prediction accuracy. With only 30 training samples used in this study, the combined RBF-ANN & P-BSHADE approach did not significantly improve the prediction accuracy.
(2) RBF-ANN is more suitable for nonlinear stochastic dynamic systems [33], whereas the relationship between AGB and environmental covariates in this study is likely a monotonically increasing function.

Conclusions
Forest AGB estimates at the plot level play a major role in connecting accurate singletree AGB measurements to relatively difficult regional AGB estimates. However, AGB estimates at the plot level are plagued by numerous uncertainties. Improving the plotlevel model of forest AGB is a key issue in producing accurate AGB maps. A variety of prediction models have been applied to make accurate AGB estimates, all of which have their own advantages and disadvantages. Different approaches complement the advantages of different models and may yield more accurate AGB estimates than would otherwise be produced by using a single method. The main goal of the current study was to determine whether combining machine learning with spatial statistics can improve plot-level AGB estimates.
This study explores the prediction performance of different AGB models, and the results show that the model combining the random forest and P-BSHADE models substantially improves the accuracy of the estimates of forest AGB. The results of this study suggest that combining machine learning with spatial statistics improves plot-level AGB estimates. The understanding gained here should help to improve AGB mapping in other regions and in different types of forests.
Supplementary Materials: The following materials are available online at https://www.mdpi.com/ article/10.3390/f12121663/s1, Figure S1: Results of spatial blocking. The 30 blue dots represent 30 sample plots, green patches represent forest patches, and the 12 red rectangles represent the results of spatial blocking. Figure S2: The one situation of spatial block cross-validation. Figure S3: Spatial autocorrelation report.