Above-Ground Biomass Estimation Based on Multi-Angular L-Band Measurements of Brightness Temperatures

There is growing interest in using passive microwave observations and vegetation optical depth (VOD) to study the above-ground biomass (AGB) and carbon stocks evolution. <italic>L</italic>-band observations, in particular, have been shown to be very sensitive to AGB. Here, thanks to the multiangle capabilities of the soil moisture and ocean salinity mission, a new approach to estimate AGB directly from multiangular <italic>L</italic>-band brightness temperatures (TBs) is proposed, thus surpassing the use of intermediate variables such as VOD. The European Space Agency (ESA) Climate Change Initiative (CCI) Biomass maps for the years 2010, 2017, and 2018 are used as the AGB reference. AGB estimates from artificial neural networks (ANN) using a purely data-driven approach explained up to 88% of AGB variability globally; even so, a decrease in retrieval performance was observed when models are applied to data from years different than the year used for their training. A new training methodology based on multiyear training sets is presented, leading to results showing more stability for temporal analyses. The best set of predictors and an optimal learning dataset configuration are proposed based on an assessment of the accuracy of the estimates. The ANN methodology using TBs is a promising alternative with respect to the common method of using a parametric function to estimate AGB from VOD. ANNs AGB estimates showed a higher correlation with CCI AGB maps (<inline-formula><tex-math notation="LaTeX">$R$</tex-math></inline-formula><sup>2</sup> <inline-formula><tex-math notation="LaTeX">$\sim$</tex-math></inline-formula>0.87 instead of <inline-formula><tex-math notation="LaTeX">$\sim$</tex-math></inline-formula>0.84) and presented a stronger agreement with their spatial structure and less differences in residual maps.


I. INTRODUCTION
D UE to its importance on the atmosphere-land carbon fluxes and land carbon stocks, above-ground biomass (AGB) has been listed as one of the essential climate variables by the Global Climate Observing System; vegetation captures atmospheric CO 2 through photosynthesis, a fraction of which is fixed and used by plants for biomass production [1], the rate of this production (which, in turn, depends on plant species, canopy structure, soil management, and climate, among others [2], [3]) minus vegetation mortality and fires, dictates the magnitude of carbon stocks magnitude in vegetation. Moreover, evidence shows that the role in carbon sequestration of globally important ecosystems such as the pantropical forests is sensitive to climate change; tropical forest could, under certain conditions, change from being a sink or a neutral agent into a net source of carbon [4], [5]. This highlights the importance of AGB estimation at a global scale, and the necessity to better understand the state and extension of vegetated areas over time.
In response to this need, various methods have been used to study AGB. The closest to a true value of vegetation biomass is obtained through in situ measurements and allometric equations [6]; however, they are local and have little temporal resolution. Remote sensing, on the other hand, has proven to be an effective tool for studying large-scale vegetation dynamics over long periods of time. The use of optical indices has been until recently, the most commonly employed remote sensing technique to follow the evolution of vegetation. Several studies have used indices such as the NDVI [7] to monitor the temporal trends of vegetation at different scales [8], [9]. Nonetheless, its use to estimate AGB is limited, as it saturates quickly over densely vegetated areas (saturation at around 50 MgC/ha [10]). More sophisticated approaches using allometric relationships, LIDAR measurements, and optical reflectances have been discussed in the literature that, at least partly, overcome that limitation [4].
Microwave remote sensing at lower frequencies has been shown to overcome these limitations. Due to their high spatial resolution, active microwave sensors such as synthetic aperture radars (SARs) have been used to estimate AGB. Early studies of airborne SAR at multiple frequencies showed the sensitivities of radar backscatter coefficients (at C-, L-, and P-band) to vegetation biomass; lower frequencies (L-and P-band) were strongly correlated to the total AGB while weaker correlations were found between total AGB and higher frequencies (Cband). Several studies have further explored this dependence and have shown higher sensitivities from low frequencies (L-and P-band, especially P) to the woody components of vegetation over densely vegetated areas [11], [12], [13]. Despite the high spatial resolution of synthetic radar observations, backscattering tends to show saturation effects with increasing biomass values, and it has even been reported that the backscattering versus biomass relationship can be inverted in case of very high biomass values [14].
Other active remote sensing instruments, such as Geoscience Laser Altimeter System on board the ICEsat [15] or the LiDAR GEDI [16] on board the ISS, have been used to estimate AGB from canopy height retrievals at high resolution, and recently Global Navigation Satellite Systems Reflectometry (GNSS-R) has also been proposed for the study of biomass [17]. From some of the previously mentioned techniques (lidar, optical reflectances, radar backscattering, and/or plot inventories), studies such as Baccini et al. [4] and Saatchi et al. [18] have produced, under different methodologies, AGB maps (representative of the period 2000-2010 and circa 2000, respectively) over the pantropical region. Recently, in the framework of the European Space Agency (ESA) Climate Change Initiative (CCI) for Biomass, three AGB maps have been produced for the years 2010, 2017, and 2018, all under a similar methodology using Cand L-band radar backscattering [19].
In the past few years, passive microwaves (PMWs) have been proposed and used as a complementary way to estimate AGB. Their uninterrupted long time-series (almost 13 years of observations for missions like the soil moisture and ocean salinity (SMOS) [20]) provide an invaluable source of information for the global study of AGB. The absorption and scattering effects of the canopy on the soil's microwave emission can be parameterized via the optical depth τ , and the single scattering albedo ω, using the so-called τ − ω model. The optical depth is related to attenuation, which is produced by scattering and absorption. Scattering also affects the albedo. Both are due to the interaction of the electromagnetic radiation with the vegetation, in particular with the water molecules contained in the vegetation. Therefore, the optical depth is frequently referred to as vegetation optical depth (VOD). VOD shows an almost linear relationship with vegetation water content (VWC) over croplands [21]. This variable has also been shown to be very sensitive to AGB [22], being lower frequencies VODs (like the SMOS L-band product) more sensitive to AGB than those from higher frequencies [23] or optical indices [24].
Liu et al. [22] explored the evolution of carbon stocks from almost two decades of PMW observations from several instruments working at C-, X-, and K-bands. However, Rodriguez-Fernandez et al. [23] showed that L-band VOD is significantly more sensitive to AGB than higher frequencies VOD. Although at a lower spatial resolution than that of optical and SAR observations, SMOS L-band VOD allows us to compute yearly maps of AGB from 2011 to 2022 and to study the temporal evolution of global AGB and carbon stocks in that period. For instance, using L-VOD, Brandt et al. [10] highlighted the importance of drylands to evaluate the global carbon stocks in Africa. L-VOD has also been used to evaluate the carbon stocks evolution in the pantropical region [25] exploring the effects that extreme events such as El Niño 2015 had on the carbon balance [26].
However, despite its demonstrated potential as a proxy for biomass, the retrieval of VOD is done through radiative transfer models, which tend to simplify physical processes and are commonly parameterized based on the characteristics of the measuring instrument. Estimating the radiation transfer through vegetation is a complex process, and even when starting from the same observations, different algorithms can compute different VOD datasets. In addition, particularly when starting from SMOSs multiangle observations, there can be some information loss when going from multiangular measurements to a single variable such as VOD. Given the wide variety of VOD products (see [27]) and the inherent complexity (or simplicity) of radiative transfer models, it is pertinent to study alternative methods to estimate AGB from L-band multiangular observations. Rodriguez-Fernandez et al. [28] showed that it is possible to find statistical relationships between L-band TB measurements and AGB. Previous studies have explored the sensitivity of variables such as SM, EVI, NDVI, or TB at different frequencies, to AGB [28], [29], [30].
However, this is the first study to propose a methodology for the estimation of vegetation biomass directly from multiangular L-band TBs. In the current study, an approach to estimate AGB directly from SMOS multiangular TBs is presented using artificial neural networks (ANN) trained on multiyear AGB maps from the ESAs CCI Biomass.
The rest of the article is organized as follows. Section II presents a description of the datasets used in this study. Section III describes the methodology applied to evaluate this new AGB retrieval method. Section IV presents the results obtained in this study, which are subsequently discussed in Section V. Finally, Section VI concludes this article.

