Assessing forest aboveground biomass dynamics over long time periods and at continental to global scale requires reference observations such as forest inventory field plots in combination with satellite Earth observation. To overcome the low availability and quality of reference data, we used the spaceborne Geoscience Laser Altimeter System (GLAS) onboard the Ice, Cloud, and land Elevation Satellite (ICESat)⁴⁸, which operated from 2003 to 2010 and acquired millions of Light Detection and Ranging (LiDAR) footprints, providing measurements of canopy height and other biophysical metrics relating to canopy structure that are highly correlated with aboveground biomass⁴⁹. The Global Ecosystem Dynamics Investigation (GEDI) LiDAR instrument²² onboard the International Space Station that commenced operation in 2019 uses similar technology to the GLAS/ICESat instrument but with smaller footprints and much denser coverage (narrower spacing between footprint locations along the orbit), providing a dense network of training and validation data on forest canopy height.

Based on a machine learning algorithm, L-band Synthetic Aperture Radar (SAR) backscatter image mosaics acquired by ALOS PALSAR-1/PALSAR-2⁴⁸ and optical multispectral Landsat-derived percent tree cover maps²⁹ were used as joint predictors to extend the canopy height obtained from GEDI to canopy height maps across the whole of Africa. Regional maps of aboveground biomass density derived from airborne LiDAR over a range of biomes in Africa were then used to derive an empirical model that estimates aboveground biomass density as a function of canopy height. The LiDAR-based aboveground biomass footprint estimates from LiDAR were used to train the machine learning model. The model was then used to produce annual Africa-wide maps of aboveground biomass density and its standard deviation (SD) at 100 m pixel spacing for the period 2007 to 2017. The maps were validated using a large independent dataset of field plot measurements across the continent. The aboveground biomass density maps and associated standard deviations were used to estimate the African aboveground woody biomass stock and its annual changes (with confidence intervals based on the uncertainty characterization described in Supplementary Materials) with the aim of contributing to improvement of carbon inventories, understanding trends, and testing whether the aboveground biomass change rate has increased, reduced or changed sign over the period 2007-2017.

Datasets

Spaceborne LiDAR from GEDI. The GEDI L2B Canopy Cover and Vertical Profile Metrics product (version 1), available from the NASA/USGS Land Processes Distributed Active Archive Center⁵⁰, was collected from April 2019 to June 2019 (Fig. 3a). The LiDAR metrics estimated by this product are representative of a 25 m diameter footprint on the ground. Version 1 of the product has a geolocation error of approximately 15-20 m, so can be difficult to use with moderate spatial resolution imagery such as Sentinel-1/-2 or Landsat (i.e. 10-30 m spatial resolution pixels).

Airborne laser scanning (ALS). This study uses gridded LiDAR-derived aboveground biomass density maps based on airborne LiDAR acquired in 2016 at 4 different sites in Gabon (Lope, Mabounie, Mondah, and Rabi)^51,52, and in 2015-16 at 2 sites in Kenya (Taita Hills and Maktau)^53,54. The canopy height to AGBD model was developed using 50% of the pixels from these datasets, while the remaining 50% were used to validate the aboveground biomass product.

Field plot data of aboveground biomass density. We collated a dataset consisting of 10,837 aboveground biomass density reference field plots across Africa (see Table S2 in supplementary material) to be used as an independent validation dataset for the 2017 map (Fig. 3). The assessment was limited to 2017 due to the lack of re-measured plot data. Field plots were measured in different years, mainly between 2000 and 2017, with the vast majority before 2010. The plot data are from different national forest inventories and research projects, so have various plot designs, sizes and shapes. Since most of the plots do not have accurate spatial coordinates due to licensing restrictions, we cannot use them directly to validate the 100 m resolution pixel maps. Instead, we followed the approach described by Santoro et al.³¹ and Araza et al.¹⁵ in which the plot values were first adjusted to minimize the temporal and areal mismatches between the plot and map estimates of aboveground biomass density. The adjustment was necessary because of the uneven spatial distribution of the reference samples, the variety of plot sizes used, variations in field survey methods, and used allometric equations to estimate biomass of the plots⁵⁵. The plots were mostly smaller than 1 ha and often represented only a small fraction of the area covered by a 1 ha biomass map pixel. To reduce the effect of random errors caused by different resolutions of the reference dataset and the biomass map, we aggregated the map and the plot data to 0.1° grid cells¹⁵. This procedure yielded 463 grid cells with reference aboveground biomass density values for validation (Fig. 3b). The standard deviation associated with each of these values was estimated by accounting for the principal plot measurement error sources (as described in Araza et al.¹⁵).

