Above-Ground Biomass Retrieval over Tropical Forests: A Novel GNSS-R Approach with CyGNSS

Abstract

An assessment of the National Aeronautics and Space Administration NASA’s Cyclone Global Navigation Satellite System (CyGNSS) mission for biomass studies is presented in this work on rain, coniferous, dry, and moist tropical forests. The main objective is to investigate the capability of Global Navigation Satellite Systems Reflectometry (GNSS-R) for biomass retrieval over dense forest canopies from a space-borne platform. The potential advantage of CyGNSS, as compared to monostatic Synthetic Aperture Radar (SAR) missions, relies on the increasing signal attenuation by the vegetation cover, which gradually reduces the coherent scattering component ${\mathsf{\sigma }}^{\mathrm{coh},0}$. This term can only be collected in a bistatic radar geometry. This point motivates the study of the relationship between several observables derived from Delay Doppler Maps (DDMs) with Above-Ground Biomass (AGB). This assessment is performed at different elevation angles ${\mathsf{\theta }}_{\mathrm{e}}$ as a function of Canopy Height (CH). The selected biomass products are obtained from data collected by the Geoscience Laser Altimeter System (GLAS) instrument on-board the Ice, Cloud, and land Elevation Satellite (ICESat-1). An analysis based on the first derivative of the experimentally derived polynomial fitting functions shows that the sensitivity requirements of the Trailing Edge $\mathrm{TE}$ and the reflectivity $\Gamma$ reduce with increasing biomass up to ~ 350 and ~ 250 ton/ha over the Congo and Amazon rainforests, respectively. The empirical relationship between $\mathrm{TE}$ and $\Gamma$ with AGB is further evaluated at optimum angular ranges using Soil Moisture Active Passive (SMAP)-derived Vegetation Optical Depth ($\mathrm{VOD}$), and the Polarization Index ($\mathrm{PI}$). Additionally, the potential influence of Soil Moisture Content (SMC) is investigated over forests with low AGB.

1. Introduction

Tropical forests are a key terrestrial ecosystem that play a leading role in the global carbon cycle. They are forested landscapes in tropical regions, i.e., land areas approximately bounded by the tropic of Cancer and Capricorn. Biomass is an Essential Climate Variable (ECV) that determines the spatial distribution of carbon in the biosphere [1]. Quantifying the spatio-temporal fluctuations of the biomass could help to understand the dynamics of tropical forests [2]. On the other hand, biomass is also an important parameter of the water cycle. Biomass affects evapotranspiration, which has a strong impact in freshwater supply in river runoff.

The Above-Ground Biomass (AGB) of tropical forests represents an important carbon pool [2]. Improving our capabilities to estimate AGB over this type of forests is essential to further understand both the carbon and the water cycles. At present, AGB estimation mainly relies on the direct estimation of the geometrical properties of the vegetation structure, such as, e.g., Canopy Height (CH) [3]. In-situ data collection over tropical regions is limited. Characterization on a regional scale of the Earth’s surface is not feasible using in-situ data because of the wide geographical extend and a complicated accessibility. Air/space-borne remote sensing instruments are the most suitable solution. The Phased Array L-band Synthetic Aperture Radar (PALSAR)/Advanced Land Observing Satellite (ALOS) was the first mission that generated systematic measurements of forest biomass. However, it saturates for moderate AGB ~ 100 ton/ha [4]. The future European Space Agency ESAs’ BIOMASS mission [5] is aimed at overcoming this limitation using a P-band fully polarimetric interferometric Synthetic Aperture Radar (InSAR) with multipass capabilities, and a ~12-m deployable reflector antenna. Alternative studies have also suggested combining data from Light Detection and Ranging (LIDAR), e.g., Geoscience Laser Altimeter System (GLAS) instrument on-board the Ice, Cloud, and land Elevation Satellite (ICESat-1) mission [6] with InSAR and/or optical systems [7]. Vegetation indices derived from optical data have been linked to Leaf Area Index (LAI); however, LAI is poorly correlated with biomass [5]. Additionally, optical sensors suffer from weather conditions and clouds, which are a predominant problem over tropical rainforests. Finally, it is worth mentioning that the use of Vegetation Optical Depth ($\mathrm{VOD}$) for biomass estimation is recently under investigation.

Global Navigation Satellite Systems Reflectometry (GNSS-R) is a sort of L-band passive multi-bistatic radar, with as many transmitters as navigation spacecrafts are in view. It provides global coverage, full temporal availability, and sampling of the Earth’s surface over multiple tracks simultaneously within a wide area up to ~ 1500 km. A potential disadvantage is the poor spatial resolution under the incoherent scattering regime. The use of GNSS radio-navigation signals for Earth remote sensing has been investigated since it was originally proposed by ESA in 1993 [8] for mesoscale ocean altimetry [8,9]. The first feasibility study from space corresponds to an experiment performed at the Jet Propulsion Laboratory (JPL) using the Space-borne Imaging Radar-C (SIR-C) on-board the Space Shuttle [10]. This experiment pushed forward the development of GNSS-R. At present, several GNSS-R space-borne studies have been performed from several platforms, including United Kingdom (U.K.) Disaster Monitoring System-1 (DMC-1) [11], U.K. Techdemosat-1 (TDS-1) [12], Soil Moisture Active Passive (SMAP) [13], and Cyclone Global Navigation Satellite System (CyGNSS) 8 microsatellites constellation [14].

Earth-reflected GNSS signals have sensitivity to a wide variety of geophysical parameters, such as, e.g., Snow Water Equivalent (SWE) [15], SMC [16,17] and $\mathrm{VOD}$ [17], flood inundation [18,19,20,21], and surface heat flux [22]. Additionally, some pioneering studies showed a promising sensitivity to forest biomass. Scattering simulations based on the Bistatic MIchigan MIcrowave Canopy Scattering (Bi-MIMICS) model at linear polarization (Horizontal H, Vertical V) suggested a better sensitivity than monostatic configurations for canopy height ~ 8 m [23]. This result was consolidated using the Tor Vergata model at circular polarization (Right Hand Circular Polarization (RHCP), Left Hand Circular Polarization (LHCP)) for higher levels of biomass up to AGB ~ 200 ton/ha [24]. Direct and multiple scattering terms were evaluated using Bi-MIMICS [25] at circular polarization (RHCP, LHCP). It was concluded that the total scattering field at both polarizations RHCP and LHCP is dominated by the scattering over the tree trunks layer. The Soil and Vegetation Reflection Simulator (SAVERS) was developed using the Tor Vergata electromagnetic model [26]. SAVERS includes capabilities to predict signal power as measured by a GNSS-R reflectometer, considering system properties, observation geometry, and applying a polarization synthesis technique to account for antenna polarization mismatch, and cross-polarization isolation [26]. Two air-borne experiments confirmed the sensitivity of the bistatic reflectivity up to high levels of AGB ~ 300 ton/ha [27]. Empirical results over boreal forests from a stratospheric balloon suggested that the coherent scattering term is roughly independent of the platform’s height [28]. The EMISVEG [29] and the Signals of opportunity Coherent Bistatic scattering model for Vegetated terrains (SCoBi-Veg) [30] simulators further analyzed this hypothesis. Lindenmayer systems [31] were used to generate fractal geometry and the Foldy’s approximation [32] was used to account for attenuation and phase change of the coherent wave propagating in the forest media. More recently, a comprehensive study [33] demonstrated a significant sensitivity of several GNSS-R observables up to AGB ~ 150 ton/ha at different elevation angles using the GLObal navigation satellite system Reflectometry Instrument (GLORI) instrument [34].

Feasibility studies for the case of GNSS-R data collected from a space-borne platform have shown a certain sensitivity to forests biomass on a global scale using the SMAP radar receiver [13] and the U.K. TDS-1 GNSS-R payload [35,36]. CH was demonstrated to be a key parameter, with a higher impact in GNSS-R signatures than Normalized Difference Vegetation Index (NDVI). In the present study, the capability of GNSS-R for biomass monitoring over tropical forests is evaluated using the CyGNSS Level 1 Science Data Record Version 2.1 [14,37,38]. The sensitivity of several scientific observables [13] derived from space-borne Delay Doppler Maps (DDMs) over dense forest canopies is here analyzed in a more quantitative manner. The selected biomass reference datasets are the integrated pantropical AGB map, as provided by Avitabile et al. [39], and ICESat-1/GLAS derived global estimates of forest CH [3,6]. Avitabile et al. AGB map was derived from ICESat-1/GLAS. As such, both biomass references were selected to be used simultaneously in this work. The following questions are addressed:

• What is the DDMs’ response to the interaction of GNSS signals with tropical forests?
• What is the GNSS signals saturation level (if any) over tropical forests?
• What is the optimum GNSS satellites elevation angle for biomass estimation over tropical forests?

This work can be understood as a follow-on study of the sensitivity analysis performed in [17]. Here, the main interest is on the sensitivity of $\Gamma$ and $\mathrm{TE}$ to AGB over tropical forests as a function of CH, $\mathrm{PI}$, and $\mathrm{VOD}$. The study in [17] was focused on the analysis of the relationship between $\Gamma$ and the normalized first Stokes parameter over a wide variety of target areas. The “tau-omega” model was found to properly fit this relationship. Several improvements have been introduced: (a) the use of raw DDMs, (b) signal calibration, and (c) the use of several observation geometries. This paper is organized as follows. Section 2 motivates the use of L-band bistatic radar for biomass studies. Section 3 describes the selected observables and the methodology. Section 4 analyses the sensitivity of CyGNSS to AGB, including a comprehensive interpretation of the experimental results. The main conclusions are included in Section 5. Information about test sites and reference data is included in Appendix A and Appendix B, respectively. Finally, Appendix C includes supplemental information for the evaluation of results in Section 4.

2. Exploiting CyGNSS for Biomass Monitoring

2.1. Background

The CyGNSS mission (down looking antenna gain ~ 14.5 dB, LHCP) was first proposed for ocean wind speed estimation over tropical cyclones using GPS L1 C/A signals of opportunity. The orbital configuration of each CyGNSS satellite is a circular Low Earth Orbit (LEO) with a height ~ 500 km and an inclination angle of ~ 35°. Each single Delay Doppler Mapping Instrument (DDMI) [12,14] aboard the 8 CyGNSS microsatellites collects forward scattered signals along four specular directions (elevation angle of the incident wave ${\mathsf{\theta }}_{\mathrm{e},\mathrm{i}}$ equals elevation angle of the reflected wave ${\mathsf{\theta }}_{\mathrm{e},\mathrm{s}}$; ${\mathsf{\theta }}_{\mathrm{e},\mathrm{i}}$ = ${\mathsf{\theta }}_{\mathrm{e},\mathrm{s}}$ = ${\mathsf{\theta }}_{\mathrm{e}}$) corresponding to four different transmitting GPS spacecrafts, simultaneously (Figure 1). As such, CyGNSS allows one to sample the Earth’s surface along 32 tracks simultaneously, within a wide range of satellites’ elevation angles ${\mathsf{\theta }}_{\mathrm{e}}$ ~ [20, 90]° over tropical latitudes ~ [−40, 40]°.

The retrieval of geophysical parameters from CyGNSS requires that the DDMI payload cross-correlates directly the collected signals against a known replica of the GPS L1 C/A code, after compensating for the Doppler shift [8,12,14]. Consequently, 32 complex reflected DDMs ${\mathrm{Y}}_{\mathrm{r}}$ are generated on-board at a time, using a short coherent integration time ${{\displaystyle \mathrm{T}}}_{\mathrm{c}}$ ~ 1 ms. Complex DDMs ${\mathrm{Y}}_{\mathrm{r}}$ are incoherently averaged along ~ 1 s to finally generate the so-called power DDMs $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$.

Power DDMs $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ can be modelled using geometrical and scattering related parameters as follows [40]:

$$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle =\frac{{{\displaystyle {{\displaystyle \mathrm{T}}}_{\mathrm{c}}^2\mathrm{P}}}_{\mathrm{T}}{{\displaystyle \mathsf{\lambda }}}^2}{{{\displaystyle (4\mathsf{\pi })}}^3}{\displaystyle \iint \frac{{\mathrm{G}}_{\mathrm{T}}{\mathrm{G}}_{\mathrm{R}}{{\displaystyle \left|\mathsf{\chi }(\mathsf{\tau }-({\mathrm{R}}_{\mathrm{T}}+{\mathrm{R}}_{\mathrm{R}})/\mathrm{c},{\mathrm{f}}_{\mathrm{D}}-{\mathrm{f}}_{\mathrm{c}})\right|}}^2{\mathsf{\sigma }}^0}{{{\displaystyle \mathrm{R}}}_{\mathrm{T}}^2{{\displaystyle \mathrm{R}}}_{\mathrm{R}}^2}}{{\displaystyle \mathrm{d}}}^2\overline{\mathsf{\rho }},$$

(1)

where ${{\displaystyle \mathrm{P}}}_{\mathrm{T}}$ is the power of the RHCP transmitted signals, $\mathsf{\lambda }$ ~ 19 cm is the wavelength of the signals at L1, ${{\displaystyle \mathrm{G}}}_{\mathrm{T}}$ and ${{\displaystyle \mathrm{G}}}_{\mathrm{R}}$ are the transmitting and receiving antenna gains, ${{\displaystyle \mathrm{R}}}_{\mathrm{T}}$, and ${{\displaystyle \mathrm{R}}}_{\mathrm{R}}$ are the ranges from the transmitter and the receiver to the specular point, respectively, $\mathsf{\chi }$ is the Woodward Ambiguity Function (WAF), $\mathsf{\tau }$ is the delay of the signal from the transmitter to the receiver, ${\mathrm{f}}_{\mathrm{D}}$ is the Doppler shift of the reflected signal, ${\mathrm{f}}_{\mathrm{c}}$ is aimed to compensate the Doppler shift of the signal, ${\mathsf{\sigma }}^0$ is the bistatic scattering coefficient at LHCP, and $\overline{\mathsf{\rho }}$ is the positioning vector of the scattering point.

The bistatic scattering coefficient is defined as follows [41]:

$${\mathsf{\sigma }}^0={\mathsf{\sigma }}^{\mathrm{coh},0}+{\mathsf{\sigma }}^{\mathrm{incoh},0},$$

(2)

where ${\mathsf{\sigma }}^{\mathrm{coh},0}$ and ${\mathsf{\sigma }}^{\mathrm{incoh},0}$ are the coherent and the incoherent scattering terms, respectively. ${\mathsf{\sigma }}^{\mathrm{incoh},0}$ describes scattering properties of the surface and the vegetation, i.e., volume scattering. ${\mathsf{\sigma }}^{\mathrm{coh},0}$ depends mainly on surface properties and it has a significantly narrower (practically, a delta function) power distribution over scattering angles than ${\mathsf{\sigma }}^{\mathrm{incoh},0}$. Consequently, DDMs consists on a sum of two terms, as follows [28,42,43]:

$$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle =\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle +\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle .$$

(3)

$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ accounts for coherent reflections, while $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ is responsible for the diffuse scattering. The relative portion of both components mainly depends on the coherent integration time ${{\displaystyle \mathrm{T}}}_{\mathrm{c}}$, the observation geometry, and both the dielectric and structural properties of the scattering media.

Several works found a transition from incoherent to coherent scattering regime over the ocean surface [43,44], where $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ is dominant for high ${\mathsf{\theta }}_{\mathrm{e}}$ [45]. Theoretical simulations have been carried out to evaluate the behavior of the incoherent reflected power as a function of ${\mathsf{\theta }}_{\mathrm{e}}$ [45]. It was stated that $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ decreases ~20% when ${\mathsf{\theta }}_{\mathrm{e}}$ decreases from ${\mathsf{\theta }}_{\mathrm{e}}$ ~90° to ${\mathsf{\theta }}_{\mathrm{e}}$ ~50°. Over land surfaces, global scale results [46] have shown that $\Gamma$ decreases with decreasing ${\mathsf{\theta }}_{\mathrm{e}}$ down to ~55°. This empirical observation was justified because of the effect of the incoherent scattering term ${\mathsf{\sigma }}^{\mathrm{incoh},0}$. On the other hand, $\Gamma$ increases for lower ${\mathsf{\theta }}_{\mathrm{e}}$ because of the increasing signal coherence.