A. L-Band Level 3 TB and VOD/SM Datasets
The SMOS satellite was launched on November 9, 2009, with the objective of measuring two key geophysical parameters: soil moisture (SM) and sea surface salinity (SSS). For this purpose, SMOS measures Earth's microwave radiation at a frequency of 1.4 GHz (21 cm) for several incidence angles (between 0 • and 60 • ) and on two orthogonal polarizations. The satellite overpasses the Equator at 6:00 (ascending) and at 18:00 (descending) following a sun-synchronous polar orbit, with a revisit time of 1-3 days depending on the latitude [20].
Level 3 TBs (L3TB) are provided in fourteen angle bins of 5 • width (excluding bin number nine). Hereafter, bin No.1 denotes incidence angles from 0 • to 5 • , bin No.2 from 5 • to 10 • and so on until bin No.14; exceptionally, the complementary bin number 9 is centered at 40 • with a 2 • width. L3TB TBs are provided in horizontal and vertical polarizations, hereafter H and V polarizations, respectively.
The Level 3 L-VOD/SM product (L3SM version 339), hereafter VOD (and SM), and TBs (version 330) are produced through the implementation of the L-band microwave emission of the biosphere (L-MEB) radiative transfer model [31] on the SMOS Level 3 retrieval algorithm [32]. Starting with version 330 of the L3 algorithm, low vegetation and forest areas are added to a common fraction of vegetated surface during the retrieval process of both VOD and SM. Therefore, it adopts the same strategy as that of the ESA Level 2 SM algorithm version 700.
Both the L3SM and L3TB datasets are provided by the CATDS (Centre Aval de Traitement de Données SMOS) from 2010 until 2022 on a daily time step and are projected on the equal-area scalable Earth (EASE) grid version 2 [33] with a grid spacing of ∼25 km. For this study, only the observations corresponding to the years of the AGB reference maps (2010, 2017, and 2018) were used (described in Section II-B). For the year 2010, however, even when the reference map is representative of the year 2010, the TBs used correspond to the year 2011, since it is the closest year to the reference with complete observations (the SMOS commission period ended in May 2010). TBs from the ascending orbits were used. To remove the dependence on seasonal variations of VWC on the signal, annual averages of each variable were calculated for the aforementioned years.
Daily filtering was performed on the TB and VOD/SM products before calculating the annual averages. The probability of the appearance of radio frequency interference (RFI) was used as a quality criterion. In the case of TBs, any observation whose probability of appearance of RFIs (number of RFIs detected/number of daily observations) by angle of incidence was greater than 0.7 was discarded. For VOD products, in addition to the RFI probability filter (> 0.2), a χ 2 coefficient (which measures the goodness of fit from TB estimates and TB observations during the minimization process of the VOD/SM) lower than 2 was selected.
Additionally, taking into account that the SMOS retrievals over coastal pixels tend to be strongly affected by the fraction of water in the observed scene (especially at high incidence angles), all observations one pixel away from the coasts (defined by the USGS land-sea mask used by the SMOS algorithm) were removed.

