Airborne Lidar Sampling Pivotal for Accurate Regional AGB Predictions from Multispectral Images in Forest-Savanna Landscapes

: Precise accounting of carbon stocks and ﬂuxes in tropical vegetation using remote sensing approaches remains a challenging exercise, as both signal saturation and ground sampling limitations contribute to inaccurate extrapolations. Airborne LiDAR Scanning (ALS) data can be used as an intermediate level to radically increase sampling and enhance model calibration. Here we tested the potential of using ALS data for upscaling vegetation aboveground biomass (AGB) from ﬁeld plots to a forest-savanna transitional landscape in the Guineo–Congolian region in Cameroon, using either a design-based approach or a model-based approach leveraging multispectral satellite imagery. Two sets of reference data were used: (1) AGB values collected from 62 0.16-ha plots distributed both in forests and savannas; and (2) an AGB map generated form ALS data. In the model-based approach, we trained Random Forest models using predictors from recent sensors of varying spectral and spatial resolutions (Spot 6 / 7, Landsat 8, and Sentinel 2), along with biophysical predictors derived after pre-processing into the Overland processing chain, following a forward variable selection procedure with a spatial 4-folds cross validation. The models calibrated with ﬁeld plots lead to a systematic overestimation in AGB density estimates and a root mean squared prediction error (RMSPE) of up to 65 Mg.ha − 1 (90%), whereas calibration with ALS lead to low bias and a drop of ~30% in RMSPE (down to 43 Mg.ha − 1 , 58%) with little e ﬀ ect of the satellite sensor used. Decomposing bias along the AGB density range, we show that multispectral images can (in some speciﬁc cases) be used for unbiased prediction at landscape scale on the basis of ALS-calibrated statistical models. However, our results also conﬁrm that, whatever the spectral indices used and attention paid to sensor quality and pre-processing, the signal is not su ﬃ cient to warrant accurate pixelwise predictions, because of large relative RMSPE, especially above (200–250 t / ha). The design-based approach, for which average AGB density values were attributed to mapped land cover classes, proved to be a simple and reliable alternative (for landscape to region level estimations), when trained with dense ALS samples.


Introduction
The vegetation in tropical Africa plays a major role in the global carbon cycle [1][2][3][4], providing valuable ecosystem services, storing vast amounts of carbon, and serving as a reservoir for climate Table 1. A selection of studies using airborne LiDAR scanning data to parameterize model-based approaches and generate wall-to-wall aboveground biomass (AGB) maps from spaceborne optical imagery. Statistics of model predictive performance (i.e., R 2 and root mean square prediction error; RMSPE) are derived from a variety of model validation strategies, including: (1) withholding of a given proportion of data for model testing, with test data selected at random (strat. 1); (2) 10-fold cross validation, with random split of data into folds (strat. 2); and (3) Monte Carlo cross-validation, with random or spatial split of data into folds (strat. 3 and 4, respectively). Details of cross-validation strategies can be found in the original studies. RMSPE provided in Mg of carbon per hectare in original studies are converted to Mg of AGB per hectare with a carbon-to-AGB conversion factor of 0.5. In contrast to this model-based approach to AGB extrapolation, a design-based approach remains well-accepted, notably recommended in the Intergovernmental Panel on Climate Change (IPCC) guidelines for green-house gas inventories [25]. In this approach, average AGB densities (referred to as emission factors in this context) are attributed to each land cover class. A precise estimation of these AGB densities requires rigorous sampling, for instance via a national forest inventory [26], which is still lacking in many tropical countries. This then imposes the need to resort to global database values for attributing AGB densities to each land cover class, resulting in lower confidences in estimates (Tier 1).

Country
In the context of a complex landscape mosaic in Cameroon, our objective was to: (i) test if intermediate ALS sampling allows improvements to large scale AGB estimations in model-based or design-based approaches; and (ii) test if the choice of multi-spectral satellite sensor or predictors can improve signal and predictive power in the case of model-based approaches.