The reflectivity $\Gamma$ accounts for contributions from both coherent and incoherent scattering terms along both angular ranges. The most relevant difference between both scattering terms is on the combination of the electromagnetic field vectors. ${\mathsf{\sigma }}^{\mathrm{coh},0}$ is the result of the coherent combination of the signals reflected on individual facets within the first Fresnel zone. ${\mathsf{\sigma }}^{\mathrm{incoh},0}$ comes from the random combination of electromagnetic signals scattered within the glistening zone that add together at the GNSS-R receiver. When the combination is fully coherent, ${\mathsf{\sigma }}^{\mathrm{coh},0}$ accounts for all the reflected power. In the occurrence of diffuse scattering, the random signs of the electric field cross-products cancel out ${\mathsf{\sigma }}^{\mathrm{coh},0}$, and the reflected power is provided by ${\mathsf{\sigma }}^{\mathrm{incoh},0}$. The changes in the angular evolution of $\Gamma$ depend on the dominant scattering mechanism.

2.2. Motivation

It is hypothesized that GNSS-R measurements have a certain sensitivity to biomass over tropical forests (Figure 2) and this sensitivity can be higher than that of current SAR missions. This hypothesis is based on the following statements: (a) The backscatter intensity measured by monostatic SAR increases with biomass up to a saturation level AGB ~ 100 ton/ha. (b) L-band GNSS signals can partially penetrate the vegetation cover because of the semi transparency of the vegetation at low frequency bands. (c) The bistatic scattering pattern of the coherent surface scattering component follows a delta function along the specular direction, while the incoherent scattered signals spread along all other directions. (d) The coherent component $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ can be higher than the incoherent one $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ even from a space-borne platform.

$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ includes coherent scattering from the surface. It can be expressed analytically as follows [47]:

$$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle =\frac{{{\displaystyle {{\displaystyle \mathrm{T}}}_{\mathrm{c}}^2\mathrm{P}}}_{\mathrm{T}}{{\displaystyle \mathsf{\lambda }}}^2{{\displaystyle \mathrm{G}}}_{\mathrm{T}}{{\displaystyle \mathrm{G}}}_{\mathrm{R}}{{\displaystyle \left|\mathsf{\chi }(\mathsf{\tau },\mathrm{f})\right|}}^2}{{{\displaystyle (4\mathsf{\pi })}}^2{({\mathrm{R}}_{\mathrm{T}}+{\mathrm{R}}_{\mathrm{R}})}^2}{\left|\mathrm{R}({\mathsf{\theta }}_{\mathrm{e}})\right|}^2\gamma \exp (-{(2\mathrm{k}\mathsf{\sigma }\sin {\mathsf{\theta }}_{\mathrm{e}})}^2),$$

(4)

where $\mathrm{R}$ is the LHCP Fresnel reflection coefficient, $\mathrm{k}$ is the signal angular wavenumber, $\mathsf{\sigma }$ is the surface height standard deviation (related to surface roughness), and $\gamma$ is the transmissivity of the vegetation:

$$\gamma =\exp (-2\mathrm{VOD}/\sin {\mathsf{\theta }}_{\mathrm{e}}),$$

(5)

where $\mathrm{VOD}$ is the Vegetation Optical Depth. $\mathrm{VOD}$ is an index of leaf and woody biomass that is used to parameterize absorption effects. $\mathrm{VOD}$ data derived at L-band are better correlated to AGB and CH than C-band data because of the higher signal penetration depth at lower frequency bands [48]. C-band is only sensitive to top canopy layers, while L-band is sensitive to the canopy from top to bottom layers [48]. Increasing levels of biomass belong to higher signal attenuation effects. On the other hand, optically derived vegetation parameters, such as, e.g., NDVI are only correlated with green leaf material. Tropical rainforests are characterized by wet biomass and the impact of signal attenuation through the vegetation cover is thus higher than scattering effects [49]. This attenuation becomes higher for decreasing ${\mathsf{\theta }}_{\mathrm{e}}$ because of the larger signal path along the vegetation cover. The impact of ${\mathsf{\theta }}_{\mathrm{e}}$ on $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ is twofold. On one hand it influences the Fresnel reflection coefficient, and on the other hand it influences the exponential decaying factor in Equation (4). For a totally flat surface, the root mean square of the surface height variation is zero, and $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ equals the square of the amplitude of the Fresnel reflection coefficient. $\mathrm{R}$ is constant for ${\mathsf{\theta }}_{\mathrm{e}}$ > 40° and it decreases significantly for lower angles [27]. However, the exponential factor increases significantly with lower ${\mathsf{\theta }}_{\mathrm{e}}$. As a consequence, $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ increases at grazing angles [46].

Finally, it is worth commenting that $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ can be assumed to be roughly independent of the platform’s height in a LEO scenario because ${\mathrm{R}}_{\mathrm{T}}\gg {\mathrm{R}}_{\mathrm{R}}$. $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ suffers lower free space losses than $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$, following the classical bistatic radar equation.

$\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ includes two main contributions: (a) surface scattering over areas with moderate to high roughness (small-scale roughness) and topography, and (b) volume scattering [29] from the vegetation cover including leaves, branches, and trunks. The impact of double bouncing scattering can be assumed to be negligible because the LHCP down looking antennas filter out the RHCP component of the total scattered electromagnetic field. The polarization of GNSS signals changes from RHCP to LHCP after direct scattering. However, they become RHCP after double reflections, first from RHCP to LHCP and then from LHCP to RHCP.

3. Methodology

3.1. Definition of the Selected CyGNSS Observables

CyGNSS Level 2.1 Science Data Record [14,37,38] was selected for this study. The original data were first filtered out using an equivalent “CyGNSS overall quality flag” over land surfaces to improve the quality of the observables. Reflected delay waveforms $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{raw}}$ were obtained from the original DDMs $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ at zero Doppler frequency:

$$\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{raw}}=\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f}=0)\right|}}^2\rangle .$$

(6)

The delay bin resolution of the original 17 lag waveforms $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{raw}}$ is ~0.2552 GPS C/A chips. A re-sampling and interpolation of the resulting 1700 lag waveforms was performed using a spline method [50] to increase the accuracy of the waveforms, before applying the algorithms to extract the observables. Several observables were considered to evaluate the sensitivity of ${\mathsf{\sigma }}^{\mathrm{coh},0}$ and the volume scattering term to forests biomass.

The trailing edge width $\mathrm{TE}$ was defined as the lag difference between the 70% power threshold of the high-resolution waveforms $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{threshold}}$ and the corresponding maximum power of the waveforms $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{peak}}$ [13]:

$$\mathrm{TE}={\mathsf{\tau }}_{\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{threshold}}}-{\mathsf{\tau }}_{\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{peak}}}.$$

(7)

Different power thresholds were tested. The 50% power threshold was found to cut off some lags when this threshold was out of the original 17 lag waveforms. On the other hand, the 90% power threshold provided a lower dynamic range as compare to the 70%. The incoherent scattering term ${\mathsf{\sigma }}^{\mathrm{incoh},0}$ is the main contribution to $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{threshold}}$, while both the coherent ${\mathsf{\sigma }}^{\mathrm{coh},0}$ and the incoherent ${\mathsf{\sigma }}^{\mathrm{incoh},0}$ scattering terms significantly contribute to the peak power of the waveforms $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{peak}}$. $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ are composed of both coherent and incoherent scattering terms over tropical forests. The width of these power waveforms depends on the scattering properties, in a differentiated manner depending on the coherent-to-incoherent scattering ratio. The width is expected to be small if the ratio is high because of the strong power contribution of the coherent term $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ to the peak of the power waveforms.

