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)48, 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 biomass49. The Global Ecosystem Dynamics Investigation (GEDI) LiDAR instrument22 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-248 and optical multispectral Landsat-derived percent tree cover maps29 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 Center50, 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.31 and Araza et al.15 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 plots55. 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 cells15. 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.15).
ALOS PALSAR/ALOS-2 PALSAR-2 radar image mosaics. JAXA’s annual mosaics of L-band SAR HH and HV polarised backscatter (γ0) were based on ALOS PALSAR from 2007 to 2010 and ALOS-2 PALSAR-2 from 2015 to 201756,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 applied22,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 in58, 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))59.
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 201729 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.9th percentile of the initial set. The Copernicus Global Land Cover dataset60 was used to exclude footprints in non-vegetated classes. We used the 11 forest classes and the shrubland class for this purpose60. 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 100 ) 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 footprints6. 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.
... continue reading