ALOS PALSAR/ALOS-2 PALSAR-2 radar image mosaics. JAXA’s annual mosaics of L-band SAR HH and HV polarised backscatter (γ⁰) were based on ALOS PALSAR from 2007 to 2010 and ALOS-2 PALSAR-2 from 2015 to 2017^56,57. No mosaics are available for 2011 to 2014. The PALSAR-2 mosaics for 2018 to 2020 are available, but the pre-processing and geolocation approaches are different from the previous mosaics, so they were excluded from this analysis. The mosaics are a calibrated, 16-look, re-projected, orthorectified and slope-corrected product with 25 m pixel spacing, to which a de-striping process has been applied^22,31,33. The PALSAR and PALSAR-2 mosaics were normalised to reduce artefacts and to ensure temporal consistency of the radar backscatter signal, allowing the same trained model to be used for the whole time series. Artefacts in the mosaics usually result from changes in moisture conditions between image acquisitions, which affects the backscatter, or appear in the pre-processing due to inadequate calibration and/or topographic corrections. We followed a similar approach to that described in⁵⁸, but instead of superpixels used a circular moving window 100 pixels in diameter (~2,000 ha). This normalises the PALSAR/PALSAR-2 imagery to a common baseline based on the mean and standard deviation of backscatter of the PALSAR mosaics (2007-2010). Implicit in this procedure is the assumption that continuous changes with scale larger than 2,000 hectares did not occur from year to year. Normalizing the images at this large scale also tends to ensure that local changes due to disturbances and vegetation growth are preserved. The analysis used both HH and HV polarisations and two additional metrics, the Cross-polarisation Ratio (CpR = HH⁄HV) and the Radar Forest Degradation Index (RFDI = (HH–HV)⁄(HH+HV))⁵⁹.

Percent tree cover data. A 30 m Landsat-based map product of percent tree canopy cover for the year 2000 and annual tree cover loss estimates for the period from 2000 to 2017²⁹ was used to generate annual percent tree cover maps for each year. For the year 2007, all pixels detected as forest cover loss were set as 0% percent tree cover, while pixels detected as having forest cover loss in previous years (i.e. from 2000 to 2006) were set to “no data”, as we have no information on regrowth after the disturbance. The canopy height and aboveground biomass predictions for these “no data” pixels were performed using only PALSAR as a predictor variable (see Modelling Framework). We repeated this process for all the years within our study period (2007-2010 and 2015-2017). PALSAR data and Landsat percent tree cover datasets were mosaicked and co-registered to generate two stacks of predictor datasets, with 50 m and 100 m pixel spacing, respectively. Woody vegetated areas were defined as pixels with equal or above 1% tree cover.

Modelling framework

