Explainable HybridEnsemble approach with golden jackal optimization for AGB estimation using multi-sensor remote sensing

Accurate spatial measurement of aboveground biomass (AGB) is essential for assessing carbon stocks in the forest ecosystem. To enhance this estimation, integrating active and passive Earth Observation data with advanced machine learning techniques offers a promising approach. This study presents an integrated HybridEnsemble model with golden jackal optimization for AGB estimation and evaluates its predictive performance against individual base-learners, including categorical boosting (CatBoost), extreme gradient boosting (XGBoost), light gradient boosting machine (LightGBM), and adaptive boosting (AdaBoost) via the synergistic application of Explainable Artificial Intelligence (XAI) and active and passive datasets. The findings revealed a clear performance ranking among these models, with the HybridEnsemble Golden Jackal Optimization (HGJO) model identified as the most effective, which yielded a correlation coefficient (R 2) of 0.821 and a Relative Root Mean Square Error (rRMSE) of 16.30%. This performance was followed by CatBoost (R 2 = 0.816, rRMSE = 16.51%), LightGBM (R 2 = 0.804, rRMSE = 17.05%), XGBoost (R 2 = 0.802, rRMSE = 17.13%), and AdaBoost (R 2 = 0.731, rRMSE = 19.97%), with all comparisons reported at 95% confidence intervals. XAI revealed that predictors from optical sensors (passive) were strongly correlated with AGB and played a significant role in predicting AGB, while features derived from SAR (synthetic aperture radar, an active sensor), less influential, provided unique backscatter and context-specific insights that enhanced the model’s performance. Forecast results indicate an increasing trend. Additionally, the analysis revealed that future AGB accumulation in the subtropical forest of Hong Kong will be highly variable and strongly dependent on initial biomass levels, with high-biomass plots likely to see the greatest gains. However, spatial uncertainty in AGB predictions varied across the study area, with higher uncertainties observed in forested areas and lower uncertainties in urban areas. Overall, this study not only enhances understanding of optimized hybrid ensemble models for biomass prediction but also offers valuable insights for forecasting forest dynamics, supporting sustainable forest management and carbon stock monitoring globally.

Forests are central to the global carbon cycle. They store substantial amounts of carbon, regulate exchanges of carbon dioxide with the atmosphere, sustain ecosystem functioning, and thereby contribute to climate change mitigation (Wang et al., 2025). Accurate, spatially explicit quantification of forest aboveground biomass (AGB) is therefore essential for constraining carbon budgets and designing effective mitigation strategies under a changing climate (Luo et al., 2024). At the national scale, consistent and comparable AGB monitoring supports climate adaptation planning and land-use decision making. It also enables sustainable forest governance by linking observed biomass trends to environmental drivers, disturbance regimes, and policy priorities (Prăvălie et al., 2023).

Historically, aboveground biomass (AGB) has been estimated from forest inventory plots and stand-level measurements, with plot values extrapolated using species or site-specific allometric equations parameterized by biometric variables such as diameter at breast height (DBH), wood density, and tree height. Although effective, plot-based approaches are labor-intensive, costly to repeat, and inherently limited in spatial coverage (Li et al., 2021; Prăvălie et al., 2023). Considerable progress has been made in developing and refining allometric models across diverse forest types. However, the growing demand for spatially continuous biomass information over large regions has intensified the need for scalable mapping frameworks (Tang et al., 2022). In response, remote sensing has emerged as a robust means of extending AGB estimation beyond plot networks. Optical imagery contributes spectral reflectance features, vegetation indices, and canopy condition proxies. LiDAR provides direct measurements of vertical forest structure and synthetic aperture radar (SAR) offers sensitivity to canopy architecture and moisture dynamics, with the advantage of all-weather, day-night acquisition (Tang et al., 2022). Together, the integration of Earth observation data with predictive algorithms enables consistent and spatially explicit AGB mapping across broad geographic extents (Li et al., 2021).

Optical sensors such as Sentinel-2 with its incorporated red-edge bands provide rich spectral information well-suited for vegetation condition and biomass assessment, while SAR contributes structural sensitivity and all-weather imaging. L-band signals from SAR sensors such as ALOS-2 PALSAR-2 penetrate deeper into canopy layers than C-band signals, but may still saturate at high biomass, as dense canopies increase scattering, reducing sensitivity to further biomass increases (Bui et al., 2024; Li et al., 2021). In essence, optical vegetation indices and texture provide details related to vegetation characteristics, whereas SAR backscatter and coherence are associated with vegetation structure information and moisture, together providing complementary drivers for AGB prediction (Bui et al., 2024; Li et al., 2021). SAR backscatter and polarimetric features respond to woody structure and dielectric properties, penetrate clouds, and operate day and night, but can also saturate at medium-high biomass, especially at shorter wavelengths. Spatial configurations that incorporate various frequency modes are often used to boost sensitivity (Li et al., 2021).

Since optical and SAR sensors capture complementary biophysical information, their fusion has typically improved AGB retrievals relative to single-sensor approaches. Multi-source integration combines spectral, textural, and structural information to mitigate signal saturation particularly in high-biomass stands and to improve model generalization across forest types and biomass density or gradients (Pham et al., 2020; Teijido-Murias et al., 2025; Wang et al., 2025). Integrating Sentinel-2 with ALOS-2 PALSAR-2 has proven to be particularly effective for AGB mapping. Results from several studies indicate that combining Sentinel-2 with ALOS-2 PALSAR-2 yields substantial accuracy and reduces residual variance compared to single-sensor models (Pham et al., 2020; Wang et al., 2025). For example, a case study by Pham et al. (2020) demonstrated that the integration of optical and SAR data substantially improves model performance, as evidenced by a marked increase in the coefficient of determination (R 2). Optical-only models achieved approximately 0.60, SAR-only models around 0.49, while the fused Sentinel-2 and ALOS-2 PALSAR-2 dataset combined with the XGBoost algorithm attains an R 2 of approximately 0.81. This significant enhancement underscores the synergistic value of combining spectral and structural information for more accurate and robust AGB estimation. Recent regional studies in northwest Spain likewise integrate Sentinel-2 with SAR-derived backscatter information and interferometric coherence, selecting XGBoost as top performers and underscoring L-band’s value for structure-sensitive biomass retrieval (Teijido-Murias et al., 2025).

Machine learning frameworks excel at modeling non-linear relations between AGB and multi-source predictors, with ensemble learners frequently ranking best (Tang et al., 2022). Within this space, gradient boosting decision tree methods such as extreme gradient boosting (XGBoost), light gradient boosting machine (LightGBM), categorical boosting (CatBoost) and related boosting families, including adaptive boosting (AdaBoost), offer strong accuracy, regularization, and feature importance diagnostics suitable for aboveground biomass prediction (Bui et al., 2024; Li et al., 2021; Tang et al., 2022). Comparative studies frequently show that boosting methods outperform other machine learning models and achieve even better results when combined through stacking or ensembling (Huang et al., 2023). This was demonstrated in a study by Borsah et al. (2025), where in-situ measurements were combined with optical and SAR images to predict AGB in a subtropical forest using machine learning models. Among the three models analyzed, the gradient boosting model demonstrated superior performance (R 2 = 0.84), surpassing both the K-nearest neighbors (KNN) model (R 2 = 0.67) and stepwise linear regression (R 2 = 0.57), while exhibiting minimal deviations in residual errors compared to the other two methods. In a related study, Li et al. (2020) employed random forest, linear regression, and XGBoost to estimate biomass in subtropical forests across China. Their analysis integrated fused imagery from Sentinel-1 and Landsat satellites. The results indicated that both random forest and XGBoost models outperformed linear regression. Furthermore, the researchers noted that parameter optimization had a more pronounced impact on the performance of XGBoost compared to random forest. Golden Jackal Optimization (GJO) offers a population-based optimization approach coupled with exploration-exploitation search strategy for automatically tuning ensemble learner hyperparameters within predefined bounds, with candidate solutions evaluated using cross-validated loss functions to identify robust configurations (Kehinde et al., 2026).

