Estimating forest aboveground carbon stock from remote sensing data with an enhanced XGBoost model

1 Introduction

China has demonstrated the most significant global increase in forest resources and maintains the world’s largest afforestation area. Furthermore, China plays a pivotal role in global greening initiatives. In particular, accurate estimation of forest aboveground carbon stock (AGC) is critical for understanding China’s carbon sequestration potential. Forests contribute approximately 90% of global atmospheric and terrestrial carbon fluxes, with an average carbon dioxide sequestration rate of 1.83 tons per cubic meter of tree growth (Yang et al., 2022). Despite covering only one-third of the global land surface, forest ecosystems store nearly 50% of terrestrial carbon stocks (Wang and Wang, 2015). Compared to alternative emission mitigation strategies, forest carbon sinks offer distinct advantages including cost-effectiveness, operational feasibility, and substantial co-benefits for both economic and ecological systems (Qian-qian et al., 2022). The total aboveground carbon stock in forest protected areas is estimated at 61.43 Gt, accounting for 26% of all terrestrial woody carbon. Compared to unprotected forests, this higher carbon stock is mainly attributed to the avoidance of emissions resulting from deforestation and forest degradation (Guo et al., 2023). The international community has established multiple treaties to address climate change challenges (Wei et al., 2022; Liang et al., 2022). Although forest soil also plays an important role in the global carbon cycle, storing approximately 73% of global soil carbon and regulating climate, the primary focus of this study is on aboveground carbon stocks (Li Y. et al., 2025). Climate change has profound impacts on global ecosystems, water security, and the frequency of extreme weather events; therefore, assessing changes in forest aboveground carbon stocks is crucial for understanding responses to climate change (Qingchun et al., 2022). Within the framework of global climate mitigation efforts, forest carbon storage research has attracted considerable scholarly attention worldwide. Accurate quantification of forest aboveground carbon stocks has emerged as a critical component in global climate change and carbon cycle studies (Siddiq et al., 2021; Xue-wei et al., 2019).

Conventional methods for estimating forest aboveground carbon storage predominantly depend on field surveys, which are not only resource-intensive in terms of manpower and materials but may also induce ecological disturbances in local forest ecosystems. Remote sensing technology has gained widespread application in regional-scale terrestrial carbon storage estimation due to its distinct advantages in spatial coverage and spatiotemporal resolution. Numerous studies have employed both passive and active remote sensing datasets for forest aboveground carbon stock estimation in specific geographical contexts (López-Serrano et al., 2020; Schumacher et al., 2019; Dong et al., 2020; Jianmin et al., 2022). These techniques have proven particularly valuable for extracting forest structural parameters owing to their information richness and superior spatial resolution. Optical remote sensing represents the longest-standing remote sensing modality for forest ecosystem carbon quantification (Cao et al., 2022). While passive optical remote sensing exhibits advantages in capturing horizontal vegetation structure, its effectiveness is constrained by forest phenology and meteorological conditions, often resulting in signal saturation that limits its utility for vertical forest aboveground carbon stock inversion (Zhang et al., 2019). SAR systems offer unique all-weather and day-night operational capabilities, with the additional advantage of vegetation penetration capacity for retrieving trunk and primary branch information, thereby enabling vertical forest structure parameter inversion. This capability provides distinct advantages for aboveground carbon storage estimation and partially mitigates the saturation effects observed in other approaches (Qiongfen et al., 2024). Data fusion of active and passive remote sensing modalities has been shown to improve the accuracy of forest aboveground carbon stock estimations (Chunhong et al., 2024). The complementary strengths of SAR and optical sensors suggest that their synergistic integration can enhance the precision of forest aboveground carbon stock assessments (Hao and Huang, 2023).

Currently, a significant number of scholars are utilizing remote sensing data and machine learning algorithms to estimate regional forest biomass and carbon stock. Hidayah et al. (2024) utilized the VV and VH polarization backscatter values from SAR Sentinel 1 imagery, along with field measurements of aboveground biomass, to estimate the carbon storage of the mangrove ecosystem in Benoa Bay, Bali, through a regression model. The results indicate that the VH polarization model (R² = 0.635) outperforms the VV polarization. He Y. et al. (2024) conducted modeling based on the normalized vegetation index model using Sentinel-2A multispectral imagery combined with community structure indicators, achieving a coefficient of determination R² = 0.70. Jiang et al. (2019) utilized Planet Labs imagery as the data source and established a forest volume inversion model using a variance rate-optimized k-NN method, achieving an R² of 0.69 and a root mean square error (RMSE) of 67.6 m³ hm⁻². Han et al., 2024 constructed a multiple linear regression model using GF-2 as the data source, achieving a coefficient of determination R² = 0.73. However, traditional regression methods alone are insufficient to meet the requirements of practical production applications and cannot fully and effectively express the complex nonlinear relationship between remote sensing parameters and forest carbon storage (Chen et al., 2024). He C. R. et al. (2024) utilized Landsat 8 OLI multispectral remote sensing imagery combined with topographic and field survey data. The Boruta algorithm was applied for feature selection, and an XGBoost model was constructed to estimate forest aboveground carbon storage in Beijing with high accuracy (R² = 0.69, RMSE = 5.18 t/hm²). The results showed that the total forest AGC in Beijing in 2014 was approximately 8.93 × 10⁶ t, with a spatial distribution pattern characterized by higher values in the northwest and lower values in the southeast. Luo et al. (2024) integrated Landsat 8 OLI and Sentinel-2A multispectral remote sensing imagery, along with topographic and climatic factors, to develop an aboveground carbon storage estimation model for Pinus densata forests using the Random Forest (RF) algorithm. The optimal accuracy (R² = 0.90, RMSE = 9.53 t/ha) was achieved after incorporating the annual precipitation factor. Rajput et al. (2025) estimated the distribution and aboveground carbon storage of non-timber forest products in Gujarat, India, using high-resolution Sentinel-2 satellite imagery combined with field survey data and geographic information analysis. The RF model demonstrated the best performance in carbon storage estimation (R² > 0.6, RMSE = 24.72 t/hm²). Li J. et al. (2025) developed a forest carbon storage inversion model for Changsha using Landsat-8 remote sensing data, based on dominant tree species classification and the RF algorithm. The model achieved a high R² ranging from 0.9339 to 0.9679, significantly outperforming MLR and support vector regression (SVR) methods.