GEDI footprint selection and clustering. We used the maximum footprint height in a footprint, provided by the GEDI L2B footprint product, as reference canopy height metric, but performed a filtering process to select only the highest quality footprints for training and validation purposes. Coverage footprints were excluded due to their lack of laser light penetration in dense forests and only footprints acquired by the full power beams were used. Only night acquisitions (solar elevation < 0°) with a beam sensitivity greater than 95% were retained. Low quality footprints, as indicated by the L2B quality assurance layers, were discarded, as were footprints with canopy height above the 99.9^th percentile of the initial set. The Copernicus Global Land Cover dataset⁶⁰ was used to exclude footprints in non-vegetated classes. We used the 11 forest classes and the shrubland class for this purpose⁶⁰. After the filtering, approximately 1.8 million footprints were available, distributed across the African continent (Fig. 3). The GEDI footprints were grouped into 4-footprint clusters along the track direction, in each of which the top canopy height values (RH₁₀₀) were averaged to correct for sampling and geolocation errors. This provided the main reference data for training and validating our canopy height model. The average by cluster increases the sampled area of our reference unit from 0.05 ha (1 footprint) to 0.2 ha (4 footprints) and has been demonstrated to increase accuracy when training models with spaceborne LiDAR footprints⁶. The larger sampled area helps to average out various errors typical of small sampling units (e.g. small inventory plots), such as sampling error and geolocation error. Only clusters with 4 consecutive footprints were used.

Canopy height modelling. We used a non-parametric machine learning Random Forests (RF) regression algorithm⁶¹, following the same framework as in⁶² and ⁶³, to generate a canopy height model (CHM) and its associated error at pixel level using PALSAR/PALSAR-2 radar backscatter and Landsat Percent Tree Cover for the same years as predictors. We used 100 trees for each RF model run. Generation of the canopy height model uses two different resolutions (Fig. 4): (i) Remote sensing signatures were extracted from the four 50 m pixels that overlap the 0.2 ha area corresponding to each cluster of four GEDI footprints and averaged; these four pixels are equivalent to the area of 1 ha. (ii) The model training and the prediction output are at 1 ha pixel size, i.e. 100 m by 100 m.

The conversion from 50 m to 100 m was needed because geolocation errors in the GEDI footprint products (15–20 m for version 1 and around 10 m for version 2) prevent small (< 50 m) pixels being adequately matched to single GEDI footprints. Additionally, as with small forest inventory plots, a large tree within the footprint biases the canopy height and makes the value unrepresentative of the canopy height within 1 ha. This two-scale scheme follows Saatchi et al.⁵ and Baccini et al.⁶, and aims to average out these errors by combining several footprints. We assume that 4 footprints located within 1 ha (four 50m x 50m pixels) are more representative of the actual canopy height than 1 or 2 footprints located within a 100 m pixel.

The GEDI canopy height dataset was randomly partitioned into 2 datasets: 10% for training the CHM model (36,944 GEDI clusters) and the remaining 90% (328,243 clusters) for independent validation. Training was performed within a jack-knife/k-fold framework in which the training reference data (10% of the original dataset) were spatially partitioned into k subsamples with k = 10. The jack-knife / k-fold allocation was applied spatially to avoid over-optimistic accuracy metrics due to the possible spatial autocorrelation of the training data⁶⁴. We used a very large dataset for independent validation (i.e., 90% of the original dataset) to avoid a large proportion of the validation dataset being autocorrelated with the training data. Hence, we divided our training dataset (36,944 GEDI clusters) over the whole of Africa into 10 regions with approximately the same number of footprint clusters. A single spatial subsample was kept aside as the validation data for testing the RF model, while the remaining k-1 spatial subsamples were used as training data. The cross-validation process was then repeated k times. Thus, 10 canopy height maps were generated. The mean value of all 10 predictions for each pixel was then used as the final canopy height estimate, and their standard deviation as the prediction error. We also propagated the GEDI footprint height measurement error, the sampling error and the error arising from the temporal difference between GEDI footprints and the satellite imagery to estimate the total error of our canopy height maps (see Supplementary Material).