Recent workflows increasingly integrate Explainable Artificial Intelligence (XAI), such as SHAP, to attribute predictions to specific spectral, textural, and radar features, thereby enhancing interpretability and facilitating variable selection and hyperparameter tuning (Bui et al., 2024; Zhao et al., 2025; Teijido-Murias et al., 2025). While XAI tools were developed for various applications, its analyses reveal adaptable practices for interpretable biomass modeling that can be directly implemented in AGB workflows. XAI has also been paired with boosting models to quantify variable contributions, aiding in model simplification and uncertainty assessment (Ali et al., 2023; Zhao et al., 2025). However, an important gap remains in subtropical AGB research. Existing studies have not systematically compared how different boosting and ensemble implementations vary in their optimization behaviour, feature selection efficiency, residual error correction capacity, and interpretability when applied to multi-sensor datasets. In particular, it remains unclear which model architecture and optimization strategies yield the most accurate and explainable AGB predictions under subtropical forest conditions. This missing optimization and model comparison gap limits both methodological selection and operational deployment. Therefore, requires a rigorous comparative evaluation of boosting and ensemble frameworks within an interpretable AGB modelling pipeline. Thus, this study aims to. Develop HybridEnsemble model with golden jackal optimization for AGB estimation

Compare its predictive performance with individual base-learners, including CatBoost, XGBoost, LightGBM, and AdaBoost

Use XAI to interpret model behaviour and predictor contributions

Apply the best performing model to forecast future biomass and assess the reliability of the resulting AGB maps through uncertainty analysis.

This study makes novel contributions by (i) introducing GJO-optimized HybridEnsemble that jointly optimizes ensemble weights and base-learner hyperparameters within a single search framework, enabling coordinated tuning of both model structure and weighting strategy, (ii) integrating SHAP-based explainability with spatial uncertainty mapping, and (iii) producing district-level AGB products and interpretations for Hong Kong forest landscape to support monitoring and management while demonstrating a transferable workflow for subtropical forest settings.

Study area

Hong Kong (Figure 1), encompassing roughly 1,100 square kilometers, sits at the northern fringe of Asia’s tropical zone between 22°09′0″-23°37′0″ N and 113°52′0″-114°30′00″ E. The territory adjoins China’s Guangdong province northward while the South China Sea encircles its remaining boundaries. Its distinctive climate melds tropical and subtropical elements featuring extended hot, humid summers followed by briefer, cooler, drier winter periods (Dudgeon and Corlett, 1994). Hong Kong’s original evergreen and semi-evergreen forests were largely cleared between 1000 and 1700 A.D., before British rule. Today, level terrain is mostly urbanized, farmed, or abandoned, whereas steeper hillsides are dominated by grassland, shrubland, and secondary forest. These secondary forests represent just about 8 percent of Hong Kong’s total area and exist primarily as fragmented patches (Ashworth et al., 1993).

Data acquisition

Remote sensing data acquisition

Aboveground biomass in the subtropical forest was predicted by integrating ALOS-2 PALSAR-2 L-band dual-polarimetric data with Sentinel-2 imagery from 2016, 2018, 2020, 2022, and 2024. Details of both datasets are provided in Table 1.

Field data sampling

To correlate the parameters from remote sensing images with biomass from the field plots, a total of 125 sample plots, encompassing 6,532 different types of tree stands, were utilized in this research. Square sample plots, each measuring 20 m by 20 m, were established to align with the image resolution of roughly 20 m. These plots were placed in uniform forested areas, avoiding steep inclines, and were located about 15 m away from features such as buildings, water bodies, roads, and other structures. The plots were set up with the T300 GNSS receiver, which is a Dual Frequency RTK Receiver. Position measurements were recorded only when Position Dilution of Precision (PDOP) was below 4 to ensure positional accuracy. Each tree in the plots was then scanned using the Green Valley LiBackpack DGC50 LiDAR System. This is a three-dimensional scanning device that is mounted on the back of an individual, produced by Beijing Green Valley Technology Co. Ltd. The main components of the LiBackpack are the VLP-16 LiDAR Sensor and the LiDAR360 software. The LiDAR360 provides modules for spatial analysis, as well as many features for processing and analyzing LiDAR point cloud data (Xie et al., 2020). For each tree in the designated square plots, the height of the tree and the Diameter at Breast Height (DBH) of the tree were measured. Finally, the biomass for each tree in the plots and the total biomass for all the trees in the plots were estimated using the allometric equation (Table 2) for the area of study.

To develop region-specific or local allometric model, destructive sampling of subtropical trees in the study area was conducted. A preliminary survey was conducted before harvesting, and surrounding vegetation was cleared to avoid contamination. After harvest, each tree was separated into leaves, branches, and trunk, and each fraction was assigned an identification number. Several samples of each component were collected, and their fresh weight, height, and diameter were measured. A total number of one hundred subtropical trees were destructively sampled. To facilitate oven drying, the tree samples were transported to the lab as quickly as possible. Each sample was kept under controlled dry conditions properly packaged and documented. Tree samples were oven-dried at 65 °C until constant weight and weighed immediately after removal to avoid moisture related mass gain. As a result, the oven-dried trunk measurements were noted as trunk sample biomass (m stem,s), the oven-dried branches as branch sample biomass (m branch,s), and the oven-dried leaves as leaf sample biomass (m leaf,s).

Prior to destructive sampling, all targeted trees were subjected to terrestrial laser scanning with the LiBackpack LiDAR system. This procedure generated a 3D point cloud of the scanned trees. After preprocessing the point clouds data, a 3-Dimensional shape model for each of the selected trees in the destructive sampling area was developed via a tree quantitative structural model. This process allowed precise measurement of tree characteristics including diameter at breast height (DBH), trunk volume, crown area, branch volume, and tree height. Individual tree biomass was then calculated using both general biomass equations and partial volume-to-weight ratio equations developed by Borsah et al. (2025) specifically for this purpose as illustrated in the allometric model development.

Allometric model development

The methodology (Figure 2) comprises allometric model development for field plot biomass estimation, image processing for feature extraction and the application of HybridEnsemble Golden Jackal Optimization (HGJO) for future aboveground biomass prediction.

The trunk, branches, and leaves were treated as separate components for the estimation of individual tree AGB. Equations 2, 3 follow from the basic biomass relation in Equation 1: where ρ is wood density (dry mass divided by volume) and V is volume. Stem and branch volumes were obtained via tree reconstruction. For foliage, Equation 4 estimates leaf biomass by scaling from sampled leaves using either crown area or, equivalently, the product of leaf area index and ground area. Hence, the AGB of each component is:

Where the stem biomass is denoted as B stem (kg), while branch and leaf biomass are represented as B branch (kg) and B leaf (kg), respectively. The sample dry masses of the stem, branch, and leaves are indicated as m stem, s, m branch, s, and m leaf, s (kg). The total volumes of the stem and branches are defined as V stem and V branch (m 3), with the sample volumes for the stem and branch indicated as v stem, s and v branch, s (m 3). For cylindrical samples, the sample volumes are calculated using the formula v stem, s = v branch, s = π r s 2 h s, where r s represents the sample radius (m) and h s is the sample height (m). Additionally, the sampled leaf area is referred to as A sample (m 2), while A crown describes the crown area (m 2). The leaf area index (LAI) is considered dimensionless, and A ground represents the ground area beneath the sampled canopy (m 2), with the relationship A crown = LAI · A ground established. Finally, total AGB per tree is computed as the sum of component biomasses through Equation 5.