This investigation concentrates on the Millennium Forest located in Hebei Prov-ince’s Xiong’an New Area, employing a synergistic approach that integrates Sentinel-1 and Sentinel-2. By comprehensively utilizing spectral characteristics, textural features, and polarization signatures derived from the remote sensing data, we implemented and compared four distinct machine learning approaches: Random Forest (RF), XGBoost, Multiple Linear Regression (MLR), and Multi-Layer Perceptron (MLP) for estimating above-ground biomass carbon stocks within the study area. Specifically, for the XGBoost model, the Optuna optimization framework is introduced to adaptively adjust its key hyperparameters. It deeply integrates the characteristics of the machine learning training process and terminates unpromising trials midway through training through an intelligent pruning mechanism. This is particularly suitable for models like XGBoost that require multiple iterations, thereby significantly improving the accuracy and generalization capability of carbon stock estimation models and providing theoretical foundations and methodological support for regional-scale forest carbon sink remote sensing monitoring.

2 Study area and data

2.1 Study area

The Millennium Forest Research Area is located within the Xiong’an New Area of Hebei Province, with geographical coordinates ranging from 116°12′to 116°20′E and latitude from 38°53′to 39°02′N (Figure 1). The terrain is relatively flat and open, with a lower elevation, which is conducive to the planting and growth of trees. The main types of vegetation consist primarily of artificially cultivated trees, including various native species such as the Chinese scholar tree (Sophora japonica), oil pine (Pinus tabuliformis), white-bark pine (Pinus bungeana), as well as introduced species suitable for local growth, forming a relatively rich vegetation com-munity. The carbon sequestration potential of the afforestation project in the Xiong’an New Area is expected to significantly increase after 2030, but it currently consists primarily of young forests, and the carbon stock is still in the early stages of accumulation. It is adjacent to several important functional areas surrounding the Xiong’an New Area. Its unique geographical position plays a critical role in improving the ecological environment of the entire Xiong’an New Area, and it also facilitates the formation of an ecologically coordinated development pattern with the surrounding regions.

2.2 Measured data

In 24 June 2024, the sample plot survey was carried out in the Millennium Forest area of Xiong’an New Area, Hebei Province. A total of 80 square survey plots, each measuring 20 m × 20 m, were established. GPS was utilized to obtain the geographical coordinates of the centers of these plots, allowing for an assessment of the spatial dis-tribution of vegetation within the study area. Tools such as diameter tape, measuring tape, and altimeters were employed to systematically measure data including the di-ameter at breast height and tree height of the trees within the plots. Using the form number method, given the average diameter at breast height, average tree height, and form number of the forest stand, one can first calculate the average volume of a single tree. This value is then multiplied by the total number of trees in the stand to obtain the stand’s volume (Equation 1) (Lin, 2018). Turner et al. (1995), a linear relationship between stand volume and aboveground carbon stocks of forests was established to obtain the above-ground carbon stocks of the forests (Equation 2). The parameters table of the sampled plots is presented (Table 1). V represents the volume in cubic meters, m³; N represents the form factor at breast height; G represents the cross-sectional area at breast height in square m m²; HF represents the tree height factor. C is the carbon storage of a certain forest type, tons of carbon, tC; A is the land area, hm²; V is the volume, m3/hm2; BCEF was the biomass transformation and expansion coefficient, t/m³; R is the rhizome ratio; Cf carbon coefficient. According to the rhizome ratio and carbon coefficients recommended by IPCC (2006), the rhizome ratio was 0.24 when the aboveground biomass was greater than 75 t/hm², 0.39 when the aboveground biomass was less than 75 t/hm², the carbon content coefficient was 0.48 for broad-leaved forest, and 0.51 for coniferous forest.

V = N × G × H F (1)

C = A × V × B C E F × 1 + R × C f (2)

TABLE 1. Statistical table of field survey data.

2.3 Remote sensing data

The imaging data sources mainly include Sentinel-1 and Sentinel-2 data from June 2024. The Sentinel-1 satellite is a Synthetic Aperture Radar (SAR) imaging satellite equipped with C-band synthetic aperture radar, capable of providing continuous im-aging and interferometric measurements. It offers four polarization combinations: HH, VV, HH HV, and VV VH. There are four imaging modes for data acquisition: Strip map mode, Interferometric Wide Swath mode, Extra Wide Swath mode, and Wave mode (Lei et al., 2020; Wang et al., 2024). This article uses Sentinel-1 imagery under the interference broad pattern, with a polarization mode of VV and VH, specifically the product S1A_IW_GRDH_1SDV_20240624T101343_20240624 acquired on 24 June 2024. The Sentinel-2 satellite is equipped with a Multi-Spectral Instrument (MSI) that covers 13 spectral bands, with a swath width of up to 290 km. Its multispectral information is rich, containing 13 bands, with a spectral resolution ranging from 15 to 180 nm, featuring a unique ‘red edge’ band ie, sensitive to vegetation (Lei et al., 2020). Remote sensing features like spectral reflectance, vegetation in-dices, and SAR backscatter coefficients enhance mapping accuracy (Yang et al., 2025). All field survey data and remote sensing data originate from the same period in June 2024, constituting a single-date dataset. Specif-ically, the field survey was conducted from June 15 to 25 June 2024, and the Sentinel-2 image used in this study is the Level-2A product S2B_MSIL2A_20240625T030529_N0510_R075, acquired on 25 June 2024.

