Estimation of forest carbon storage with small samples based on Bayesian hierarchical model

1 Introduction

According to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, the global climate is experiencing a clear warming trend, which has become one of the major environmental challenges worldwide. Forest ecosystems are an important part of terrestrial ecosystems and play a vital role in regulating the global carbon cycle and mitigating climate change through carbon sequestration (Li et al., 2025). Improving forest carbon sequestration capacity is considered an effective approach for achieving carbon neutrality and addressing climate change. Therefore, accurate estimation of forest carbon storage and carbon sequestration potential is essential to support international climate governance and national carbon peaking and carbon neutrality goals. Against the backdrop of ongoing climate change, the assessment of forest carbon storage and sequestration potential has become an important research topic in ecology and climate change.

Field survey data serve as a critical source for estimating forest biomass and carbon storage. When samples are evenly distributed across space, the sample size is generally sufficient to support reliable prediction of aboveground biomass and carbon storage for a given tree species (Frank et al., 2019). However, acquiring field data in forest ecosystems remains challenging in practice. Factors such as large forest distribution areas, complex terrain, limited accessibility, and difficult survey conditions often lead to uneven spatial distribution of samples and insufficient sample sizes. This problem is particularly evident in remote mountainous regions and densely forested areas, where sample scarcity is more common. Under these conditions, traditional statistical models that rely on large-sample assumptions may show limited applicability and reduced reliability in forest carbon stock estimation.

Small area estimation (SAE) is a statistical method that generates reliable estimates for regions or domains with insufficient or even zero sample data, while enabling quantitative evaluation of estimation accuracy (Pfeffermann, 2013). In practical applications, conventional field survey- and empirical model-based approaches are frequently constrained by limited sample sizes, thereby lowering the accuracy and reliability of forest carbon estimates. By incorporating auxiliary information, SAE can improve estimation performance under sample-scarce conditions. SAE methods are generally divided into design-based and model-based approaches. The latter are commonly developed within either the frequentist or Bayesian framework, and in some cases combine both, which is referred to as empirical Bayes in SAE studies (Pfeffermann, 2013). As a classic empirical Bayes technique, the Bayesian Hierarchical Model integrates Bayesian inference with classical statistical modeling and produces more robust and reliable estimates than conventional survey-based methods, especially under small-sample conditions.

Although machine learning (ML) algorithms can capture complex nonlinear relationships between predictors and response variables and have been widely applied in forest biomass and carbon stock estimation, their performance largely depends on sufficient representative training samples. Limited sample sizes render ML models susceptible to overfitting, which impairs model stability and prediction reliability (Jiang X. et al., 2020; Lu and Jiang, 2024). In addition, traditional regression and mixed-effects models are generally based on deterministic parameter estimation frameworks, making it difficult to adequately characterize the uncertainty associated with model parameters and predictions (Zhang et al., 2014). To address the instability of model estimation and the difficulty in uncertainty quantification under limited sample conditions, Bayesian hierarchical models provide an effective probabilistic modeling framework. Due to its excellent applicability under limited sample conditions and its effective ability to quantify model uncertainties, it has received extensive attention in recent years (Junttila and Laine, 2017; Ver Planck et al., 2018). By constructing hierarchical structures, this method can effectively estimate complex parameter relationships and achieve good model performance with relatively small training datasets and highly correlated multidimensional data (Junttila et al., 2015). Meanwhile, Bayesian approaches can characterize parameter uncertainty through posterior probability distributions and directly account for uncertainty in small area estimation rather than providing only point estimates (Ver Planck et al., 2018). Furthermore, the Bayesian framework allows prior information, such as historical data, expert knowledge, or physical models, to be incorporated into the probabilistic inference process and combined with new observations for parameter updating, thereby providing a more comprehensive uncertainty assessment for forest carbon stock estimation under complex environmental conditions (Jia et al., 2025).

The Three Parallel Rivers Region refers to the area where the Jinsha, Lancang, and Nujiang Rivers flow in parallel (Song et al., 2025). Located at the core of the Hengduan Mountains, this region features highly complex alpine terrain with steep slopes, deep valleys, fragmented landscapes, and substantial elevation gradients. It exhibits a subtropical monsoon climate, with cold and humid conditions at high elevations and warm, relatively dry conditions in river valleys. These complex topographic and climatic conditions limit field accessibility, resulting in sparse and unevenly distributed survey samples that cannot fully represent the study area. Based on this context, this study uses the 2022 continuous forest inventory data of Yunnan Province and multi-source remote sensing data to construct and compare three models, including the Bayesian Hierarchical Model, Ordinary Least Squares model, and Random Forest model. The objectives are to evaluate the performance of the Bayesian Hierarchical Model under small-sample conditions, improve the accuracy of forest carbon stock estimation, and identify the most suitable remote sensing data source. Finally, the optimal model is applied to map forest carbon stocks and analyze their spatial patterns, with further exploration of the roles of topography and vegetation in shaping carbon distribution.

2 Study area