Drawing on prior research related to subtropical trees (Lin et al., 2017; Williams et al., 2005), several regression-based model formulations were tested using our dataset, with AGB as the response variable and DBH and height (H) as predictor variables. Log-transformed values of DBH and AGB were utilized, and the best model fit (Table 2) was determined by assessing goodness-of-fit statistics, including Coefficient of Determination (R 2), Root Mean Square Error (RMSE), Akaike Information Criterion (AIC), and Residual Standard Error (RSE). The coefficients’ p-values (P < 0.001) demonstrated the strong statistical significance of the model. Furthermore, biases that arose during the back-transformation of the natural logarithm equation were mitigated by applying a correction factor (CF), calculated using the inverse transformation method (Picard et al., 2012). This CF, expressed as exp(σ 2 /2) with σ equal to the RSE of 0.352, was subsequently implemented.

Image pre-processing and feature extraction

The first phase of the process involved preparing the ALOS-2 PALSAR-2 and Sentinel-2 images captured in 2016, 2018, 2020, 2022, and 2024. This preparation included radiometric calibration to ensure that the images were radiometric consistent and accurate followed by resampling and feature extraction. For the Sentinel-2 images, the preprocessing was performed using version 9.0.0 of Sentinel Application Platform (SNAP) software sourced from the European Space Agency. The images were radiometrically corrected via the Sen2cor package. After this correction, the surface reflectance images were resampled to 20 m spatial resolution to conform with the field plot measurements. From this refined Sentinel-2 imagery, vegetation indices and Tasseled Cap transformed features (Table 3) were extracted for modeling and predicting aboveground biomass. These indices and transformed features offer valuable insights into vegetation health and other characteristics that are essential for assessing AGB.

Additionally, Japan Aerospace Exploration Agency (JAXA) data from the ALOS-2 PALSAR-2 satellite, featuring dual polarization in the L-band with a resolution of 25 m, was utilized for this study. Images captured under both Horizontal transmit/Horizontal receive (HH) and Horizontal transmit/Vertical receive (HV) polarizations, contributed to the estimation of aboveground biomass (AGB). The PALSAR-2 data was preprocessed by JAXA via mosaicking algorithm prepared by Shimada and Ohtaki (2010). After geocoding, resampling, and sub-setting, the backscatter signals of these radiometrically and geometrically corrected images were converted into γ° in decibel values via Equation 6. where γ ° represents the backscatter coefficients, Â represents the image digital numbers, and X is the calibration factor, valued at −83, for both HV and HH polarizations (Shimada et al., 2009).

Spectral bands, vegetation indices, and Tasseled Cap transformed features from optical Sentinel-2 data and backscatter polarization band (viz: HH, HV) from SAR data were used as predictor variables for AGB predictions. Backscatter band math such as HH + HV, HH-HV, HH2+HV2, and HH2-HV2 were computed. Overall, six backscatter predictors from SAR data and twenty-three spectral predictors from Optical data were utilized in this study (Table 3).

Model training and prediction

AGB prediction process

In AGB mapping, no single algorithm consistently outperforms others across different regions, feature sets, or validation frameworks. The HybridEnsemble approach for AGB estimation strategically integrates complementary base-learners using techniques such as performance-weighted averaging or stacking with a simple meta-learner trained on out-of-fold predictions. This method enhances prediction stability across varying forest types and spatial clusters while reducing biases associated with individual models (Luo et al., 2024). When integrating SAR and optical features such as Sentinel-1 backscatter and texture with Sentinel-2 or Landsat red-edge indices, these hybrid models exhibit superior accuracy and robustness, particularly under spatial cross-validation schemes that address spatial autocorrelation (Bui et al., 2024). In the present study, AGB was modelled using a weighted HybridEnsemble that combined four gradient boosting regressors: CatBoostRegressor, XGBRegressor, LGBMRegressor, and an AdaBoostRegressor with a shallow decision-tree base-learner. The dataset was then randomly partitioned into training and testing subsets using an 80/20 split. The ensemble prediction was formulated as a weighted average of the four base-learners via a voting regressor, where the ensemble weights were constrained to be non-negative and were normalized to sum to one. We optimized the HybridEnsemble configuration using Golden Jackal Optimization (GJO) over a continuous search space with a 16-dimensional decision vector encoding: (i) four ensemble weights and (ii) twelve model hyperparameters (three for CatBoost, three for XGBoost, four for LightGBM, and two for AdaBoost). That is a continuous decision vector x ∈ R D with D = 16 variables. The vector in Equation 7 encoded both ensemble weights and base-learner hyperparameters as:

The first four elements x 1, …, x 4 represented raw, nonnegative weights for the four base-learners (CatBoost, XGBoost, LightGBM, AdaBoost). These were transformed into simplex‐constrained weight vector via Equation 8:

Where w ∼ j is the nonnegativity transform applied to raw weight components to prevent negative weights, w ∼ q is the intermediate, nonnegative (unnormalized) weight for base-learner q, computed from the raw decision variable before normalization, x j is the j -th decision variable in x, and w j is the normalized weight assigned to the j -th base-learner in the ensemble.

The remaining elements were used to parameterize each learner as follows: CatBoost was configured with depth, learning rate, and iterations corresponding to x 5, x 6, and x 7, respectively; XGBoost was assigned max depth, learning rate, and number of estimators based on x 8, x 9, and x 10; LightGBM incorporated x 11, x 12, x 13, x 14 for maximum depth, learning rate, number of estimators, and the number of leaves; and AdaBoost utilized x 15, x 16 for the number of estimators and learning rate. Integer-valued parameters were obtained by rounding the respective variables: depth was set as round x 5, iterations as round x 7, max depth as round x 8, number of estimators as round x 10, LightGBM’s max depth as round x 11, LightGBM’s number of estimators as round x 13, number of leaves as round x 14, and AdaBoost’s number of estimators as round x 15. All decision variables were constrained to user-defined bounds L ≤ x ≤ U, matching the bounds list used in the implementation. Where L and U represent the lower and upper bound vectors defining feasible ranges for each decision variable.

GJO minimized the mean K -fold cross-validation root-mean-square error (CV-RMSE) of the weighted VotingRegressor HybridEnsemble. For fold k with n k validation samples, RMSE was calculated through Equation 9: and the optimization objective function was computed through Equation 10.

Where, k∈{1, …,K} is the index of the CV fold, y ^ i k x denotes predictions from the HybridEnsemble instantiated from x after decoding and weight normalization, y i k is the observed/true target value (AGB) for the i -th sample in the validation split of fold k, and n k is the number of validation samples in fold k.

The sixteen decision variables were defined according to the specified bounds. Variables one through four represented the raw ensemble weights for CatBoost, XGBoost, LightGBM, and AdaBoost. Variables five through seven corresponded to CatBoost parameters, including depth, learning rate, and iterations (Supplementary Figure S2). Variables eight through ten denoted XGBoost parameters, comprising max depth, learning rate, and number of estimators. Variables eleven through fourteen specified LightGBM parameters, consisting of maximum depth, learning rate, number of estimators, and number of leaves. Variables fifteen and sixteen indicated AdaBoost parameters, specifically the number of estimators and learning rate. This encoding facilitated the joint optimization of both model capacity and ensemble composition within a single continuous metaheuristic search. Supplementary Table S1 in the Supplementary Material specifies the hyperparameter configuration used for the analysis.

Sensitivity analysis was conducted using a held-out test set and the finalized HGJO model, with the training data used to define perturbation magnitudes and imputation constants. The held-out test set represents 20% of the dataset and is distinct from the training (calibration) dataset used to optimize the HybridEnsemble. First, baseline predictive performance was established by computing the HybridEnsemble root mean square error (RMSE) and coefficient of determination (R 2) on the unmodified test set. Sensitivity was subsequently evaluated through four different methods: one-at-a-time (OAT) feature perturbation, noise robustness, missingness robustness, and ensemble weight sensitivity.