The preprocessing of Sentinel-1 data was performed using the SNAP platform following a standardized workflow to mitigate system errors and environmental interference, thereby producing geometrically accurate backscatter coefficient images with clear physical significance. This study selects interferometric wide-swath mode imaging, processed into Level-1 geocoded products, including VV and VH dual-polarization channels with high spatial resolution. First, precise orbit correction was applied, which replaced the initial orbit information with accurate ephemeris data, reducing geolocation errors from tens of meters to the centimeter level and establishing a foundation for subsequent radiometric and geometric processing. Next, thermal noise removal and radiometric calibration were conducted to convert the raw digital numbers into the Sigma0 backscatter coefficient, which quantifies surface scattering characteristics. To suppress the inherent coherent speckle noise in SAR imagery, multi-looking combined with a Refined Lee filter was employed. This approach significantly enhanced the radiometric quality while preserving edge and texture details. Finally, terrain radiometric correction was implemented using the SRTM 1 arc-second digital elevation model. The Range-Doppler terrain correction method was applied to eliminate radiometric distortions and geometric displacements caused by topography. All data were uniformly resampled to a 10-m spatial resolution, projected onto the WGS84 coordinate system, and ultimately output as a backscatter coefficient map in decibels.

Sentinel-2 is equipped with a Multispectral Instrument (MSI) covering 13 spectral bands from the visible to the shortwave infrared range. This study uses Sentinel-2 L2A product-level satellite data as the source of optical remote sensing data. The data were obtained from the official ESA website, with the download link at https://scihub.copernicus.eu/. The L2A level products are atmospheric bottom-of-atmosphere reflectance products generated on the basis of geometrically corrected L1C-level top-of-atmosphere reflectance using specialized processors such as Sen2Cor for atmospheric correction. Its core feature is that it eliminates the interference of atmospheric components such as aerosols and water vapor, directly providing the true reflectance data of surface objects, and by default includes practical scene classification layers such as masks for clouds, cloud shadows, snow, and water bodies. This preprocessing enables the direct use of L2A data for quantitative inversion, allowing for the calculation of key biophysical parameters such as NDVI and LAI without the need for additional atmospheric correction. Therefore, L2A products are the most commonly used and recommended data level for fields requiring quantitative analysis, such as precision agriculture monitoring, forest resource surveys, and land cover change detection. A distinct advantage of this sensor is its inclusion of three red-edge bands, which exhibit high sensitivity to vegetation chlorophyll content, internal leaf structure, and water status. These bands provide critical spectral information for monitoring vegetation physiological conditions. In this study, Level-2A surface reflectance products were directly utilized. These products have been processed by the Sen2Cor algorithm for atmospheric correction, radiometric calibration, and geometric precision correction. They also include auxiliary data layers such as cloud masks, cloud shadow detection, and land cover classification. This preprocessing enables the direct calculation of vegetation indices and the inversion of biophysical parameters, significantly improving both the reliability and computational efficiency of the analysis.

3 Methods

Based on Sentinel-1 and Sentinel-2 remote sensing images as well as field data, the technical route for inverting forest aboveground carbon storage is illustrated (Figure 2). Firstly, polarization data is extracted from Sentinel-1, and spectral information, veg-etation data, and texture features are extracted from Sentinel-2. Then, feature selection is performed using the Pearson and Boruta algorithms. Subsequently, RF, XGBoost, MLR, and MLP models are constructed using Python. Based on the model comparison results, this study adopts a Bayesian optimization approach to automatically tune the hyperparameters of the most accurate XGBoost model, thereby developing a higher-precision forest aboveground carbon stock estimation model and providing reliable technical support for regional-scale forest carbon stock assessment.

3.1 Feature variable extraction and filtering

Due to an excessive number of features leading to data redundancy, which in turn reduces the inversion accuracy of the model, it is necessary to select advantageous feature variables that are beneficial for carbon storage inversion in order to enhance the model’s accuracy. Feature selection can significantly reduce the computation time of the algorithm by eliminating irrelevant and redundant data, thereby preventing the model from becoming overly complex and avoiding the phenomenon of overfitting (Moghimi et al., 2024). The remote sensing factors involved in the feature variable selection of this study in-clude the following: two types of polarization data (VV, VH) from Sentinel-1 imagery, totaling two remote sensing factors; twelve individual band reflectances from Senti-nel-2 imagery along with their corresponding eight second-order texture feature values and seven commonly used vegetation indices; as well as the band combinations b243, b4712, b843, b348, and b83 (Table 2).

TABLE 2. Band combination variables and calculation formulas.

B 2 , B 3 , B 4 , B 8 represent the original bands 2 to 8.

The seven vegetation indices include: Normalized Difference Vegetation Index (NDVI), Green Normalized Difference Vegetation Index (GNDVI), Soil Adjusted Vegetation Index (SAVI), Difference Vegetation Index (DVI), Ratio Vegetation Index (RVI), Enhanced Vegetation Index (EVI), and Optimized Soil Adjusted Vegetation Index (OSAVI) (Table 3).

TABLE 3. Vegetation index and calculation formula.

L is the soil conditioning coefficient, L = 0.5.

In texture analysis, this study employs a typical representative method: the Gray Level Co-occurrence Matrix (GLCM) (Kadiyala and Kumar, 2018), with a window size of 5 × 5 (Lu and Batistella, 2005). The eight second-order texture variables include: mean, variance, correlation, contrast, dissimilarity, homogeneity, angular second moment, and entropy (Table 4).

TABLE 4. Texture eigen index and calculation formula.

i , j: The grayscale value in the gray-level co-occurrence matrix at the i-th row and j-th column. P i , j: The probability of the gray value at (i, j). N: The number of rows in the image. μ: The mean of P. σ: standard deviation.

This study employs the Pearson correlation coefficient method to conduct a correlation analysis between the extracted remote sensing factors and forest aboveground carbon storage. Weakly correlated influencing factors (|r|<0.1) were excluded, resulting in a total of 18 variables identified as candidate factors for modeling (Figure 3). Subsequently, the Boruta algorithm was applied to further refine the selected 18 influencing factors, ultimately yielding 10 influencing factors to be used as modeling factors. Different types of feature variables can have varying impacts on the modeling of above-ground carbon storage in forests.

In order to investigate the impact of different feature variables and their combina-tions on the model accuracy, this study utilized Sentinel-1 data in conjunction with Sentinel-2 remote sensing data, and established six groups of variable combination scenarios (Table 5). By comparing the predicted results of different scenarios, we aim to identify the optimal combination of feature variables that can enhance the accuracy of estimating aboveground carbon storage in forests.

TABLE 5. Combination scenarios of different quantities of features.

3.2 Remote sensing inversion model