The reflectivity $\Gamma$ is estimated as the ratio of the reflected $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{peak}}$ and the direct $\mathrm{W}{\mathrm{F}}_{\mathrm{d},\mathrm{peak}}$ power waveforms peaks, after compensating for the noise power floor and the antennas’ gain patterns as a function of the elevation angle ${\mathsf{\theta }}_{\mathrm{e}}$ as:

$$\Gamma =\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{peak}}/\mathrm{W}{\mathrm{F}}_{\mathrm{d},\mathrm{peak}}.$$

(8)

The antennas’ gain patterns were compensated versus the gain at the corresponding boresight direction (down looking gain ∼14.5 dB, ${\mathsf{\theta }}_{\mathrm{e}}$ ~62° and up looking gain ∼4.7 dB, ${\mathsf{\theta }}_{\mathrm{e}}$ ~90°), and the difference of both gains at boresight. The compensation of the antennas’ gain was implemented as a function of ${\mathsf{\theta }}_{\mathrm{e}}$. This is relevant for an accurate estimation of $\Gamma$ because the transmitted signal power depends on ${\mathsf{\theta }}_{\mathrm{e}}$, and because both gain patterns have a different relationship with this variable. The CyGNSS Level 1 Science Data Record variables used in this procedure are highlighted here: raw DDM and Zenith Signal-to-Noise Ratio (SNR) for the estimation of $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{peak}}$ and $\mathrm{W}{\mathrm{F}}_{\mathrm{d},\mathrm{peak}}$, while Specular point Rx antenna gain for the information of the down looking antenna gain in the direction of the specular point. The up looking antenna is omnidirectional with a gain ∼4.7 dB at boresight and a half power beam width of ∼57°. The noise power floor was estimated as the mean value of the DDM subset where the signal is not present [51]. The delay separation between the DDM subset and DDM peak was at least 0.75 chips, and estimation of $\Gamma$ was discarded if this distance was shorter than 0.75 chips.

The half power beam width (−3 dB) of each down-looking CyGNSS antenna is ~30°, and the −6 dB field-of-view is ~45°. On the other hand, each antenna pointed to the Earth’s surface with an elevation angle of ${\mathsf{\theta }}_{\mathrm{e}}$ ~62° (antenna boresight). As such, GNSS signals reflected at low angles ${\mathsf{\theta }}_{\mathrm{e}}$ down to ~20° can be accurately collected through the main lobe of each antenna.

The nominal mission lifetime is 2 years. This milestone was achieved on March 2019. Experience with CyGNSS wind measurements shows that absolute power measurements is currently one of the challenges of this mission [52,53]. In this work, the direct signal was used to partially improve the estimation of $\Gamma$. At present, the mission is operating nominally, and the Automatic Gain Control (AGC) for the Zenith channel was disabled from August 2018. This fact helps to improve the accuracy in the estimation of $\Gamma$ using Equation (8). On the other hand, $\Gamma$ could also be inverted from scattering models. However, this approach relies on the assumption that the scattering is totally coherent [51] or totally incoherent [54]. This was the main motivation to use Equation (8) in the estimation of the reflectivity over forests regions, where $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ are composed of both coherent $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ and incoherent $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ scattering terms.

3.2. Auxiliary Forests Vegetation Parameters

The SMAP-Enhanced L3 Radiometer Global Daily 9 km Level L3 SPL3SMP_E Version 1.0 product was used to generate auxiliary vegetation parameters. This product was derived from SMAP radiometer (6 A.M. – 6 P.M. data in different arrays) and ancillary data over the global 9-km Equal Area Scalable Earth (EASE 2-0) grid [55,56,57]. The SMAP L-band mission provides a better signal penetration depth over the vegetation as compared to higher frequency bands. As such, it is an adequate candidate to evaluate the relationship between CyGNSS derived $\mathrm{TE}$ and $\Gamma$ with AGB, as a function of $\mathrm{VOD}$, $\mathrm{PI}$, and SMC.

(1) Vegetation Optical Depth ($\mathrm{VOD}$): The retrieval is based on a priori NDVI information obtained from visible-near infrared reflectance data from the NPP/JPSS VIIRS instrument, and land cover type assumptions [55,56,57]. It is used to retrieve the transmissivity $\gamma$ of the vegetation as an estimation of the attenuation of the electromagnetic signal through the vegetation layer. SMOS-derived $\mathrm{VOD}$ has been also used in, e.g., [58].

(2) Polarization Index ($\mathrm{PI}$): It is defined as the difference between the microwave radiometric brightness temperatures at V-pol ${\mathrm{T}}_{\mathrm{BV}}$ and H-pol ${\mathrm{T}}_{\mathrm{BH}}$, normalized to their average value as follows [59,60,61,62]:

$$\mathrm{PI}=\frac{{\mathrm{T}}_{\mathrm{BV}}-{\mathrm{T}}_{\mathrm{BH}}}{({\mathrm{T}}_{\mathrm{BV}}+{\mathrm{T}}_{\mathrm{BH}})/2}.$$

(9)

$\mathrm{PI}$ normalizes the measurements of the brightness temperatures and it is independent of the physical temperature of the surface. The potential sensitivity of $\mathrm{PI}$ to dense forest canopies is based on the assumption that V-polarized emission from soil is attenuated by vegetation, which, on the other hand, emits much less polarized radiation depending on the amount of biomass. This parameter was found to decrease when the biomass increases, independently of the crop type. The performance over tropical forests is evaluated here.

(3) Soil Moisture Content (SMC): The retrieval of SMC is based on the use of the V-pol single channel algorithm [55,56,57] only when favorable surface conditions appear at a given grid cell. Corrections for surface roughness, effective soil temperature, and Vegetation Water Content (VWC) were also applied.

3.3. Strategy

$\mathrm{TE}$ and $\Gamma$ were classified into different groups according to different ranges of the satellites’ elevation angles ${\mathsf{\theta }}_{\mathrm{e}}$ in steps of ~10° from ${\mathsf{\theta }}_{\mathrm{e}}$ ~90° to ${\mathsf{\theta }}_{\mathrm{e}}$ ~20°. ${\mathsf{\theta }}_{\mathrm{e}}$ is an important parameter that remarkably affects the ratio of the coherent $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ to incoherent $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ scattering components. The contribution of the topography to the incoherent scattering term $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ was partially filtered out to evaluate only the effect of the volume scattering term. In so doing, only pixels over the surface with a Terrain Ruggedness Index (TRI) lower than ~ 15 were considered [63,64]. This threshold was found, empirically, to limit the effect of topography. Several TRI thresholds were tested. It was found that levels higher than ~ 15 significantly increase the spreading of the $\mathrm{TE}$. On the other hand, there is a remaining incoherent scattering contribution due to the small-scale surface roughness. However, the coherent scattering term $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ largely overpass the incoherent scattering term $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ associated to the small-scale roughness over land surfaces [28,42,47]. Topography is an important factor that disturbs ${\mathsf{\sigma }}^0$ because the local surface slopes modify the size of the scattering area [65].

Quantifying and analyzing the relationship between $\mathrm{TE}$ and $\Gamma$ against AGB and CH is the main focus of this work. In so doing, SMAP-derived $\mathrm{VOD}$, $\mathrm{PI}$, and SMC products were also used to further analyze the results. A common grid was thus required, because of the variety of products with different footprint sizes and spatio-temporal sampling characteristics. A 0.1° × 0.1° latitude/longitude grid was selected. All products were averaged using a moving window of 0.2° in steps of 0.1°. The spatial resolution is ~ 20 km × 20 km at equatorial latitudes, which is wide enough to account for the across-track spreading ~20 km of the DDMs due to $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$. The ultimate resolution at a cell center is ~ 20 km × 20 km because multiple satellites overpasses with different orientations were considered. This grid also allows us to account for $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$. This term is limited by the along-track resolution ~7.6 km, related with the orbital motion and the incoherent integration time ~1 s. Finally, it is worth pointing out that this strategy provides a filtering of potential fluctuations of $\Gamma$ and $\mathrm{TE}$ due to remaining (after quality flags removal) noise sources, such as Attitude Determination and Control Subsystem (ADCS), the geolocation of the nominal specular reflection points, and the use of eight different DDMs sources (eight distributed microsatellites) in the moving averaging filter.