OAT feature perturbation was performed by independently increasing each predictor in the test set by a small additive increment equal to 10% of that predictor’s training-set standard deviation (σ). For each perturbed predictor, the HybridEnsemble was rerun and sensitivity was quantified as (a) the change in RMSE relative to baseline (ΔRMSE) and (b) the mean absolute change in predictions (mean |Δprediction|) compared with the baseline predictions. Predictors were ranked by these sensitivity scores, and the most influential predictors were visualized using horizontal bar charts.

Noise robustness was evaluated through a fast Monte Carlo scheme in which zero-mean Gaussian noise was injected into all predictors in the test set simultaneously. Noise levels were set as fixed fractions of the training set σ (0, 0.05, 0.10, and 0.20), applied per feature to preserve scale differences across predictors. For each noise level, multiple stochastic runs were generated using controlled random seeds, the HybridEnsemble predictions were recomputed, and the resulting RMSE distributions were summarized across runs. Robustness patterns were visualized using both violin plots (with quartile overlays), with the baseline RMSE indicated as a reference line.

Missingness robustness was examined by randomly masking test-set predictor values to simulate incomplete inputs at missing rates of 0, 0.10, 0.20, and 0.30. Missing values were introduced under a random mask and then imputed using the corresponding training-set median for each predictor, reflecting a simple but stable replacement strategy. As with the noise experiment, multiple Monte Carlo runs were executed per missingness level, RMSE was computed for each run, and the distributional shifts in error were visualized using violin plots relative to the baseline RMSE.

Ensemble weight sensitivity was quantified by perturbing each optimized base learner weight (CatBoost, XGBoost, LightGBM, and AdaBoost) by ±20% while holding the other weights fixed. After each perturbation, weights were constrained to be non-negative and renormalized to sum to one before recomputing the weighted ensemble prediction from the stacked base-prediction matrix. The impact of each weight perturbation was summarized as ΔRMSE relative to baseline, and results were presented as a tornado-style plot (Figure 3) to highlight which base learner’s weight uncertainty most strongly affected ensemble accuracy.

Predictive uncertainty was estimated using a nonparametric bootstrap on the HybridEnsemble model. Following the 80/20 train–test partition, the training data were resampled with replacement 100 times, preserving the original training-sample size in each replicate. The HybridEnsemble was refitted to each bootstrap sample and used to predict AGB for the fixed test set. For each test observation, the bootstrap distribution of predictions was summarized by its mean and standard deviation, and 95% prediction intervals were obtained from the 2.5th and 97.5th percentiles. Empirical interval coverage was computed as the proportion of observed AGB values contained within these intervals. Uncertainty was visualized using observed-versus-predicted plots with bootstrap derived error bars, together with the distributions of predictive standard deviation and interval width.

AGB projection process

To forecast the future AGB, annual AGB multi-temporal stacked rasters for 2016, 2018, 2020, 2022, and 2024 were first aligned to a common template grid. For each valid non-zero AGB pixel, a linear model, AGB t = β 0 + β 1 t, was fitted using the observed years as the temporal predictor, where β 1 represents the annual rate of biomass change and β 0 the intercept. The ‘trend_predict_2036’ function in Python 3.10.0 version was designed to analyze time series imagery dataset across specified years and predict future values for a target year, specifically 2036. The function accepted several parameters, including a list of integers representing the years for which data was available and subsequently maps each year to its corresponding data file path. Initially, a mask was used to identify forest areas from the AGB raster of the available years, followed by the calculation of mean biomass values for each year, where valid data existed. Deviations from these mean values were then used to compute the slope and intercept of the linear trend, allowing predictions for the target year to be calculated. Additionally, fitted values were computed for all years using the linear model, and residuals were utilized to calculate the sum of square errors and total sum of squares. The R 2 value was computed to assess how well the model explained the variability in the data. Predicted values, slopes, intercepts, and R 2 values were then written to their respective output files. Moreover, the ‘mean_agb_rasters’ function was designed to compute the mean AGB from a series of raster datasets representing different years. The function iterates over the input files, reads the corresponding raster data for each year, and appends this data to a stack file. Furthermore, it computes the mean AGB across the stack, and finally, the mean values were written to the output raster. This projection assumes that (i) the temporal AGB trajectory at each valid pixel was adequately approximated by a linear trend over 2016–2024, (ii) the annual AGB maps were comparable through time because they were aligned to a common grid and derived from consistently preprocessed inputs, (iii) forecastable pixels were represented by valid non-zero AGB cells, which operationally defined the forested analysis domain, and (iv) abrupt future disturbances, management changes, or non-linear regime shifts were not explicitly modeled.

Model accuracy assessment and SHAP method

Shapley additive exPlanations (SHAP) offers a novel way to examine the interrelationships between individual variables for predicted scenarios and the intricacies of predictive model outcomes (Ghafarian et al., 2022). A helpful technique for classifying effects and dissecting predictions into distinct feature importance (Nohara et al., 2022; Toharudin et al., 2023). It offers additive, game-theoretic feature attributions that are both local, explaining individual pixel or plot predictions, and global, summarizing importance across the entire study area (Hou et al., 2025; Ankah et al., 2025). SHAP attributes any single prediction to input features by computing Shapley values, a game-theoretic allocation of the model’s output relative to a baseline. According to Lundberg and Lee, (2017), the Sharley value for feature i in Equation 11 is expressed as:

This averages the marginal contribution of feature i over all subsets S of features. M is the total number of features, ! is factorial, |S| is the size of subset S, f x (S) denotes the model output when only the features in subset S are known, and f x ∪ i represents the model output when the features in subset S plus feature i are known. The local explanation satisfies additivity as presented in Equation 12: where ∅ 0 represents the expected value of the model output over all input features, ∑ i ∅ i and represents the sum of all contributions of SHAP value for feature i (Lundberg and Lee, 2017; Hou et al., 2025). In practice, local explanations provided by per plot SHAP values quantify how specific features, such as HH, HV, or red-edge reflectance, influence predictions either upward or downward. For global importance, the mean (|SHAP|) ranks features, while SHAP summary (“beeswarm”) and dependence plots uncover nonlinear thresholds and interaction patterns commonly observed in vegetation structure. Remote-sensing studies utilize SHAP in conjunction with boosted trees to interpret the contributions of different SAR and optical indices in driving predictions. This is illustrated through SHAP summary and dependence style analyses in related SAR estimation tasks.

To improve model interpretability, SHAP was implemented after model calibration to quantify the contribution of each predictor to AGB estimates. The SHAP analysis was conducted for the CatBoost, XGBoost, LightGBM, AdaBoost, and optimized HGJO models using the predictor variables listed in Table 3, comprising optical Sentinel-2 bands, vegetation indices, Tasseled Cap features, and ALOS-2 PALSAR-2 backscatter metrics. Global importance was summarized as the mean absolute SHAP value across samples, which ranked predictors according to their overall influence on model output. Local explanations were obtained from per-sample SHAP values, where positive values indicated an increase in predicted AGB relative to the model baseline and negative values indicated a decrease. These local contributions were visualized using SHAP decision plots to illustrate how optical and SAR predictors cumulatively affected individual predictions. This implementation enabled direct comparison of feature-attribution patterns among the component learners and the optimized ensemble, thereby strengthening the interpretability of the modelling framework and helping identify ecologically meaningful drivers of AGB variation.

The performance of the ensemble models for predicting aboveground biomass was evaluated using metrics, including R 2, rRMSE, and rBias, defined in Equations 13–16. where e i and é i are the measured and predicted AGB whiles ê represent the mean of measured AGB and n represent the sample size. To provide inferential support for model comparison, point estimates of predictive performance were supplemented with statistical uncertainty analysis and paired model-comparison tests. For each model, 95% confidence intervals for R 2, rRMSE, and rBias were estimated by bootstrap resampling of the held-out test predictions. In addition, because all models were evaluated on the same test plots, pairwise comparisons between the HGJO model and each base-learner were conducted using paired error statistics computed at plot level. Paired tests were applied to absolute prediction errors, with normality of paired differences assessed prior to test selection. A paired t -test was used when the normality assumption was satisfied, otherwise the Wilcoxon signed-rank test was applied Multiple comparison adjustment was performed using the Holm procedure.