In order to explore the application of different types of algorithm models in the inversion of forest above-ground carbon storage, the most suitable algorithm model for carbon storage inversion research in the study area is selected. This study uses 80% of the samples as the training set and 20% of the samples as the testing set to construct four forest AGC models: RF, XGBoost, MLR, and MLP, and to validate the estimation effectiveness of different models for forest aboveground carbon storage in the study area.

To ensure unbiased model evaluation, the entire dataset was first randomly stratified into a training set and an independent test set based on forest type. All feature selection procedures including Pearson correlation analysis for multicollinearity screening and the Boruta algorithm for variable importance ranking were conducted exclusively on the training set. The test set was sequestered from the entire model development process and used only once for the final accuracy assessment. Additionally, cross-validation was applied within the training set to optimize model parameters and prevent overfitting.

All experiments were implemented in Python 3.8. Key packages and their versions included: scikit-learn 1.0.1, XGBoost 1.5.0, Optuna 2.10.0, pandas 1.3.4, numpy 1.21.3, joblib 1.1.0, pyarrow 6.0.0, FastAPI 0.70.0, uvicorn 0.15.0, and loguru 0.5.3.

3.2.1 Random forest

RF employs decision trees as base learners to construct Bagging ensembles, pri-marily encompassing training and classification processes (Fan et al., 2025). The overall computation of RF is formally defined in Equations 3–6. Key advantages include high accuracy, robustness to overfitting, and the ability to handle high-dimensional data with missing values. In this study, RF was implemented using default parameters (n_estimators = 100, max_depth = 5, random_state = 42).

The final prediction of the regression problem is represented by the following formula:

y ^ R F = 1 T ∑ t = 1 T y ^ t (3)

T Is the total number of decision trees in the forest, y ^ t is the prediction of the t-th tree.

Final prediction of the classification problem:

y ^ R F = m o d e y ^ 1 , y ^ 2 , … , y ^ T (4)

m o d e Represents the mode.

Common criteria for node splitting: Within the feature subset F, select the splitting feature f and the splitting point s such that the MSE of the two resulting subsets, left set L and right set R, is reduced the most: M S E p a r e n t is the mean squared error of the parent node, |L| and |R| are the number of samples in the child subsets, N is the total number of samples in the parent node, and ȳ is the mean of the samples within the node. The algorithm selects the split that maximizes the Gain.

G a i n = M S E p a r e n t − L N . M S E L + R N . M S E R (5)

M S E = 1 N ∑ i = 1 N y i − y ¯ 2 (6)

3.2.2 eXtreme gradient boosting

XGBoost is based on the gradient boosting framework, which gradually improves the model’s performance by iteratively constructing decision trees (Katongtung et al., 2026). The mathematical formulation of XGBoost is given in Equations 7–11. It incorporates regularization terms to prevent overfitting and supports parallel computation. In this study, baseline XGBoost used default parameters (num_folds = 4, num_trials = 100 learning_rate = 0.2, max_depth = 6), while the optimized version used parameters derived from Optuna tuning.

In the additive model, the predictive function assumes there are K trees, and the prediction for sample i is the sum of the outputs of all the trees:

y ^ i = ∑ k = 1 K f k x i , f k ϵ F (7)

f k Is the k-th decision tree, and F is the function space of all possible CART trees.

The regularization objective function consists of two parts: the loss function L and the regularization term Ω.Normalization is the item Ω f = γ T + 1 2 λ ω 2。 T:Number of leaf nodes in the tree. ω:the weight vector of the leaf nodes. γ and λ:hyperparameters that control the model complexity.

O b j θ = ∑ i = 1 n L y i , y ^ i + ∑ k = 1 K Ω f k (8)

Optimization goal for the t-round iteration: Assuming that in the t-round iteration, we are training the f t of the t-tree, the predicted value is y ^ i t = y ^ i t − 1 + f t x i. The objective function becomes:

O b j t = ∑ i = 1 n L y i , y ^ i t − 1 + f t x i + Ω f t + C (9)

Second-order Taylor expansion approximation: Perform a second-order Taylor expansion of the loss function at y ^ i t − 1:

O b j t ≈ ∑ i = 1 n L y i , y ^ i t − 1 + g i f t x i + 1 2 h i f t 2 x i + Ω f t (10)

Splitting Gain Criterion: When constructing a tree, it is necessary to select the optimal feature and split point. By comparing the change in the target value Obj before and after the split, the gain is defined as:

G a i n = 1 2 G L 2 H L + λ + G R 2 H R + λ − G L + G R 2 H L + H R + λ − γ (11)

L and R represent the left and right child nodes after the split, respectively. Choose the split scheme that maximizes Gain. If the maximum Gain is negative, stop splitting.

3.2.3 Multiple linear regression

MLR models the linear relationship between independent variables and the dependent variable (Attua, 2008) (num_folds = 4). The functional form of this linear relationship is given in Equation 12.

For a dataset with n observations and p independent variables, the model can be expressed as: y i: The dependent variable for the i-th sample. β 0: The intercept term, representing the baseline predicted value when all independent variables are zero. β 1 , β 2 , … , β p:The regression coefficients for each independent variable, indicating the average change in the dependent variable for a one-unit change in the respective independent variable, holding all other variables constant. x i 1 , … , x i p:The values of the independent variables for the i-th sample. ϵ i:The random error term, representing the random fluctuations that the model cannot explain, typically assumed to follow a normal distribution with a mean of zero and constant variance.

y i = β 0 + β 1 x i 1 + β 2 x i 2 + … + β p x i p + ϵ i , i = 1 , 2 , … , n (12)

3.2.4 Multi-layer perceptron

MLP is a feedforward artificial neural network consisting of an input layer, hidden layers, and an output layer (Pujadas et al., 2025). It uses backpropagation to minimize prediction errors and can approximate complex nonlinear functions. In this study, MLP was configured with two hidden layers, ReLU activation, and a maximum of 200 training iterations.

3.2.5 XGBoost model optimization