The selected temporal window was 6 months (01/08/2018 to 31/01/2019). This temporal length is large enough to provide statistically significant results. This specific time period was selected because the AGC was disabled on July 2018, so as to improve the performance in the estimation of $\Gamma$ using Equation (8). The temporal length was selected as large as possible at the time of starting this work because this study is not focused on the analysis of potential temporal fluctuations of different geophysical parameters. The temporal variations of AGB and CH during this period were assumed to be negligible.

4. Performance Analysis

4.1. Introduction

CyGNSS-derived $\mathrm{TE}$ and $\Gamma$ were used in this study to improve our understanding in the performance of GNSS-R for biomass studies (Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8). The analysis was performed over different types of forests, such as Congo rain (Figure 3a–d, Figure 4a–d, and Figure 5a–d), Amazon rain (Figure 3e-h, Figure 4e-h, and Figure 5e-h), coniferous (Figure 6a–d, Figure 7a–d, and Figure 8a–d), dry (Figure 6e-h, Figure 7e-h, and Figure 8e-h), and moist (Figure 6i–l, Figure 7i–l, and Figure 8i–l) tropical forests. These target areas were selected so as to cover a wide variety of forests on a pantropical scale.

AGB and CH data were obtained from the Avitabile et al. product [39] and the ICESat-1/GLAS space-borne instrument [3,6], respectively. Both biomass references were displayed over the specific target areas, enabling a first qualitative evaluation (Figure 3 and Figure 6). A certain correlation of $\mathrm{TE}$ and $\Gamma$ with AGB and CH was found over the selected study sites (

Table 1. Statistics (Correlation (Spearman for rain, Pearson for other forests), root mean square error (RMSE), order polynomial “O” fit, and slope for linear regressions) of the functional correlation between trailing edge width $\mathrm{TE}$ and reflectivity $\Gamma$ with Above-Ground Biomass (AGB) over Congo (CG) and Amazon (AM) rain, coniferous (CF), dry (DR), and moist (MT) tropical forests at different elevation angles ${\mathsf{\theta }}_{\mathrm{e}}$ ~[80, 90]°, ${\mathsf{\theta }}_{\mathrm{e}}$ ~[60, 70]°, ${\mathsf{\theta }}_{\mathrm{e}}$ ~[40, 50]°, and ${\mathsf{\theta }}_{\mathrm{e}}$ ~[20, 30]°.

	${\mathsf{\theta }}_{\mathrm{e}}$ ~[80, 90]°								${\mathsf{\theta }}_{\mathrm{e}}$ ~[60, 70]°
	$\mathrm{TE}$				$\Gamma$				$\mathrm{TE}$				$\Gamma$
	r	RMSE [ton/ha]	O	Slope [(ton/ha)/m]	r	RMSE [ton/ha]	O	Slope [(ton/ha)/dB]	r	RMSE [ton/ha]	O	Slope [(ton/ha)/m]	r	RMSE [ton/ha]	O	Slope [(ton/ha)/dB]
CG	0.73	40.6	3	-	−0.78	26.8	4	-	0.91	31.5	3	-	−0.92	29.8	4	-
AM	0.75	33.7	4	-	−0.93	13.5	4	-	0.68	26.6	6	-	−0.88	13.6	6	-
CF	0.36	12.2	1	0.0131	−0.69	6.8	1	−1.1201	0.34	9.2	1	0.0112	−0.67	8.4	1	−1.0972
DR	0.38	11.4	1	0.0164	−0.62	7.6	1	−1.3260	0.53	12.3	1	0.0176	−0.67	5.6	1	−1.3675
MT	0.25	5	1	0.0087	−0.26	5.3	1	−0.2802	0.10	3.9	1	0.0035	−0.43	4.2	1	−0.4402
	${\mathsf{\theta }}_{\mathrm{e}}$ ~[40, 50]°								${\mathsf{\theta }}_{\mathrm{e}}$ ~[20, 30]°
	$\mathrm{TE}$				$\Gamma$				$\mathrm{TE}$				$\Gamma$
	r	RMSE [ton/ha]	O	Slope [(ton/ha)/m]	r	RMSE [ton/ha]	O	Slope [(ton/ha)/dB]	r	RMSE [ton/ha]	O	Slope [(ton/ha)/m]	r	RMSE [ton/ha]	O	Slope [(ton/ha)/dB]
CG	0.89	33.4	3	-	−0.65	24.1	4	-	0.89	40.1	5	-	−0.66	41.2	4	-
AM	0.84	24.6	4	-	−0.82	16.6	4	-	0.90	18.4	4	-	−0.93	16.9	4	-
CG	0.25	8.6	1	0.0085	−0.64	6.7	1	−1.0194	0.11	5.9	1	0.0026	0.34	9.9	1	−1.0666
DR	0.47	7.5	1	0.0173	−0.66	6.3	1	−1.1877	0.20	7.1	1	0.0080	0	6.7	1	0.0452
MT	0.40	4.3	1	0.0108	−0.55	4.7	1	−0.6470	0.36	3.8	1	0.0078	−0.53	3.7	1	−0.6118

). This correlation improves over more dense forest canopies. $\Gamma$ decreases with increasing biomass levels. On the other hand, the spreading of $\mathrm{TE}$ increases with increasing levels of biomass. Avitabile’s map and CyGNSS sampling properties show a good performance, while GLAS-derived CH data are provided with less spatial density. This observation motivates the selection of AGB from Avitabile et al. as the main biomass reference. CH was used as an auxiliary product to assess the potential impact of this term in the interaction of GNSS signals with the vegetation.

4.2. Sensitivity Analysis as a Function of the Elevation Angle

Figure 4, Figure 5, Figure 7 and Figure 8 show the scatter plots of $\mathrm{TE}$ and $\Gamma$ with AGB over Congo rain (Figure 4a–d and Figure 5a–d), Amazon rain (Figure 4e–h and Figure 5e–h), coniferous (Figure 7a–d and Figure 8a–d), dry (Figure 7e–h and Figure 8e–h), and moist (Figure 7i–l and Figure 8i–l) tropical forests. These plots are depicted as a function of CH. They were computed using the mean values of the biomass reference data in steps of $\mathrm{TE}$ ~ 1 m and $\Gamma$ ~ 0.05 dB. The values of these steps were found empirically. This strategy allowed us to filter noise at a pixel level, such as speckle in ${\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})$, and potential errors in GLAS-derived AGB. On the other hand, it provided enough sampling density to compute the correlation coefficients, i.e., Spearman for rainforests, and Pearson for other forests (Table 1). Different steps were also tested. Decreasing values of these steps belonged to a noisier behavior in the scatter plots, while increasing values reduced sampling density. The extension of the target areas is large, providing enough data. The errors were assumed to be randomly distributed. Thus, noise could be reduced after averaging, so as to find the underlying relationship between CyGNSS observables and AGB. This work is focused on the analysis of this empirical relationship on a regional scale. The correlation between the mean values of both CyGNSS observables and AGB is clearly visible. Increasing levels of AGB belongs to decreasing $\Gamma$ and increasing $\mathrm{TE}$.

Speckle is a source of noise that involves diffuse scattering. If the scattering medium is rough with respect to $\mathsf{\lambda }$, the different heights of the scatterers over the glistening zone will randomly shift the phases of the GNSS signals. Some of these paths will destructively interfere and others will constructively interfere. Thus, the signal power level as measured by CyGNSS will change randomly. From this point, it can be stated that the SNR of the reflected signals is affected by tropical forests. The strategy described previously was applied to mitigate this effect. At the same time, the AGB reference data could potentially have some errors because this dataset is derived from CH data and assumed continental allometric relationships. Nonetheless the validation of Avitabile’s map [39] shows a lower root mean square error (RMSE) (reduction down to ~74%) and bias (reduction down to ~153%) than both original datasets [66,67] over all the continents. An interpretation of the results is provided for each study site:

(1) Tropical rainforests (Figure 4 and Figure 5): The empirical relationships of $\mathrm{TE}$ and $\Gamma$ with AGB were fitted by empirically derived polynomial functions (Table 1) at different angular ranges ${\mathsf{\theta }}_{\mathrm{e}}$ ~[80, 90]°, ~[60, 70]°, ~[40, 50]°, and ~[20, 30]° over Congo (Figure 4a–d and Figure 5a–d) and Amazon (Figure 4e–h and Figure 5e–h). The first derivative of these functions was computed at different mean levels of AGB in steps of ~50 ton/ha, for different angular ranges (

Table A1. First-derivative [(ton/ha)/m] polynomial fitting functions of CyGNSS derived trailing edge width $\mathrm{TE}$ over Congo rainforests.

AGB [ton/ha]	50	100	150	200	250	300	350
[20, 30]°	1.53	1.21	0.90	0.55	0.35	0.33	0.30
[40, 50]°	-	1.05	0.93	0.77	0.58	0.37	0.17
[60, 70]°	-	1	0.87	0.71	0.53	0.35	0.19
[80, 90]°	-	-	-	0.62	0.45	0.22	-

,

Table A2. First-derivative [(ton/ha)/m] polynomial fitting functions of CyGNSS derived trailing edge width $\mathrm{TE}$ over Amazon rainforests.

AGB [ton/ha]	50	100	150	200	250	300	350
[20, 30]°	-	1.07	0.75	0.40	0.10	-	-
[40, 50]°	-	1.15	0.85	0.53	0.18	-	-
[60, 70]°	-	-	0.70	0.60	0.20	-	-
[80, 90]°	-	-	0.70	0.42	0.12	-	-

,

Table A3. First derivative [(ton/ha)/dB] of the polynomial fitting functions of CyGNSS derived reflectivity $\Gamma$ over Congo rainforests.

AGB [ton/ha]	50	100	150	200	250	300	350
[20, 30]°	-	-	-	-	17.1	15.5	-
[40, 50]°	-	-	-	19.3	25.2	27.4	23.1
[60, 70]°	-	-	10.5	21.2	26	26.1	17.5
[80, 90]°	-	-	-	-	16.2	21.3	16.2

and

Table A4. First derivative [(ton/ha)/dB] of the polynomial fitting functions of CyGNSS derived reflectivity $\Gamma$ over Amazon rainforests.

AGB [ton/ha]	50	100	150	200	250	300	350
[20, 30]°	-	-	15.8	16	9.5	-	-
[40, 50]°	-	17.2	18	16	11.5	9	-
[60, 70]°	-	16	17.8	16	10.4	-	-
[80, 90]°	-	-	15.8	14	9.5	-	-

). The values of the derivative reduce with increasing AGB levels. Uncertainties in AGB estimation were computed as the values of the first derivatives times the uncertainties in the GNSS-R observables. This point is translated into smaller uncertainties in the estimation of AGB with increasing biomass levels. A budget (ignoring noise) is provided assuming uncertainties of, e.g., $\mathrm{TE}$ ~ 1 m and $\Gamma$ ~ 0.1 dB (

Table A5. AGB (ton/ha) uncertainty (ignoring noise) over Congo rainforests assuming an uncertainty of trailing edge width $\mathrm{TE}$ ~ 1 m.

AGB [ton/ha]	50	100	150	200	250	300	350
[20, 30]°	1.53	1.21	0.90	0.55	0.35	0.33	0.30
[40, 50]°	-	1.05	0.93	0.77	0.58	0.37	0.17
[60, 70]°	-	1	0.87	0.71	0.53	0.35	0.19
[80, 90]°	-	-	-	0.62	0.45	0.22	-

,

Table A6. AGB (ton/ha) uncertainty (ignoring noise) over Amazon rainforests assuming an uncertainty of trailing edge width $\mathrm{TE}$ ~ 1 m.

AGB [ton/ha]	50	100	150	200	250	300	350
[20, 30]°	-	1.07	0.75	0.40	0.10	-	-
[40, 50]°	-	1.15	0.85	0.53	0.18	-	-
[60, 70]°	-	-	0.70	0.60	0.20	-	-
[80, 90]°	-	-	0.70	0.42	0.12	-	-

,

Table A7. AGB (ton/ha) uncertainty (ignoring noise) over Congo rainforests assuming an uncertainty of reflectivity $\Gamma$ ~ 0.1 dB.

AGB [ton/ha]	50	100	150	200	250	300	350
[20, 30]°	-	-	-	-	1.71	1.55	-
[40, 50]°	-	-	-	1.93	2.52	2.74	2.31
[60, 70]°	-	-	1.05	2.12	2.6	2.61	1.75
[80, 90]°	-	-	-	-	1.62	2.13	1.62

and

Table A8. AGB (ton/ha) uncertainty (ignoring noise) over Amazon rainforests assuming an uncertainty of reflectivity $\Gamma$ ~ 0.1 dB.

AGB [ton/ha]	50	100	150	200	250	300	350
[20, 30]°	-	-	1.58	1.6	0.95	-	-
[40, 50]°	-	1.72	1.8	1.6	1.15	0.9	-
[60, 70]°	-	1.6	1.78	1.6	1.04	-	-
[80, 90]°	-	-	1.58	1.4	0.95	-	-

).

In this work, the main interest is in the sensitivity, ignoring noise in the measurements when observations deviate from the fits [68]. As such, at a first order, the required GNSS-R sensitivity and the AGB uncertainties are related by the first derivative of the fitting functions [69]. These results suggest that different types of forests have a different response. As such, different models should be applied to retrieve biomass at each type of forest.

Figure 4 and Figure 5 show that $\mathrm{TE}$ increases with increasing AGB up to $\mathrm{TE}$ ~ 800 m over Congo (Figure 4d) and $\mathrm{TE}$ ~ 600 m over Amazon (Figure 4e–h). This experimental evidence is justified because, for increasing vegetation levels (Figure A1a,c), (a) a higher volume scattering term in ${\mathsf{\sigma }}^{\mathrm{incoh},0}$ increases the tail of $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$, and (b) the coherent scattering term ${\mathsf{\sigma }}^{\mathrm{coh},0}$ is gradually more attenuated. On the other hand, the measurable AGB dynamic range increases with decreasing ${\mathsf{\theta }}_{\mathrm{e}}$, from ~140 ton/ha at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[80, 90]° to 300 ton/ha at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[20, 30]° over the Congo target area, while it increases from ~ 140 ton/ha at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[80, 90]° to 170 ton/ha a ${\mathsf{\theta }}_{\mathrm{e}}$ ~[20, 30]° over the Amazon site. The measurable AGB range is a function of ${\mathsf{\theta }}_{\mathrm{e}}$ because of the impact of the different scattering mechanisms at different ${\mathsf{\theta }}_{\mathrm{e}}$. The coherent scattering increases at grazing angles, belonging to a higher measurable AGB dynamic range, as compared to a Nadir looking configuration. The improved measurable dynamic range at lower angles is explained because of the higher coherent reflectivity at this geometry Equation (4). This scattering property belongs to an increment of the SNR dynamic range, which in turn improves the sensitivity of $\mathrm{TE}$ to the attenuating cover. Coherence effects could also appear after scattering over the vegetation cover if the coherent integration time would set to be long enough, e.g., ${{\displaystyle \mathrm{T}}}_c$ ~ 20 ms [27,28,33], so as to filter out the noise and the volume scattering term. However, the volume scattering term is high for lower integration times, e.g., ${{\displaystyle \mathrm{T}}}_c$ ~ 1 ms.