The Three Parallel Rivers Region is located in northwestern Yunnan Province, China, covering Lijiang City, Diqing Tibetan Autonomous Prefecture, and Nujiang Lisu Autonomous Prefecture. As a core region of the Hengduan Mountains, this area is characterized by complex mountainous terrain dominated by high mountains and deep valleys, with strong topographic relief and highly heterogeneous landscapes. The elevation varies greatly across the region, ranging from approximately 1,000 m in river valleys to 6,740 m at Kawagebo Peak of Meili Snow Mountain. The steep terrain, deeply incised valleys, and large elevation differences have resulted in pronounced vertical environmental gradients, with some areas having slopes exceeding 30°. The region has a subtropical monsoon climate, and the substantial elevation gradients generate distinct variations in temperature and moisture conditions, which support diverse forest ecosystems. The main forest types considered in this study include coniferous forests, broad-leaved forests, and mixed forests. Coniferous forests are mainly composed of Pinus yunnanensis, Pinus densata, Picea spp., and Abies spp. Among them, Picea and Abies are mainly distributed in cold and humid high-elevation areas, whereas Pinus yunnanensis commonly occurs in relatively dry and sunny mid–low elevation mountains. Broad-leaved forests are mainly composed of Quercus spp., Cyclobalanopsis Oerst., and other evergreen broad-leaved species, occurring across diverse habitats from warm valleys to mountainous areas. Mixed forests consist of co-occurring coniferous and broad-leaved species and are mainly distributed in transitional zones along elevation gradients. Human activities in the region are mainly concentrated in accessible valleys and low-elevation areas, including settlement expansion, agricultural activities, and infrastructure development. In contrast, higher-elevation mountainous forests are generally less affected by direct human disturbance. The location of the study area is shown in Figure 1.

3 Materials and methods

3.1 Data source

3.1.1 Remote sensing data

Considering the different characteristics of these satellite sensors, Landsat 8 OLI and Sentinel-2 data were integrated to obtain complementary information for forest condition characterization. Landsat 8 OLI provides continuous observations and reliable spectral information, particularly in the shortwave infrared bands, which are sensitive to vegetation moisture and canopy structural characteristics. Nevertheless, its relatively coarse spatial resolution and absence of red-edge bands limit its ability to capture fine-scale forest structural variations. In contrast, Sentinel-2 provides finer spatial resolution (10–20 m) and additional red-edge bands, which are more sensitive to vegetation structure and biochemical properties. Furthermore, the shorter revisit interval of Sentinel-2 complements the temporal continuity of Landsat observations. However, Sentinel-2 suffers from a shorter operational archive compared with the long-term Landsat record. Therefore, integrating Landsat 8 OLI and Sentinel-2 data enables a more comprehensive characterization of forest structural variability and vegetation conditions, which partly explains the superior performance of the combined dataset for forest carbon stock estimation (Claverie et al., 2018; Roy et al., 2018).

All remote sensing data were obtained from the Google Earth Engine platform, including Landsat 8 OLI, Sentinel-2 imagery, and digital elevation model (DEM) data. A total of 95 Landsat 8 OLI images covering the Three Parallel Rivers Region from June to December 2022 were used. Cloud and cloud-shadow contamination were removed before analysis, and a median composite was generated to produce a representative image. Six spectral bands (Bands 2–7) were selected for analysis, with a spatial resolution of 30 m.

For Sentinel-2 data, all images covering the study area from June to December 2022 were selected. Cloud and cloud-shadow effects were removed, and a median composite was generated, resulting in 12 spectral bands. The Sentinel-2 data were then resampled to a unified spatial resolution of 30 m and integrated with other datasets. The digital elevation model (DEM) was obtained from the 30 m SRTM DEM dataset. Based on the DEM, terrain analysis was conducted to derive slope and aspect variables for the study area.

3.1.2 Climatic data

The monthly precipitation data (Peng et al., 2017) (1901–2024) and the monthly average temperature data (Peng et al., 2018) (1901–2024) obtained from the National Qinghai-Tibet Plateau Data Center, both with a spatial resolution of 0.0083333° (approximately 1 km), have monthly precipitation units of 0.1 mm and monthly average temperature units of 0.1 °C. Based on the data, we processed the data by creating NetCDF raster layers and pixel statistics in ArcMap platform, calculated the annual average temperature and precipitation, and converted the data units to °C and mm.

3.1.3 Spatial distribution data of vegetation types

The vegetation type data were obtained from the 30 m global fine land-cover dynamic monitoring product provided by the Data Sharing Service System for Earth Big Data Science Engineering, Chinese Academy of Sciences (CAS)¹. This dataset covers the period from 1985 to 2022 with a spatial resolution of 30 m (Zhang et al., 2021). In this study, the vegetation type data for 2022 were selected for analysis. The downloaded vegetation data were processed and extracted in ArcGIS software and subsequently reclassified into three categories: coniferous forest, broadleaved forest, and mixed forest. Before extracting raster values for sample points, all datasets were transformed into a unified coordinate system (WGS 84, EPSG: 4326), and the nearest-neighbor interpolation method was applied to resample the datasets to a spatial resolution of 30 m. This procedure was conducted to ensure consistency in spatial resolution and to minimize error propagation during data registration and integration (Robinson, 2009).

3.1.4 Factor extraction and selection

Spectral information is a crucial source for identifying ground objects in remote sensing imagery and can reflect the basic characteristics of surface targets. Vegetation indices are constructed through mathematical operations on different spectral bands, are widely used to characterize vegetation growth conditions and estimate forest carbon stocks. Texture features of remote sensing imagery can further describe the spatial variation of pixel gray levels and reflect differences in vegetation structure. In this study, texture features were extracted in ENVI software using the Gray-Level Co-occurrence Matrix (GLCM) method with a 3 × 3 moving windows. The extracted texture metrics included mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment, and correlation, which provided auxiliary information for forest carbon stock estimation.

Based on the Sentinel-2, Landsat 8 OLI and meteorological data, the following variables were extracted, respectively, (as shown in Table 1).

Table 1. Feature factors.

3.2 Sample plot data

The sample data were derived from the 2022 continuous forest resources survey data of Yunnan Province. A total of 87 sample data were selected for modeling covering diverse forest types (as shown in Table 2). The distribution of inventory plots among different forest types is presented in Table 3.