Bayesian optimization is a global optimization method based on probabilistic surrogate models, which approximates the optimal configuration with relatively few iterations by mapping hyperparameters to model performance. Optuna is an efficient engineering implementation of this theoretical framework, specifically designed for machine learning hyperparameter optimization scenarios. Compared with traditional algorithmic-level implementations of Bayesian optimization, Optuna enhances practicality in several key aspects: its intelligent pruning mechanism can terminate unpromising trials early during training, which is particularly suitable for iterative models such as XGBoost; its Define-by-Run approach enables dynamic definition of parameter spaces, facilitating the management of complex conditional dependencies among hyperparameters; furthermore, its asynchronous parallel architecture and visualization diagnostic tools transform hyperparameter optimization from a black-box process into an interpretable and manageable engineering practice, providing a more convenient and reliable solution for model tuning.

Bayesian optimization is an intelligent hyperparameter tuning method based on probabilistic models (Haixia et al., 2024). Its core idea is to gradually approach the optimal solution by establishing a probabilistic mapping relationship between the objective function and the hyperparameters. The mathematical framework of Bayesian optimization, including the surrogate model and acquisition function, is formally defined in Equations 13–18. The Bayesian optimization algorithm can identify the optimal hyperparameter combination for XGBoost with a limited number of evaluations, thereby enhancing the model’s predictive accuracy (Figure 4).

Assuming the set of hyperparameters is θ, the training objective of XGBoost given θ is to minimize its regularized objective function:

L t θ = ∑ i = 1 n l y i , y ^ i t − 1 + f t x i ; θ + Ω f t ; θ (13)

Ω includes a complexity penalty term controlled by θ. The average performance of the model under k-fold cross-validation is the objective that Bayesian optimization seeks to optimize. Θ represents the search space of hyperparameters.

P e r f θ = − 1 k ∑ f o l d = 1 k L v a l i d f o l d θ (14)

θ * = arg max θ ∈ Θ P e r f θ (15)

The core of Bayesian optimization is addressed through two key components: the surrogate model Gaussian process, which assumes that the objective function f(θ) follows a Gaussian process prior: m θ: Mean function, usually set as a constant. k θ , θ ′:A covariance function (kernel function) that measures the similarity of the hyperparameter θ to θ^'. Matten kernel or square exponential kernel is commonly used.

f θ ∼ GP m θ , k θ , θ ′ (16)

The acquisition function α(θ) is used to balance exploration and exploitation. Expected improvement is the most commonly used acquisition function. Let the best function value observed so far be defined as f * = max ⁡ _ i ≤ t f θ _ i。The improvement at a point θ is given by: I θ = max ⁡ θ , f θ − f *。Since f(θ) is a random variable, the expected improvement is the expected value of the improvement I(θ):

E I θ = E I θ = μ t θ − f * − δ Φ Z + σ t θ ϕ Z , i f σ t θ > 0 0 , i f σ t θ = 0 (17)

The next step in Bayesian optimization is to select the point θ t + 1 that maximizes the acquisition function E I θ as the next evaluation point.

θ t + 1 = arg max θ ∈ Θ E I θ (18)

To systematically evaluate and enhance the model’s stability and generalization capability, this study employed a robust validation strategy during hyperparameter optimization. Specifically, we conducted 100 independent hyperparameter search trials using the Optuna framework. The implementation of this objective function within the Optuna framework, including the trial definition and cross-validation aggregation, is formally expressed in Equations 19–21. In each trial, a 4-fold cross-validation was applied to assess candidate parameter sets, with the objective function being the mean R² across the 4-fold validation (n_estimators = 100, max_depth = 6, learning_rate = 0.08, subsample = 0.85, colsample_bytree = 0.9, gamma = 0.1)θ represents the set of hyperparameters to be optimized; Ω(θ) is the early stopping penalty term in Optuna.

L * θ = 1 N ∑ i = 1 N L y i , y ^ i + Ω θ (19)

Optuna’s TPE Sampling Strategy:

p θ | L = p L | θ p θ p L (20)

Construct two density distributions by dividing the experiments into superior and inferior categories: p θ | L < τ , p θ | L > τ.

Select the parameter with the maximum Expected Improvement:

θ n e w = arg max θ E I θ = arg max θ ∫ − ∞ τ τ − L p L | θ d L (21)

3.3 Model evaluation and accuracy testing

To assess the accuracy of the inversion of forest aboveground carbon storage, this study intends to use the coefficient of determination R², root mean square error (RMSE), and mean absolute error (MAE) to validate the model’s precision. R², also known as the coefficient of determination, is used to quantify the degree of fit of a model to the data. It measures the strength of the relationship between the independent variables and the dependent variable in the model, reflecting the proportion of variance explained by the model. Its value ranges from 0 to 1, with a value closer to 1 indicating a better fit of the model, meaning that the independent variables can explain a greater proportion of the variance in the dependent variable. RMSE is used to measure the deviation between predicted values and actual values, with its units consistent with the original data, making it easier to understand the magnitude of the model’s prediction errors. A smaller RMSE value indicates higher predictive accuracy of the model. MAE, or Mean Absolute Error, is a metric used to evaluate the performance of predictive models. It assesses the accuracy of the model by calculating the average of the absolute errors between predicted values and actual values. The mathematical formulations of R², RMSE, and MAE are given in Equations 22, 23. y i, y ¯ i, y ^ i are the measured aboveground carbon storage values, the average measured aboveground carbon storage values, and the model-predicted aboveground carbon storage values, respectively and n is the verification sample size.

R 2 = 1 − ∑ i = 1 n y i − y ^ i 2 ∑ i = 1 n y i − y ¯ i 2 (22)

R M S E = ∑ i = 1 n y i − y ^ i 2 n (23)

M A E = 1 n ∑ i = 1 n y i − y ^ i (24)

4 Results

4.1 Optimal feature variables

The core logic of the classification scenario design is to systematically control input variables to isolate the contribution of different feature types to the model’s predictive accuracy. Six inversion scenarios were designed, encompassing four categories of variables: spectral bands, texture features, vegetation indices, and backscatter coefficients, along with their various combinations (Figure 5). This systematic approach was employed to investigate the influence of feature combinations on the estimation of forest above-ground carbon storage and to evaluate the performance of the XGBoost model in this task. The results indicate that variable selection and combination strategy is a key determinant of model accuracy, with the coefficient of determination R² ranging from 0.13 to 0.60 across different scenarios, demonstrating significant performance disparities (Table 6).