We trained 2 different random forest (RF) models to predict canopy height. The first one (main model) used annual percent tree cover and SAR backscatter data, while the second one only used SAR backscatter data. Our main canopy height model uses L-band PALSAR/PALSAR-2 radar imagery together with the Landsat percent tree cover product. This model was applied for all pixels from the year 2000 to the given year that were undisturbed according to Hansen et al.²⁹. However, for the pixels detected as disturbed (i.e., forest cover loss), a second canopy height model based only on PALSAR/PALSAR-2 was fitted and used for the years after the disturbance, as this is better suited to detecting any recovery (the Landsat Percent Tree Cover product assumes no regrowth after a forest loss event).

Conversion of canopy height to aboveground woody biomass density. We developed a single empirical linear model to estimate the square root of aboveground biomass density as a function of canopy height. The model was based on six airborne LiDAR aboveground biomass density maps covering closed-canopy moist tropical forest, mangrove, montane forest, drylands and open savannas (see supplementary material).

The square root of AGBD is derived to reduce heteroscedasticity⁶⁵. Conversion to AGBD requires correction of the inherent negative bias in the square root transformation⁶⁶ using the ratio estimator proposed by Snowdon et al.⁶⁷:

Boveground biomass density change analysis. For a given pixel, we identified a significant biomass loss between times t1 and t2 if AGBD_t1–SD_t1 > AGBD_t2 + SD_t2, while there is a significant biomass gain if AGBD_t1 + SD_t1 < AGBD_t2–SD_t2, where the AGBD values at t1 and t2 are AGBD_t1 and AGBD_t2 respectively, and their corresponding SD values are SD_t1 and SD_t2. Changes between consecutive years were considered significant when the error bounds on the estimates of AGBD at the two times did not overlap, otherwise the measured change was considered insignificant. For pixels showing a significant annual gain or loss, we calculated the cumulative significant change over the time series. If this was negative at the end of the period (i.e. 2017), we accepted it as significant biomass loss, and assign this loss value to the year with the highest loss (when significant loss is detected in two or more consecutive pairs of images). Biomass losses usually result from forest cover loss events, but could also arise from other causes, such as natural or anthropogenic fires. A similar approach was taken to biomass gains. However, small gains and losses are hard to detect using this method, and biomass gains were very scarce in forest areas with high aboveground biomass density (i.e. mature forest). Either change is negligible in such areas due to a balance between mortality and growth, or it simply cannot be detected due to insufficient sensitivity of the Earth observation signal when biomass is high.

For pixels that were undisturbed (i.e. with positive cumulative significant biomass change or non-significant biomass change), the slope of the temporal linear regression for the given time period was used to represent the average rate of vegetation growth or progressive loss. We calculated the long-term Sen’s slope⁶⁸ and only accepted significant slopes (p < 0.05). Additionally, we followed Xu et al.²⁴ in assuming that biomass gain cannot be properly measured once aboveground biomass density exceeds a certain level, as evidenced by the steep increase in the standard deviation of our estimates for high aboveground biomass density (Fig. S5). Thus, for undisturbed pixels with aboveground biomass density above 50 Mg ha^-1 we used the IPCC aboveground biomass density growth factors for mature forest⁶⁹.

Once all the significant changes across the time-series were identified, we generated a new aboveground biomass density dataset using the aboveground biomass density map from 2017 as reference (the year used to train the model, equation¹), and then rolled only the significant changes back to 2007 to generate a consistent temporal dataset.

Biome-level AGB quantification: We used terrestrial biomes as spatial units to calculate regional AGB statistics. Biomes are not strictly defined by the physical land cover types (e.g. forests, shrublands etc.), but rather a complex combination of climate and vegetation conditions⁷⁰. While biome extents may be dynamic, due to gains and losses in woody vegetation across landscapes, alongside changes in regional climatic conditions, to our knowledge, there is no data product that accounts for land cover dynamics to delineate biomes on an annual basis. We therefore used the most up-to-date static product available, the Ecoregions2017 Resolve biome map¹⁴ and assumed that biome extents have not significantly changed over the study period. Quantifying woody aboveground biomass losses / gains is itself an indicator of biomass changes in forests and shrublands.