To complement the primary accuracy metrics, two additional diagnostic visualizations were generated from the held-out test set. First, an UpSet plot was used to summarize overlap among models in identifying well-predicted test samples. A test sample was classified as well predicted for a given model when its absolute prediction error was less than 1.0 t ha −1. The UpSet intersections therefore represent the number of held-out test samples that satisfied this criterion for unique or shared combinations of models. Second, a Taylor diagram was constructed using the same held-out test predictions to jointly summarize correlation with observed AGB, the standard deviation of predicted AGB, and RMSE. The Taylor diagram was used as a compact multimetric diagnostic and not as an independent validation analysis.

Golden jackal optimization (GJO)

The GJO exhibited a clear and stable convergence pattern (Figure 4) in terms of cross-validated error RMSE. The best RMSE remained essentially constant across the iterations, indicating that a near-optimal HybridEnsemble configuration was identified at the beginning of the search and was preserved throughout the run. In contrast, the mean RMSE of the candidate population continued to decreased from approximately 47.41 t ha −1 at iteration 0 to about 47.17 t ha −1 by iteration 4, an absolute improvement of approximately 0.24 t ha −1 of RMSE, which was about 0.51%, with a short plateau around iterations 1 to 2 followed by a more pronounced drop between iterations 3 and 4. The progressive narrowing of the gap between mean and best RMSE suggested that the population had increasingly concentrated around the best performing region of the search space, reflecting improved solution consistency and reduced dispersion as optimization proceeded.

Noise, missingness, ensemble weight and OAT feature perturbation sensitivity analysis

Figure 3 demonstrates the impact of noise and missingness on the HGJO model’s predictive accuracy. Noise introduces uncertainty and variability into observed feature values, while missingness reduces the available information for prediction and necessitate imputation which could introduce further errors. When these two sources of perturbation occur, a model faces both corrupt and incomplete inputs, amplifying the challenge of accurate prediction. The sensitivity analysis results indicated that the HGJO model maintained its best predictive accuracy under unperturbed conditions with RMSE approximately 32 t ha −1, but its performance deteriorated systematically as input quality was degraded. When additive Gaussian noise was introduced scaled to the training-set standard deviation, RMSE continued to increase from the baseline to median values of approximately 36 t ha −1 (0.05σ), 38 t ha −1 (0.10σ), and 42 t ha −1 (0.20σ), and the error distributions broadened at higher noise levels, indicating increasing variability in prediction reliability. A stronger sensitivity was observed under missing data perturbations. As the missing rate increased, median RMSE rose more steeply to approximately 37 t ha −1 (0.10), 43 t ha −1 (0.20), and 49 t ha −1 (0.30), accompanied by wider distributions and longer tails, suggesting substantial instability when AGB predictor information content was reduced. Overall, these results showed that the HGJO model was more vulnerable to incomplete predictor information than to moderate measurement error. These results indicate that data completeness and effective handling of gaps are critical determinants of biomass map reliability and temporal comparability in biomass estimation and justify treating missing data effects as a major uncertainty source alongside sensor noise in the modelling process.

The ensemble weight sensitivity analysis (Figure 5) showed that prediction error responded differently to ±20% perturbations depending on which base-learner weight was altered. Overall, the HGJO model performance remained centered near the baseline for all perturbations, but the spread of ΔRMSE revealed which weights the ensemble depended on most. The largest sensitivity was observed for AdaBoost (w_AdaBoost), where ±20% changes produced the widest ΔRMSE range of approximately −0.19 t ha −1 to +0.20 t ha −1. This pattern indicated that small-to-moderate perturbation of the AdaBoost weight could either slightly improve performance (negative ΔRMSE) or noticeably degrade it (positive ΔRMSE), implying that this component exerted the strongest leverage on the ensemble’s error. CatBoost (w_CatBoost) showed the second-highest sensitivity, with ΔRMSE spanning roughly −0.11 t ha −1 to +0.13 t ha −1, suggesting that the ensemble was also meaningfully dependent on the CatBoost contribution, although less strongly than on AdaBoost. In contrast, LightGBM (w_LightGBM) and XGBoost (w_XGBoost) exhibited comparatively narrow ranges both near ±0.05, with XGBoost slightly smaller, indicating that the ensemble error was relatively insensitive to moderate perturbations of these weights. When considered collectively, Figure 6 suggested that the optimized weighting scheme was most critical for the AdaBoost and CatBoost components, whereas the LightGBM and XGBoost weights behaved more like secondary tuning parameters within the tested perturbation bounds. Figure 6 slight asymmetry further suggested that increasing or decreasing a given weight did not produce perfectly mirrored effects, consistent with non-linear interactions among base-learner predictions in the weighted aggregation.

The one-at-a-time (OAT) feature perturbation analysis showed that the HGJO model responded unevenly to independent +10% standard deviation (σ) increases in individual predictors. Two complementary outcomes were examined: the change in RMSE, which indicated the influence of each predictor on overall predictive accuracy, and the mean absolute change in predicted AGB, which reflected the magnitude of output sensitivity (Figure 6). ND45 produced the largest increase in RMSE, followed by TCA, TCW, and GNDVI, indicating that model accuracy depended most strongly on vegetation sensitive optical predictors, particularly tasseled-cap components and greenness indices. Moderate error increases were associated with NDVI, TCB, and TCG, whereas SAVI, EVI, selected red-edge and SWIR bands, and most remaining predictors showed comparatively smaller effects. In contrast, prediction shift rankings identified TCA as the dominant driver of changes in predicted AGB, followed by ND45, GNDVI, TCW, and NDVI. SAR backscatter variables and polarization derived combinations showed measurable but generally weaker effects on both RMSE and prediction shifts. The analysis indicated that optical vegetation related predictors exerted the strongest control on HGJO model performance and output variability.

The combined effect (Figure 7) identifies TCA as the most significant positive contributor across all the models, suggesting that TCA was a robust predictor. In contrast, DVI showed a strong negative impact across models, particularly affecting HybridEnsemble model, LightGBM model, XGBoost model, and CatBoost model. Vegetation indices such as NDVI, RDVI, SAVI, NDI45, and GNDVI contributed positively, although their importance varied among different models. Additionally, several spectral bands, including B04, B06, and B07, also provided positive contributions, with their significance depending on the specific model used. While SAR features had relatively minor contributions compared to optical indices, they still played a role, especially in XGBoost, LightGBM, and HybridEnsemble. The stronger contribution of optical predictors suggests that optical reflectance, red-edge information, and vegetation indices were highly responsive to canopy greenness, chlorophyll content, leaf-area development, and photosynthetic condition, all of which are closely associated with biomass variation in subtropical secondary forests. Although structurally informative, L-band SAR-AGB relationships are often weaker in complex subtropical forests because of signal saturation at medium-to-high biomass, sensitivity to canopy moisture, terrain effects, and mixed scattering from heterogeneous land cover. Ultimately, the combined SHAP values demonstrated that both SAR and optical features contributed to AGB prediction.

Correlation analysis of AGB and predictor variables

