Forest Aboveground Biomass Estimation and Response to Climate Change Based on Remote Sensing Data

: As the largest and most important natural terrestrial ecosystem, forest plays a crucial role in reducing the concentrations of greenhouse gases in the atmosphere, mitigating global warming, maintaining the global ecological balance, and promoting global biological evolution and community succession. The accurate and rapid assessment of forest biomass is highly signiﬁcant for estimating the regional carbon budget and monitoring forest change. In this study, Landsat images and China’s National Forest Continuous Inventory data of 1999, 2004, 2009, and 2014 were used to establish extreme gradient boosting (XGBoost) models for forest aboveground biomass (AGB) estimation based on forest type in the Xiangjiang River Basin, Hunan Province, China. Kriging interpolation of the AGB residuals was used to correct the error of AGB estimation. Then, a new XGBoost model was established using the ﬁnal corrected AGB maps and climate data to estimate the AGB under different climate scenarios during the 2050s and 2070s. The results indicated that AGB estimation using the XGBoost model with correction via Kriging interpolation of the AGB residuals can signiﬁcantly improve the accuracy of AGB estimation. The total AGB of the study area increased over time from 1999 to 2014, indicating that the forest quality improved in the study area. Under the different climate scenarios, the total AGB during the 2050s and 2070s was predicted to decline continuously with increasing of greenhouse gas emissions, indicating that greenhouse gas emissions have a negative impact on forest growth. The results of this study can provide data support for evaluating the ecological function and value of forest ecosystems, and for formulating reasonable forest management measures to mitigate the effects of climate change.


Introduction
The Fifth Assessment Report (AR) of the Intergovernmental Panel on Climate Change (IPCC) confirms that the global climate has experienced significant changes characterized by global warming over the past century [1].The sixth AR of the IPCC provides stronger evidence to further confirm the objective facts of global warming in the last 100 years and a clearer signal of the impact of human activities on climate warming [2].Global climate change has a significant impact on natural ecosystems and the social economy, which will threaten human survival and pose long-term and severe challenges to the sustainable development of human societies; thus, global climate change has evolved from a scientific issue to a political, economic, and environmental issue of global concern [3][4][5].As an important part of terrestrial ecosystems, the forest ecosystem has an important impact on the terrestrial biosphere and other surface processes, and plays a crucial role in promoting global biological evolution and community succession and maintaining global ecological balance [6][7][8][9].
Forest biomass, which is the most basic quantitative feature of the forest ecosystem and an important indicator of its carbon sources and sinks, can reflect the complex relationship between forest material cycles and energy flows with the environment [10][11][12][13].Changes of forest aboveground biomass (AGB) can reflect the quality and status of the forest ecosystem, as well as the effects of natural disturbances, human activities, and climate change [14][15][16].Therefore, regional and global estimates of AGB under a changing climate not only provide a theoretical basis for studies of the terrestrial ecosystem carbon cycle and global climate change, which play an important role in understanding and monitoring the response of forest ecosystems to greenhouse gas emissions, but also provide strategic guidelines for sustainable forest management, which is crucial for the rational use of forest resources and improving the forest ecological environment [17][18][19][20].
A rapid and accurate estimation of AGB on a regional scale has always been a key issue in forest research.The main approaches used for AGB estimation in previous studies can be classified into (1) field measurement and (2) combined use of remote sensing data and field measured plot data [21,22].The field measurement method is suitable for a small forest stand or forest sample plot; it cannot be used to study the distribution and variation of AGB at the regional scale because it is too costly, time consuming, labor intensive, and damaging to the forest [23,24].It should be emphasized that the field measurement plot data is essential for AGB estimation in combination with remote sensing data.Remote sensing, as an emerging Earth observation technology, not only provides multi-scale data covering a large area for the rapid and quantitative evaluation of regional AGB but also enables real-time monitoring of surface parameters for analyzing and evaluating the spatiotemporal patterns and changes in trends of AGB [25][26][27].Numerous studies showed that remote sensing data was highly correlated with AGB, effectively predicting and monitoring AGB at a regional scale [28][29][30][31].Therefore, regional AGB estimation can be predominantly conducted using remote sensing data.The AGB estimation based on remote sensing data mainly involves two aspects: remote sensing data source and estimation model.Various remote sensing data (mainly including optical remote sensing, microwave radar, and LiDAR) have been used to estimate the regional forest AGB [32][33][34].Among all the available satellites, Landsat has certain advantages: it provides cross-calibrated data spanning more than 40 years, has high spectral characteristics and spatial fidelity for global coverage, and provides free and open access to details [34][35][36].The estimation model includes an empirical model, a physical model, a mechanism model (also called a process model), and a comprehensive model [8,[37][38][39]].The empirical model, which establishes the relationship between the measured plot AGB and the remote sensing variables based on statistical analysis, is the most commonly used method for AGB estimation based on remote sensing data [21,40].
Xiangjiang River, which is an important tributary of the Yangtze River, is regarded as the mother river in Hunan Province, China, with the Xiangjiang River Basin forming the core economic region and development zone in the Hunan Province.Therefore, rapid and accurate estimates of the status and dynamic changes of AGB in the Xiangjiang River Basin are required to evaluate the forest ecological service function of the basin and help formulate scientific policies and development models for coping with climate change, protecting ecological security, and promoting sustainable forest development in the basin.In this study, we use National Forest Continuous Inventory (NFCI) data and Landsat data in combination with the extreme gradient boosting (XGBoost) model to establish an AGB estimation model and draw AGB maps for the study area.On the basis of these results, we analyze the response of AGB to climate change, and estimate the distribution of AGB under different future climate change scenarios.The accuracy of XGBoost for AGB estimation based on remote sensing data was significantly improved compared with the accuracy of linear regression [22].Furthermore, the methods used in this study such as XGBoost, spatial analysis, and Kriging interpolation, are common methods, which are convenient for other study areas.

Study Area
The Xiangjiang River (9.47 × 10 4 km 2 , 24 • 31 N-29 • 01 N, 110 • 30 E-114 • 15 E) is the largest tributary of the Dongting Lake in the Yangtze River Basin (Figure 1).The Xiangjiang River Basin, which is located within both the Yangtze River Economic Belt and South China Economic Circle, is characterized by densely distributed towns, a concentrated population, convenient transportation, and a developed economy; thus, it is the core area of socioeconomic development in Hunan Province.The Xiangjiang River Basin is a long basin with complex topography; the east, south, and west parts are surrounded by mountains, the north is relatively flat, and the center is predominantly hills, with the terrain tilted from south to north in a horseshoe shape (Figure 1).The entire basin is located in a humid, subtropical monsoon climate zone with abundant sunlight, water, and heat, and obvious seasonality.It is hot in summer and cold in winter, with an annual average temperature of 17.5 • C. Annual precipitation is unevenly distributed throughout the year, i.e., concentrated in spring and summer, and the annual average precipitation and evaporation are 1400 mm and 1200 mm, respectively.The typical vegetation is subtropical evergreen forest, with a forest cover percentage of 54.4% [41].