Table 2. Plot statistics.

Table 3. Distribution of inventory plots among different forest types.

3.3 Sample area carbon storage calculation

Carbon storage can be calculated by converting forest volume into biomass and then applying the corresponding carbon content coefficient. The biomass of forest vegetation was estimated using the biomass–volume relationship model developed by Fang et al. (1996), and the calculation formula is as follows:

B = aV + b

In the formula, B represents the biomass per unit area (unit:t / h m 2); V represents the volume per unit area (unit:m 3 / h m 2); a and b are model parameters. For species-specific calculations, parameters from closely related tree species were used for species lacking specific biomass–volume conversion models. The parameters of the biomass–volume regression models for dominant tree species are presented in Table 4. Most studies have estimated forest vegetation carbon storage by applying a constant carbon fraction of 0.5 or 0.45. However, because carbon content coefficients vary considerably among tree species or species groups, the use of a single fixed value may reduce the accuracy of carbon storage estimates. Therefore, this study followed the “Guidelines for Carbon Sink Measurement and Detection of Afforestation Projects” (LY/T2253-2014) previous studies (Sun et al., 2004; Yang and Guan, 2008; Yang, 2021) for carbon storage estimation. For tree species (groups) without measured carbon content rates, the values of closely related species were adopted as substitutes.

Table 4. Biomass models and carbon content coefficients for dominant tree species.

3.4 Carbon storage estimation model

3.4.1 Random forest

Random Forest (RF) is a machine learning method based on an ensemble of multiple decision trees. The model constructs multiple regression trees through bootstrap sampling and random feature selection, and the final prediction is obtained by averaging the outputs of all trees (Biau and Scornet, 2016). Owing to its ensemble structure, RF can effectively reduce model variance and improve prediction stability. In addition, RF can effectively capture nonlinear relationships and handle high-dimensional data, while showing relatively strong robustness against noise and overfitting. Therefore, Random Forest has been widely applied in remote sensing and forest carbon stock estimation studies.

In this study, the Random Forest model was further optimized to ensure a fair comparison with the proposed Bayesian hierarchical model. A grid search strategy combined with five-fold cross-validation was adopted to determine the optimal hyperparameter combination. The optimized parameters included the number of trees (n_estimators), maximum tree depth (max_depth), minimum samples required for leaf nodes (min_samples_leaf), minimum samples required for node splitting (min_samples_split), and the number or proportion of features considered at each split (max_features). The optimized Random Forest model was subsequently evaluated using the same predictor variables and leave-one-out cross-validation (LOOCV) procedure as the Bayesian hierarchical model. The parameter information from different data sources is shown in Table 5.

Table 5. Parameter configurations of the random forest models for different data sources.

3.4.2 Ordinary least square

Ordinary Least Squares (OLS) is a widely used statistical approach for estimating parameters in linear regression models. The method estimates regression coefficients by minimizing the sum of squared residuals between observed and predicted values, thereby obtaining the best-fitting linear relationship between predictors and response variables (Georganos et al., 2021).

3.4.3 Bayesian hierarchical model

The Bayesian Hierarchical Model (BHM) is a statistical modeling approach that can effectively describe hierarchical structures and variability in complex datasets, and is particularly suitable for small-sample studies (Ver Planck et al., 2018). In this study, the model was constructed using three levels, including the data layer, process layer, and parameter layer.

Before model construction, all continuous predictor variables and the response variable were standardized to improve numerical stability during Bayesian inference and to place parameters on comparable scales for prior specification.

At the data level, the field-measured forest carbon storage was treated as the response variable, and the probabilistic relationship between observations and the model-predicted mean was formulated. Considering that the observation errors may deviate from normality, the Student’s t-distribution was adopted as the likelihood function to account for potential heavy-tailed errors.

y i ~ Studen t t ( v , μ i , σ g i )

where y i denotes the observed forest carbon storage at the i-th plot, μ idenotes the predicted mean carbon storage estimated by the hierarchical model, σ g i represents the observation error scale parameter corresponding to the vegetation type of the i-th plot, and νis the degrees-of-freedom parameter of the Student’s t-distribution controlling the tail heaviness of the distribution.

The model linear predictor was formulated by integrating fixed effects with hierarchical random effects associated with different forest vegetation types:

μ i = α g i + X i β + Z i b g i

where α g i represents the random intercept for the g i-th forest vegetation type, X i is the matrix of fixed-effect variables covariates describing population-level effects, and Z i represents the design matrix for hierarchical random effects, linking observations to vegetation-type-specific deviations from the population-level relationship. The parameters β and b g irepresent the fixed-effect coefficients and the random-effect coefficients, respectively.

At the process level, a hierarchical prior structure was introduced to model parameter variations among different forest vegetation types. The random intercept was modeled using a hierarchical structure centered around the overall mean:

α g = α 0 + δ g

δ g ~ N ( 0 , σ α 2 )

where α 0 denotes the population-level mean of the random intercepts across all forest vegetation types; δ g represents the deviation of the g-th forest vegetation type from the overall mean; and σ α represents the variability scale of the random intercept.

The random slope parameters associated with selected hierarchical predictors were modeled using a similar hierarchical structure:

b gj = b 0 j + γ gj

γ gj ∼ N ( 0 , σ b j 2 )

where b 0 j denotes the population-leve mean random slope for the j-th predictor variable; γ gj represents the deviation of the g-th forest vegetation type from the overall slope; and σ b j represents the variability scale of the random slope parameter for the j-th predictor variable.

This hierarchical structure follows the Bayesian partial pooling framework, which allows information sharing among different vegetation types while preserving their individual characteristics. Compared with models that estimate parameters independently for each group, partial pooling can effectively utilize shared information among vegetation types, reduce parameter uncertainty for groups with limited samples, and improve model generalization capability.