$\Gamma$ gradually decreases down to ~ −25 dB with increasing AGB because of the increasing attenuation of the GNSS signals along the vegetation cover (Figure 5). The measurable AGB dynamic range first increases from ${\mathsf{\theta }}_{\mathrm{e}}$ ~[80, 90]° to ~[60, 70]°, and later it reduces with decreasing angles from ${\mathsf{\theta }}_{\mathrm{e}}$ ~[60, 70]° to ~[20, 30]°. On the other hand, $\Gamma$ is significantly stronger at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[20, 30]° than at higher observation angles. This is a symptom that the coherent scattering term ${\mathsf{\sigma }}^{\mathrm{coh},0}$ is dominant at grazing geometries, even over tropical rainforests [46]. The coherent scattering term ${\mathsf{\sigma }}^{\mathrm{coh},0}$ is mainly associated with surface scattering. As such, the sensitivity of $\Gamma$ to AGB reduces when ${\mathsf{\sigma }}^{\mathrm{coh},0}$ is dominant. Indeed, it is found that the RMSE of the scatter plots increases with increasing $\Gamma$ levels. Potential fluctuations of surface properties have a stronger impact at this geometry (Figure 5). With this respect, the performance of $\mathrm{TE}$ is better than $\Gamma$. The information in $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{peak}}$ is poorly related to biomass when ${\mathsf{\sigma }}^{\mathrm{coh},0}$ is dominant. However, the higher SNR dynamic range is linked to a higher $\mathrm{TE}$ dynamic range, which in turn increases the sensitivity of $\mathrm{TE}$ to AGB, although the sensitivity of $\Gamma$ is poor. These observations suggest that the sensitivity of $\Gamma$ to forests canopies improves when the incoherent scattering term is dominant over the coherent one ${\mathsf{\sigma }}^{\mathrm{coh},0}$ [46]. Indeed, the highest AGB dynamic range appears at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[60, 70]°. The interpretation is two-fold: (a) a larger signal propagation path through the vegetation, as compared to that corresponding to higher ${\mathsf{\theta }}_{\mathrm{e}}$, improves the sensitivity to biomass, and (b) the coherent scattering term ${\mathsf{\sigma }}^{\mathrm{coh},0}$ starts to be dominant for angles lower than ${\mathsf{\theta }}_{\mathrm{e}}$ ~ 60° [46]. A complementary description of the impact of ${\mathsf{\theta }}_{\mathrm{e}}$ can be found in [46].

Overall, it is also found that the dispersion in the scatter plots is very high over the strongest levels of $\Gamma$, where surface reflectivity starts to play an important role in the overall reflectivity. This significant dispersion is related to the fluctuations of certain surface properties, depending on the surface area under consideration. This dispersion increases for grazing angles, which suggests that surface properties such as rivers start to play an important role. As such, $\Gamma$ appears quite high even for high AGB levels.

The sensitivity over Congo rainforests shows a higher angular variability than over the Amazon. This experimental evidence could be linked to different structural properties of the vegetation cover. The AGB distribution over the Amazon is a bimodal one with the two peaks at ~200 and ~250 ton/ha (Figure A1c). The re-radiation pattern could have isotropic properties because of the small separation between both modes. This characteristic could explain the small angular variability of the sensitivity. On the other hand, a bimodal distribution is also found over the Congo (Figure A1a). The two modes are much more differentiated, with peaks at ~300 and ~500 ton/ha. It is hypothesized that this characteristic could belong to a much more complex structure of the canopy layer, and a higher angular variability.

(2) Coniferous, dry, and moist tropical forests (Figure 7 and Figure 8): The functional correlations of AGB with $\mathrm{TE}$ (Figure 7) and $\Gamma$ (Figure 8) can be fitted by empirically derived 1st order polynomial functions. Table 1 summarizes the statistics of these fits at different angular ranges ${\mathsf{\theta }}_{\mathrm{e}}$ ~[80, 90]°, ~[60, 70]°, ~[40, 50]°, and ~[20, 30]°, over coniferous (Figure 7a–d and Figure 8a–d), dry (Figure 7e–h and Figure 8e-h), and moist (Figure 7i–l and Figure 8i–l) forests. A distinguished characteristic of these forests is the low level of AGB, combined with moderate to high levels of CH (Figure A1e–j). These forests are presented following a decreasing CH level from coniferous (Figure A1f) to moist forests (Figure A1j).

The assessment using $\mathrm{TE}$ and $\Gamma$ shows that the maximum correlation coefficients and slopes of the linear regressions (Table 1) appear gradually at decreasing ${\mathsf{\theta }}_{\mathrm{e}}$ with decreasing CH levels: ${\mathsf{\theta }}_{\mathrm{e}}$ ~[80, 90]° over coniferous, ${\mathsf{\theta }}_{\mathrm{e}}$ ~[60, 70]° over dry, and ${\mathsf{\theta }}_{\mathrm{e}}$ ~[40, 50]° over moist forests (Table 1). The angular differences are small. The RMSE is high even for low AGB levels. A reasonable justification is provided, despite the fact that these aspects complicate the interpretation of the results.

Larger signal propagation paths through the vegetation cover are required to increase the sensitivity with decreasing CH levels. On the other hand, the surface contribution partially masks the vegetation effects at grazing geometries. This latter aspect can be observed more clearly at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[20, 30]°, where the correlation coefficients with AGB are quite low (Table 1). This behavior is different to that over rainforests, because low AGB levels belong to a stronger contribution of surface scattering in $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{peak}}$. Finally, it is worth pointing out that the correlation of AGB with $\Gamma$ is much higher than with $\mathrm{TE}$. This observation is also justified because of the dominant effect of surface scattering (Figure 2). On the other hand, the influence of volume scattering is low because of the low AGB levels. This point is translated into lower correlation coefficients with $\mathrm{TE}$, because the sensitivity of $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{threshold}}$ to AGB is almost negligible in this scenario. Over forests with low AGB ~ < 40 ton/ha, the soil surface contribution to the case of GNSS-R is dominant. This explains the lower spreading of the waveforms, as compared to rainforests (Figure 2). In the latter scenario, the volume scattering is very strong.

4.3. Evaluation with SMAP-derived VOD, PI, and SMC

At present, most of biomass monitoring studies are based on the use of visible and infrared indices. However, these types of measurements are only sensitive to the upper canopy layers and they additionally suffer from weather conditions and clouds. On the other hand, this study is based on measurable biomass parameters (AGB, CH). These considerations motivate to complement this work with microwave radiometry-derived parameters, such as $\mathrm{VOD}$ and $\mathrm{PI}$. There is an increasing interest on the use of both parameters to estimate some vegetation properties, such as VWC and AGB [58,59,60,61,62]. Additionally, SMC is considered over coniferous, dry, and moist forests, because of the potential impact of this parameter in the performance of CyGNSS over forests with low AGB levels (Figure A1). Over tropical rainforests, the impact of SMC is assumed to be negligible [24]. As a matter of fact, SMC data are not fully reliable over very dense forest canopies.

The functional relationship of AGB with both $\mathrm{TE}$ and $\Gamma$ is further evaluated at the specific angular ranges where the optimum performance is found in terms of correlation coefficients and AGB dynamic ranges (Figure 4, Figure 5, Figure 7 and Figure 8). This is a case by case study that depends on the characteristics of each type of forest: Congo rain ($\mathrm{TE}$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[20, 30]° and $\Gamma$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[60, 70]°), Amazon rain ($\mathrm{TE}$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[20, 30]° and $\Gamma$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[60, 70]°), coniferous ($\mathrm{TE}$ and $\Gamma$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[80, 90]°), dry ($\mathrm{TE}$ and $\Gamma$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[60, 70]°), and moist ($\mathrm{TE}$ and $\Gamma$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[40, 50]°) forests. An interpretation of the results is provided for each target area:

(1) Tropical rainforests (Figure 9 and Figure 10): $\mathrm{VOD}$ increases with increasing AGB, so as to provide an estimation of the attenuation by the vegetation cover. The sensitivity of $\mathrm{VOD}$ to AGB is high up to AGB ~350 ton/ha ($\mathrm{TE}$ ~ 800 m and $\Gamma$ ~ −25 dB) over Congo (Figure 9a and Figure 10a), and up to AGB ~250 ton/ha ($\mathrm{TE}$ ~ 600 m & $\Gamma$ ~ −25 dB) over Amazon (Figure 9b and Figure 10b). This sensitivity is qualitatively similar to that of CH (Figure 4d,h and Figure 5b,f). Recently, a good correlation of Soil Moisture Ocean Salinity (SMOS)-derived $\mathrm{VOD}$ with CH has also been reported [48,70]. On the other hand, SMAP-derived $\mathrm{VOD}$ is linearly related to VWC trough a parameter that accounts for the structural effects of the vegetation cover [55]. As such, it seems reasonable to find this dependence with CH, although VWC is indirectly inferred from NDVI.