The relationship between AGB and multi-source remote sensing predictors was characterized using Pearson’s r coefficients, revealing distinct spectral and polarimetric sensitivities (Figure 8). Vegetation indices derived from the optical dataset emerged as the most robust correlates, with SAVI, NDVI, and TNDVI exhibiting strong positive associations (r = 0.56–0.58, p < 0.001), confirming their direct relationship with vegetation greenness and biomass. Moderate positive correlations were also found for RDVI (r = 0.46), TCG (r = 0.40), and DVI and EVI (both r = 0.33). In contrast, several visible and shortwave bands including B02 (r = −0.58), B3 (r = −0.52), B04 (r = −0.65), and B11 (r = −0.31) showed strong inverse correlation, consistent with spectral darkening in the visible and shortwave range as biomass accumulates. In contrast, SAR metrics showed only weak associations with AGB (r = 0.06 to 0.17), confirming their lower direct sensitivity in this dataset. The predictor matrix also revealed substantial collinearity among optical variables, especially among the red-edge and near-infrared bands (B06, B07, B08, and B08A: r = 0.92–0.99) and among vegetation indices derived from similar spectral information, with DVI and EVI perfectly correlated (r = 1.00) and NDVI, SAVI, and TNDVI also nearly identical. Tasseled-cap variables were more differentiated, with TCB and TCD closely linked to optical reflectance structure, whereas TCW showed contrasting relationships consistent with canopy moisture effects. SAR variables were strongly intercorrelated within their own group but generally weakly related to optical predictors, indicating that they contributed complementary structural information rather than redundant spectral information. Collectively, the heatmap highlighted the synergistic potential of using both SAR and optical features together. While optical data provided strong, direct signals for biomass, SAR data contributed unique, orthogonal information that could enhance the biomass prediction.

Model explainability

SHAP decision plots (Figure 9) were utilized to visualize how each feature contributes to the model’s output for individual predictions. This XAI demonstrated the cumulative effect of features on the predicted AGB from top to bottom. In CatBoost, the top influencing features included B04, TCA, B02, GNDVI, TNVDI, and B06, with a wide spread of SHAP values indicating strong influence and variability among samples. This suggested that the most influential features were primarily optical bands, while SAR features, such as HH2+HV2 and HH-HV, ranked lower. Consequently, the model relied heavily on optical and vegetation indices for AGB prediction. Similarly, XGBoost identified important features like B04, TCA, NDVI, GNDVI, B02, B06, and NDI45, showing a broader spread in SHAP values for the top three features, which implied a strong and variable effect on predictions. Although polarimetric SAR features like HH + HV, HH-HV, and HH2-HV2 were present, they were less impactful compared to spectral indices, reinforcing the model’s prioritization of optical and vegetation features. LightGBM followed suit, with a slightly different order of feature importance but still highlighting major contributors such as B04, TCA, GNDVI, NDVI, and B02, where optical bands and vegetation indices dominated. Notably, GNDVI ranked higher than in CatBoost and XGBoost, while polarimetric SAR features remained consistently lower but moderately influential. Similar to LightGBM, HybridEnsemble identified B04, TCA, GNDVI, and NDVI as important features showing a broader spread in SHAP values. In contrast, AdaBoost exhibited a more concentrated and narrower feature impact, with top features like B04, TCA, and NDI45 dominating, while polarimetric SAR features were among the least influential. Overall, across all models, spectral features related to specific bands (viz: B04, B02), vegetation indices (viz: NDVI, GNDVI, SAVI), and the Tasseled Cap transform (TCA) were prominent in defining AGB predictions, while polarimetric SAR features played a supporting but less influential role. This consistent ranking across models strengthened the reliability of these features for AGB estimation, and the transparency provided by SHAP plots confirmed that the models made decisions based on ecologically meaningful variables, enhancing their value in XAI.

The UpSet plot (Figure 10a) displays the unique contribution of each model, specifically the predictions made by each model. The horizontal bars represent the total number of correct predictions for each individual model, while the vertical bars indicate the size of each intersection, showing the number of samples accurately predicted by specific combinations of models. The HGJO model uniquely predicted seven cases, followed by CatBoost with six predictions and XGBoost with two predictions. In contrast, AdaBoost and LightGBM each accounted for only one unique prediction. Model overlap was limited, with a single shared intersection evident between HybridEnsemble and XGBoost. The analysis highlighted that most predictions were uniquely captured by the HGJO model, underscoring the importances of hybrid ensemble optimization approach for achieving comprehensive coverage.

The Taylor diagram (Figure 10b) compares the performance of five models against observed data using three key metrics (viz: correlation coefficient, standard deviation, and RMSE). The correlation coefficient indicates how well the predictions match the observations. The standard deviation reflects the spread of predictions, and the RMSE represents the distance from the observed values. In a model-by-model analysis, HGJO achieved the best performance, characterized by the highest correlation with the observed data, the best match to observed variability, and the lowest RMSE. The HybridEnsemble point lay closest to the reference, indicating the most favorable joint performance in terms of agreement with observations, variance representation, and lower RMSE relative to the individual gradient-boosting learners. Additionally, CatBoost, XGBoost, and LightGBM emerged as strong performers, positioned close to HybridEnsemble. However, AdaBoost underperformed, as it recorded the lowest correlation, the lowest standard deviation, and the highest RMSE. The diagnostics indicated that HGJO significantly enhanced predictive performance relative to single-model implementations, as it yielded predictions that were not only more consistent with observed patterns but also exhibited reduced error rates.

AGB model prediction and accuracy assessment

The analysis revealed varying performances among the individual gradient boosting models and the HGJO model for predicting AGB (Figures 11, 12; Supplementary Figure S1). For CatBoost, the points generally followed the 1:1 line, indicating a reasonable correlation between predicted and measured AGB, with an R 2 of 0.816, which explained approximately 81.6% of the variance, reflecting a good fit, though not the best among the models. Its rRMSE was 16.51%, suggesting a moderate level of prediction error, while the rBias of −3.39% indicated a slight tendency to underestimate AGB. Additionally, the residuals spread above and below zero exhibiting a slight upward trend at higher AGB values. In contrast, XGBoost showed points widely scattered around the 1:1 line, demonstrating less predictive performance, especially at higher AGB values. This model achieved an R 2 of 0.802, explaining 80.2% of the variance, coupled with a lower rRMSE of 17.13%, indicating less accurate predictions and underestimation reflected in its rBias of −3.84. The residuals were evenly distributed around zero, showing minimal systematic error and a narrower spread, suggesting consistent predictions. LightGBM’s performance was slightly higher than XGBoost, with points closely following the 1:1 line, though it exhibited slightly less scatter, particularly at the extremes. It produced an R 2 of 0.804, explaining 80.4% of the variance, and its rRMSE of 17.05% was slightly smaller than that of XGBoost, while its rBias of −3.61% indicated minimal bias, with evenly distributed residuals. AdaBoost, however, displayed more scatter and deviation from the 1:1 line, particularly at higher AGB values, achieving the lowest R 2 of 0.731, which explained only 69.2% of the variance, alongside the highest rRMSE of 19.9%, indicating less accurate predictions and a slight underestimation with an rBias of −3.78, while its residuals showed more spread and a slight upward trend at higher AGB. Finally, the HGJO model exhibited points very tightly clustered around the 1:1 line, indicating the best fit among all models, with minimal scatter even at higher AGB values, achieving the highest R 2 of 0.821, which explained 82.1% of the variance. Its rRMSE was 16.30%, reflecting very low error, and its rBias of −2.50 indicated minimal bias, while the residuals were evenly distributed around zero with no discernible trend, showcasing the narrowest spread and highly consistent predictions across the entire AGB range. Bootstrap analysis showed that its R 2, rRMSE, and rBias estimates were associated with 95% confidence intervals that were narrower and/or more favorable than those of the base-learner models (Figure 13a).