Study Site
The field site is crossed by the Sanaga river and falls within a forest-savanna mosaic of the Guineo-Congolian region with gallery forests along the river courses [27] and covers an area of 216 km 2 (4 • 00'-4 • 30' N and 11 • 30'-12 • 00' E; Figure 1a).  [24] In contrast to this model-based approach to AGB extrapolation, a design-based approach remains well-accepted, notably recommended in the Intergovernmental Panel on Climate Change (IPCC) guidelines for green-house gas inventories [25]. In this approach, average AGB densities (referred to as emission factors in this context) are attributed to each land cover class. A precise estimation of these AGB densities requires rigorous sampling, for instance via a national forest inventory [26], which is still lacking in many tropical countries. This then imposes the need to resort to global database values for attributing AGB densities to each land cover class, resulting in lower confidences in estimates (Tier 1).
In the context of a complex landscape mosaic in Cameroon, our objective was to: (i) test if intermediate ALS sampling allows improvements to large scale AGB estimations in model-based or design-based approaches; and (ii) test if the choice of multi-spectral satellite sensor or predictors can improve signal and predictive power in the case of model-based approaches.

Study Site
The field site is crossed by the Sanaga river and falls within a forest-savanna mosaic of the Guineo-Congolian region with gallery forests along the river courses [27] and covers an area of 216 km 2 (4° 00'-4° 30' N and 11° 30'-12° 00' E; Figure 1a).  The area is under the influence of an equatorial climate of Guinean type [28], which is hot and humid with an average annual temperature of 25 • C. Mean annual rainfall is 1,500 mm, with a rainfall distribution characterized by a dry season lasting over three months (December-March), during which Remote Sens. 2020, 12, 1637 4 of 20 the monthly rainfall is less than 70 mm. The average altitude is 600 m. Dominating soil types in the area are ferralitic and hydromorphic soils [29]. The area is a target for small-and large-scale agriculture (artisanal cocoa farms, palm plantations) and industrial activities (e.g. hydroelectric dam construction), leading to a strong anthropogenic pressure on, for instance, the vegetation structure [30].
The vegetation is dominated by savanna formations (dominated by pyrophylous savanna species Terminalia glaucescens Planch. ex Benth; Annona senegalensis Pers.; Bridelia ferruginea Benth.) interspersed by gallery forests or semi-deciduous forests where Triplochiton scleroxylon K. Schum. and Terminalia superba Engl. & Diels dominate. Land cover classes comprising a woody component of significant AGB are agroforests (AF), degraded secondary forests (DF), old-growth secondary forests (OF), shrubby savanna (SS), and woody savanna (WS) (Figure 1c). Their proportions vary across the study area, with a higher proportion of AF, DF, and OF in the East (folds 3 and 4) and a higher proportion of SS and WS in the Western part (folds 1 and 2; Figure 1b,c).

Data Acquisition and Processing
The general methodology used in in upscaling aboveground biomass from plot scale to satellite scale with model-based and design-based approaches is presented in Figure 2.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 22 The area is under the influence of an equatorial climate of Guinean type [28], which is hot and humid with an average annual temperature of 25 °C. Mean annual rainfall is 1,500 mm, with a rainfall distribution characterized by a dry season lasting over three months (December-March), during which the monthly rainfall is less than 70 mm. The average altitude is 600 m. Dominating soil types in the area are ferralitic and hydromorphic soils [29]. The area is a target for small-and large-scale agriculture (artisanal cocoa farms, palm plantations) and industrial activities (e.g. hydroelectric dam construction), leading to a strong anthropogenic pressure on, for instance, the vegetation structure [30].
The vegetation is dominated by savanna formations (dominated by pyrophylous savanna species Terminalia glaucescens Planch. ex Benth; Annona senegalensis Pers.; Bridelia ferruginea Benth.) interspersed by gallery forests or semi-deciduous forests where Triplochiton scleroxylon K. Schum. and Terminalia superba Engl. & Diels dominate. Land cover classes comprising a woody component of significant AGB are agroforests (AF), degraded secondary forests (DF), old-growth secondary forests (OF), shrubby savanna (SS), and woody savanna (WS) (Figure 1c). Their proportions vary across the study area, with a higher proportion of AF, DF, and OF in the East (folds 3 and 4) and a higher proportion of SS and WS in the Western part (folds 1 and 2; Figure 1b,c).