At the parameter level, prior distributions were specified for model parameters within the Bayesian framework. Weakly informative normal priors were assigned to the fixed-effect regression coefficients to provide regularization and constrain parameter estimates. Hierarchical priors were specified for the population-level parameters and variability scales of the random intercepts and random slopes. In addition, independent observation error scale parameters were assigned to different forest vegetation types and constrained using half-normal priors. The degrees-of-freedom parameter of the Student’s t-distribution was assigned an exponential prior. By integrating prior information with plot-level observations, the model can quantify parameter uncertainty and improve the stability and reliability of forest carbon storage estimation under limited-sample conditions. The specific prior distributions were specified as follows:

β ~ N ( 0 , 0.4 2 )

α 0 ~ N ( 0 , 0.5 2 )

b 0 j ~ N ( 0 , 0.4 2 )

σ α ~ HalfNormal ( 0 , 0.5 )

σ b j ~ HalfNormal ( 0 , 0.6 )

σ g ~ HalfNormal ( 0 , 0.6 )

v ~ Exponential ( 1 / 1.5 ) + 2

where β represents the fixed-effect coefficient; the HalfNormal distribution is used to constrain scale parameters to positive values; and the degrees-of-freedom parameter is shifted by adding a constant of 2 to ensure that the Student’s t distribution has a finite variance. These prior specifications provide moderate parameter regularization under limited-sample conditions while avoiding the introduction of overly strong prior assumptions that may affect model inference.

Model parameters were estimated using the Markov Chain Monte Carlo (MCMC) algorithm (Lye et al., 2022). Specifically, the No-U-Turn Sampler (NUTS), an improved algorithm based on Hamiltonian Monte Carlo (HMC), was employed for posterior sampling to enhance sampling efficiency and parameter space exploration capability in complex hierarchical models. The model was implemented using the PyMC 5.0 framework in a Python 3.9 environment. Four independent chains were run, with each chain consisting of 1,000 warm-up iterations and 1,000 posterior sampling iterations to obtain stable posterior parameter distributions. The target acceptance rate was set to 0.99 to improve sampling stability. No thinning was applied, as the NUTS sampler generally produces posterior samples with low-autocorrelation posterior samples, making thinning unnecessary.

Through iterative sampling, posterior distributions of model parameters were obtained based on Bayes’ theorem:

p ( θ ∣ y , X ) ∝ p ( θ ) × p ( y ∣ X , θ )

where p ( θ ) represents the prior distribution of the parameters, and p ( y ∣ X , θ ) represents the likelihood function. The posterior distribution combines the prior information and the information from the observed data, providing a comprehensive estimation of the model parameters.

For a new predictor variable x 0, the predictive distribution of the target variable can be expressed as:

p ( y 0 ∣ y , X , x 0 ) = ∫ p ( y 0 ∣ θ , x 0 ) p ( θ ∣ y , X ) dθ

This predictive framework incorporates both parameter uncertainty and observed information, enabling more reliable prediction for new samples.

3.5 Model accuracy

Cross-validation is an effective approach for evaluating and comparing the performance of different predictive models, particularly when dealing with small sample datasets. Leave-one-out cross-validation (LOOCV) can provide nearly unbiased estimates of generalization error and can be used as a reliable basis for parameter selection (Ding et al., 2019). Based on the results obtained from cross-validation, model performance is evaluated using the Coefficient of Determination (R 2), the Root Mean Squared Error (RMSE), and the Relative Root Mean Square Error (rRMSE) after cross-validation. Its calculation formula is as follows:

R fit 2 = 1 − S S fit S S tot

R fit 2 represents the modeling accuracy, S S fit is the sum of squares of the residuals of the fitted model, S S tot is the sum of squares of the differences between each observation and the overall average of the observations.

R cv 2 = 1 − S S cv S S tot

R cv 2represents the accuracy of the model validation results, while S S cv represents the sum of the squares of the differences between all the observed values and predicted values obtained through cross-validation.

RMSE = ∑ i = 1 n ( pr e i n − ob s i ) 2 n

rRMSE = RMSE obs ¯

4 Results

4.1 Model factor selection

Predictor variables were selected using a stepwise regression approach. Specifically, a stepwise regression equation was first established using all training samples to identify global predictors representing overall variation patterns. Subsequently, the training samples were grouped according to forest types, and stepwise regression models were independently developed within each group to extract hierarchical predictors. Finally, the global and hierarchical predictors were integrated to form the final predictor set (Jian et al., 2025). Considering that different predictor selection strategies may introduce potential bias into model comparison, the identical predictor set was used to construct the Bayesian Hierarchical Model, Ordinary Least Squares, and Random Forest models, ensuring a fair comparison among different model approaches under the same input information conditions. The screening results are shown in Table 6. Under different data sources, the selected variables differ in type and combination, mainly including remote sensing spectral and texture features, and climatic factors.

Table 6. Modeling factors.

X 30 m refers to data resampled to 30 m resolution; X 10 m indicates a resolution of 10 m; ASM stands for angular second moment.

4.2 Model construction and accuracy evaluation

The performance metrics of the three models constructed based on different data sources are shown in Table 7. The Bayesian Hierarchical Model showed superior fitting performance across all data sources: in the cases of Landsat 8, Sentinel-2, and multi-source fusion data, its coefficients of determination R 2 reaches 0.80, 0.82, and 0.83 respectively, all remaining above 0.80; in contrast, the maximum R 2 of the ordinary least squares model was only 0.17 (in the multi-source data scenario), and the R 2 of the RF is consistently within the range of 0.47 to 0.57, with the predictive performance of these two models being significantly lower than that of the Bayesian Hierarchical Model. Moreover, by comparing the fitting effects of the same model under different data sources, it was found that the combined data of Landsat 8 and Sentinel-2 achieved better model performance than Landsat 8 or Sentinel-2 data alone. Among them, the fitting effect using the BHM in the multi-source data scenario is the most ideal.