B. AGB Reference Maps
Three AGB maps from ESAs CCI biomass project with global coverage were used as vegetation biomass references for the years 2010, 2017, and 2018 (datasets version 3.0) [19], from now on referred to as CCI2010, CCI2017, and CCI2018, respectively. Reference maps are provided with a spatial resolution of 100 m, and represent biomass density in Mg/ha. The retrieval algorithm used for the production of these datasets is based on a weighted linear combination of biomass estimates from C-and L-band SAR observations. In densely vegetated areas (such as tropical forests) more weight is given to L-band observations, in other regions to C-band observations. The CCI2010 is produced from ALOS-1 (PALSAR-1) and Envisat ASAR datasets, the CCI2017-18 are produced from ALOS-2 (PALSAR-2) and Sentinel-1 observations of SAR backscatter (additional datasets were used during the calibration of models) [34].
According to the validation protocol applied by the CCI biomass project (comparing the CCI maps with true AGB observations from forest field inventories), the CCI2017 and CCI2018 were in better agreement with AGB references than the one from 2010. Although the CCI2010 showed a lower agreement with field estimations, it still offers an approximation of the biomass for that reference year, and therefore, it was used in this study. In spite of their possible uncertainties, the CCI AGB maps will be referred to as "references," assuming that they are a good approximation to the actual state of vegetation biomass globally. For this study, all AGB reference maps have been regridded to the common grid of SMOS L3 products (EASE grid v2) by averaging the original 100-m pixels within the ∼25-km EASE grid cells.

A. AGB Estimation From Multiangular L-Band TBs
The proposed methodology for estimating AGB maps from ANN is based on a purely statistical relationship between inputs, in this case, multiangular TBs, and targets, the CCI reference maps. Yet, each polarization of the TBs is grouped into 1-14 bins (according to their angle of incidence-see Section II-A). As a result, the number of possible predictors can be as large as 28. In order to determine the optimal set of predictors to retrieve biomass from ANNs, an evaluation of the estimated AGBs was performed. Taking advantage of the multiangle capabilities of the SMOS mission, the quality of the produced biomass maps was evaluated in two ways. First, the influence of the TBs angle of incidence on the quality of the AGB inversions was assessed. Then, the influence of different predictors was analyzed, being TBs H and V the main ones. Other variables such as the polarization ratio [PR, (1)] and soil temperatures (TSs) were used as complementary predictors. TS was also used to estimate emissivities (E) from the TBs, calculated as E = TB/TS; henceforth, E in horizontal and vertical polarizations will be referred to as E h and E v , respectively. The polarization ratio is defined as Scattering events within the vegetation strongly depolarizes the ground emission. Therefore, low PR values are associated with dense vegetation and higher values with sparsely vegetated soils [35], [36], [37].
In addition, a hybrid approach for the estimation of AGB was evaluated. In this approach, the outputs of a physically based radiative transfer model (VOD or SM) were used along with the TBs as complementary predictors in the ANN input vector. To summarize, AGB estimates will be presented from two ANN approaches: The first referred to as 1) the data-driven approach (ANN D ), which is based solely on TBs (and their derived variables), and the second referred to as 2) the hybrid approach, which uses both TBs and VOD (or SM) as inputs for the networks (ANN H ), see Fig 1. An additional case using parametric functions based solely on the VOD will be used and discussed later in the document (VOD function).
The learning sets used by the networks are divided in three subsets. The first one, the training subset, is used to train the network and adjust its parameters (weights and biases) in order to reduce the difference between estimates and targets. The validation subset is used to monitor the performances of the ANN at each step of the training process and to detect overtraining if the performances continue to increase on the training subset but start decreasing on the validation subset. The last subset, the test one, is used to evaluate the ability of the ANN to generalize the patterns learned during the training stage to a set of new data. The training, validation, and test subsets were built from a random selection of the global observations contained in the learning set, and represent respectively 50 (∼58 000 pixels/year), 20 (∼23 000 pixels/year), and 30% (∼35 000 pixels/year) of the original set.
The optimization of the ANN parameters is done using the Levenberg-Marquardt algorithm; the performance of the network is calculated using the mean square error (MSE). A simple architecture is used: a single hidden layer with 10 neurons with a logistic-sigmoidal activation function, together with an output layer with a linear activation function. Initially, the training is carried out on input-target pairs of the same year. However, thanks to the three available years, temporal cross-validation is applied. In other words, once a model is trained on data from a specific year, it is then applied to data from another year in order to compare the resulting inversion to the reference map of that period, e.g., an ANN trained with the input-target pair: TBs 2017-CCI2017 is subsequently applied to TBs of 2018 to compare their inversion to the CCI2018 map, and so on for all other possible combinations. An alternative training scheme, based on a multiyear learning-set (a concatenation of inputs and references from the three available years) was explored to assess the flexibility of such an ANN for temporal cross-validation; the same ANN architecture mentioned before is used, the construction of the learning subsets (training, validation, and test) is carried out from a random selection of the three-year concatenated dataset.

B. AGB Estimation From Parametric Functions
As mentioned in Section I, currently, the common approach to estimate AGB from PMWs is based on determining a relationship between VOD and AGB using a functional fit. Previous studies have explored this methodology and have proposed different parametric functions to estimate AGB from VOD [22], [23]. Here, the one proposed by Rodríguez-Fernández et al. [23] was used. To do so, VOD datasets are sliced into bins of 0.05 units of width. For the corresponding AGB values of each VOD bin, the mean and the 5th and 95th quantiles are found. The resulting AGB-VOD curves are then fitted using (2), in which a, b, c, and d are the set of adjustable parameters used to fit the function to the AGB dataset (2)