Data acquisition and Processing
The general methodology used in in upscaling aboveground biomass from plot scale to satellite scale with model-based and design-based approaches is presented in Figure 2.

Field Inventory Data
Field data acquisition campaigns were conducted from February to December 2018. We established eighteen 0.16-ha plots in savannas and eleven 1-ha plots in forests. Field plots in savanna were set Remote Sens. 2020, 12, 1637 5 of 20 along a tree density gradient, from open grassy savanna to closed woody savanna [31] (Figure 1c). Savanna plots size (i.e., 40 × 40 m) was set so as to sample relatively homogenous vegetation structure within each plot, given the high spatial heterogeneity of the vegetation in this ecosystem. Forest plots had a standard 1-ha (100 m × 100 m) area, and were distributed among the main types of closed-canopy vegetations found in the area (see Figure 1b and Section 2.1). Plots geolocation was recorded with a Trimble Geo7X Global Navigation Satellite System (GNSS) receiver, by collecting a point every 20 m along each plot contour in order to increase geolocation accuracy [7]. In each plot, the diameter at breast height (DBH; in cm) was measured for all trees with DBH ≥ 10 cm. Trees were identified in the field by expert botanists, and herbarium specimens were collected on each species for cross-identification at the herbarium of Université Libre de Bruxelles (BRLU). In savanna plots, the height (H) of all trees with DBH ≥ 10 cm was measured using a graduated pole for short trees (H ≤ 7 m) and a laser range finder device (Trupulse 360"R) for taller trees (> 7 m). In forest plots, tree height was measured on a subsample of trees per plot (approximately 50 sampled along the DBH range) with their crown apex visible using the laser range finder device. We subsequently used the BIOMASS R package [32] to fit plot-level height-diameter allometric models via a three-parameter Weibull function and predict the height of unmeasured trees. In total, we censused 4309 trees in the forests, belonging to 150 species, and 3487 trees in the savanna from 43 species. In the following analyses, we used the savanna plot size (40 × 40 m) as our minimum mapping unit, and thus split each 1-ha forest plots in four subplots of 40 * 40 m (selecting each time the subplots located on the external edge of the 1-ha plots; Appendix A). The field dataset thus consisted in 62 field plots (18 and 44 plots in savanna and forests, respectively). Table 2 presents a summary statistics of the 62 plots data installed in forest and savanna sites.

Airborne LiDAR Data
Airborne LiDAR Scanning data were acquired in 2012 with a Riegl LMQ-560 sensor mounted on an airplane of type Pilatus PC6 with a flight height of 600 m above ground level and an average ground speed of 167 km.h −1 . The scan angle was 60 • with a band swath of 690 m and 50% of overlap among adjacent flight lines which resulted in an average point density of 8.4 points.m −2 . The ALS point error was 10 cm vertically and 15 cm horizontally. We used the lidR R package (v2.0.2; [33]) to process the ALS point cloud and generated a 2-m resolution canopy height model (CHM).

Spaceborne Data
We considered three multispectral sensors, namely Spot 6-7, Landsat 8, and Sentinel 2. Among cloud-free images in sensors archives, we looked for dry-season images acquired: (i) at (approximately) the same date (to mitigate cross-sensor differences associated to change in land cover or vegetation phenology); and (ii) as close as possible to the ALS acquisition date. The year 2015 matched our criteria and we collected Level-1C images for Spot 6/7 (acquired on January 9th, Row/Column 4912/3514), Landsat 8 (acquired on 1 December; Path/Row 185/057) and Sentinel 2B (acquired on December 19th mosaic of T32NQK & T32NRL).
We processed spaceborne optical data using the Overland algorithms [34]. Overland is a satellite image processing chain developed by AIRBUS DS Geo which aims to produce cloud and shadow masks and perform image atmospheric corrections, especially for areas with a high degree of cloudiness like western Central Africa. It is primarily coded in the IDL language (Harris) for image processing algorithms, with a core scene model and model inversion engine that has been developed in Matlab (MatWorks) Overland uses look-up tables from LOWTRAN and performs an inversion of a coupled atmospheric scene model [35] to estimate atmospheric parameters and discard influences from sky, aerosols, and clouds on the surface reflectance. Another feature of Overland is the ability to partition the reflectance of individual pixels into respective contributions of soil, photosynthetic vegetation (green matter), and the non-photosynthetic matter (dead wood), and characterize the self-cast shadows of the rough vegetation canopies. For this, it implements a vegetation model by combining PROSPECT [36], SAIL [37], and a soil model. By inverting this vegetation model, we can notably derive the fractional cover of green vegetation (fCover), the canopy shade factor (CSF), and the leaf area index (LAI).
To predict vegetation AGB, we considered the three variables provided by Overland, as well as images spectral bands and vegetation indices from corrected images summarized in Table 3. The dataset of spaceborne optical variables thus consisted of 11 variables for Spot 6/7, 13 variables for Landsat 8 and 24 variables for Sentinel 2.

Field-Based AGB Estimates
We used the R BIOMASS package [32], which relies on the Global Wood Density Database [38,39], to attribute to each tree a wood density (WD) value based on tree taxonomy. For trees identified at the species or genus level, the average WD of the respective taxonomic level was used. For trees identified at the family level, or unidentified trees, the plot average WD was used. We then computed tree AGB using allometric models based on DBH, H and WD. For forest trees, we used the pantropical AGB model of Chave et al. [40]. For trees from savanna plots, we used an AGB model developed for semi-arid savanna by Colgan et al. [41]. Lastly, we computed plot AGB density (AGB FIELD ) as the sum of individual tree AGBs over the plot area (expressed in Mg.ha −1 ).

LiDAR-Based Reference AGB Map and Land-Cover Classification
We extracted vegetation structural metrics from the CHM (2 -m resolution) within square windows of 40 × 40 pixels so as to be an integer divisor of the size of our field plots (40 × 40 m). We computed several vegetation metrics (quantiles of pixels height distribution, canopy gap fraction, statistics of leaf area vertical distribution). These metrics allowed us to derive a simple land cover classification map (Figure 1b) via supervised classification (maximum likelihood) in Envi 5.0.
To predict vegetation AGB density from the LiDAR canopy height model, we evaluated the predictive power of a set of models using a leave-one-out cross-validation (LOO-CV) procedure. LOO-CV consists of iteratively training the model on N-1 plots (with N the total number of plots), each time withholding a different plot for testing. To account for potential autocorrelation between 0.16-ha forest subplots extracted from the same 1-ha plot (which would violate the independence hypothesis between training and test sets of the CV procedure and result in overly optimistic CV statistics), we discarded 0.16-ha forest subplots of the same 1-ha plot from the training set each time a forest subplot was tested. The vector of independent AGB predictions was then used to compute CV statistics, namely the squared correlation between AGB predictions and AGB FIELD (henceforth R 2 ; Equation (1)) and the root mean squared prediction error (RMSPE), with: where n is the number of plots, y i is the AGB FIELD estimate for plot i, andŷ i is its AGB prediction. We found that a simple linear model (Equation (2)) based on vegetation median height (MCH) yielded the best results as in [16] with an R 2 of 0.81 and a RMSPE of 52.7 Mg.ha −1 (Appendix A): where AGB FIELD is the plot AGB density (Mg.ha −1 ) derived from field data, MCH is the vegetation median height (m) over the plot on the 2 m CHM obtained in Section 2.3.2. This best LiDAR-based model was used to predict vegetation AGB over the entire study area (AGB ALS map; res = 40 m) that constituted our intermediary scale for the cal/val of spaceborne-based and design-based AGB models (see next section).

Design-Based AGB Estimates
As a reference AGB prediction method, we followed the recommendation of the IPCC [25,49] by averaging AGB density values (either AGB ALS or AGB FIELD ) per (woody vegetation) land cover class. In the case of ALS, the area was fully characterized, and hence not sampled stricto sensu.

Model-Based AGB Estimates
Imagery products from the different satellite sensors were co-aligned to the AGB ALS raster and aggregated to 40 m resolution. We used the Random Forest (RF) algorithm [50] implemented in the randomForest R package [51] to model vegetation AGB density from spaceborne optical variables. RF is a popular machine learning technique in remote sensing studies due to its ability to handle high-dimensional datasets, to account for non-linear relationships between response and predictor variables, and to its relative robustness to multicollinearity, model overparameterization, and overfitting [52]. Here, we built several RF models, considering each time a different set of spaceborne optical variables (i.e., from Spot 6/7, Landsat 8, or Sentinel 2) and using either field-based Remote Sens. 2020, 12, 1637 8 of 20 AGB estimations (strategy 1, henceforth RF FIELD ) or the LiDAR-based AGB ALS map (strategy 2, henceforth RF ALS ) for model specification and training. We then evaluated the predictive power of the RF models by comparing RF models' predictions to independent reference AGB estimations.
Since RF have shown to overfit training data when predictors are correlated (e.g. [32,33]) or spatially auto-correlated [53], both properties being expected in our sets of optical variables, we performed a spatial forward variable selection procedure (see Appendix B). This procedure starts with no variable in the model, computes the decrease in model's relative RMSPE (calculated by dividing the RMSPE by the mean of "observed" AGB values) that the addition of each candidate explanatory variable would lead to, and adds to the model the variable leading to the largest relative RMSPE decrease. The procedure is iterated as long as adding a supplementary variable in the model leads to relative RMSPE decrease larger than 1%. In the case of RF FIELD , the RMSPE at each iteration of the procedure was computed using a LOO-CV over the 62 AGB FIELD estimates (as in Section 2.3.2). In the case of RF ALS , the number of AGB ALS estimates (i.e., 117,415 pixels) makes the LOO-CV computationally prohibitive. We thus used a four-fold block CV for variable selection.

Four-Fold Cross-Validation
We evaluated the ability of design-and model-based approaches to predict vegetation AGB outside training areas, using the LiDAR-based AGB ALS map as target AGB density values, and the four folds (blocks) defined in Figure 1c.
For design-based approaches, average AGB densities per class were assessed in the three training folds, either using ALS or field data from these folds. In the remaining fold, the mean fold AGB density was computed by multiplying the AGB density of each land cover class by its respective area, relative to the total woody vegetated area of the fold. Validation could only be performed (qualitatively) at the scale of each land cover class or the whole fold.
For model-based approaches, in the case of RF ALS , we generated AGB density predictions outside training areas using the four-fold block CV. In the case of RF FIELD , to circumvent the relatively small plot number, all 62 AGB FIELD estimates were used to train the model (see LOO-CV procedure in Section 2.3.2) and generate AGB density predictions in each of the four folds. The vector of independent AGB predictions was used to compute the CV statistics as in Section 2.3.2 (i.e., R 2 and RMSPE). We also computed the mean signed deviation (MSD; Equation (3)) as an indicator of model bias: where n is the number of pixels,ŷ i is the AGB prediction of pixel i by RF FIELD or RF ALS and y i is its AGB ALS value.

Results
In model-based approaches, the variables selected varied depending on the training data (calibration on field AGB; RF FIELD or calibration on the intermediate AGB ALS map; RF ALS ; see Appendix B) and satellite sensor used for model calibration. Table 4 shows the variables selected (RF FIELD & RF ALS ) and the model performances in cross-validation for the different models. Independently of the satellite sensor, RF FIELD models gave the poorest performances, with R 2 values around 0.6, and RMSPE of up to 65 Mg.ha −1 (i.e., 90%), whereas RF ALS models greatly improved the prediction accuracy (drop of~30% in RMSPE and relative RMSPE). RF ALS models based on Landsat 8 or Sentinel 2 predictors lead to a decrease of 10% both in RMSPE and relative RMSPE (R 2 = 0.7; RMSPE of c. 43 Mg.ha −1 and relative RMSPE of c. 60%) compared to RF ALS models based on Spot 6/7 predictors (R 2 = 0.6; RMSPE = 48 Mg.ha −1 and relative RMSPE = 66.5%) As model-based approaches provide pixelwise predictions, we can have a detailed look at the scatterplots of observed vs. predicted AGB density values in Figure 3 (for the Sentinel 2 sensor), and Appendix C & Appendix D (Spot 6-7 and Landsat 8 respectively). Concordance between predictions and observations was greatly improved, i.e., closer to the 1:1 line, for models calibrated on ALS data relative to those calibrated with field data, whereas the spaceborne sensor used seemed to make little difference. Predictions were capped around 250 Mg.ha −1 for all RF ALS models. At first sight, RF FIELD models seemed capable of predictions over slightly higher AGB ranges than RF ALS (point cloud extending beyond predicted AGB values > 250 Mg.ha −1 in Figure 3a), but predictions were in fact much less accurate overall and strongly biased in this upper range.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 22 RMSPE of c. 60%) compared to RFALS models based on Spot 6/7 predictors (R 2 = 0.6; RMSPE = 48 Mg.ha −1 and relative RMSPE = 66.5%) As model-based approaches provide pixelwise predictions, we can have a detailed look at the scatterplots of observed vs. predicted AGB density values in Figure 3 (for the Sentinel 2 sensor), and Appendix C & D (Spot 6-7 and Landsat 8 respectively). Concordance between predictions and observations was greatly improved, i.e., closer to the 1:1 line, for models calibrated on ALS data relative to those calibrated with field data, whereas the spaceborne sensor used seemed to make little difference. Predictions were capped around 250 Mg.ha −1 for all RFALS models. At first sight, RFFIELD models seemed capable of predictions over slightly higher AGB ranges than RFALS (point cloud extending beyond predicted AGB values > 250 Mg.ha −1 in Figure 3a), but predictions were in fact much less accurate overall and strongly biased in this upper range.   To better characterize the uncertainties of RF models estimates and possible dependence to AGB, we grouped the mean values of the RMSPE and the MSD (all in Mg.ha −1 ) into equally sized bins of 25 Mg.ha −1 . We decomposed error and bias along the axes of both predicted and observed AGB. From a prediction perspective, RF FIELD models lead to a systematic overestimation of AGB ALS (MSD > 0, hence predictions are higher than observations) across the whole predicted AGB range, with both increasing bias and error towards large AGB values (Figure 3a-c; the diagrams display model performances for Sentinel 2 predictors). This means that in any predicted AGB bin, average predictions were both largely inaccurate and imprecise when calibrating extrapolation models with field data, and all the more so at high AGB values. When binning errors along the observed AGB axis, the story was a little different, with a hinge point around 200 Mg.ha −1 below which predictions overestimated observed AGBs (MSD > 0) and above which predictions increasingly underestimated observed AGBs. RF ALS models allowed a near complete bias reduction in predicted AGB values with a MSD close to zero (Figure 3d,e) along the whole range of predicted AGB values except for the highest values. In other words, despite considerable scatter, predictions were on average unbiased at the landscape level across all AGB bins up to about 200 Mg.ha −1 . Along the observed AGB axis, error and bias were distributed similarly to RF FIELD , with a slightly lower hinge point.
To have a better understanding of the implication of the different AGB prediction approaches outside of training areas, we had a closer look at the distribution of predicted AGB density in the four blocks used in the 4-fold CV (Figure 4). We first investigated ALS-trained approaches. Using RF ALS , (a model-based approach, based on Sentinel 2 imagery data), we could use the continuous distribution of biomass predictions. In all folds, but more markedly in the 3rd and 4th folds, the shape of the density curve diverged from the reference (ALS) curve over the higher range of AGB density values, with the expected drop around 250 Mg.ha −1 . Folds 3 and 4 comprise a much higher proportion of forest vegetation than the two others (Figure 1c), which explains the difficulty of the model-based approach to faithfully reproduce de AGB density distribution. At the fold scale however (colored arrows beneath density plots), the mean predicted AGB density was close (below +/-10%) to the reference ALS value, except for the 4th fold (+13% difference). This illustrates that a well-trained statistical model, even with poor per-pixel predictive power, can provide unbiased prediction at the landscape level, at least as long as the landscape matrix is not too different from the training conditions. The design-based approach, which does not require any biomass prediction model from optical satellite predictors and 'blindly' applies an average LiDAR-based AGB density value to each land cover class, appeared to perform equally well, with mean AGB densities between -7.2% and +12% of the ALS reference value.
When focusing on specific land cover classes (histogram insets in Figure 4) the model-based approach, as could be expected, systematically underestimated AGB density relative to LiDAR-based estimates in the old-growth secondary forest land cover type. The expected opposite trend (overestimation) at the other end of the biomass gradient, was more subtle in low biomass vegetations such as woody savanna, with higher overestimations on the 4th fold. Design-based predictors did not present a systematic bias tendency in any vegetation type.
A similar analysis with models calibrated with AGB FIELD is presented in Appendix E. Here model-based predictions showed clearly aberrant density curves, and produced fold-level mean AGB density predictions comprised between +36.6% and +35.4% of the ALS reference. The design-based approach performed better in three folds, with meant overall AGB density below 12%, but showed a 30.9% bias in fold 3. These variations can be explained by the poor sampling rate (and design) in the training folds due to the small number of available plots, which results in some large errors in the estimation of the AGB density of some of the land cover classes (up to 100%).  To have a better understanding of the implication of the different AGB prediction approaches outside of training areas, we had a closer look at the distribution of predicted AGB density in the four blocks used in the 4-fold CV (Figure 4). We first investigated ALS-trained approaches. Using RFALS, (a model-based approach, based on Sentinel 2 imagery data), we could use the continuous distribution of biomass predictions. In all folds, but more markedly in the 3 rd and 4 th folds, the shape of the density curve diverged from the reference (ALS) curve over the higher range of AGB density values, with the expected drop around 250 Mg.ha −1 . Folds 3 and 4 comprise a much higher proportion of forest vegetation than the two others (Figure 1c), which explains the difficulty of the model-based approach to faithfully reproduce de AGB density distribution. At the fold scale however (colored arrows beneath density plots), the mean predicted AGB density was close (below +/-10%) to the reference ALS value, except for the 4 th fold (+13% difference). This illustrates that a well-trained statistical model, even with poor per-pixel predictive power, can provide unbiased prediction at the landscape level, at least as long as the landscape matrix is not too different from the training conditions. The design-based approach,