Table 7. Model accuracy.

In terms of the model’s generalization performance, the validation R 2 of the BHM ranges from 0.73 to 0.76. Even when using Landsat 8 data alone, it still maintains a relatively high level of 0.73. The maximum RMSE of the model was 8.93 t / h m 2, and the rRMSE is within the range of 42.4 to 45.3%. On the other hand, the validation R 2 of the ordinary least squares model is only between 0.01 and 0.07. The validation R 2 of the random forest model is between 0.12 and 0.21. Both types of models have an RMSE exceeding 15 t / h m 2, and the rRMSE is significantly higher.

Given the superior performance achieved using multi-source fusion data compared with individual data sources, subsequent analyses were conducted based on the multi-source dataset.

4.2.1 Spatial autocorrelation analysis of model residuals

We used Moran’s I (Moran, 1950) to test for spatial autocorrelation in model residuals. The results showed that the Moran’s I value of the residuals was −0.0299, with a p-value of 0.766, which was greater than the significance level of 0.05. Therefore, the null hypothesis of spatial randomness of the residuals could not be rejected. Although Moran’s I exhibited a slightly negative value, the statistically insignificant result does not provide evidence of a systematic spatial pattern in the residuals (Beale et al., 2010). No significant spatial dependence was detected in the model residuals, suggesting that the assumption of spatial independence of observation errors adopted in this study was reasonable. Therefore, considering the current spatial scale and sampling conditions, the non-spatial Bayesian hierarchical model framework was retained for forest carbon stock estimation.

4.2.2 Robustness assessment under different data partitioning schemes

To further assess the robustness of the Bayesian hierarchical model, five-fold cross-validation was conducted as a supplementary analysis. Unlike LOOCV, which was used for the primary performance evaluation, five-fold cross-validation was performed to assess whether model performance remained stable under different training-validation partition schemes. The results are shown in Figure 2. The results indicate that the validation R 2 remained within the range of 0.70–0.75 across different random partitions, with limited variation. This suggests that the predictive performance of the model is relatively insensitive to the training sample partition and does not depend on a specific training set, indicating good robustness and reducing the risk of overfitting.

4.3 Bayesian hierarchical model convergence diagnostics and posterior parameter analysis

The convergence of the Bayesian hierarchical model was comprehensively evaluated using trace plots, the Gelman–Rubin statistic (R-hat), and the effective sample size (ESS). The results are presented in Figure 3 and Table 8. As shown in Figure 3, the four MCMC chains of the representative parameters (a 0, β, and v) exhibited good mixing and stable sampling behavior, without apparent trends or chain separation, indicating satisfactory convergence of the posterior sampling process for these key parameters. According to the convergence diagnostics of all model parameters (Table 8), the R-hat values ranged from 0.999 to 1.003, all below the recommended threshold of 1.05, indicating satisfactory convergence among different MCMC chains. Meanwhile, the ESS values of all parameters exceeded 1700, which is substantially higher than the commonly accepted criterion of 200, indicating that the posterior samples contained sufficient effective samples for reliable inference.

Table 8. Parameter estimates and convergence diagnostics.

Further examination of the posterior distributions of the key parameters (a 0, β, and v) revealed highly consistent distribution patterns among different MCMC chains. The posterior mean of the overall intercept parameter α₀ was 0.095, indicating that the model effectively captured the overall response characteristics of forest carbon storage. The posterior mean of the degrees-of-freedom parameter ν of the Student’s t distribution was 4.013, suggesting that the observation errors exhibited certain heavy-tailed characteristics. This finding supports the use of the Student’s t likelihood function, which improves the model’s robustness against potential outliers. In addition, the variance parameters of the random effects (σ a and σ b) obtained stable posterior estimates, indicating the existence of inter-group heterogeneity among different vegetation types. The hierarchical structure successfully captured vegetation-type-specific variations and enhanced the stability of model estimation under small-sample conditions through information sharing.

4.4 Comparison of model accuracy under different sample sizes

To further evaluate the performance of the Bayesian hierarchical model under limited sample conditions, different sample sizes (25, 30, 50, and 80 samples) were randomly selected for modeling. The prediction performance of each model was evaluated using leave-one-out cross-validation (LOOCV), and the results are presented in Table 9. The results demonstrate that the Bayesian hierarchical model remained effective even with limited sample sizes. Even with only 25 samples, the proposed model achieved relatively high and stable prediction performance (R cv 2 ≈ 0.65, with a lower RMSE). In contrast, OLS exhibited limited predictive capability, with R cv 2 values ranging from −0.29 to 0.08 across different sample sizes, indicating weak explanatory and predictive performance. Although RF achieved improved prediction compared with OLS, with R cv 2 values ranging from −0.09 to 0.28, its performance remained substantially inferior to that of the BHM.

Table 9. Model accuracy at different sample sizes.

4.5 The parameter estimation results and convergence performance under different sample sizes