C. Quality Evaluation of AGB Estimates
The quality of the AGB inversions produced by ANNs was evaluated under various criteria. First, with the coefficient of determination (R 2 ), to assess how well the variance of a target, in this case, the AGB references can be explained by a set of predictors (measured between 0 and 1, higher values represent better approximations). Second, by calculating the dispersion of the predictions with respect to the reference map using the root mean squared error (RMSE); expressed as percent error by computing: RMSE/mean(AGB reference) × 100. Finally, a qualitative evaluation was done under two conditions: 1) the spatial distribution of residuals (difference between the inversion and the reference map) and 2) an analysis of the spatial patterns of the AGB maps produced.

A. Quality Evaluation and Predictors Selection
Variations in the R 2 scores between AGB references and estimates, with respect to the angle of incidence and different predictors, are shown in Fig. 2. Several combinations of ANN predictors were tested, using single or grouped variables. Higher R 2 values were found in retrievals from higher incidence angles (bin > 8, incidence angles > 35 • ). In single predictors, the H-polarization and PR presented the highest R 2 values (around 0.72 and 0.80, respectively). V and emissivities alone presented nonsignificant values, around 0.1 for the most part and 0.45 as the highest value (for E h between 60 to 65 • [Bin No. 14]). The highest R 2 was obtained when multiple predictors were used at the same time. The most performant combinations were Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.  the cases j, k, and l in Fig. 2, all with similar values close to 0.85; the use of E v and V did not significantly increase the coefficient' scores. The case f, which did not use TBs but only variables derived from TB and TS, and case i, which used only TBs and PR, showed similar scores for angles on incidence > 35 • (R 2 ∼0.82). In inversions with multiple predictors, the less performant approach was that of case g, with slightly lower R 2 scores, at around 0.81.
Thus, from two of the most performant groups of predictors: case j, with the highest values of R 2 , and case i with relatively lower results but using exclusively TBs, the quality of AGB estimates (in terms of R 2 ) was evaluated for different ranges of angles of incidence on a global and continental scale (see  . Three additional sets of predictors were studied, two for the hybrid approach (using VOD and SM as complementary predictors to H and PR) and one using TS due to its strong influence in TB observations (inputs f, g, and h in Fig. 3).
Estimates from a single angle of incidence had the lowest R 2 values among all combinations. However, using one angle with the predictors (H, PR, E h ) presented similar scores to those produced from angles 35 •  A similar analysis to that of Figs. 2 and 3 was done for 2010 and 2018 (at a global scale and for Africa and South America), and similar results were obtained. Consequently, henceforth the following groups of predictors will be used: 1) the optimal combination of predictors (H, PR, E h ), as it gives the best compromise between angles of incidence and higher R 2 scores of the estimates, 2) the combination (H, V, PR), which despite having relatively low values for the coefficient of determination, exhibits the results that can be expected from multiangular TBs alone, and 3) both hybrid cases, as they presented considerable R 2 values without the use of emissivities; despite having similar results to the latter, the combination that includes TS as complementary predictor will not be used, since it is considered that the effect of the surface temperature on the AGB estimates was better captured from the emissivities (derived from TB and TS).
Model inversions from the selected cases were evaluated in both their ability to explain the variability of the reference AGB maps (R 2 ) and their dispersion with respect to them (RMSE). The first three matrices (from top to bottom) of Fig. 4, present the R 2 values obtained for the AGB inversions produced for the three reference years. In Fig. 4, the title of each matrix represents the year of the dataset that was used for the training of Models; therefore, the second matrix titled "2017" contains the R 2 of all the models trained using input variables from 2017 as predictors and the CCI2017 as reference. The columns, on the other hand, represent the year of the dataset [reference map] to which the trained model was applied to [compared to], e.g., the value of the upper right corner of the "2017 matrix" (so 0.87), represents the R 2 between: 1) the CCI2018 reference map and 2) the AGB estimate from the ANN D (H, PR, E h ) trained with 2017 data and applied to datasets from 2018.
On average, when the AGB estimates produced by ANN D (H, V, PR) were evaluated with respect to their respective reference map, R 2 of around 0.85 units was obtained; this approach presented slightly lower R 2 values than the other cases. When looking at the remaining ANNs and the temporal extrapolation of models, two patterns were observed. First, both ANN H (H, PR, VOD) and ANN D (H, PR, E h ) presented similar performances (R 2 ∼0.87) for each year. However, the use of E h or VOD as a complementary predictor yielded relatively (but not significantly) higher values than the use of SM; the inclusion of horizontal polarized emissivities as an additional predictor of ANNs slightly increased the R 2 value (in ∼0.03 units) with respect to the ANN D (H, V, PR). Second, the temporal cross-validation of models showed a decrease in R 2 values. Models trained in 2010 and applied to 2018 reduced their R 2 by about 0.2 units, when applied to 2017 datasets the difference was not significant. A similar decrease was observed when applying models trained in 2018 to data from 2010 (∼0.2). The difference in R 2 between 2017 and 2018 models was not significant.
Similar patterns were observed in the RMSE values (not shown here), i.e., 1) an improvement from the inclusion of E h (and similar results from the VOD hybrid case) and 2) a deterioration when applying models to datasets from years different than those used during training. AGB estimates from ANNs trained using (H, V, PR)  Finally, ANNs were trained with learning-sets composed of the three reference years (multiyear learning-set mentioned at the end of Section III-A); the temporal cross-validation of models trained under these conditions is shown in the "3 years" matrix in Fig. 4 (bottom matrix). The evaluation is done using only test subset (see Section III-A). The ANNs trained under this configuration, as in the previously mentioned cases, showed slightly higher results from both the ANN D (H, PR, E h ) and the ANN H (H, PR, VOD) approach (an increase of ∼0.03 units) than from the ANN D (H, V, PR). Additionally, a significant improvement in the cross-validation of models trained with multiple years was observed. The aforementioned interannual differences in R 2 (0.2 for models trained in 2010 [2018] and applied to 2018 [2010] datasets) were not significant under the 3-years training configuration (<0.01 units). When comparing the R 2 (and RMSE) from multiyear trained models with the results from single-year models (trained and applied on data from the same year), a slight (but not significant) decrease in R 2 was observed in some cases. However, to carry out the evaluation of AGB estimates, the models trained with a multiyear learning set (concatenation of three years) were chosen, since they are the ones that showed greater stability for interannual analyses.