TABLE 6. Performance of XGBoost models with different numbers of features.

Scenario d, which delivered optimal performance, is based on ten key variables selected through feature optimization. This set includes Sentinel-2 spectral bands, the NDVI and SAVI vegetation indices, five texture features, and the Sentinel-1 VH-polarized backscatter coefficient. This combination achieved the highest accuracy on the test set, with an R² of 0.60. This suggests that a moderate, optimized feature set effectively balances the retention of valid information with the avoidance of overfitting, thereby yielding the best generalization capability. These multi-source, multi-dimensional features collectively elucidate the underlying mechanisms driving forest carbon storage, significantly enhancing the model’s explanatory and predictive power.

In contrast, Scenario b, which utilized only radar backscatter coefficients, performed the worst, with a test set R² of only 0.13. This indicates that in dense forest areas, microwave signals are susceptible to interference from factors such as moisture and surface roughness, making them insufficient for high-precision inversion when used in isolation. Scenario c, which relied solely on Sentinel-2 optical and texture features, also yielded limited accuracy. This limitation underscores the model’s inability to fully capture complex biophysical processes without the support of multi-dimensional information. Furthermore, the scenario involving an unfiltered combination of all variables performed worse than the optimized combination Scenario d, implying that feature redundancy may introduce multicollinearity and noise, ultimately degrading generalization performance.

Scenario a, representing a preliminary combination of original multi-source data, achieved moderate accuracy test set R² = 0.52. This result demonstrates that even without refined feature selection, the integration of key spectral and vegetation information can establish a model with reasonable explanatory power. However, its accuracy remains notably lower than that of Scenario d, highlighting the critical role of feature optimization in extracting complementary information and mitigating saturation effects.

Through multi-scenario comparison, this study demonstrates that while the XGBoost model holds potential for forest carbon stock estimation, its performance is highly contingent upon the physical and ecological relevance of the input features and their compatibility with regional conditions. The optimal results were achieved using a diversified feature system that integrates spectral, textural, and microwave features. Intelligent feature selection within this system helped counteract some of the interference caused by surface heterogeneity, thereby enhancing the model’s stability and interpretability in complex environments.

It is important to note that the accuracy of remote sensing inversion is also profoundly influenced by regional environmental heterogeneity. Variations in climate, vegetation type, and soil conditions across different study areas directly affect the quality of data acquisition and spectral response characteristics. Additionally, terrain undulations lead to spatial variations in surface illumination and vegetation growth patterns, which alter reflectance properties and introduce uncertainty, thereby reducing inversion accuracy. Consequently, when constructing a feature system, it is imperative not only to use intelligent selection to integrate multi-dimensional information but also to fully account for the influence of region-specific environmental factors on feature representativeness.

4.2 Model evaluation of AGC

The random forest model demonstrated better generalization characteristics, achieving a test set R² of 0.42 and exhibiting a narrower gap between its training and testing performance compared to other baseline methods. In contrast, the XGBoost algorithm achieved a substantial performance leap, with the test set R² increasing to 0.60—an improvement of 0.18 over random forest (Figure 6). This enhancement stems from the synergistic effect of its gradient boosting framework and integrated regularization mechanisms: the former sequentially fits residuals to enhance the model’s expressive power, while the latter effectively controls model complexity through structural constraints (Table 7).

TABLE 7. Comparison results of six models.

Following Bayesian optimization of its hyperparameters, the XGBoost model’s performance was further improved, achieving a test set R² of 0.65 and RMSE of 27.63 Mg C/ha. The optimization process constructed a probabilistic mapping between the hyperparameter space and model performance, enabling an efficient search through a balanced strategy of exploration and exploitation. Notably, the performance disparity between the training and test sets was further minimized after optimization, indicating that hyperparameter tuning not only improves predictive accuracy but also enhances the model’s generalization capability.

The adoption of the Optuna framework for optimization yielded a more substantial breakthrough. The Optuna-tuned XGBoost model attained a test set R² of 0.72, representing an improvement of 0.12 over the baseline XGBoost model and 0.07 over the Bayesian-optimized version. These gains are attributed to three key innovative mechanisms within the Optuna framework: (1) dynamic pruning, which early terminates trials with low potential; (2) an asynchronous parallel architecture, which improves computational resource utilization; and (3) the TPE, which optimizes sampling efficiency in mixed hyperparameter spaces.

Regarding error metrics, multiple linear regression showed consistent performance in both RMSE and mean absolute error (MAE), aligning with its least squares optimization objective. The advantage of random forest in MAE reflects the inherent robustness of tree-based ensemble methods against outliers. Conversely, the XGBoost series of models led across all error metrics, confirming the strong adaptability of the gradient boosting framework for remote sensing inversion tasks.

From an applied science perspective, this performance evolution carries clear ecological significance. The increase in the coefficient of determination (R²) from 0.38 to 0.72 indicates that the model’s capacity to explain the spatial variability of carbon stocks nearly doubled. Concurrently, the RMSE decreased from 44.69 to 24.48 Mg C/ha, signifying a substantial enhancement in per-hectare carbon stock estimation accuracy. At a regional scale, such improvements in accuracy can be directly translated into a more reliable scientific foundation for policy-making and ecological management.

The results of this study indicate that the integration of ensemble learning and intelligent optimization techniques provides a robust methodological framework for the remote sensing inversion of forest structural parameters. The Optuna-enhanced XGBoost model achieves an optimal balance of accuracy, efficiency, and stability under current technological conditions, with a test set R² of 0.72 and an RMSE of 24.48 Mg C/ha, thereby establishing a reliable technical benchmark for small-scale forest carbon monitoring. This advancement not only confirms the significant potential of machine learning in ecological remote sensing but also demonstrates a viable pathway to addressing complex geoscientific problems through algorithmic innovation. With the continuous advancement of Earth observation technologies and computational methods, such intelligent inversion models are expected to propel carbon cycle research toward a more quantitative and refined stage across larger spatial and temporal scales.

4.3 Spatial distribution of carbon storage in forest soils