Figure 4 presents the posterior distributions of the global intercept ( a 0 ), fixed effect coefficient β 1, and random effect standard deviation ( s d b ) in the Bayesian hierarchical model under different sample sizes (n = 80, 50, 30, 25), as well as the MCMC convergence diagnostic results. Overall, the posterior distributions of each parameter exhibit a single peak and continuous characteristics, and the center positions of the distributions are basically the same across different sample sizes, indicating that the model parameters have good identifiability. As the sample size decreases, the width of the parameter posterior distributions increases, reflecting a reasonable expansion of uncertainty, but no significant deviation or abnormal fluctuations are observed. The MCMC trace plots show that the sampling process of each parameter is stable and there is no obvious trend drift. The sampling chains are well mixed. The results of the autocorrelation analysis indicate that the parameter autocorrelation rapidly decays after a small number of lags, indicating a high effective sample size and reliable sampling results. In summary, the Bayesian hierarchical model constructed shows good numerical stability and convergence across different sample sizes.

Figure 5 shows the prediction results of different models under different sample sizes. From Figure 5, it can be seen that the Bayesian hierarchical model can simultaneously provide the predicted mean and the corresponding uncertainty interval. In cases where the sample size is small or the observed values fluctuate greatly, the prediction intervals become wider to reflect the estimated uncertainty, thereby avoiding over-interpreting the prediction results. This uncertainty expression is something that traditional OLS and random forest methods cannot directly provide.

4.6 Fitting of different forest types under multiple data sources

As shown in Figure 6, the prediction errors of the three models varied considerably among different forest types. Overall, BHM consistently achieved lower RMSE and rRMSE values across all three forest types, whereas OLS and RF showed larger differences in error magnitude among forest types. In terms of RMSE, BHM achieved the lowest RMSE values for all three forest types. For mixed forests, RF showed the highest RMSE, followed by OLS, while BHM exhibited the lowest error. Among coniferous forests, RF showed the largest RMSE among the three models. For broad-leaved forests, OLS produced a substantially higher RMSE than RF and BHM. In terms of rRMSE, the differences among models were most pronounced in mixed forests, where the rRMSE values followed the order of RF > OLS > BHM. The differences in rRMSE among the three models became less pronounced in coniferous forests; however, BHM still maintained the lowest relative error. For broad-leaved forests, OLS showed a considerably higher rRMSE, whereas RF and BHM exhibited similar and substantially lower values.

Considering both error indicators, BHM demonstrated the most consistent and accurate estimation performance across the three forest types. RF showed moderate performance, with relatively better results for broad-leaved forests but larger absolute errors for coniferous forests. In contrast, OLS exhibited weaker overall performance, particularly in broad-leaved forests, where prediction errors were substantially higher.

The above results fully demonstrate that the Bayesian Hierarchical Model has higher accuracy and generalization ability than the ordinary least squares method and RF when the sample size is small, making it highly effective in estimating multiple carbon stocks. Since the Bayesian Hierarchical Model outperforms the OLS and RF in terms of reliability, prediction precision, and stability, the subsequent analysis will be based on the modeling results of this model.

4.7 Spatial distribution characteristics of forest carbon storage

The results of the estimation of the Three Parallel Rivers Region based on the Bayesian Hierarchical Model and multi-source data are shown in Figure 7. Figure 7A is the distribution map of carbon storage in broad-leaved forests, Figure 7B is the distribution map of carbon storage in coniferous forests, Figure 7C is the distribution map of carbon storage in mixed forests, and Figure 7D is the distribution map of overall forest carbon storage.

As shown in Figure 7, the carbon storage of broadleaved forests in the study area is mainly at the median level, mainly concentrated in the western and southern local areas of the study region. The low and high value areas are relatively scattered. The carbon storage of mixed forests has a relatively smaller distribution range, with a point-like or small patch-like distribution, mainly distributed in Nujiang Lisu Autonomous Prefecture, and the carbon storage is relatively low. Overall, in the forest vegetation of the study area, the carbon sequestration effect of broadleaved forests is more significant than that of coniferous forests. As can be seen from Figure 7D, the carbon storage of the forests in the study area is widely distributed, with the carbon storage ranging from 0 to 245.86 t / h m 2, presenting a spatial distribution pattern of higher values in the north and middle parts and lower values in the south and east parts. This may be due to the fact that the northern area is located in the Hengduan Mountainous Area, with complex terrain, less human activity interference, and intact original forests; while the eastern area is because the tourism industry is developed in this area, and interference from human activities has potential impacts on the forest ecosystem.

Spatial characterization of uncertainty in forest carbon stock estimation was conducted based on the width of the 95% credible interval (CI width) derived from the Bayesian hierarchical model predictions, the results are shown in Figure 8. The results showed that the uncertainty of forest carbon stock estimates exhibited evident spatial heterogeneity, while the overall uncertainty level remained relatively low. The 95% credible interval width of most forest areas in the study region was mainly concentrated within the range of 0–30 t / h m 2, indicating that the model predictions were relatively stable and that the posterior predictive distributions were highly concentrated. In particular, forest areas with continuous spatial distribution generally exhibited smaller CI widths, suggesting stronger relationships between remote sensing variables and forest carbon stocks, as well as better sample representation, allowing the model to effectively capture spatial variations in forest carbon storage. From the perspective of spatial distribution patterns, areas with low uncertainty were widely distributed across the study region, whereas areas with higher uncertainty showed a scattered distribution pattern. Compared with other regions, the western, northern, and parts of the eastern areas exhibited relatively higher uncertainty levels, with some regions showing CI widths exceeding 50 t / h m 2, and a small proportion of areas exceeding 70t / h m 2. These high-uncertainty areas were mainly distributed in regions characterized by complex terrain conditions and strong ecological heterogeneity. The increased uncertainty may be associated with variations in topographic conditions, differences in vegetation composition, and limited field sample availability, which could increase parameter uncertainty and spatial extrapolation uncertainty during model prediction.