Study Area
The Xiangjiang River (9.47 × 10 4 km 2 , 24°31′ N-29°01′ N, 110°30′ E-114°15′ E) is the largest tributary of the Dongting Lake in the Yangtze River Basin (Figure 1).The Xiangjiang River Basin, which is located within both the Yangtze River Economic Belt and South China Economic Circle, is characterized by densely distributed towns, a concentrated population, convenient transportation, and a developed economy; thus, it is the core area of socioeconomic development in Hunan Province.The Xiangjiang River Basin is a long basin with complex topography; the east, south, and west parts are surrounded by mountains, the north is relatively flat, and the center is predominantly hills, with the terrain tilted from south to north in a horseshoe shape (Figure 1).The entire basin is located in a humid, subtropical monsoon climate zone with abundant sunlight, water, and heat, and obvious seasonality.It is hot in summer and cold in winter, with an annual average temperature of 17.5 °C .Annual precipitation is unevenly distributed throughout the year, i.e., concentrated in spring and summer, and the annual average precipitation and evaporation are 1400 mm and 1200 mm, respectively.The typical vegetation is subtropical evergreen forest, with a forest cover percentage of 54.4% [41].

Forest Inventory Data
Forest inventory data are important for AGB estimation based on remote sensing.In this study, the NFCI data of the Xiangjiang River Basin in Hunan Province from 1999, 2004, 2009, and 2014 were used.The NFCI data are currently the highest-quality data available at the provincial scale in China [42,43].The size of each square plot is 25.82 m × 25.82 m (approximately 0.0667 ha), and the plots were systematically allocated based on a grid of 4 km × 8 km [34].The data of the plots situated on non-forestry land (such as water areas, cropland, urban land, and bare land) or covered by clouds in the remote sensing images were eliminated.Ultimately, 1685 plots were used for modeling in this study (Figure 2).

Remote Sensing Data
Landsat Surface Reflectance products, which were derived from Landsat 5 TM (Thematic Mapper), Landsat 7 ETM+ (Enhanced Thematic Mapper), and Landsat 8 OLI (Operational Land Imager) satellite images from 1999, 2004, 2009, and 2014, were used in this study.Li et al. analyzed Landsat images taken in various seasons to examine the accuracy of AGB estimation and demonstrated the highest levels of accuracy in images taken in autumn [45].Therefore, Landsat images taken in autumn were considered for use in this study.The images were downloaded from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/,accessed on 1 May 2021).There were 30 scene images (Figure 1 and Table S1 in the Supplementary Information section).The images were resampled to a pixel size of 25.82 m × 25.82 m, which was the same size as the inventory plots.Lu et al. explored the relationship between the Landsat 5 TM image textures and forest AGB in Rondônia, the Brazilian Amazon, in 2005.The results indicated the correlation between the textures and forest AGB, which might improve the AGB estimation [46].Therefore, the textures were considered for use in this study.The texture images were calculated using a grey-level co-occurrence matrix algorithm with 3 × 3-, 5 × 5-, and 7 × 7pixel windows.In addition, 20 vegetation indices were generated for this study (Table 3; Table S2, in the Supplementary Information section).Landsat data were processed by En- The AGB of a plot is the sum of the AGBs of all the trees recorded in the plot.The AGB of a tree was calculated by using the general one-variable aboveground biomass model, which can be expressed as follows [44]: where M a (kg) is the AGB of a tree, D (cm) is the diameter of the tree at breast height, and p (g/cm 3 ) is the basic wood density (Table A1).The plot AGB was converted to biomass per hectare (Mg/ha).
Based on the species standing volume according to China's technical regulation for forest continuous inventory, the plots were classified into three types, i.e., broadleaved, coniferous, and mixed forest (Table 1 and Figure 2).The statistical results of AGB are reported in Table 2.

Remote Sensing Data
Landsat Surface Reflectance products, which were derived from Landsat 5 TM (Thematic Mapper), Landsat 7 ETM+ (Enhanced Thematic Mapper), and Landsat 8 OLI (Operational Land Imager) satellite images from 1999, 2004, 2009, and 2014, were used in this study.Li et al. analyzed Landsat images taken in various seasons to examine the accuracy of AGB estimation and demonstrated the highest levels of accuracy in images taken in autumn [45].Therefore, Landsat images taken in autumn were considered for use in this study.The images were downloaded from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/,accessed on 1 May 2021).There were 30 scene images (Figure 1 and Table S1 in the Supplementary Information section).The images were resampled to a pixel size of 25.82 m × 25.82 m, which was the same size as the inventory plots.Lu et al. explored the relationship between the Landsat 5 TM image textures and forest AGB in Rondônia, the Brazilian Amazon, in 2005.The results indicated the correlation between the textures and forest AGB, which might improve the AGB estimation [46].Therefore, the textures were considered for use in this study.The texture images were calculated using a grey-level co-occurrence matrix algorithm with 3 × 3-, 5 × 5-, and 7 × 7-pixel windows.In addition, 20 vegetation indices were generated for this study (Table 3; Table S2, in the Supplementary Information section).Landsat data were processed by Environment for Visualizing Images (ENVI) software (version 5.3.1;Boulder, CO, USA).