Discussion
Transitional landscapes cover broad extents in Central Africa, and thus significantly affect the carbon budget of the continent [54]. Achieving unbiased estimation of carbon stocks and fluxes in these highly dynamic environments, characterized by mosaics of very different land uses and covers, is a critical challenge. Major climate change mitigations strategies, such as national greenhouse gas inventories (MRV-REDD+) or high carbon stock (HCS) approaches, depend on our ability to meet this challenge.
Despite known limitations, multispectral spaceborne data remain widely used for AGB extrapolations. A widely held assumption is that Airborne LiDAR scanning (ALS) data ensure better model calibration, and hence partly compensates signal limitations [14,16,18,55]. Our results indeed show that model fit and error can be drastically improved, with an R 2 of 0.7 and a RMSE decrease of 30 %, when using a Random Forest model calibrated with AGB ALS reference data (RF ALS ) instead of field plots (RF FIELD ). Predictions of the latter moreover proved highly biased (inaccurate) and imprecise along the whole range of AGB densities. Looking closely at the AGB predictions of RF ALS , we can however see that the signal of multispectral data, whatever the sensor, does not in fact allow accurate AGB predictions in low and high ranges of actual vegetation biomass (i.e., AGB ALS ). Indeed, low biomass ranges are systematically overestimated, and high biomass ranges are underestimated. This behavior has been already evidenced across the tropics [23,56] in similar studies, highlighting the fact that errors in AGB estimations when using optical signals are unfortunately still both unavoidable and crippling for large scale wall-to-wall AGB mapping [7]. It is thus important to better benchmark the values and limitations of optical signals in varying contexts through landscape-scale studies integrating highly informative ALS data. Here we showed that improved calibration of spaceborne models by ALS data did ensure unbiased estimation of AGB overall. In other words, although the signal does not allow accurate predictions of AGB in a given forest or savanna location, RF ALS models still provide accurate prediction of average AGB levels across the landscape. This is only true, however, inasmuch as the balance between land cover types in the predicted landscape is comparable to that of the training area. Notably, the underestimation of total AGB will plummet with the share of land harboring high AGB forests, above the hinge point of saturation of about 225-250 Mg.ha −1 . This threshold ought to be kept in mind in any further applications.
Regarding possible effects of the compromise between spatial and spectral resolution allowed by different spaceborne sensors, RF ALS models based on lower spatial resolution and narrow-wavelength spaceborne images (i.e., Landsat 8 & Sentinel 2) seemed to perform slightly better than models based on higher spatial resolution broad-wavelength imagery (Spot 6/7) with a minor decrease of 10% in RMSPE. Similarly, the variability in the spectral predictors selected by our forward model selection procedure does not allow being conclusive regarding the relative interest of any given spectral index over others.
The well-accepted design-based approach, on the other end, provided a simple and accurate alternative for landscape-level AGB estimation, when trained on a dense sample of ALS data. As a single AGB value is attributed to each land cover class, this approach does not provide a detailed intra-class variation map. However, if the model-based approach does provide such a map, it is unreliable anyways, as we have shown, as long as the available predictors remain poorly correlated to AGB. It might be better not to lure the user with the pretense of a high-resolution product, when the estimates are only valid at large scale.