By superimposing the DEM data on the forest above-ground carbon storage distribution in the study area, different scenarios of carbon storage at different elevation gradients, slopes, and aspects can be obtained, as shown in Figure 9. Vertically, the above-ground carbon storage of broadleaved forests is predominantly concentrated at elevations between 2,500 and 3,000 m, representing 32.13% of the total above-ground carbon in broadleaved forests within the Three Parallel Rivers region. In contrast, coniferous forests store most of their carbon above 3,000 m, whereas mixed forests are mainly distributed at 2,000–2,500 m, accounting for 31.97% of the above-ground carbon in mixed forests across the region. Overall, the forest carbon storage in the Three Parallel Rivers Region is mainly concentrated in the area above 2000 m. From the aspect of slope, the forest carbon storage in the Three Parallel Rivers Region is mainly distributed within the range of 20–40°. Areas with a slope exceeding 40° have a lower proportion of forest carbon storage due to the larger slope. From the aspect of slope direction, the area with no specific slope direction has the smallest proportion of carbon storage, while the sunny slope and semi-sunny slope have the highest proportion, which may be mainly due to the sufficient light conditions, which are conducive to the photosynthesis of forest vegetation and the accumulation of more biomass, thereby increasing the carbon storage.

5 Discussion and conclusion

5.1 Discussion

Estimation of forest carbon storage is prone to data sparsity and sample heterogeneity. In reality, a small sample size often leads to poor modeling performance and prediction accuracy, thereby increasing the uncertainty of the estimation results. This is because the established model lacks robustness and representativeness (Lu et al., 2016; Ploton et al., 2020). This study is the first to apply the BHM to the estimation of carbon storage in the Three Parallel Rivers Region under small-sample conditions. The results show that compared with traditional modeling methods (such as OLS and RF), the BHM has the most ideal effect. Moreover, even under limited sample size conditions, it can still to some extent alleviate the uncertainty problem caused by data sparsity in the estimation, demonstrating good numerical stability and convergence. This conclusion is consistent with the studies conducted by Zapata-Cuartas et al. (2012) and Wang et al. (2019). It should be noted that the inferior performance of traditional models became particularly evident under extremely limited sample conditions. Negative R cv 2 values were observed for both the ordinary least squares (OLS) regression and random forest (RF) models. When only 25 samples were available, the OLS and RF models achieved R² values of −0.29 and −0.09, respectively, indicating that their predictive performance was worse than simply using the mean value of observations as the prediction. This phenomenon can be attributed to the insufficient sample size for reliable parameter estimation and model training. Conventional statistical and machine learning models generally require adequate observations to capture underlying relationships and maintain model stability. When the sample size is insufficient, parameter uncertainty increases substantially, resulting in poor generalization ability and even negative R cv 2 values. In contrast, the Bayesian hierarchical model maintained superior predictive performance under the same small-sample conditions, achieving an R² value of 0.65 with only 25 samples. These results demonstrate that the BHM can effectively address the uncertainty and instability associated with limited sample availability. These results indicate that the Bayesian hierarchical model can effectively alleviate these limitations and maintain robust predictive performance under the same small sample conditions.

This improved performance can be mainly attributed to the hierarchical structure and parameter modeling strategy of the BHM. The advantage of the BHM lies in its hierarchical structure and parameter modeling approach. The model allows different parameters to be specified at different levels and adopts hierarchical prior information, which helps to reduce the excessive influence of prior information on the model estimation results. By comprehensively considering various factors at multiple levels and treating the parameters in the model as random variables (White et al., 2021), it can reveal the heterogeneity of sample data between different groups and the dissimilarity of different sample data within the same group, thereby improving the flexibility and accuracy of the modeling. In particular, the model incorporates factors such as forest types, which have random variation characteristics, into the hierarchical structure, explicitly including the differences among different types in the model fitting process (Babcock et al., 2018; Veenman et al., 2024). At the same time, it systematically characterizes the uncertainty by considering both fixed effects and random effects (Browne and Draper, 2006). In addition, by incorporating weak informative priors to impose regularization constraints on the parameters, the model can effectively integrate overall and hierarchical information even when the sample information is limited; this enables it to alleviate the uncertainty caused by insufficient sample sizes and improve model accuracy (van de Schoot et al., 2014). The results from the verification of forest type accuracy also reveal significant differences in the estimation accuracy of different forest types, providing a direction for further improving the precision of forest carbon storage estimation. Unlike traditional statistical methods that only provide point estimates, the Bayesian hierarchical model quantifies uncertainty and quantifies the prediction uncertainty of each plot through the posterior distribution and confidence intervals (Li et al., 2012). This enables the research results to not only be directly compared but also provide decision-makers with an assessment on the credibility of the prediction, thereby enhancing the interpretability and application value of carbon storage estimation. At the same time, Bayesian models use Markov chain Monte Carlo (MCMC) methods to assess parameter convergence and obtain posterior distributions, allowing uncertainty to be directly quantified. In contrast, frequentist approaches generally rely on sampling distributions and asymptotic theory for inference, where parameters are treated as fixed but unknown quantities (Bayarri and Berger, 2004; Gelman and Rubin, 1992).