After the optimization and screening using the Pearson and Boruta algorithms, ten critical feature factors were selected as variables for the OP-XGBoost model to invert the forest aboveground carbon storage for 2024. The performance comparison confirms that the OP-enhanced XGBoost model effectively improves estimation accuracy. Furthermore, an inversion map of forest aboveground carbon storage for 2024 was created using ArcGIS. As illustrated in Figure 7, the minimum estimated aboveground carbon storage in the study area is 32 Mg C/ha, while the maximum estimated value reaches 284 Mg C/ha. Notably, the spatial distribution pattern of carbon storage shows higher values in the middle to upper zones, exhibiting significant variability.

5 Discussion and conclusion

5.1 Discussion and outlook

This study presents an exploratory approach that aims to estimate the AGC in the Millennial Forest of the Xiong’an New Area. The feasibility of a model integrating Sentinel-1 and Sentinel-2 remote sensing data with an Optuna-optimized XGBoost algorithm was investigated. Experimental results demonstrate that the proposed model achieves high accuracy and stability. The study region is characterized by a warm temperate semi-humid continental monsoon climate and flat terrain. The forests consist primarily of intensively managed, relatively young artificial plantations, dominated by species such as poplar and ash, exhibiting a uniform stand structure. This artificial afforestation pattern is typical of plain regions and is widely representative of agricultural areas in North and Northeast China, as well as new urban greening projects. Consequently, the model developed in this study serves as a direct reference for carbon monitoring in other managed plain forests under similar climatic conditions.

Our model provides an effective tool for dynamic forest carbon monitoring, supporting carbon sink assessment and management while offering a low-cost, high-efficiency solution for similar plain afforestation regions. Nevertheless, several limitations exist: scale mismatch between plots and pixels may introduce aggregation errors, particularly in heterogeneous stands, and model generalizability varies between mixed conifer-broadleaf and pure plantations due to differential responses of remote sensing features to carbon stocks. We recommend stratified sampling and locally calibrated models for specific forest types. Future work should integrate UAV LiDAR for three-dimensional structure and incorporate time-series features to capture phenological effects on carbon accumulation.

Contribution of different variable groups to AGC estimation. Based on the scenario comparisons in Section 4.1, the optimal feature combination (Scenario d) included spectral bands (B3, B4, B8), vegetation indices (NDVI, SAVI), texture features (Correlation, Entropy, Homogeneity, Mean, Second Moment), and the Sentinel-1 VH backscatter coefficient. Among these, the spectral bands and vegetation indices provided the foundational spectral information for AGC estimation, as they directly capture vegetation greenness and canopy coverage. The texture features contributed substantially by characterizing canopy roughness, uniformity, and spatial heterogeneity—information ie, not available from spectral bands alone. This is evidenced by the performance improvement from Scenario a (R² = 0.52, without texture) to Scenario d (R² = 0.60, with texture). The VH polarization from Sentinel-1 provided complementary information on vertical canopy structure and moisture content, although its independent contribution was limited (Scenario b achieved only R² = 0.13). Collectively, these results indicate that the synergistic integration of multi-source features—rather than any single variable group—was the key driver of model performance. Future work should use SHAP analysis or ablation studies to quantify the exact contribution of each variable.

From a remote sensing mechanism perspective, the improved performance is attributed to the synergistic and complementary effects of multi-source features. Sentinel-1 C-band backscatter is sensitive to forest vertical structure but tends to saturate under high biomass conditions (Forkuor et al., 2020). In contrast, Sentinel-2 red-edge bands are strong indicators of chlorophyll content. The synergistic use of optical and SAR data has been widely demonstrated to improve AGB estimation accuracy, with recent studies showing that XGBoost outperforms traditional models in such tasks (Jia et al., 2024). Specifically, Zhang et al. (2024) achieved an R² of 0.814 for AGC estimation in coniferous forests by combining optical and SAR time-series data, while Ozdemir and Abdikan (Ozdemir and Abdikan, 2025) reported an R² of 0.74 using multi-frequency SAR with optical imagery. In our study area of planted young forests, the Optuna-optimized XGBoost achieved an R² of 0.72, which is comparable to these benchmarks considering differences in forest type and sample size.

Nevertheless, several uncertainties and limitations in this study warrant further discussion. First, as noted in recent systematic reviews (Ullah et al., 2025), the scale mismatch between field plots and remote sensing pixels remains a common challenge in regional AGC estimation. Second, while our feature selection approach (Pearson + Boruta) improved model performance, recent studies have adopted more rigorous methods such as SHAP analysis for variable importance quantification, approaches that should be considered in future work (Xu et al., 2025). Third, the generalizability of our model to other forest types or regions requires further validation, as Zheng et al. (2025) demonstrated that model performance varies considerably across different ecological zones.

Limitations regarding carbon stock calculation methods. This study employed a standardized volume-to-carbon conversion approach based on the form number method and IPCC default coefficients (BCEF, root-to-shoot ratio, and carbon fraction), rather than species-specific allometric equations. While this approach is widely used for regional-scale forest carbon estimation and is appropriate for our study area—dominated by young planted pure forests with relatively uniform stand structure—we acknowledge its limitations. In mixed-species or structurally complex natural forests, species-specific equations would likely yield more accurate estimates. Future studies should develop locally calibrated allometric equations for dominant tree species to reduce estimation bias.

Limitations regarding feature contribution analysis. In this study, we evaluated feature importance through six predefined variable combinations rather than conducting an ablation study on the final selected feature set. While our scenario-based approach effectively identified the relative importance of different feature types (e.g., spectral, texture, SAR), it does not precisely quantify the contribution of each individual variable. A more rigorous ablation study or SHAP analysis would provide finer-grained insights into variable-level contributions. Future research should adopt such methods to better interpret model mechanisms.