From the analysis, HGJO model achieved the highest R 2 values of 0.821 demonstrating its ability to explain most variance in the measured AGB. CatBoost with R 2 of 0.816 followed closely behind, while LightGBM (R 2 = 0.804), XGBoost (R 2 = 0.802) and AdaBoost (R 2 = 0.731) ranked lower. Additionally, HybridEnsemble exhibited the lowest rRMSE of 16.30%, signifying the most accurate predictions, whereas AdaBoost had the highest rRMSE of 19.97%, reflecting less accurate predictions overall. Pairwise statistical comparison based on plot-level absolute errors further showed that HGJO accuracy did not differ significantly from that of CatBoost (p adj = 0.471), XGBoost (p adj = 0.333), or LightGBM (p adj = 0.471), whereas it significantly outperformed AdaBoost (p adj = 0.000; Holm-adjusted p<0.05). Thus, the inferential analysis suggests that the numerical advantage of HGJO was statistically confirmed only relative to AdaBoost, while differences among HGJO, CatBoost, XGBoost, and LightGBM remained comparatively small (Figure 13b). Although all models displayed low rBias, HGJO model was the closest to zero, indicating minimal systematic error. Overall, the HGJO model outperformed all individual models in predicting aboveground biomass, and its superior performance was attributed to the optimized ensemble approach that combined the strengths of multiple gradient boosting algorithms, including CatBoost, XGBoost, LightGBM, and AdaBoost. This resulted in higher accuracy, lower error, and minimal bias, making it the most suitable model for AGB prediction in the study area.

The AGB quartile plots (Supplementary Figure S3) evaluated the ability of the five models to reproduce measured biomass across the full AGB gradient. Comparatively, all models predicted the lowest and highest AGB quartiles more accurately than the intermediate quartiles, indicating stronger performance at the distribution extremes than in the mid-range. Among these models, HGJO showed a better balance across quartiles, particularly by reducing misclassification in the middle biomass ranges relative to the base-learner individual-boosting models. In contrast, AdaBoost performed the least consistently, with particularly low accuracy in the intermediate quartiles, indicating a stronger bias toward lower biomass classes. Although all models captured extreme AGB conditions reasonably well, the analysis showed that the HGJO model provided the most balanced quartile-level predictive performance.

Aboveground biomass mapping

The AGB maps shown in Figure 14 represent model predictions for 2024, the most recent observed year in the study series. It shows the spatial distribution maps of predicted AGB in tons per hectare across Hong Kong, employing the five different models. CatBoost (Figure 14a) demonstrated a moderate spread of high AGB values especially in the central and northern districts. In contrast, XGBoost (Figure 14b) predicted higher AGB values across several districts, especially in the northern and eastern regions. This model appeared more sensitive to high AGB zones, possibly identifying biomass-rich areas. LightGBM (Figure 15c) exhibited a similar pattern to XGBoost but with slightly less intensity in high AGB zones, resulting in a smoother distribution with fewer extreme values. AdaBoost (Figure 14d) predicted the lowest AGB overall, with most areas depicted in green shades less than 146 t ha −1, indicating sparse and less pronounced high AGB zones, which suggests that it underestimated AGB, particularly in known biomass-rich areas. Lastly, the HybridEnsemble map (Figure 14e) integrated the strengths of individual models, resulting in balanced predictions with clear spatial differentiation. It effectively captured high AGB areas, while urban districts were mapped within lower classes. This hybrid ensemble approach with GJO reduced instances of over and under prediction, thereby offering a more realistic representation of spatial distribution.

The histogram compares the mean AGB in tons per hectare predicted by each model across Hong Kong’s districts. The HGJO model consistently predicted the highest mean AGB across most districts, particularly in rural and forested areas such as North, Tai Po, Yuen Long, Sai Kung, and the Islands. This finding suggested that the HGJO model effectively captured more biomass information by leveraging the strength of multiple algorithms. In contrast, AdaBoost predicted the lowest mean AGB in nearly every district, thereby confirming its conservative nature and tendency to underpredict. Meanwhile, CatBoost, XGBoost, and LightGBM exhibited intermediate AGB values, Notably, rural and forested districts demonstrated higher mean AGB across all models, with the HybridEnsemble model amplifying these differences. Conversely, urban districts, including Central and Western, Kowloon City, Sham Shui Po, and Yau Tsim Mong, exhibited lower mean AGB, with AdaBoost recording the lowest values. Ultimately, the Hybrid Ensemble provided the most differentiated and realistic predictions.

Future biomass prediction

The HybridEnsemble model with GJO, identified as the best prediction model, was employed to predict future aboveground biomass. Figure 15a displays the prediction trend of mean AGB in tons per hectare from 2016 to 2036. The mean AGB began at approximately 140 tons per hectare in 2016. A slight dip occurred in 2018, followed by a steady increase through 2020, 2022, and 2024. By 2024, the mean AGB had risen to approximately 160 t ha −1. The forecast from 2024 to 2036 projected a continued upward trend, with mean AGB expected to reach about 210 tons per hectare by 2036. This projection was based on a second-order polynomial fit which is a statistical extrapolation of the observed AGB time series and should not be interpreted mechanistically. It accounted for both linear and nonlinear growth patterns. The polynomial trend suggested that the rate of increase might slightly accelerate during the forecast period. This projected increase is an extrapolation of the observed trend under the assumption that recent dynamics persist. The standard deviation was relatively wide, particularly in the forecast period of 2036, which implied increasing uncertainty in future predictions. This uncertainty could stem from climate variability, land use changes, model limitations, or unforeseen disturbances. Despite these uncertainties, the upward trend suggested ongoing forest growth in the subtropical forest, with the increasing AGB implying greater carbon storage potential.

The violin plots in Figure 15b display the AGB distribution values for each year. The distributions were generally unimodal but exhibited some skewness and spread, especially in 2036. The range of AGB values increased over time, with the widest distribution appearing in 2036, indicating greater heterogeneity in the predicted biomass. Additionally, the median and interquartile range shifted upward over time, consistent with the increasing mean AGB observed in Figure 15a. The 2036 distribution also showed more extreme values, both low and high, reflecting increased model uncertainty or possible changes in forest dynamics.

Figure 15c visualizes the distribution of annual changes in AGB for each time interval and the forecasted period from 2024 to 2036. From the analysis, the median biomass change (ΔAGB) from all intervals was found to be close to zero, which indicated that most plots experienced modest changes in biomass from year to year. Specifically, from 2016 to 2024, the interquartile range (IQR) was relatively narrow with most changes occurring between −50 and +50 t ha −1. Although there were some positive and negative outliers, most plots clustered near zero. In contrast, the forecast period from 2024 to 2036 revealed a much wider distribution, as the IQR expanded significantly, featuring more extreme positive changes of up to +500 t ha −1 and some negative changes down to −200 t ha −1. The forecast period exhibited positive skewness, with a greater number of plots experiencing large biomass gains compared to losses. This suggests that future years would likely witness greater variability in biomass changes with some plots accumulating much more biomass than others.

Uncertainty analysis

To quantify uncertainty, bootstrapping was applied to all error components expressed as standard deviations, and the mean values were considered. Figure 16a compares observed AGB with the HybridEnsemble predictions and their bootstrap derived from 95% prediction intervals. The points clustered closely around the 1:1 line, indicating that the HybridEnsemble model reproduced variability in plot-level AGB reasonably well across the full range of biomass. Most observations fell within relatively narrow 95% prediction intervals, which suggested that the model generated stable predictions for the majority of plots. Slight deviations from the 1:1 line at the upper end of the AGB range implied a tendency to under or overestimate biomass in a few high biomass stands, but no strong systematic bias was visually apparent.