From the perspective of spatial distribution characteristics, forest carbon storage exhibited substantial spatial heterogeneity under different terrain conditions. More than 80% of the total forest carbon storage in the study area was distributed in areas above an elevation of 2,500 m, with areas above 3,000 m accounting for the largest proportion (63.16%). The western and northern regions of the study area showed relatively higher carbon storage, whereas the eastern region exhibited comparatively lower values. The results of this study are consistent with those estimated by Cheng et al. (2024) and Liang and Li (2023). This spatial pattern may be related to the geographical characteristics of northwestern Yunnan, where complex terrain, high elevation, and relatively low human disturbance contribute to the conservation of natural vegetation and the maintenance of forest carbon stocks (Tang and Jianlin, 2019). In addition, differences among forest types also contributed to the spatial variation in forest carbon storage, reflecting the potential influence of vegetation composition and forest structural characteristics on carbon accumulation. From a management perspective, priority should be given to the conservation of existing high-carbon-storage forests in the western and northern regions to maintain their current carbon sequestration functions. In contrast, areas with relatively lower carbon storage may have greater potential for improvement through appropriate forest management and restoration measures. Given the strong spatial heterogeneity of forest carbon storage, future management strategies should consider local topographic conditions, vegetation characteristics, and ecological constraints rather than adopting uniform approaches across the region.

Although this study has achieved certain results, it also has some limitations. The hierarchical structure model usually requires certain simplifying assumptions about the error distribution and hierarchical relationships. In complex ecosystems, these assumptions may not be fully met (Gelman et al., 1995). Furthermore, the BHM relies on numerical methods such as MCMC for parameter estimation. For complex models or large-scale training data, calculating the posterior distribution can be extremely time-consuming, may have slow convergence rates, and involve high computational costs (Zhang et al., 2025). These issues may limit the practical application of this model in large-scale regional or long-term sequence studies of forest carbon storage. Xie et al. (2020) demonstrated in their study on Bayesian statistics that the reasonable selection of prior information directly affects the fitting results. In this study, the weak informative prior adopted has little impact on the parameters. In the future, the construction of prior information can be optimized by combining machine learning algorithms (such as random forests, neural networks). The construction and selection of characteristic variables are also important factors influencing the uncertainty of forest carbon storage estimation. Stepwise regression has been widely applied in forest remote sensing modeling studies due to its advantages of automated and data-driven variable selection, which enables the identification of key variables from numerous candidate predictors and the removal of redundant variables with limited contributions (Jiang et al., 2021; Jiang F. et al., 2020; Tariq et al., 2021). However, this method relies on significance-test-based variable selection, which may result in biased parameter estimates, increased model selection uncertainty, and reduced reliability of subsequent statistical inference (Chatfield, 1995). From the perspective of multimodel inference, a single optimal variable subset obtained through automated selection may not fully capture the uncertainty among different candidate models (Burnham and Anderson, 2002). Therefore, future studies could integrate multimodel inference strategies with more flexible variable selection methods, such as Bayesian variable selection and shrinkage estimation approaches, to further optimize the variable selection process and improve model stability and uncertainty assessment. Furthermore, the formation of forest carbon storage is influenced by various factors such as geographical environment, climate change, stand structure, and human activities, as well as their interactions. This study only selected some key variables for modeling and has not fully covered information on forest structure and habitat conditions. Future research can, based on the availability of samples and data, introduce environmental variables such as stand structure parameters and soil attributes, or combine laser radar and other technical means to obtain more accurate information on forest height and vertical structure, in order to further reduce estimation uncertainty and improve model accuracy.

Overall, this study demonstrates that the BHM, through the organic integration of a hierarchical structure and prior information, provides an effective method for estimating forest carbon storage under small sample conditions, which takes into account both accuracy and uncertainty characterization. With the continuous enrichment of multi-source remote sensing data and auxiliary information, the application potential of this method in regional forest carbon sink assessment and refined management is expected to further increase.

5.2 Conclusion

Sentinel-2 imagery showed slightly better performance than Landsat 8 OLI imagery, while the fusion of both datasets achieved the highest estimation accuracy. Based on the experimental results of this study, the following conclusions are drawn:

• By comparing the forest carbon stock estimation results obtained from individual Landsat 8 OLI data, Sentinel-2 data, and their combination, it was found that the fused dataset achieved the best estimation performance. The three models constructed using the combined Landsat 8 OLI and Sentinel-2 data exhibited relatively higher fitting accuracy and lower RMSE values, indicating that the complementary information from the two satellite datasets can improve the accuracy of remote sensing-based forest carbon stock estimation compared with using either dataset alone.
• The comparison of the three remote sensing estimation models indicates that the Bayesian Hierarchical Model (BHM) performs the best overall in estimating forest carbon stocks using remote sensing variables. This model achieved the highest fitting accuracy, with an R 2 value of 0.83. Compared with traditional modeling methods such as ordinary least squares and random forests, the Bayesian Hierarchical Model can maintain relatively stable predictive performance even when the sample size is reduced. BHM can effectively address the issues of unstable parameter estimation and decreased predictive performance in traditional models under very small sample conditions, providing a more robust solution for forest carbon stock estimation with limited data.
• The spatial analysis of forest carbon stock in the Three Parallel Rivers Region based on model inversion revealed distinct spatial variations. Forest carbon stock was generally higher in the northern and central regions and lower in the southern and eastern regions. Areas with higher carbon stock were mainly distributed in regions with higher elevations and relatively intact forest ecosystems. Significant differences in carbon stock were observed among different terrain conditions and vegetation types, indicating that terrain conditions and vegetation types are important factors associated with the spatial variation of forest carbon stock.
• The proposed Bayesian hierarchical modeling framework provides a potential approach for improving forest carbon stock estimation under limited sample conditions. By integrating multi-source remote sensing data and uncertainty quantification, this framework enables more reliable carbon stock estimation and provides useful information for forest carbon monitoring and assessment in mountainous regions.

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: Matthias Peichl, Swedish University of Agricultural Sciences, Sweden

Reviewed by: Zhou Huoyan, Yunnan University, China

Tsikai Solomon Chinembiri, University of Cape Town, South Africa

Footnotes

https://data.casearth.cn/