The improved performance from multi-source feature fusion can be explained by sensor physics and their response to forest structure. Sentinel-1 C-band VH cross-polarization is sensitive to planted forest structure—low crown base height, continuous canopy, and distinct vertical stratification enable VH signals to capture volume scattering from trunks and branches, establishing a strong correlation with aboveground carbon stocks, though radar signals tend to saturate under high-biomass conditions. Sentinel-2 optical data provide complementary physiological information. NDVI and SAVI perform prominently in young forests, where the partially open canopy makes soil background influence significant. SAVI reduces this interference via a soil adjustment factor, while NDVI accurately reflects young tree growth vigor through its sensitivity to absorbed photosynthetically active radiation. As stands mature, red-edge bands—indicating chlorophyll content and canopy structure—become more critical, partially compensating for radar saturation. Texture features derived from the Gray-Level Co-occurrence Matrix quantify canopy spatial heterogeneity, indirectly capturing variations in crown size, shape, and arrangement through spatial dependencies of pixel gray values. The synergy lies in complementary integration: the Optuna framework allocates feature weights so that radar structure information dominates under moderate biomass, while optical physiological and textural information supplements under high-biomass or complex conditions, jointly mitigating saturation and underestimation from single data sources.

Sentinel-1 and Sentinel-2 offer significant complementary advantages for forest carbon monitoring, yet challenges remain including data interference, limited cross-ecosystem reliability, and difficulties in long-term prediction. Future efforts should focus on deep integration of multi-source data from hyperspectral, LiDAR, and UAV platforms, optimized by advanced AI, alongside standardized observation networks and international collaboration to enhance monitoring accuracy in support of global carbon neutrality goals.

The trade-off between model complexity and interpretability remains a pressing concern. The shift from parameter-transparent models like MLR to powerful but opaque ensemble learning methods often improves predictive performance at the cost of mechanistic interpretability. Future research should prioritize developing interpretability frameworks tailored for ecological applications. Techniques such as feature importance analysis and local dependency decomposition can elucidate the physical links between remote sensing features and forest carbon cycle processes. Promoting the integration of explainable AI (XAI) with remote sensing science is a critical direction for quantitative remote sensing.

It is noteworthy that climate change exerts a significant influence on forest growth and the carbon cycle, with temperature and humidity serving as key factors regulating tree growth (Guo et al., 2024a). Therefore, in studies that use remote sensing and machine learning techniques to invert aboveground forest carbon storage, the incorporation of key climatic drivers such as temperature and precipitation is expected to markedly enhance the model’s explanatory power and predictive capability. Guo pointed out that climatic variables play an important role in vegetation dynamics and ecosystem processes, thereby regulating carbon accumulation and turnover (Guo et al., 2024b). Similarly, He found that during the rapid urbanization stage, natural factors are the primary drivers of ecological and environmental changes; whereas in the stage of urbanization quality improvement, the interaction between socioeconomic factors and natural factors becomes significantly stronger (He et al., 2025). This further suggests that, in constructing carbon storage models, in addition to remote sensing and climatic variables, the comprehensive consideration of socioeconomic and other anthropogenic driving factors may have positive implications for improving the model’s applicability within complex human–natural coupled systems.

5.2 Conclusion

Taking the Millennium Forest in Xiong’an New Area, Hebei Province as the study area, this research utilized Sentinel-1 microwave and Sentinel-2 optical data to estimate forest aboveground carbon stocks. Optimal predictive variables were selected by combining Pearson correlation analysis with the Boruta algorithm, and multiple machine learning models were developed and evaluated. Bayesian optimization was applied to refine the hyperparameters of the XGBoost model. Based on the optimized model, a spatial distribution map of aboveground carbon stocks for 2024 was generated. The main conclusions are as follows:

• Pearson correlation analysis combined with the Boruta algorithm effectively identified the optimal feature combination by eliminating multicollinearity and reducing information redundancy. Among the original 34 features, 10 core features were selected, resulting in an increase in R² of 0.18, a decrease in RMSE of 8.96 Mg C/ha, and a reduction in MAE of 7.9%. The selected Sentinel-1 VH polarization features demonstrated sensitivity to the three-dimensional structural parameters of planted forests, while Sentinel-2 vegetation indices, spectral bands, and texture features also contributed substantially to carbon stock inversion.
• Among the evaluated models, the OP-XGBoost model optimized under the Optuna framework achieved the highest performance, with a test set R² of 0.72, RMSE of 24.48 Mg C/ha, and MAE of 18.56%, significantly outperforming random forest, multilayer perceptron, multiple linear regression, and the unoptimized XGBoost model. This model achieved a favorable balance between fitting capability and generalization performance.
• The model developed in this study is applicable to planted timber forests and ecological public-welfare forests in the North China Plain, demonstrating stable performance in even-aged pure plantations or structurally homogeneous stands with moderate canopy density (0.4–0.7) and stand ages of 5–15 years. Several limitations should be acknowledged. The model was established based on single-temporal remote sensing data (June 2024), without incorporating time-series features, which limits its ability to capture seasonal dynamics and interannual trends in forest carbon stocks. The study area is predominantly composed of planted pure forests, and the applicability of the model to natural forests, mixed conifer-broadleaf forests, or structurally complex uneven-aged stands remains to be validated. Field plots covered only a portion of forest types, and the limited sample size may affect estimation stability under extreme conditions. Future research should extend to broader forest types and geographical regions, and incorporate time-series remote sensing data as well as LiDAR-derived three-dimensional structural information to further enhance model generalizability and estimation accuracy.

The comparison between training and testing accuracy across all scenarios suggests no significant overfitting. Scenarios a and c exhibited the smallest R² differences between training and test sets (both 0.01), indicating the best model stability. Scenario d showed the largest R² difference (0.15), yet still achieved the highest test R² (0.60) and the lowest RMSE (31.76 Mg C/ha) among all scenarios, confirming that this parsimonious feature combination, despite slight generalization variability, delivered the best overall predictive performance. Scenario e exhibited an R² difference of 0.17, potentially attributable to over-adaptation caused by excessive variable dimensionality (34 features). Scenario b performed poorly on both training and test sets (R² ≤ 0.16), corroborating that radar-only data are inadequate for forest carbon stock estimation. Collectively, the training-test accuracy discrepancies across all models remain within acceptable limits, affirming reliable model generalization.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Weijing M. A., Lanzhou University, China

Reviewed by: Guangping Qie, Southern Illinois University Carbondale, United States

Cenliang Zhao, Sun Yat-sen University, China