Conclusions
We showed that airborne LiDAR-based AGB data can significantly improve the calibration of prediction models from spaceborne multispectral data, with an error reduction of~30% compared to a field-based AGB calibration. It is however crucial to acknowledge that, due to signal limitations, irrespective of the multispectral sensor and mix of spectral indices used, predictions are only unbiased at the landscape or regional level, and for land cover conditions similar to the training area. Data from upcoming new radar-based (BIOMASS, NISAR) and LiDAR-based space missions (GEDI, ICESat-2, MOLI) are expected to improve our extrapolation capacities by providing both global coverage and signals with better relationships to vegetation structure [7,57]. However, some of those valuable data will need spatial interpolation (GEDI or Ice) from optical data or will not be available for long-term monitoring. Therefore improving the use of optical data will remain an issue. Attributing average AGB density values to broad land-cover classes (referred to as design-based sampling here) will continue to remain a valid alternative to obtain regional unbiased AGB estimations.
Author Contributions: L.B.T.S., P.P. and N.B. conceived and designed the experiments. H.P. helped with image processing using the Overland processing chain. L.B.T.S. performed the experiments and analyzed the results, with the help of P.P. and N.B. L.B.T.S. wrote the first drafts of the paper and N.B., P.C., P.P., H.P. and B.S. contributed edits and advice on content. All authors have read and agreed to the published version of the manuscript.
Funding: This work is financially supported by the Nachtigal Hydropower Company (Contract n • C006C007-DES-2017) under the environmental impact study associated to the construction of a hydroelectric dam on the Sanaga River in the Central region of Cameroon. We also benefitted from the "Allocations de recherche pour une thèse au Sud (ARTS)" grant program to fund research stays in France.