The spatial map (Figure 16b) reveals significant spread of uncertainty. It utilizes color coding to represent different levels of uncertainty. Dark blue indicates the lowest uncertainty (0–5 t ha −1), blue represents low to moderate uncertainty (5–10 t ha −1), yellow denotes moderate to high uncertainty (10–15 t ha −1), and green highlights the highest uncertainty (15–20 t ha −1). Key observations revealed that urban districts, such as Kowloon City, Kwun Tong, and Yau Tsim Mong, predominantly displayed dark blue and blue areas, signifying lower uncertainty in AGB predictions. This lower uncertainty is likely attributed to more homogeneous land cover, consisting of urban and built-up environments, which offered less vegetation and more consistent remote sensing signals. Conversely, forested districts like Tai Po, North, Sai Kung, Yuen Long, and Islands showed a noticeable presence of green and yellow, indicating higher uncertainty concentrated in these areas due to the complexity and heterogeneity of vegetation, topography, and mixed land cover types. Areas such as Sha Tin and Tsuen Wan exhibited mixed uncertainty levels, with patches of moderate to high uncertainty interspersed with lower uncertainty, reflecting a mosaic of urban, suburban, and natural landscapes that increased prediction variability. Overall, the uncertainty in AGB predictions is spatially variable across Hong Kong, with the highest values in forested and rural districts and the lowest in urban areas.

Understanding future biomass change dynamics

Across Hong Kong’s administrative districts, the mean AGB values in t ha −1 derived from the HybridEnsemble predictions spans from moderate-to-high gradient that mirrors the city’s forest structure and urbanization pattern (Figure 14f). These district levels mean AGB values generally cluster between 100 t ha −1 and 165 t ha −1, with a few highly urbanized districts nearer 70–90 t ha −1, a pattern that is consistent with independent measurements for Hong Kong’s forest types. Specifically, the territory’s secondary forests AGB have been reported to be around 164 t ha −1, with remnant old-growth Fung Shui Woods reaching approximately 346 t ha −1 (The Conservancy Association, 2025). This indicates that districts with more mature, protected stands occupied the upper tail of our distribution, while built-up districts fell lower. Recent handheld mobile laser scanning further differentiated an average AGB of 305.7 t ha −1 in old-growth plots from secondary plots, reinforcing the interpretation that our district means reflect a landscape dominated by secondary and mixed-age stands rather than extensive old-growth stands (Mak et al., 2025). When benchmarked against old-growth subtropical forests in South China, such as Dinghushan’s evergreen broadleaf forest, our AGB values were below their 223 t ha −1 average (Lin et al., 2012). At a broader scale, China’s latest 30 m national product indicated a national mean AGB density of 123 t ha −1, placing many Hong Kong districts at or above the national average, feasible given the city’s productive evergreen broadleaf hillslopes, orographic rainfall, and protected country park network (Cai et al., 2025).

In this study, we employed agreement and bias assessment from Bland-Altman graph, and relationship mapping from Hexbin graph, to gain a multi-faceted understanding of the future biomass change dynamics in the study are. The Bland-Altman graph (Figure 17a) compares the difference in AGB between 2024 and 2036 against the mean AGB of those 2 years for each plot. The blue line represents the mean difference of +49.99 t ha −1, while the red dashed lines denote the limits of agreement (LoA), which are −76.65 and +176.63 t ha- 1. On average, plots were predicted to gain approximately 50 t ha −1 in AGB between 2024 and 2036. While most data points fell within the limits of agreement (LoA), some outliers were observed, particularly at higher mean AGB values. The difference in AGB increased with higher mean AGB, indicating that plots with greater initial biomass would accumulate more biomass over time. Furthermore, the spread of differences also increased with higher mean AGB, suggesting that variability in biomass change was greater in high-biomass plots. The analysis revealed that biomass gains were not uniform across the landscape, with plots possessing higher baseline AGB likely to experience greater increases, possibly due to favorable site conditions.

The Hexbin graph (Figure 17b) shows the relationship between baseline AGB in 2024 and the change in AGB from 2024 to 2036. The color intensity represents the log10 (count) of data points in each hexbin, highlighting areas of high density. Clear positive relationship was observed, as plots with higher baseline AGB in 2024 showed larger increases in AGB by 2036. The highest density of points occurred around a baseline AGB of 100–200 t ha −1 and a change of 0–100 t ha −1, indicating that most plots started with moderate biomass and experienced moderate gains. However, some plots with low baseline AGB exhibited negative changes, reflecting biomass loss, while others with high baseline AGB recorded substantial gains of up to +400 t ha −1. Additionally, a secondary cluster emerged at higher baseline AGB values (greater than 250 t ha −1), suggesting that certain forest patches may have been particularly productive. Although initial biomass proved to be a strong predictor of future biomass gain, considerable variability was observed, especially at higher baseline values, with some plots potentially experiencing losses.

The distribution of year-to-year changes in AGB indicated that variability in biomass changes increased significantly during the forecast period from 2024 to 2036, suggesting that some forest plots were expected to experience much larger gains or losses in biomass compared to previous years. The Bland-Altman graph demonstrated a positive association between the difference in AGB between 2036 and 2024 and the mean AGB of those years, showing that plots with higher average biomass tend to have larger increases, and the spread of differences grows with increasing mean AGB. Similarly, the Hexbin plot supported this finding by showing a strong positive relationship between baseline AGB in 2024 and the change in AGB by 2036. Most plots with higher initial biomass were predicted to have greater gains, although considerable variability remained, including some plots experiencing biomass loss. This combined effect highlighted that not all forest areas would respond similarly over time, as some would see substantial biomass gains while others might stagnate or even lose biomass, particularly during the forecast period when model uncertainty and ecological variability were greatest. Given that both the Bland-Altman and Hexbin graphs indicated that high-biomass plots were likely to contribute most to future biomass increases, management strategies should have prioritized these areas for conservation, monitoring, and carbon sequestration efforts. Overall, the analysis revealed that future biomass accumulation in the forest would be highly heterogeneous, with plots of higher initial biomass likely to experience greater increases, yet the range of possible changes remained wide, especially in the forecast period, underscoring the importance of understanding forest dynamics for informed management decisions.

Research limitations and future work

The proposed framework demonstrated strong predictive performance within the sampled study domain. However, future work should consider implementing spatial block cross-validation and, where feasible, external validation across independent forest regions to assess spatial transferability and generate more conservative, operational estimates for wall-to-wall mapping. Furthermore, integrating mechanistic constraints, such as seasonal composites and multi-year forecasting, could enhance predictive accuracy and provide valuable insights into the future biomass change dynamics. Finally, while uncertainty was quantified and visualized, the sources of uncertainty were not fully delineated into categories, such as field measurement error, model uncertainty, and preprocessing uncertainty. Future work could explore these areas to enable attribution of uncertainty to specific components of the workflow.

This study presents an integrated HybridEnsemble model with golden jackal optimization for AGB estimation and forecasting, and seeks to evaluate its predictive performance against individual base-learners, including CatBoost, XGBoost, LightGBM, and AdaBoost via the synergistic application of XAI and active and passive datasets. It applied uncertainty analysis to the predicted AGB maps to assess the reliability of biomass estimates. Key findings include a performance ranking, with HGJO model (R 2 = 0.821) emerging as the best model, followed by CatBoost (R 2 = 0.816), LightGBM (R 2 = 0.804), XGBoost (R 2 = 0.802), and AdaBoost (R 2 = 0.731), with HGJO model showing the most balanced residuals pattern due to its ability to integrate the strengths of the individual gradient boosting algorithms (CatBoost, XGBoost, LightGBM, AdaBoost) to enhance prediction. All models demonstrated high accuracy at the lowest and highest AGB quartiles, while HGJO model showed improved balance in the middle quartiles. XAI analysis revealed that optical predictors were strongly correlated with AGB, consistently yielding high SHAP values, while SAR features provided additional context-specific insights. A positive trend in AGB over time was observed in the subtropical forest, with spatial uncertainty varying across the study area. Higher uncertainty was noted in forested regions compared to urban districts. The analysis also revealed that future AGB accumulation in the subtropical forest will be highly variable and strongly dependent on initial biomass conditions. High-biomass plots are likely to see the greatest gains, but the range of possible outcomes is wide, especially in the forecast period. These findings underscore the capabilities of explainable HybrideEnsemble modeling approach coupled with golden jackal optimization for AGB prediction, providing valuable insights for future biomass mapping and monitoring efforts.