Microwave radiometric brightness temperatures, ${\mathrm{T}}_{\mathrm{BH}}$ and ${\mathrm{T}}_{\mathrm{BV}}$, depend on frequency, polarization, and the geometrical and dielectric properties of the scattering media. Different polarization characteristics of bare soil and vegetation suggest the potential use of $\mathrm{PI}$ for biomass studies [59,60,61,62]. $\mathrm{PI}$ should gradually decrease with increasing AGB levels because (a) V-polarized emissivity from soil is attenuated by the upwelling vegetation cover, and (b) the radiation component due to forest canopies is roughly unpolarized. Figure 9 and Figure 10 confirm these theoretical expectations. The polarization degree decreases with increasing AGB, without an apparent saturation. Additionally, it is found that $\mathrm{PI}$ takes different values for the Congo and Amazon for a given AGB level. This empirical observation suggests that $\mathrm{PI}$ slightly depends on CH additionally to AGB, at least over tropical rainforests. A similar dependence with CH is also found for $\mathrm{VOD}$.

(2) Coniferous, dry, and moist tropical forests (Figure 11 and Figure 12): A significant sensitivity of $\mathrm{VOD}$ is also found for low levels of AGB, over coniferous (Figure 11a and Figure 12a), dry (Figure 11b and Figure 12b), and moist (Figure 11c and Figure 12c) forests. The sensitivity improves with increasing CH levels, from moist to coniferous forests, despite the roughly similar AGB levels (Figure A1e,g,i). On the other hand, $\mathrm{PI}$ also performs well over this type of CH-driven forests (Figure 11d–f and Figure 12d–f). The $\mathrm{PI}$ dynamic range over dry forests ~[0.05 0.075] is slightly lower than over coniferous forests ~[0.05 0.1] because of the lower CH dynamic range (Figure A). At the same time, it is found that $\mathrm{PI}$ is higher over moist forests ~[0.06 0.085] than over dry forests, because of the slightly lower CH levels (Figure A).

The results are independent of SMC over coniferous and dry forests. Over coniferous forests, increasing levels of SMC up to ~0.2 m³/m³ are found for gradually increasing $\mathrm{TE}$ (Figure 11g) and decreasing $\Gamma$ (Figure 12g). This observation means that CH (Figure A1f) has a higher impact than SMC in both observables, $\mathrm{TE}$ and $\Gamma$. An increment in SMC would reduce $\mathrm{TE}$ and increase $\Gamma$, in the case that SMC would be dominant. Over dry forests, SMC ~0.18 m³/m³ is quite homogenous along both $\mathrm{TE}$ (Figure 11h) and $\Gamma$ (Figure 12h) dynamic ranges. As such, changes on $\mathrm{VOD}$ and $\mathrm{PI}$ are assumed to be mainly related to the effects of the vegetation cover. Finally, in the case of moist forests, SMC influences $\mathrm{TE}$ (Figure 11i) and $\Gamma$ (Figure 12i). Increasing levels of SMC up to ~0.25 m³/m³ are found for decreasing $\mathrm{TE}$ (Figure 11i) and increasing $\Gamma$ (Figure 12i). Even in this scenario, $\mathrm{VOD}$ (Figure 11c and Figure 12c) and $\mathrm{PI}$ (Figure 11f and Figure 12f) are found to have a certain degree of sensitivity to canopy forests.

4.4. First Evidences of Biomass Retrieval Using an Empirical Approach

Finally, the ability to retrieve biomass using an empirical approach is evaluated (Figure 13). To do so, the empirically derived polynomial fitting functions corresponding to the relationship between $\mathrm{TE}$ and AGB at low elevation angles ${\mathsf{\theta }}_{\mathrm{e}}$ ~[20, 30]° are used to provide a preliminary AGB estimation at the two selected test sites over Congo and Amazon tropical rainforests. The percentage of data used as training and validation subsets ranged from 30% and 70%, respectively. On the other hand, the AGB map provided by Avitabile et al. is used as the AGB reference. Small absolute errors in AGB estimation are found except in areas with AGB > ~ 350 ton/ha and in areas with little biomass and the presence of inland water bodies (Figure 3). These promising preliminary results should be further investigated in future works also using ground truth data for validation.

5. Conclusions

Earth-reflected GNSS signals collected in a bistatic radar configuration by CyGNSS, show a significant sensitivity to biomass over rain, coniferous, dry, and moist tropical forests. Future studies with the U.K. TDS-1 mission would help to also cover boreal forests. The sensitivity has been evaluated using several parameters derived from DDMs, as a function of ${\mathsf{\theta }}_{\mathrm{e}}$. The optimum results in terms of correlation coefficients and AGB dynamic ranges were found at different angles: Congo rain ($\mathrm{TE}$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[20, 30]° and $\Gamma$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[60, 70]°), Amazon rain ($\mathrm{TE}$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[20, 30]° and $\Gamma$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[60, 70]°), coniferous ($\mathrm{TE}$ and $\Gamma$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[80, 90]°), dry ($\mathrm{TE}$ and $\Gamma$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[60, 70]°), and moist ($\mathrm{TE}$ and $\Gamma$ at ${\mathsf{\theta }}_{\mathrm{e}}$ ~[40, 50]°). The SNR reduces at low elevation angles for both direct and reflected signals. However, we use $\Gamma$. This observable is more related with geophysical parameters than the SNR. In particular, $\Gamma$ increases at lower angles because of the lower effective surface roughness at this geometry.

Over rainforests, the optimum geometry for the use of $\mathrm{TE}$ appears at grazing angles because: (a) a higher volume scattering term increases the tail of $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$, and (b) a higher coherent scattering term $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{coh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ belongs to an improved sensitivity of $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{peak}}$ to the attenuating cover. On the other hand, the optimum geometry for $\Gamma$ is found at moderate ${\mathsf{\theta }}_{\mathrm{e}}$, when the incoherent scattering term $\langle {{\displaystyle \left|{\mathrm{Y}}_{\mathrm{r},\mathrm{incoh}}(\mathsf{\tau },\mathrm{f})\right|}}^2\rangle$ is still dominant and the signal propagation path along the vegetation cover is higher than at a Nadir-looking configuration. Then, the first derivative of the corresponding polynomial fitting functions of the relationship between AGB with $\mathrm{TE}$ and $\Gamma$ was computed. It was found that there is a certain sensitivity to AGB up to ~350 and 250 ton/ha, respectively, over Congo and Amazon target areas, without an apparent saturation. The first maps of absolute error in AGB retrieval from space using GNSS-R were included, showing a small error. These findings suggest the possibility of accurate AGB retrieval using GNSS-R on-board small satellites such as, e.g., CyGNSS.

Over forests with low AGB ~ < 40 ton/ha, the optimum geometry for the use of both parameters $\mathrm{TE}$ and $\Gamma$ appears at decreasing ${\mathsf{\theta }}_{\mathrm{e}}$ with gradually decreasing CH levels. This point is explained because an increasing signal propagation path through the vegetation is required to improve the sensitivity with decreasing CH levels. The correlation coefficients of AGB with $\Gamma$ are higher than those with $\mathrm{TE}$. The sensitivity of $\mathrm{TE}$ to AGB is low because of the negligible volume scattering term over these types of forests. In this scenario, the spreading of ${\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})$ is low, which reduces the sensitivity of $\mathrm{W}{\mathrm{F}}_{\mathrm{r},\mathrm{threshold}}$ to biomass.

Finally, it is worth pointing out that the application of a moving averaging filter minimizes potential errors in the estimation of $\Gamma$ due to inaccuracies in the signal calibration. The uncertainty in establishing the relationship between CyGNSS measurables and forest parameters is due to (a) potential errors in CH data derived from GLAS and (b) noise at the pixel level, such as speckle in ${\mathrm{Y}}_{\mathrm{r}}(\mathsf{\tau },\mathrm{f})$. Overall, RMSE in the scatter plots increases with increasing biomass levels, while the correlation coefficients increase.