B. Evaluation of AGB Estimates From ANNs
Once the quality metrics for the different cases were evaluated, and the optimal training configuration defined, the spatial distribution of residuals was explored (that is, the differences between references and estimates). This analysis was performed exclusively on three cases of interest, more precisely:  different reference years (2010, 2017, and 2018), the obtained AGB estimates were, in turn, compared to their respective reference maps. Given that ANNs were trained using a three-year concatenation, the evaluation was carried out only on the fraction of the original dataset (from each reference year) that was not used during the training stage of the models (test-set in Section III-A). Thus, the ability of the networks to generalize the patterns learned during training (on a set of data unknown to them) was evaluated. The residual maps are presented in Fig. 5. A close-up of the African continent and South America is added.
In general, similar patterns were observed in the different residual maps, in northern latitudes (>60.0 • N), estimated biomass was larger than the references, mainly in the northernmost regions of Russia, China, and Canada. At lower latitudes (Latitude: 40.0 • N -50.0 • N) estimates were lower than the reference maps. Over the pan-tropical region (Latitude: 30.0 • S -30.0 • N), characterized by tropical forests, AGB inversions were lower and higher than the CCI maps; pixels located on the highly humid Pacific coast of Colombia, as well as a large part of the Amazon rainforest, presented larger estimates for the same locations on the AGB references. However, vegetation biomass inversions over the fraction of the Amazon rainforest within the Peruvian border were systematically and in all cases lower than the CCI maps. The line that delimits these two areas of estimates larger-than and lower-than the reference was strictly marked by the Andes mountain range. Two points where AGB inversion surpassed the CCI maps were identified, the first between the border between Bolivia and Paraguay, and the second in the triple border between Colombia, Brazil, and Venezuela.
Overall, vegetation biomass over the Congo rainforest was above and below the values from the reference AGB maps. Australia's southern border is, in most cases, characterized by AGB inversions lower than the references. Finally, biomass estimates over the Malay Archipelago, although represented mostly by tropical forests, are generally larger than the CCI in most inversions in the western islands and lower in the eastern islands.
Finally, the evaluation of the spatial patterns of AGB estimates produced from the selected ANN D and ANN H over tropical Africa and America are shown in Fig. 6. In order to represent the spatial distribution of vegetation biomass, the AGB was regrouped into bins of 35 Mg/ha of width. In the left column, and in descending order, are the AGB maps from: CCI2018, the data-driven cases ANN D (H, V, PR) and ANN D (H, PR, E h ), and the hybrid case ANN H (H, PR, VOD) (the final panel will be discussed in Section IV-C). In the right column the count per AGB bin. All inversions are from models trained with the multiyear dataset and applied to 2018 observations. Just like it was shown by the other quality metrics (  [38]. Similarly, even if the main structural features of reference maps were captured by ANNs, e.g., lower values in northern and southern tropical Africa, and higher values in the more equatorial region of both America (Amazon rainforest) and Africa (Congo forest), subtle spatial variations like the transition from densely to sparsely vegetated areas in the Amazon rainforest (region between Lat: [−10, 0], Lon: [−70, 60]) were not reproduced accurately. Nevertheless, the cluster of high AGB values from the reference map, located near the Atlantic Ocean on the tropical America, was more closely reproduced by ANN D (H, PR, E h ) and ANN D (H, V, PR), than by the hybrid approach. Regardless, ANNs tended to correctly replicate the main features of the distribution of the reference map (histograms on the right-hand side of Fig. 6).

C. Comparison of Estimates From ANN and VOD-Fitted Functions
In view of the demonstrated potential of the methodology here proposed to estimate AGB, a comparison between AGB estimates produced from ANNs with those produced from the fitted VOD function was finally explored. The same evaluation criteria applied to the ANNs were applied to the VOD-based estimates. To ensure a fair comparison, the training of models and the comparison of estimates were performed on the same datasets for both methods.
The fifth row of each matrix in Fig. 4 (R 2 results) and the lower panel of Fig. 5 (distribution of residuals) give a comparison of the performance of the VOD function with respect to ANNs. The R 2 scores between references and AGB estimates from the VOD functions presented similar results independently of year used for the training. Despite this, estimates from the VOD parametric function exhibited larger and lower differences with respect to the CCI maps than those from ANNs in the residual maps. R 2 values from ANN D (H, PR, E h ) are generally ∼0.03 units higher than those from the VOD parametric function; a similar difference was observed when compared to ANN H (H, PR, VOD). The RMSE values associated with the inversion of the VOD function were around 2.8 Mg/ha larger than those on the ANN H (H, PR, VOD), the difference with respect to the ANN D (H, PR, E h ) case was on the order of ∼2 Mg/ha.
The lower panels of Fig. 6 show the spatial patterns and histogram of AGB estimates produced from VOD. Estimates of the parametric function tended to simplify the spatial variations of biomass. In regions with higher AGB densities, the estimates of the parametric function deviated from the references more strongly than those from ANNs, resulting in inversions larger than the reference AGB in the Amazon rainforest while lower than the AGB reference in tropical Africa. Furthermore, while it is true that both methods (the ANNs and the parametric functions) underestimated the highest AGB values (see Fig. 6), it is also true that the ANNs were closer to the latter than the VOD function. In general, the AGB reference distribution was reproduced more accurately by the ANNs (right-hand-side figures in Fig. 6). For example, the Amazon rainforest was mostly represented by values between 315 and 385 Mg/ha under the VOD-based models while according to the estimates of the reference, it can encompass densities between 245 and > 500 Mg/ha.
To get further insight into the differences between estimates and references displayed in Fig. 5, global heatmaps between reference and estimated maps (not limited to Africa and South America tropical region) were computed for all models (including the VOD function) (see Fig. 7). Overall similar distributions were observed; two major clusters of points are visibly separated, one at AGB values below 200 Mg/ha and the other at values around 340 Mg/ha. The scatterplots for ANN D (H, PR, E h ) and ANN H (H, PR, VOD) showed a slightly more linear shape. The estimations with the VOD function showed the highest dispersion.
One might think that a more complex approach such as using an ANN with VOD as only input instead of a parametric function could give better performances; however, an ANN with VOD as input did not change significantly the results obtained with the parametric function. Fig. 8 shows scatterplots of AGB versus VOD (in gray dots) and the AGB predictions from VOD using both a parametric function or an ANN. The marked differences, with respect to the references, from estimates produced through the VOD function are explained by the dispersion of the highest and lowest reference AGB values; VODs ≥1 correspond to AGBs around 300 and 450 Mg/ha. Hence, when the fit of the VOD to AGB relationship (solid black line in Fig. 8) is used to estimate biomass, a value of ∼375 Mg/ha would be predicted from a VOD ≥ 1, which is both larger and lower than the reference AGB values. Similar behavior was observed in the ANN(VOD) (solid red line), which was also limited to describing a rough approximation of the mean distribution of the AGB-VOD relationship, producing a threshold similar to the one mentioned above but with a clear saturation at ∼340 Mg/ha.

D. Uncertainty Estimation
A modified version of the methodology proposed by Aires et al. [39] to estimate the uncertainties of temperature profiles produced by ANN was applied here to AGB estimates. However, Fig. 7. Heatmap between reference AGB (x-axis) and estimated AGB (yaxis); inversions are produced by models trained with a multiyear dataset and applied to 2018 observation, the reference is the CCI2018 map. as noted by the latter, there is currently no simple way to estimate the uncertainties of ANN inversions; thus, the methodology here proposed does not provide uncertainties in the strict sense, but rather an estimation of the errors of the inversions with respect to the reference maps used. To do so, AGB inversions were regrouped in 35 Mg/ha bins from 0 to 490 Mg/ha (estimates greater than 490 Mg/ha were grouped under the same bin); as an example, within the bin with limits 140-175 Mg/ha, all the pixels whose estimates are included between these two values are found. For each estimated value (which corresponds to a pixel with a specific geographic location), there is necessarily a pixel on the reference map with the same position; therefore, there is an estimate-reference pair for each location (red dots in Fig. 9). The uncertainties calculated here correspond to the standard deviation (STD) of the differences between all the estimate-reference pairs contained in each bin. They, the STDs, are centered around the median of the observations/references of the corresponding bin (green dots in Fig. 9). Thus, a metric of the reliability of the inversions was established by intervals of AGB, which describes to what extent the estimated values deviate from the "true" value (the reference maps). Fig. 9 presents the uncertainties for the estimates shown in Fig. 7.
The calculated uncertainties varied according to the AGB interval. On average, across all evaluated cases, the uncertainties' values started at around ±14 Mg/ha for AGBs < 35 Mg/ha and increased to a maximum of ±83 Mg/ha between ∼210 and 245 Mg/ha AGBs. Above this value (245 Mg/ha), they remained in a range between ±50 to ±70 Mg/ha. Among all AGB estimates, the largest uncertainties were produced by the  Fig. 9). The lower uncertainties were those of the data-driven case (ANN D (H, PR, E h )) and the hybrid case (ANN H (H, PR, VOD)).

V. DISCUSSION
According to the analysis presented in Section IV, for a purely data-driven approach, the optimal combination of predictors for the estimation of AGB from ANNs includes the use of TBs in H polarization, E h and PR; higher synergies were observed for incidence angles > 35 • . The use of V or E v as unique (or complementary) predictors had a negligible influence on the accuracy of the estimates (see Fig. 2). It is well known that the behavior of TBs' H and V polarizations (with respect to incidence angle) is close to the Fresnel law in the case of bare soil; same values for H and V are expected at 0 • , just as an increasing difference between them as the incidence angle increases. In contrast, in the presence of vegetation, H and V exhibit closer values to each other at higher and lower angles of incidence [40].
Hence, despite their low performances (when used as only predictors) to estimate biomass, V-polarization TBs carry latent information essential for the study of AGB. This information is effectively extracted by PR, which probably explains the multicollinearity effects of using V (or E v ) with PR. Previous studies have highlighted the correlations of TBs captured at high incidence angles to the presence of vegetation [31], [41]; here, this relationship has been demonstrated from a purely statistical approach (ANN) using the multiangular TBs from SMOS; results showed that the use of TBs from various angles of incidence (8:14 in Fig. 2) can produce more accurate estimates than if a single angle is used (10 in Fig. 2). Moreover, perhaps the similarity between both polarizations at high angles of incidence (over vegetated areas), or the fact that the emission from the ground [canopy] is less [more] predominant on the observed area as the observation angle increases, can explain the higher synergies at high angles of incidence observed in this article.
The use of complementary predictors such as the VOD presented slightly higher quality statistics (R 2 , RMSE) than the pure data-driven model (ANN D (H, PR, E h )). However, the spatial structure of the estimates and the residual maps showed that the hybrid approach does not provide significantly more precise estimates than those produced from the exclusive use of TBs and emissivities. On the other hand, the use of horizontally polarized emissivities (E h ) proved to be useful for a more accurate reproduction of the AGB. This, probably explained by the fact that emissivities depend directly (among other things) on the physical temperature (TS) of the observed object, which in densely vegetated warm-biomes such as tropical forests, can play an important role in the estimation of plant biomass. Interestingly, the use of TS as an external complementary predictor [so ANN D (H, PR, TS)-not shown in this document] slightly reduced R 2 and increased uncertainties with respect to the ANN D (H, PR, E h ); the effect of temperature on the estimation of AGB, seems to be better captured by the angle-dependent E h than by the TS.
Additionally, and thanks to the three years of information available, a temporal cross-validation of AGB estimates was performed (see Fig. 4). Temporally extrapolating models (trained on single-year datasets) showed a clear deterioration in the quality of the estimates when the year of training differs from the year of application of the model. The alternative proposed here, the multiyear training of models, minimizes the dependence to a single reference map and reduces the possible inaccuracies that may be transmitted by one of the selected references. Encouraging preliminary results were found, as this training scheme presents an improvement in the quality of AGB estimates when the latter are compared with different references over time; a potential advantage with respect to the methods used by other studies to estimate the temporal evolution of carbon/AGB stocks, which presume that a spatial relationship found at a given time between VOD and AGB can be used to extrapolate to other periods. This is a strong assumption that can be questioned [42]. However, a thorough study of the influence that reference maps may have on the temporal extrapolation of models is necessary. In any case, the robustness of training the ANNs with multiyear datasets is a significant step in that direction.
In relation to the influence of datasets on AGB inversions, when observing the R 2 values in Fig. 3, it can be seen how estimates produced from ANNs trained on continental learning sets (Africa or South America) tend to have higher correlations with the reference maps than the global ones. Probably, this difference in performance is due to the seasonality of higher latitude pixels and their influence on retrievals done on global observations, as remarked in [43]. It is true that when computing annual averages at a global scale (as is done here), the seasonality of pixels in the northernmost latitudes (where ground temperatures can pass the freezing point) could have an effect on TBs, and therefore, on estimates. However, an additional filter to the TBs before the calculation of annual means was made (pixels with temperatures lower than 275 K or with a snow depth greater than 10 −4 mm were not taken into account), the estimated AGB values over the boreal forest from datasets filtered this way (not shown here) did not diverge significantly from the results obtained with the filters used here (see Section III-A). A more detailed work on the influence of seasonality on the input data and the AGB estimates should be considered in future studies.
Regarding the spatial structure of the AGB estimates presented in Fig. 6, the use of multiangular TBs proved to be useful to reproduce the spatial gradients of biomass in tropical areas, and the employment of E h reduced the discrepancies between references and estimates in boreal forests (see Fig. 5). Finally, the method discussed in the current study presents several potential advantages with respect to the method used by all the published papers studying the AGB and carbon stock evolution from L-VOD (e.g., [10], [22], [25]). First, estimating VOD from PMW is complex and requires estimating the radiation transfer through the vegetation layer, which is far from an easy exercise. In the case of a sensor with multi-incidence angle and full-polarization capabilities, such as SMOS, this implies moving from a high dimensional space (two polarizations times the number of incidence angle bins) to a one-dimensional quantity. Therefore, there is a potential information loss when the final goal is not VOD but AGB. This is probably the origin of the observed degeneration of the VOD to AGB relationship (for a given VOD several AGB values may correspond), which can be partially compensated by using multidimensional data as input. An example is the large dispersion found on AGB estimations from VOD, in particular for large AGB values (see Fig. 8). In other words, a VOD value of 1 can correspond to an AGB from 300 to 450 Mg/ha. When using just the mean value, there are large differences with respect to some of the reference AGB values. To our knowledge, most of the published papers using this method only take into account the mean distribution of these parametric functions, and therefore, possible conclusions in carbon stocks evolution could be within the uncertainties of AGB estimations from VOD. In contrast, when starting from multiangular TBs and a statistical approach without the intermediate step of VOD, it is possible to break that degeneracy in the VOD to AGB relationship, which turns out in significantly lower differences in estimates with the AGB references (see Fig. 5) and lower uncertainties (see Fig. 9). As shown in Fig. 6, these effects can also affect the spatial patterns of AGB. The AGB(VOD) map for the equatorial forest in the Amazonian region shows a flat pattern with very little structure due to the underestimation in some areas and overestimation in other areas, in contrast to the actual structure seen in the reference AGB map, which is much better reproduced by the ANN D (H, PR, E h ).
Therefore, starting from a multidimensional TBs space and a pure data-driven approach [ANN D (H, PR, E h )] to estimate AGB seems to present some advantages with respect to using just VOD. On the other hand, it has been shown that VOD observations (produced from radiation transfer models) can be used as a sophisticated feature extraction approach, thus providing VOD as a predictor to ANNs (in addition to TBs) improves the performance of AGB estimates. In this context, the method becomes hybrid as it makes use of both machine learning and physical modeling. In view of the results presented in this study, the proposed methodology can be considered as an alternative to the methods currently used for the estimation of biomass from PMW, however, more work is needed to perform a more robust estimation of AGB variations with time from VOD or TBs.

VI. CONCLUSION
This article presents an alternative methodology for the estimation of AGB from L-band PMW observations, which directly exploits the statistical relationships between biomass and multiangle TB observations using ANNs. The results show that the estimation of AGB from the proposed method presented certain advantages with respect to the method using parametric functions used by most of the works published to date on the subject. The estimated AGB maps, presented a slightly higher correlation than the current methods, with the reference AGB maps (R 2 ∼0.87 instead of ∼0.84), a stronger agreement with the structural patterns of the reference maps over tropical regions, residual maps with lower absolute values and smaller uncertainties when compared to the state-of-the-art AGB maps (see Fig. 9). The results suggest that accurate biomass estimates can be produced directly from SMOS L-Band multiangular TBs. Being the optimal set of predictors for the estimation of AGB directly from TBs: the TBs in H polarization, the PR, and the E h (for angles of incidence > 35 • ).
In addition, it has been shown that using a single year to establish the relationship between the input variables and AGB introduces significant uncertainties when this relationship is used to study the evolution of AGB/Carbon Stocks. An alternative training scheme was proposed based on multiyear learning sets, which accounted for the differences mentioned above and produced more stable estimates for different years. Despite all of the above, a more detailed study should be carried out on the possible synergies that the observations of different instruments (frequencies) could have with other canopy components, and on the possible applications that this new methodology may have in the future, like the estimation of carbon stocks or AGB temporal evolution. For the last 25 years, he has worked for various geophysical laboratories, putting to stress advanced computer science and applied mathematics paradigms against real problems, particularly in the context of inverse problems applied on remote sensing observations. He is currently with Centre d'Etudes Spatiales de la BIOsphère, Toulouse. His research interests include signal processing, nonlinear modeling, and inverse problem, particularly using artificial neural networks such as for real-time signal processing controller of a radio receiver dedicated to solar wind plasma line tracking onboard the WIND/WAVES spacecraft, or for direct-inverse modeling of ocean surface wind from ERS 1-2 C-band/NSCAT Ku-band scatterometer or biophysical parameters, LAI, chlorophyll, etc., from POLDER optical directional multispectral reflectances, or using more traditional iterative minimization approaches like for soil moisture retrieval from SMOS L-band brightness temperatures, its validation, and future improvements he has been working on since 2005. Stéphane Mermoz received the Ph.D. degree in hydrology and radar remote sensing from the National Institute of Scientific Research, Quebec, QC, Canada, and the University of Rennes 1, Rennes, France, in 2010, with a dissertation on remote sensing of river ice in Canada using polarimetric radar data.
In 2011, he received a postdoctoral grant from the National Centre for Space Studies (CNES) to carry out research works at the Center for the Study of the Biosphere from Space (CESBIO), where he stayed until 2018, in the frame of supporting activities for the Biomass mission as a Postdoctoral Researcher. He created GlobEO in 2018 and is currently working on the development of methodologies to monitor forest and agriculture using remote sensing data. He is particularly involved in forest above-ground biomass estimation and forest loss detection using radar data. He participated in more than 25 projects mainly founded by ESA, CNES, and the European Commission. He also authored or coauthored more than 30 scientific papers in peer-reviewed journals.
Yann H. Kerr [1987][1988]. He was with CESBIO as the Deputy Director from 1995 to 1999 and the Director from 2007 to 2016. He is involved in many space missions. He was an EOS Principal Investigator (interdisciplinary investigations), and PI and precursor of the use of the SCAT over land. In 1989, he started to work on the interferometric concept applied to passive microwave Earth observation and was subsequently the Science Lead on the MIRAS project for ESA with MMS and OMP. His research interests include theory and techniques for microwave and thermal infrared remote sensing of the Earth, with an emphasis on hydrology, water resources management, and vegetation monitoring.
Dr. Kerr was also a Coinvestigator on IRIS, OSIRIS, and HYDROS for NASA. He was a Science Advisor for MIMR and Co I on AMSR. He is a member of the SMAP Science Team. In 1997, he first proposed the natural outcome of the previous MIRAS work with what was to become the SMOS Mission to CNES, proposal that was selected by ESA in 1999, with him as the SMOS mission Lead Investigator and the Chair of the Science Advisory Group. He is also in charge of the SMOS science activities coordination in France. He has organized all the SMOS workshops and was the Guest Editor on three IEEE Special issues and one RSE. He is currently involved in the exploitation of SMOS data, in the Cal Val activities and related level-2 soil moisture and level-3 and -4 development and SMOS Aquarius SMAP synergistic uses and on the soil moisture essential climate variable. He is also working on the SMOS-Next SMOS-HR concepts and is involved in both the Aquarius and SMAP missions. He was nominated as She is currently a Senior Researcher with the Centre d'Etudes Spatiales de la Biosphère (CESBIO), the P.I. of the BIOMASS mission, and the Co-Chair of the ESA BIOMASS Mission Advisory Group. Her research activity has been in the area of microwave remote sensing for land applications, including experimental and modeling studies and SAR image analysis, with a focus on forest and agriculture. Her current research interest focuses on the use of SAR data to monitor changes in the land surfaces under climatic and human impacts.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. He designed the neural network algorithm of the ESA SMOS soil moisture Near Real-Time product and the one that is used for operational data assimilation of SMOS soil moisture at ECMWF. He is also working in aperture synthesis algorithms and new satellite mission design (SMOS-High Resolution, FRESCH Fine Resolution Explorer for Salinity Carbon and Hydrology). His current research concerns microwave remote sensing for Earth observation, in particular machine learning approaches for the estimation of geophysical parameters such as soil moisture and biomass and their assimilation into numerical weather prediction and carbon cycle models.