Land Cover Data
The European Space Agency (ESA) Climate Change Initiative (CCI) project consistently delivered annual global land cover (LC) maps at a 300 m spatial resolution from 1992 to 2015 [47].The CCI-LC maps for 1999, 2004, 2009, and 2014 were downloaded from the ESA website (http://maps.elie.ucl.ac.be/CCI/viewer/index.php, accessed on 1 May 2021).The CCI-LC maps were then consolidated into seven types based on the typology of the CCI-LC, i.e., broadleaved, coniferous, and mixed forests, as well as cropland, urban land, water, and others (including bare land, grassland, etc.) (Figure 1).
To validate the accuracy and consistency between the classifications of the NFCI and CCI-LC, the attributes of the CCI-LC maps were extracted by the NFCI plot center.The CCI-LC maps were highly accurate, and can therefore satisfy the research needs of this study.For example, the producer accuracies of the CCI-LC map of broadleaved, coniferous, and mixed forests in 2014 were 0.91, 0.88, and 0.82, respectively, and the user accuracies were 0.93, 0.91, and 0.92, respectively.Furthermore, the overall accuracy and kappa coefficient of broadleaved, coniferous, and mixed forests were 0.92 and 0.88, respectively.The CCI-LC maps were then converted from raster to polygon to clip the AGB maps based on forest type.

Climate Data
The climate data used in this study can be divided into two types, namely historical climate data and future climate scenario data.The climate scenarios for the 2050s and 2070s were derived from the downscaled global climate model (GCM) of the IPCC AR5.The Beijing Climate Center Climate System Model version 1.1 (BCC-CSM1-1), which was developed by the Beijing Climate Center of the China Meteorological Administration, provided the data on China for the IPCC AR5 [48].The BCC-CSM1-1 included data on four climate change scenarios under different representative concentration pathways (RCPs), i.e., RCP2.6,RCP4.5, RCP6.0, and RCP8.5, which represent possible ranges of total radiative force (the cumulative measure of human emissions of greenhouse gases from all sources expressed) values are respectively 2.6, 4.5, 6.0, and 8.5 W/m 2 for the year 2100 [48,49].The climate scenario data used in this study were downloaded from the WorldClim Global Climate Data website (http://www.worldclim.org/,accessed on 1 June 2016).The data include 19 climate variables derived from monthly temperature and precipitation values to generate more biologically meaningful data that reflect a range of temperature and precipitation summaries (e.g., trends, seasonality, and extremes) (Table 4).The spatial resolution of the data was 30" (approximately 1 km × 1 km).

Methods
This study first has accurately estimated the AGB in 1999, 2004, 2009, and 2014, and then based on this result has estimated AGB under different climate scenarios for the future.

AGB Estimation Combined with Spatial Interpolation
The non-parametric XGBoost model was used to estimate the AGB in this study.XG-Boost was proposed by Chen et al. [50] in Kaggle machine learning competitions.XGBoost performs second-order Taylor expansion for the objective function and uses the second derivative to accelerate the convergence speed of the model during training.Furthermore, a regularization term is added to the objective function to control the complexity of the tree, thereby obtaining a simpler model and avoiding overfitting [50].XGBoost is a highly scalable tree-structured enhanced model that can effectively handle sparse and missing data and greatly improve the speed of the algorithm and compress the computational memory in large-scale data training [50].XGBoost exhibits excellent performance and is very popular in data mining and machine learning activities globally [34].
In this study, XGBoost was used to establish AGB models for different forest types.The model errors were then corrected by spatial interpolation to improve the accuracy of AGB estimation.The procedure comprised the following three steps (Figure 3).

AGB Estimation Combined with Spatial Interpolation
The non-parametric XGBoost model was used to estimate the AGB in this study.XGBoost was proposed by Chen et al. [50] in Kaggle machine learning competitions.XGBoost performs second-order Taylor expansion for the objective function and uses the second derivative to accelerate the convergence speed of the model during training.Furthermore, a regularization term is added to the objective function to control the complexity of the tree, thereby obtaining a simpler model and avoiding overfitting [50].XGBoost is a highly scalable tree-structured enhanced model that can effectively handle sparse and missing data and greatly improve the speed of the algorithm and compress the computational memory in large-scale data training [50].XGBoost exhibits excellent performance and is very popular in data mining and machine learning activities globally [34].
In this study, XGBoost was used to establish AGB models for different forest types.The model errors were then corrected by spatial interpolation to improve the accuracy of AGB estimation.The procedure comprised the following three steps (Figure 3).(1) AGB modeling and estimation.The XGBoost models based on forest type were established and the optimal predictor variables were selected based on variable importance.The AGB and AGB residuals were then calculated for the entire study area.XGBoost defines two measures for variable importance that can be used to rank variables: the first measure is calculated by the fractional contribution (Gain) of each feature to the model based on the total gain of this variable's splits; the second measure is calculated by the relative number (Frequency) of times a feature be used in the trees [51].More important predictor variables have higher Gain and Frequency percentages.The method of variable selection was derived from the study of Li et al. [34].(2) Spatial analysis and spatial interpolation of AGB residuals.First, the Global Moran's Index was used to evaluate the spatial autocorrelation of AGB residuals.Then, if the spatial autocorrelation of AGB residuals was determined, the spatial interpolation of AGB residuals was performed by Kriging interpolation.The Global Moran's Index can measure the spatial distribution pattern of AGB residuals, i.e., clustered, dispersed, or random [52].This index and both the z-score and p-value can be used to evaluate the significance of that spatial distribution pattern.In general, the Global Moran's Index value is bounded by −1.0 and 1.0, i.e., a positive correlation exists if (1) AGB modeling and estimation.The XGBoost models based on forest type were established and the optimal predictor variables were selected based on variable importance.The AGB and AGB residuals were then calculated for the entire study area.XGBoost defines two measures for variable importance that can be used to rank variables: the first measure is calculated by the fractional contribution (Gain) of each feature to the model based on the total gain of this variable's splits; the second measure is calculated by the relative number (Frequency) of times a feature be used in the trees [51].More important predictor variables have higher Gain and Frequency percentages.The method of variable selection was derived from the study of Li et al. [34].(2) Spatial analysis and spatial interpolation of AGB residuals.First, the Global Moran's Index was used to evaluate the spatial autocorrelation of AGB residuals.Then, if the spatial autocorrelation of AGB residuals was determined, the spatial interpolation of AGB residuals was performed by Kriging interpolation.The Global Moran's Index can measure the spatial distribution pattern of AGB residuals, i.e., clustered, dispersed, or random [52].This index and both the z-score and p-value can be used to evaluate the significance of that spatial distribution pattern.In general, the Global Moran's Index value is bounded by −1.0 and 1.0, i.e., a positive correlation exists if the value is greater than zero, a negative correlation exists if the value is less than zero, and no correlation exists if the value equals zero.Kriging interpolation, which is a geostatistical method of interpolation, assumes that the distance or direction between sample points reflects a spatial correlation that can be used to explain variation in the surface; it fits a mathematical function to a specified number of measured points, or all measured points within a specified radius, to determine the value for each point [53,54].Moreover, Kriging interpolation has the capacity to produce a prediction surface as well as some measure of the certainty or accuracy of the predictions [54].
Ordinary Kriging, which is the most common and extensively used Kriging method, was used to predict the AGB residual maps for different forest types in this study.(3) Correction and evaluation of the AGB map.The final AGB maps were obtained by correcting the original predicted AGB maps with the AGB residual interpolation maps, then evaluating the accuracy of the corrected maps.

Response of AGB to Climate Change
In this study, the final predicted AGB values for 1999, 2004, and 2009 were used as the response variables, whereas the climate variables in the corresponding years were used as the predictor variables to establish the XGBoost model for AGB estimation.Then, the AGB result for 2014 was used as validation data to evaluate the performance of the model.Finally, the AGB maps were predicted under different climate change scenarios in the 2050s and 2070s using the developed model (Figure 3).

Evaluation Criteria
In this study, the coefficient of determination (R 2 ), the root-mean-square error (RMSE) and the percentage root-mean-square error (RMSE%) were also used to evaluate the performance of the models [34]: where y i is the observed AGB value, ŷi is the predicted AGB value based on models, y is the arithmetic mean of all the observed AGB values, and n is the sample number.In general, a higher R 2 value and lower RMSE and RMSE% values indicate a better estimation performance of the model.The method of cross validation was also used to verify the performance of the model and evaluate the generalization ability of the model.The recommended 10-fold cross validation procedure was selected.A standard deviation map, which can analyze the stability of cross validation and the uncertainty of AGB estimation, was generated using the estimation results of 10-fold cross validation.A higher standard deviation indicates higher uncertainty, whereas a lower standard deviation indicates lower uncertainty.Figure 5 shows the final selected predictor variables according to variable importance for the different forest types in the four years.Despite the fact that the predictor variables differed for the different forest types in the four years, the textures of bands 3 (Red band), 4 (NIR band), and 6 (SWIR2 band) were extremely important in the models.The correlation and mean values of these bands were included in all models, which indicated that these variables contained abundant information that enhances the performance of AGB estimation models.In addition, the spectral variables, particularly B3 (Red band), were also very important.Moreover, according to the value and shape of the graph of variable importance, XGBoost models tended to centralize the importance at a single variable.For example, B3 (Red band) was significantly correlated with B3T7Mea at a significance level of 0.01 with a value of 0.854, whereas the importance of B3 (Red band) was significantly higher than that of B3T7Mea for the broadleaved forest model in 2014.Figure 5 shows the final selected predictor variables according to variable importance for the different forest types in the four years.Despite the fact that the predictor variables differed for the different forest types in the four years, the textures of bands 3 (Red band), 4 (NIR band), and 6 (SWIR2 band) were extremely important in the models.The correlation and mean values of these bands were included in all models, which indicated that these variables contained abundant information that enhances the performance of AGB estimation models.In addition, the spectral variables, particularly B3 (Red band), were also very important.Moreover, according to the value and shape of the graph of variable importance, XGBoost models tended to centralize the importance at a single variable.For example, B3 (Red band) was significantly correlated with B3T7Mea at a significance level of 0.01 with a value of 0.854, whereas the importance of B3 (Red band) was significantly higher than that of B3T7Mea for the broadleaved forest model in 2014.Figure 5 shows the final selected predictor variables according to variable importance for the different forest types in the four years.Despite the fact that the predictor variables differed for the different forest types in the four years, the textures of bands 3 (Red band), 4 (NIR band), and 6 (SWIR2 band) were extremely important in the models.The correlation and mean values of these bands were included in all models, which indicated that these variables contained abundant information that enhances the performance of AGB estimation models.In addition, the spectral variables, particularly B3 (Red band), were also very important.Moreover, according to the value and shape of the graph of variable importance, XGBoost models tended to centralize the importance at a single variable.For example, B3 (Red band) was significantly correlated with B3T7Mea at a significance level of 0.01 with a value of 0.854, whereas the importance of B3 (Red band) was significantly higher than that of B3T7Mea for the broadleaved forest model in 2014.

Evaluation of XGBoost Models
Model performance was expressed by scatterplots, which show the relationship between the predicted and observed AGB values (Figure 6).The results showed that the broadleaved forest model had the highest accuracy, followed by coniferous and mixed forests models.Moreover, the combined accuracy of broadleaved, coniferous, and mixed forests models was higher than that of the model for all forest plots with no classification of forest type.This indicates that AGB estimation is more accurate when performed according to forest type.The models had the highest accuracy in 2004, followed by 2009, 2014, and 1999.According to the statistical results, the RMSE% values were relatively high when AGB was low (<30 Mg/ha), but very low otherwise (less than 50%) which indicates that the accuracy of AGB estimation was relatively high (Table A2).However, problems of over-and underestimation of AGB based on remote sensing data using the XGBoost model, which have been noted in previous studies, were observed in all models, although model performance was significantly improved.

Evaluation of XGBoost Models
Model performance was expressed by scatterplots, which show the relationship between the predicted and observed AGB values (Figure 6).The results showed that the broadleaved forest model had the highest accuracy, followed by coniferous and mixed forests models.Moreover, the combined accuracy of broadleaved, coniferous, and mixed forests models was higher than that of the model for all forest plots with no classification of forest type.This indicates that AGB estimation is more accurate when performed according to forest type.The models had the highest accuracy in 2004, followed by 2009, 2014, and 1999.According to the statistical results, the RMSE% values were relatively high when AGB was low (<30 Mg/ha), but very low otherwise (less than 50%) which indicates that the accuracy of AGB estimation was relatively high (Table A2).However, problems of over-and underestimation of AGB based on remote sensing data using the XGBoost model, which have been noted in previous studies, were observed in all models, although model performance was significantly improved.

Spatial Interpolation of AGB Residuals
The spatial correlation of AGB residuals was analyzed after the AGB residuals were calculated for all models.Table 5 shows that the Global Moran's Index for different forest types in the four years.The Global Moran's Indices of AGB residuals were above zero for all models, indicating that the AGB residuals were positively spatially correlated.In addition, the z-scores and p-values of AGB residuals indicated a statistically significant aggregation rather than a random distribution pattern.Therefore, Kriging interpolation was deemed appropriate for the interpolation of AGB residuals (Tables S3-S5 and Figure S1, in the Supplementary Information section).

Spatial Interpolation of AGB Residuals
The spatial correlation of AGB residuals was analyzed after the AGB residuals were calculated for all models.Table 5 shows that the Global Moran's Index for different forest types in the four years.The Global Moran's Indices of AGB residuals were above zero for all models, indicating that the AGB residuals were positively spatially correlated.In addition, the z-scores and p-values of AGB residuals indicated a statistically significant aggregation rather than a random distribution pattern.Therefore, Kriging interpolation was deemed appropriate for the interpolation of AGB residuals (Tables S3-S5 and Figure S1, in the Supplementary Information section).Figure 7 presents scatterplots between the AGB residuals and predicted AGB residuals for different forest types.The accuracy was highest for broadleaved forest, followed by coniferous and mixed forest.Furthermore, the problem of under-and overestimation also occurred using Ordinary Kriging interpolation of the AGB residuals; thus, the variability of the AGB residuals of the sample points was not sufficiently estimated using Ordinary Kriging.Although the high and low values were similar between AGB values and the AGB residuals, the error derived from Ordinary Kriging interpolation of AGB residuals will be transferred to the AGB map when using the Kriging interpolation map to correct the errors in the AGB model; meaning that the problems of under-and overestimation of AGB were not completely eliminated.Figure 7 presents scatterplots between the AGB residuals and predicted AGB residuals for different forest types.The accuracy was highest for broadleaved forest, followed by coniferous and mixed forest.Furthermore, the problem of under-and overestimation also occurred using Ordinary Kriging interpolation of the AGB residuals; thus, the variability of the AGB residuals of the sample points was not sufficiently estimated using Ordinary Kriging.Although the high and low values were similar between AGB values and the AGB residuals, the error derived from Ordinary Kriging interpolation of AGB residuals will be transferred to the AGB map when using the Kriging interpolation map to correct the errors in the AGB model; meaning that the problems of under-and overestimation of AGB were not completely eliminated.8 and 9 show the predicted AGB residual maps and the standard error maps for different forest types, respectively.The predicted AGB residual map of the combined forest model revealed a higher range of AGB residuals and a higher standard deviation than the model for all plots with no classification of forest type, which indicates that the Ordinary Kriging interpolation of AGB residuals based on forest type improved the estimation accuracy of AGB residuals.Moreover, the combined model exhibited a lower average standard error and a larger area of low value.Furthermore, the standard errors were lower than the AGB residuals, which also indicates that the Ordinary Kriging method achieves more accurate estimation of AGB residuals.
Sustainability 2022, 14, x FOR PEER REVIEW 13 of 27 Figures 8 and 9 show the predicted AGB residual maps and the standard error maps for different forest types, respectively.The predicted AGB residual map of the combined forest model revealed a higher range of AGB residuals and a higher standard deviation than the model for all plots with no classification of forest type, which indicates that the Ordinary Kriging interpolation of AGB residuals based on forest type improved the estimation accuracy of AGB residuals.Moreover, the combined model exhibited a lower average standard error and a larger area of low value.Furthermore, the standard errors were lower than the AGB residuals, which also indicates that the Ordinary Kriging method achieves more accurate estimation of AGB residuals.

AGB Mapping
The final AGB maps were obtained after correction using the spatial interpolation map of AGB residuals.Figure 10 shows the estimation accuracy of the final AGB maps for different forest types in the four years.The estimation accuracy was highest for broadleaved forest, followed by coniferous and mixed forests.Moreover, the accuracy of the combination of broadleaved, coniferous, and mixed forests was greatly improved from that of all forest plots with no classification of forest type, and the accuracy of the final corrected AGB maps was greatly improved from that of the original AGB map created by XGBoost models.The RMSE and RMSE% values of the final corrected AGB maps were lower than the original AGB maps, which further verifies that AGB maps corrected by AGB residual maps can effectively improve the accuracy of AGB estimation based on remote sensing (Table A3).
The F-test was used to determine any significant differences between the results of original and corrected AGB maps based on forest type and those based on all plots with no classification (Figure 11).The confidence level was set to 95%.The F-test results showed significant differences in the predicted AGB maps at a confidence level of 95%, although the p-values and t-values differed.The F-test results also indicated the need to improve the accuracy of AGB estimation corrected by AGB residuals.

AGB Mapping
The final AGB maps were obtained after correction using the spatial interpolation map of AGB residuals.Figure 10 shows the estimation accuracy of the final AGB maps for different forest types in the four years.The estimation accuracy was highest for broadleaved forest, followed by coniferous and mixed forests.Moreover, the accuracy of the combination of broadleaved, coniferous, and mixed forests was greatly improved from that of all forest plots with no classification of forest type, and the accuracy of the final corrected AGB maps was greatly improved from that of the original AGB map created by XGBoost models.The RMSE and RMSE% values of the final corrected AGB maps were lower than the original AGB maps, which further verifies that AGB maps corrected by AGB residual maps can effectively improve the accuracy of AGB estimation based on remote sensing (Table A3).
The F-test was used to determine any significant differences between the results of original and corrected AGB maps based on forest type and those based on all plots with no classification (Figure 11).The confidence level was set to 95%.The F-test results showed significant differences in the predicted AGB maps at a confidence level of 95%, although the p-values and t-values differed.The F-test results also indicated the need to improve the accuracy of AGB estimation corrected by AGB residuals.

Model Construction
The Pearson correlation coefficients between the climate variables and AGB indicated that all climate variables had a significance level of 0.01 with AGB; the two variables with the highest correlation coefficient were Bio06 (minimum temperature of the coldest month) and Bio01 (annual mean temperature), with a value of −0.62 and −0.50, respectively (Table S6, in the Supplementary Information section).
Figure 13 shows the validation result and variable importance of the XGBoost model based on climate variables.The predicted AGB values of the XGBoost model were mainly distributed in the middle range (50-150 Mg/ha).The accuracy of the XGBoost model based on climate variables for high AGB values (>30 Mg/ha) with a lower RMSE% was greater than that of the XGBoost model based on Landsat variables because there were more samples available to establish the decision tree in this range; however, the model had relatively low accuracy for low AGB values (<30 Mg/ha) with a higher RMSE% (Table 6).Similarly, the problem of under-and overestimation also existed in this XGBoost model.In terms of the importance of variables, the two most important variables for Gain were Bio06 (minimum temperature of the coldest month) and Bio14 (precipitation of the driest month), whereas the two most important variables for Frequency were Bio01 (annual mean temperature) and Bio02 (mean diurnal range).

Model Construction
The Pearson correlation coefficients between the climate variables and AGB indicated that all climate variables had a significance level of 0.01 with AGB; the two variables with the highest correlation coefficient were Bio06 (minimum temperature of the coldest month) and Bio01 (annual mean temperature), with a value of −0.62 and −0.50, respectively (Table S6, in the Supplementary Information section).
Figure 13 shows the validation result and variable importance of the XGBoost model based on climate variables.The predicted AGB values of the XGBoost model were mainly distributed in the middle range (50-150 Mg/ha).The accuracy of the XGBoost model based on climate variables for high AGB values (>30 Mg/ha) with a lower RMSE% was greater than that of the XGBoost model based on Landsat variables because there were more samples available to establish the decision tree in this range; however, the model had relatively low accuracy for low AGB values (<30 Mg/ha) with a higher RMSE% (Table 6).Similarly, the problem of under-and overestimation also existed in this XGBoost model.In terms of the importance of variables, the two most important variables for Gain were Bio06 (minimum temperature of the coldest month) and Bio14 (precipitation of the driest month), whereas the two most important variables for Frequency were Bio01 (annual mean temperature) and Bio02 (mean diurnal range).Figure 14 shows the distribution of AGB under four different climate scenarios predicted for the 2050s and 2070s.The overall spatial distribution trends of AGB for climate scenarios in the 2050s and 2070s were similar to the current distribution of AGB, that is, high AGB was predominantly distributed in areas of high altitudes and less human activity, whereas low AGB was predominantly distributed in the valleys and plains with low altitude and frequent human activities in the central areas of the study area.However, the spatial distribution of AGB exhibited significant differences between the four different climate scenarios.The mean AGB for the study area in the 2050s was 207.15, 192.30, 180.24,

AGB Estimation under Climate Scenarios
Figure 14 shows the distribution of AGB under four different climate scenarios predicted for the 2050s and 2070s.The overall spatial distribution trends of AGB for climate scenarios in the 2050s and 2070s were similar to the current distribution of AGB, that is, high AGB was predominantly distributed in areas of high altitudes and less human activity, whereas low AGB was predominantly distributed in the valleys and plains with low altitude and frequent human activities in the central areas of the study area.However, the spatial distribution of AGB exhibited significant differences between the four different climate scenarios.The mean AGB for the study area in the 2050s was 207.15, 192.30, 180.24, and 176.26Mg/ha for RCP2.6,RCP4.5, RCP6.0, and RCP8.5, respectively, and the mean AGB in the 2070s was 196.55, 191.80, 175.08, and 170.58 Mg/ha, respectively.Thus, the total AGB consistently decreased in the following order: RCP2.6 > RCP4.5 > RCP6.0 > RCP8.5, indicating that greenhouse gas emissions have a negative effect on forest growth and quality improvement.

Methods of Improving AGB Estimation
In this study, 170 predictor variables were extracted from Landsat images to establish the XGBoost model of AGB estimation.However, not all predictor variables can be used for modeling and estimation due to their high correlations and high number, which affect the performance and operational efficiency of the XGBoost model.Variable selection, which is one of the most important processes of non-parametric models, can reduce the memory space and data dimensions, increase the speed of calculation, and improve the performance and interpretability of the model [55].The results of this study revealed that the performance of XGBoost models was significantly improved through variable selection, which reduced the number of predictor variables from hundreds to merely several (Figures 4 and 5).Moreover, texture and spectral variables were very important in this study.The role of spectral variables and texture variables differed depending on the forest structure.The multiple species and canopy layers of broadleaved and mixed forest were sufficiently expressed by abundant spectral information for AGB estimation; thus, the models tended to choose spectral variables.Conversely, the more limited species composition of coniferous forest meant that the spectral characteristics were not significantly different, whereas the textural information effectively explained the AGB estimation; thus, the models tended to select texture variables [21,56,57].In addition, the decision tree

Methods of Improving AGB Estimation
In this study, 170 predictor variables were extracted from Landsat images to establish the XGBoost model of AGB estimation.However, not all predictor variables can be used for modeling and estimation due to their high correlations and high number, which affect the performance and operational efficiency of the XGBoost model.Variable selection, which is one of the most important processes of non-parametric models, can reduce the memory space and data dimensions, increase the speed of calculation, and improve the performance and interpretability of the model [55].The results of this study revealed that the performance of XGBoost models was significantly improved through variable selection, which reduced the number of predictor variables from hundreds to merely several (Figures 4 and 5).Moreover, texture and spectral variables were very important in this study.The role of spectral variables and texture variables differed depending on the forest structure.The multiple species and canopy layers of broadleaved and mixed forest were sufficiently expressed by abundant spectral information for AGB estimation; thus, the models tended to choose spectral variables.Conversely, the more limited species composition of coniferous forest meant that the spectral characteristics were not significantly different, whereas the textural information effectively explained the AGB estimation; thus, the models tended to select texture variables [21,56,57].In addition, the decision tree of XGBoost was generated based on the previous tree to correct the residuals error; therefore, if a predictor variable was selected by the previous tree, the importance of other highly correlated predictor variables would be greatly reduced in the subsequent decision tree [34,58].
Traditional statistical theory assumes that the studied variable is a random variable.In fact, many factors of the forest ecosystem, such as volume, biomass, and soil physiochemical properties, are not completely independent and randomly distributed due to the influence of the natural environment and human disturbances; instead, they are autocorrelated within a certain spatial range, showing structural unity and randomness [59].Geostatistical methods are based on the regionalized variable theory and take the empirical variogram as a basic tool to study phenomena exhibiting either random and structural spatial distributions, or spatial autocorrelation and dependence [60][61][62].Because of its significant advantage for the spatial calculation and analysis of uncertain factors, the geostatistical method has become the preferred method of studying the spatial randomness and structure of forest ecosystem factors.
As an important method of geostatistics, Kriging can provide interpolated surfaces as well as a measure of uncertainty, indicating the accuracy of the predictions [63].The measurement of uncertainty is critical to informed decision making, as it provides information on the possible values for each location rather than just one interpolated value [63].The results of this study showed that the predicted standard error of Kriging interpolation of AGB residuals was relatively low; thus, the highly accurate Kriging interpolation maps of AGB can be used for AGB correction (Figure 7).However, the results also showed that the variability of AGB residuals was not adequately estimated by Kriging; thus, problems of under-and overestimation also occurred in the Kriging interpolation results, which led to errors in the correction of AGB maps using Kriging interpolation maps (Figure 10).
The non-parametric XGBoost model was used for AGB estimation in combination with the Kriging interpolation of AGB residuals in this study.The results indicated that the problems of under-and overestimation of AGB based on remote sensing data were efficiently reduced and the accuracy of AGB estimation was significantly improved, although this error still existed after correction by AGB residuals.Most early AGB estimation studies were based on parametric models, such as linear regression, multiple regression, and non-linear regression [22].Parametric models have ideal assumptions for data distribution on which the model is established, e.g., the data obeys a normal distribution [64].However, in reality, it is difficult to estimate the distribution of data or there are no obvious features due to the complex relationship between the remote sensing factors and AGB.In contrast, non-parametric models are a statistical analysis method that makes no assumptions about the overall distribution of the sample and directly analyzes the sample [65].Non-parametric models have a strong fitting ability and can obtain a good prediction result, although needing more data and more time for training, and there is a risk of overfitting the training data [65].
The problem of under-and overestimation, which also existed in previous studies, was the main factor influencing the accuracy of AGB estimation based on remote sensing data in this study [66].This problem is partly determined by the algorithm itself.The decision trees, which are key components of the XGBoost model, were built on the training data; thus, the accuracy is reduced when extrapolations are performed outside the training data [34].However, the under-and overestimation problem is also related to the data.For remote sensing data, saturation in multispectral data is the main reason for underestimation of plots with high AGB values [21,67,68].The pixels in Landsat images with relatively low spatial resolution are mixed, so cannot accurately express the spectral information of land cover for plots with low AGB values.Furthermore, NFCI data contain fewer plots with high or low AGB values because of the large proportion of secondary forest with an uneven age distribution, which is composed primarily of young and middle-age trees and few mature and over-mature trees; thus, the XGBoost model is unable to obtain sufficient information to establish the decision trees.
To solve the problem of under-and overestimation and further improve the accuracy of AGB estimation, the following method is suggested for future research.For remote sensing data, high spatial and radiometric resolution data, such as hyperspectral data, synthetic aperture radar data, and LiDAR data, can be solely or synergistically used for AGB estimation.Recent research has proved that this approach can effectively reduce the problem of mixed pixels and multispectral saturation, thereby improving the accuracy of multi-sensor data synergy for AGB estimation.As such, this approach has become a new trend in AGB estimation based on remote sensing data [69][70][71].In addition, as well as the variables used in this study, other factors such as soil, climate, topography, and phenology can be as the predictor variables for AGB estimation.For NFCI data, it is recommended to establish and survey as many temporary plots (especially plots with high or low AGB values) as possible in accordance with the technical regulations of the continuous NFCI of China to ensure the accuracy and consistency of NFCI data and increase the amount of training samples.In addition, existing plots should be screened to select representative plots and eliminate unrepresentative plots to ensure a more reasonable representation of the distribution of AGB in the study area.

Impact of Climate Change on Forests
Numerous studies have revealed the current and future impacts of climate on plant phenology, forest succession, the composition and distribution of forest structure, biodiversity, forest productivity, and the forest carbon sink function of forest ecosystems [5,[72][73][74].Close material and energy exchanges exist between the climate and forests, which affect and restrict each other [75][76][77][78].Therefore, the forest ecosystem, as the largest carbon source and sink in the terrestrial ecosystem, will respond to climate change and exhibit feedbacks and regulation effects in relation to climate change [79][80][81][82].
The direction and magnitude of climate change impacts on forests vary according to environmental factors and forest types [38,83].Previous studies have shown that the productivity of forest ecosystems typically increases due to the growth extension benefit caused by rising temperatures and the fertilization benefit caused by rising CO 2 concentrations; however, some studies have reached the opposite conclusion, that is, climate change has a negative impact on the productivity of forest ecosystems [81, [84][85][86].The results of this study indicated that climate change has a negative impact on AGB.
As an important aspect of global change, climate change is predominantly manifested as temperature increases, precipitation pattern changes, and greater climate extremes [48].In the forest ecosystem, stomatal conductance between the leaves and canopy will decrease significantly, leading to inhibited plant growth when plants are stressed by high temperatures and CO 2 concentrations [87][88][89].The greenhouse effect results in a higher rate of temperature increase at night than during the day; thus, the dark respiration of plants cancels out the organic matter of forests accumulated during the day [90].Meanwhile, the greenhouse effect increases soil temperatures, which stimulates the growth and activity of edaphon and increases the mineralization rate of soil organic matter; however, the carbon nitrogen ratio of soil will also increase with an increase of soil carbon, which will inhibit the respiration of edaphon, slowing the litter decomposition rate and preventing the prompt replacement of nutrients, thereby hindering normal plant growth [91,92].
Water availability also has a significant impact on plant height, leaf area, branch number, photosynthesis, and growth [93].The greenhouse effect accelerates the transpiration of plants and increases soil evaporation; thus, plants are subjected to water stress, generating physiological drought, which will affect the photosynthesis and growth of plants when the soil moisture content is below the minimum threshold required by plants [72,94].Climate change not only causes temperature changes in different regions, but also changes in precipitation.An increase in precipitation during the driest month will result in a longer growth season, thereby facilitating plant growth [95].Furthermore, climate extremes have important impacts on the diversity, distribution, and growth of plants [96,97].For example, a decrease in the minimum temperature of the coldest month will result in premature freez-ing injuries to plants, and long-term low temperatures will lead to the death of plants [98].Conversely, an increase in the maximum temperature of the warmest month will destroy the water balance of plants and promote the coagulation of proteins and the internal mechanism of harmful metabolites, thereby hindering plant growth [99].The results of this study indicated that variables related to climate extremes were very important for AGB estimation based on climate variables (Figure 13).
The temporal changes in AGB were similar under all climate scenarios in this study.The differences in AGB in different climate scenarios were slightly greater than the change in AGB over time (Figure 14 and Table 7).We hypothesized that the AGB dynamics may be driven by climate change.We also hypothesized that the spatiotemporal changes in AGB in response to climatic changes were equivalent, but the impact of climate change on AGB actually differed among forest types, stand ages, and regions [76,93,100], hence lacking a consensus on the impact of climate change on AGB [38,[101][102][103].The impact of climate change on AGB is a relatively long-term process [104,105], and long-term field measurement is necessary to accurately predict future changes in AGB under changing climate regimes.In addition, human activity, which is another major factor affecting forest ecosystems, has more direct, rapid, and drastic impacts [106,107].The IPCC has revealed the unprecedented speed and scale of the impact of human activities on environmental change over the last 100 years, and shown that climate change is closely related to human activities [1].Humans can implement a variety of forest management measures that will either exacerbate or slow down the impact of climate change caused by human activities [108,109].The estimation of AGB distribution under the influence of human activities and climate change, as well as the scientific and rapid assessment of future forest changes, are crucial for formulating targeted management plans.In this study, climate data were used to establish an AGB estimation model and predict AGB under four different climate scenarios in the 2050s and 2070s; however, we did not consider the impact of human activities on the forest ecosystem or climate.This was predominantly because the comprehensive impact of human activities is too complex to quantitatively evaluate, and there is no reasonable index for quantifying the impact of human activities or predicting future human activities.In the future, basic research on the response of forest ecosystems to climate change should be conducted at different levels (e.g., typical tree species, typical forest communities, and forest landscapes), along with forest management planning with a focus on climate change adaptation, and quantitative scientific research on the impact of human activities on forest ecosystems and climate change.Furthermore, it will be important to explore the comprehensive impacts of human activities and climate change on forest ecosystems from the perspective of human-natural system coupling and the need for multidisciplinary development.
It is important to notice that the results and conclusions of this study cannot be directly used at different spatial and temporal scales, because forest stands with different AGB and different geography conditions have different characteristics reflected in remote sensing images.The methods adopted in this study, such as XGBoost model, semivariance analysis, Kriging interpolation, and evaluation criteria, have been frequently used in previous studies as well and are easy to implement in the spatial analysis of AGB estimation.

Conclusions
In this study, the variable selection method based on variable importance was used to select optimal variables from Landsat predictor variables, then establish an XGB model for AGB estimation according to forest type.Then, Kriging interpolation of the AGB residuals was employed to correct the errors of the XGBoost model.Subsequently, a new XGBoost model was established based on climate variables using the final corrected AGB maps to estimate the AGB under four different climate scenarios during the 2050s and 2070s.The major conclusions were as follows.(1) Variable selection significantly improved the performance of the XGBoost model.(2) Kriging interpolation of AGB residuals based on spatial autocorrelation improved the accuracy of AGB estimation.(3) The total AGB increased and the forest quality improved over time in the study area.(4) The total AGB under different climate scenarios decreased in the following order: RCP2.6 > RCP4.5 > RCP6.0 > RCP8.5, indicating that greenhouse gas emissions have a negative effect on forest growth and quality.The observed negative impacts of climate change on forest growth and productivity have prompted humans to implement a variety of measures to protect forests.However, the impact of forest management on ecosystems and the evaluation of its effectiveness are long-term processes; therefore, the impact of climate change must be considered when formulating long-term forest management plans.Determining the dynamic changes of AGB and the main driving variables under different climate scenarios is highly significant for formulating appropriate forest management plans.

Figure 1 .
Figure 1.Study area details: (a) location of the study area in China, (b) elevation, Landsat scene image numbers (P: path, R: row), and important cities (red circles), and (c) land cover map of 2014.

Figure 1 .
Figure 1.Study area details: (a) location of the study area in China, (b) elevation, Landsat scene image numbers (P: path, R: row), and important cities (red circles), and (c) land cover map of 2014.

Figure 2 .
Figure 2. Spatial distribution of the sample plots, including the forest types and observed AGB from 1999.

Figure 2 .
Figure 2. Spatial distribution of the sample plots, including the forest types and observed AGB from 1999.

Figure 3 .
Figure 3. Schematic of the workflow: AGB estimation with spatial interpolation (orange part), and AGB estimation under different climate scenarios for the future (blue part).

Figure 3 .
Figure 3. Schematic of the workflow: AGB estimation with spatial interpolation (orange part), and AGB estimation under different climate scenarios for the future (blue part).
Figure 4  illustrates how the R 2 values changed with the count of variables selected for XGBoost models in the different years.Each line represents an independent model, and the different colors indicate the different forest types.The results showed that the R 2 values of models increased as the number of predictor variables decreased, which indicates that variable selection can effectively reduce the data dimensions and improve the performance and interpretability of models.

Figure 4 .
Figure 4. Accuracy XGBoost models with the selected variable count changing based on different forest types."All", results of all forest plots with no classification of forest type.

Figure 5 .
Figure 5. Variable importance of XGBoost models of different forest types.The variable importance of each model was scaled to sum to 100."All", results of all forest plots with no classification of forest type.

Figure 4 .
Figure 4. Accuracy XGBoost models with the selected variable count changing based on different forest types."All", results of all forest plots with no classification of forest type.

Figure 4 .
Figure 4. Accuracy XGBoost models with the selected variable count changing based on different forest types."All", results of all forest plots with no classification of forest type.

Figure 5 .
Figure 5. Variable importance of XGBoost models of different forest types.The variable importance of each model was scaled to sum to 100."All", results of all forest plots with no classification of forest type.

Figure 5 .
Figure 5. Variable importance of XGBoost models of different forest types.The variable importance of each model was scaled to sum to 100."All", results of all forest plots with no classification of forest type.

Figure 6 .
Figure 6.Scatter plots of predicted and observed AGB for XGBoost models of different forest types in different years."All", results of all forest plots with no classification of forest type."Combined", merged with the respective results of broadleaved, coniferous, and mixed forest.

Figure 6 .
Figure 6.Scatter plots of predicted and observed AGB for XGBoost models of different forest types in different years."All", results of all forest plots with no classification of forest type."Combined", merged with the respective results of broadleaved, coniferous, and mixed forest.

Figure 7 .
Figure 7. Evaluation of cross validation of the Ordinary Kriging interpolation of AGB residuals for different forest types."All", results of all forest plots with no classification of forest type."Combined", merging the respective results of broadleaved, coniferous, and mixed forest.

Figure 7 .
Figure 7. Evaluation of cross validation of the Ordinary Kriging interpolation of AGB residuals for different forest types."All", results of all forest plots with no classification of forest type."Combined", merging the respective results of broadleaved, coniferous, and mixed forest.

Figures
Figures8 and 9show the predicted AGB residual maps and the standard error maps for different forest types, respectively.The predicted AGB residual map of the combined forest model revealed a higher range of AGB residuals and a higher standard deviation than the model for all plots with no classification of forest type, which indicates that the Ordinary Kriging interpolation of AGB residuals based on forest type improved the estimation accuracy of AGB residuals.Moreover, the combined model exhibited a lower average standard error and a larger area of low value.Furthermore, the standard errors were lower than the AGB residuals, which also indicates that the Ordinary Kriging method achieves more accurate estimation of AGB residuals.

Figure 8 .
Figure 8. AGB residual predicted using Ordinary Kriging interpolation for all plots with no classification of forest type (upper panels) and combined forest types (lower panels)."All", results of all forest plots with no classification of forest type."Combined", merging the respective results of broadleaved, coniferous, and mixed forest.

Figure 8 .
Figure 8. AGB residual predicted using Ordinary Kriging interpolation for all plots with no classification of forest type (upper panels) and combined forest types (lower panels)."All", results of all forest plots with no classification of forest type."Combined", merging the respective results of broadleaved, coniferous, and mixed forest.

Figure 8 .
Figure8.AGB residual predicted using Ordinary Kriging interpolation for all plots with no classification of forest type (upper panels) and combined forest types (lower panels)."All", results of all forest plots with no classification of forest type."Combined", merging the respective results of broadleaved, coniferous, and mixed forest.

Figure 9 .
Figure 9.Standard error of Ordinary Kriging interpolation of AGB residuals for all plots with no classification of forest type (upper panels) and combined forest types (lower panels)."All", results of all forest plots with no classification of forest type."Combined", merging the respective results of broadleaved, coniferous, and mixed forest.

Sustainability 2022 , 27 Figure 9 .
Figure 9.Standard error of Ordinary Kriging interpolation of AGB residuals for all plots with no classification of forest type (upper panels) and combined forest types (lower panels)."All", results of all forest plots with no classification of forest type."Combined", merging the respective results of broadleaved, coniferous, and mixed forest.

Figure 10 .
Figure 10.Scatterplots of the predicted AGB from the final corrected AGB maps and the observed AGB for different forest types in different years."All", results of all forest plots with no classification of forest type."Combined", merging the respective results of broadleaved, coniferous, and mixed forest.

Figure 10 .
Figure 10.Scatterplots of the predicted AGB from the final corrected AGB maps and the observed AGB for different forest types in different years."All", results of all forest plots with no classification of forest type."Combined", merging the respective results of broadleaved, coniferous, and mixed forest.

Figure 11 .
Figure 11.Comparisons (F-test) of original and corrected AGB maps in different years.Regarding the row and column labels, if the first letter is O, it refers to the original AGB; if it is C, it refers to the corrected AGB.If the second letter is A, it refers to the predicted result without classification of forest types; if it is C, it refers to the predicted result for the combination of broadleaved, coniferous, and mixed forests.Green cells are t-values and blue cells are p-values.

Figure 12
Figure12shows the final predicted AGB maps based on forest type and the standard deviation maps using 10-fold cross validation in 1999, 2004, 2009, and 2014.The mean of the AGB map was 54.623, 66.242, 68.079, and 77.579 Mg/ha, respectively, indicating an increase in total AGB and improved forest quality over time in the study area.The mean of the standard deviation map was 1.881, 3.109, 3.036, and 2.914 Mg/ha, respectively, indicating that the predicted result of 10-fold cross validation was stable with only slight changes.Regarding the spatial distribution of AGB, high AGB values were predominantly distributed in the southeastern and southern regions, which are characterized by a high altitude, steep slopes, high vegetation coverage, low population density, less human interference, and poor economic conditions.Conversely, low AGB value were predominantly distributed in the low hills and valleys of cropland, villages, and towns, which are characterized by a low altitude, gentle slopes, and more human interference.

Figure 12 .
Figure 12.Final predicted AGB maps based on forest type and standard deviation maps using 10fold cross validation.

Figure 11 .
Figure 11.Comparisons (F-test) of original and corrected AGB maps in different years.Regarding the row and column labels, if the first letter is O, it refers to the original AGB; if it is C, it refers to the corrected AGB.If the second letter is A, it refers to the predicted result without classification of forest types; if it is C, it refers to the predicted result for the combination of broadleaved, coniferous, and mixed forests.Green cells are t-values and blue cells are p-values.

Figure 12
Figure12shows the final predicted AGB maps based on forest type and the standard deviation maps using 10-fold cross validation in 1999, 2004, 2009, and 2014.The mean of the AGB map was 54.623, 66.242, 68.079, and 77.579 Mg/ha, respectively, indicating an increase in total AGB and improved forest quality over time in the study area.The mean of the standard deviation map was 1.881, 3.109, 3.036, and 2.914 Mg/ha, respectively, indicating that the predicted result of 10-fold cross validation was stable with only slight changes.Regarding the spatial distribution of AGB, high AGB values were predominantly distributed in the southeastern and southern regions, which are characterized by a high altitude, steep slopes, high vegetation coverage, low population density, less human interference, and poor economic conditions.Conversely, low AGB value were predominantly distributed in the low hills and valleys of cropland, villages, and towns, which are characterized by a low altitude, gentle slopes, and more human interference.

Figure 11 .
Figure 11.Comparisons (F-test) of original and corrected AGB maps in different years.Regarding the row and column labels, if the first letter is O, it refers to the original AGB; if it is C, it refers to the corrected AGB.If the second letter is A, it refers to the predicted result without classification of forest types; if it is C, it refers to the predicted result for the combination of broadleaved, coniferous, and mixed forests.Green cells are t-values and blue cells are p-values.

Figure 12
Figure12shows the final predicted AGB maps based on forest type and the standard deviation maps using 10-fold cross validation in 1999, 2004, 2009, and 2014.The mean of the AGB map was 54.623, 66.242, 68.079, and 77.579 Mg/ha, respectively, indicating an increase in total AGB and improved forest quality over time in the study area.The mean of the standard deviation map was 1.881, 3.109, 3.036, and 2.914 Mg/ha, respectively, indicating that the predicted result of 10-fold cross validation was stable with only slight changes.Regarding the spatial distribution of AGB, high AGB values were predominantly distributed in the southeastern and southern regions, which are characterized by a high altitude, steep slopes, high vegetation coverage, low population density, less human interference, and poor economic conditions.Conversely, low AGB value were predominantly distributed in the low hills and valleys of cropland, villages, and towns, which are characterized by a low altitude, gentle slopes, and more human interference.

Figure 12 .
Figure 12.Final predicted AGB maps based on forest type and standard deviation maps using 10fold cross validation.

Figure 12 .
Figure 12.Final predicted AGB maps based on forest type and standard deviation maps using 10-fold cross validation.

Figure 13 .
Figure 13.Scatterplot of validation result (a) and variable importance for the XGBoost model based on climate variables (b).The variable importance was scaled to sum to 100.

Figure 13 .
Figure 13.Scatterplot of validation result (a) and variable importance for the XGBoost model based on climate variables (b).The variable importance was scaled to sum to 100.

Figure 14 .
Figure 14.AGB predicted with the established XGBoost model based on climate variables for the 2050s and 2070s under four different climate scenarios.RCPs: Representative Concentration Pathway data defined by the possible range of radiative forcing values (2.6, 4.5, 6.0, and 8.5 W/m 2 , respectively) in the year 2100.

Figure 14 .
Figure 14.AGB predicted with the established XGBoost model based on climate variables for the 2050s and 2070s under four different climate scenarios.RCPs: Representative Concentration Pathway data defined by the possible range of radiative forcing values (2.6, 4.5, 6.0, and 8.5 W/m 2 , respectively) in the year 2100.

Table 1 .
Classification standard of forest types.

Table 2 .
Distribution of the plot AGB values (Mg/ha) of the different forest types.
Note: "All", results of all forest plots with no classification of forest type.

Table 3 .
Summary of predictor variables including Landsat Surface Reflectance band images, vegetation indices, and texture images for AGB estimation.

Table 4 .
Climate variables used for AGB estimation.Historical climatic data were derived from the Dataset of Monthly Surface Observation Values in Individual Years in China of the China Meteorological Data Service Centre (http://data.cma.cn/,accessed on 1 September 2022).The data include 23 climate elements from 75 surface meteorological observation stations in and around the study area.The same climate variables in 1999, 2004, 2009, and 2014 were produced by the historical climate data according to the calculation method of the climate scenario data using ANUSPLIN software (version 4.3; Canberra, ACT, Australia).

Table 5 .
Summary of the Global Moran's Index for different forest models.

Table 5 .
Summary of the Global Moran's Index for different forest models.

Table 6 .
Performance of the XGBoost model based on climate variables.

Table 6 .
Performance of the XGBoost model based on climate variables.

Table 7 .
Summary of AGB estimation for the 2050s and 2070s under four different climate scenarios.

Table A3 .
Summary of AGB maps corrected by the residual AGB interpolation maps.