An Improved Method for Pan-Tropical Above-Ground Biomass and Canopy Height Retrieval Using CYGNSS

Abstract

An improved method for retrieving Above-ground Biomass (AGB) and Canopy Height (CH) based on an observable from Cyclone Global Navigation Satellite System (CYGNSS), soil moisture from Soil Moisture Active Passive (SMAP) and location is proposed. The observable derived from CYGNSS is more sensitive to vegetation. The CYGNSS observable, soil moisture and the location are used as the input features of an Artificial Neural Network (ANN) to retrieve AGB and CH. The sensitivity analysis of the CYGNSS observable to target parameters shows that the proposed observable is more sensitive to AGB/CH than the conventional observable. The AGB/CH retrievals of the improved method show that it has better performance than that of the traditional method, especially in the areas with AGB in the range of 0 to100 Mg/ha and CH in the range of 0 to10 m. For AGB retrievals, the root mean square error (RMSE) and correlation coefficient are 64.84 Mg/ha and 0.80 in the range of 0 to 550 Mg/ha. Compared with the traditional method, the RMSE is decreased by 11.63%, while the correlation coefficient is increased by 5.26%. For CH retrievals, the RMSE and correlation coefficient are 5.97 m and 0.83 in the range of 0 to 45 m. The RMSE is decreased by 12.59%, while the correlation coefficient is increased by 5.06%. The analysis of the improved method in different areas shows that the performance of the improved method over the area with high vegetation is better than the area with low vegetation. The results obtained here further strengthens the capability of GNSS-R for global AGB/CH retrievals as well as different land cover areas.

1. Introduction

Mapping Above-ground Biomass (AGB) is crucial for quantifying the global carbon cycle, which contributes to our understanding of the global carbon balance and climate change [1,2,3,4]. Due to tropical deforestation, 5–19% of global anthropogenic CO₂ is emitted into the atmosphere [5,6]. Estimating these emissions needs accurate spatial and temporal estimation of AGB in global forests, particularly for pan-tropical areas with high AGB density [7,8].

The AGB retrievals can be routinely estimated via empirical allometric models relating field tree heights and diameter at breast heights (DBHs) [9]. As important methods for validation purposes, traditional field tree heights and DBH measurement methods have a higher accuracy in the local area. However, they require a lot of manpower and material resources, which are very difficult, costly, and time-consuming [10,11,12]. Moreover, the traditional methods are suitable for small areas with little vegetation change, as their accuracy of assessing larger region AGB distribution may be reduced [12,13,14,15]. Therefore, the traditional methods are not suitable for rapid assessment of global AGB distribution. Remote sensing instruments can be used for rapid assessment of global AGB distribution by transmitting signals and receiving reflected signals from earth’s surface [16]. The acquired remote sensing data can be analyzed to estimate forest structure, such as Canopy Height (CH), which is further correlated with AGB [17]. The first satellite mission for AGB retrievals is the Synthetic Aperture Radar (PALSAR)/Advanced Land Observing Satellite (ALOS) [18]. Due to signal saturation when AGB is larger than 60–100 Mg/ha of PALSAR/ALOS, interferometric Synthetic Aperture Radar (InSAR) and a deployable reflector antenna are used in the BIOMASS mission operated by the European Space Agency [19]. Furthermore, the information of species type and composition can be accurately retrieved by airborne hyperspectral remote sensing instruments [20] and Light Detection and Ranging (LiDAR) can provide accurate CH and volume [21]. After the tree structure information is obtained, the AGB is estimated by empirical methods. However, due to the limitation of hyperspectral resolution, different forest types often have similar spectral characteristics. The platforms of LiDAR used for AGB estimation are mostly based on an airborne platform, which has high associated data volumes and costs [22]. Therefore, they are more suitable for investigating forests over small, localized areas. As for spaceborne platforms, LiDAR measurements from the Ice, Cloud, and land Elevation Satellite (ICESat) were further used to estimate CH and AGB with RMSEs between 4.85 m and 12.66 m for CH and 58.3 Mg/ha for AGB [23]. In addition, the Global Ecosystem Dynamics Investigation (GEDI) also provides global, high-resolution observations of forest vertical structure using LiDAR [24].

Global Navigation Satellite Systems Reflectometry (GNSS-R) exploits GNSS signals reflected from the earth surface in a forward bistatic radar configuration, which has developed as a mature remote sensing technique for earth surface geophysical parameter retrieval [25,26]. GNSS-R operating L band signals have a good ability to penetrate the atmosphere, vegetation and rain, helping to sense, for example, water, soil under vegetation and clouds [27]. Furthermore, thanks to the GNSS signals, the GNSS-R receiver does not require an additional transmitter, which contributes to continuous high spatial and temporal coverages of GNSS-R observations, as well as low-volume, low-cost and low-power of the GNSS-R receiver. With its successful use in oceanographic applications, e.g., ocean wind field, sea surface height, sea roughness, [28,29,30,31], GNSS-R has recently been shown to have capacity for land surface geophysical parameter retrieval, e.g., soil moisture, wetland, flood inundation, and roughness [32,33,34,35,36]. In the application of AGB/CH retrieval, only a few early studies and experiments have been carried out based on ground-based and airborne platforms [37,38,39], which has the limitation of global coverage. The implementation of low earth orbit (LEO) satellite missions for GNSS-R, e.g., the United Kingdom-Disaster Monitoring Constellation (UK-DMC), TechDemoSat-1 (TDS-1), and Cyclone Global Navigation Satellite System (CYGNSS), has demonstrated the potential of spaceborne GNSS-R for global coverage. Recently, the focus has been on AGB/CH retrieval using TDS-1 and CYGNSS data [40,41]. In general, the algorithms of these studies are based on GNSS-R observables, such as signal-to-noise ratio (SNR) and trailing edge (TE) derived from delay doppler maps (DDMs). These data can be further used to estimate the incoherent and coherent components of scattered signals. It is important to determine GNSS-R measurements that are dominated by coherent or incoherent effects for Earth’s land physical parameter retrievals [42,43]. In most cases, scattering is neither coherent nor incoherent, but partially coherent. According to the nature of scattering, the relationship between the GNSS-R observables and the geophysical parameters is established. To estimate one geophysical parameter, the other geophysical parameters are obtained first. Therefore, to estimate AGB/CH, the information of other parameters, e.g., surface roughness or soil moisture, is needed. For AGB/CH retrievals, Santi et al. successfully applied ANN to forest biomass and achieved encouraging results with a RMSE of 76.4 Mg/ha for AGB estimation and a RMSE of 6.5 m for CH estimation [41]. However, the traditional observable method is heavily influenced by other parameters, especially soil moisture. Due to the heavy rainfall in tropical areas, soil moisture will increase and GNSS reflected signals strength will be enhanced. On the other hand, broadleaved evergreen trees are mainly distributed in the tropics, which weakens GNSS reflected signals. Therefore, it is crucial to exclude the effect of soil moisture on GNSS reflected signals for pan-tropical AGB/CH retrieval. On the other hand, soil moisture also affects the growth of vegetation. Thus, it is worth investigating soil moisture for retrieving AGB/CH. To further develop the potential of GNSS-R in the forest biomass, this paper investigates the effects of the soil moisture on AGB/CH retrievals. In addition, the effects of topography on spaceborne GNSS-R signals have been modeled and verified in previous works [44,45]. However, the purpose of this paper was to consider and remove the effects of soil moisture on AGB/CH retrievals. Since SMAP can provide very high global soil moisture accuracy (4%), to some extent, the soil moisture provided by SMAP is of potential value for AGB/CH retrieval using CYGNSS.

The focus of this paper is proposing an improved method for CYGNSS AGB/CH retrievals using soil moisture derived from SMAP. The structure of the remainder of the paper is organized as follows. Section 2 presents the introductions of CYGNSS, Soil Moisture Active Passive (SMAP), and reference AGB/CH datasets. Section 3 describes the proposed AGB/CH retrieval methods. Section 4 shows the results and discussions. Finally, Section 5 gives concluding remarks.

2. Datasets

2.1. CYGNSS

The primary purpose of the CYGNSS NASA mission is to monitor hurricanes by using a revision of the GNSS Receiver Remote Sensing Instrument (SGR-ReSI). With exploiting GNSS signals reflected from the earth surface in a forward bistatic radar configuration and using eight micro-satellites, CYGNSS increases the capacity of global coverage within approximately ±38° latitudes and also reduces the revisit time by several hours. Most recently, in addition to ocean applications, CYGNSS shows the potential for retrieving Earth’s geophysical parameters. Therefore, CYGNSS has the advantage of spatial coverage and high temporal resolution for biomass monitoring.

CYGNSS provides four different processing data levels. The version 2.1 of the CYGNSS Level 1 data in NetCDF format, from Jan to Jun, 2019, is used for further analysis in this study. The data are available at https://podaac.jpl.nasa.gov (accessed on 5 January 2021). The observables derived from Level 1 data are DDMs, latitude, longitude, etc. Some of the observables (shown in Section 3) are used to compute the surface reflectivity. After correction, the surface reflectivity is mapped into a 0.05° × 0.05° grid. Because there is variable vegetation or vegetation with different growth patterns in the pixel, different azimuth angles of CYGNSS will bring different reflectivities. To reduce error, if there are repeated values in each pixel, the average values are used for further analysis. In addition, the other observables (shown in Section 3) also are mapped into a 0.05° × 0.05° grid and averaged. Finally, the mean observables, as well as the corresponding latitude and longitude, serve as the algorithm input features.

2.2. SMAP

The SMAP mission was launched in January 2015 by NASA. By taking advantage of both active (radar) and passive (radiometer) microwaves, SMAP can provide the accuracy, resolution, and global soil moisture data for land within ±45° latitudes, every 2–3 days [46]. The data is available at https://nsidc.org/data/smap/smap-data.html (accessed on 7 January 2021). This dataset is posted to the 9 km Equal-Area Scalable Earth Grid, Version 2.0 (EASE-Grid 2.0) in a global cylindrical projection [47]. To match the CYGNSS observables, the soil moisture from Jan to Jun 2019 is mapped into a 0.05° × 0.05° grid and averaged as depicted in Figure 1. However, the blank pixels can be clearly seen in Figure 1. The reason is that SMAP provides the soil moisture product with a resolution of 9 Km, which is not sufficient for the 0.05° × 0.05° grid. Therefore, in order to make up for this, the blank pixels are interpolated using the nearest neighbor value interpolation. The results are depicted in Figure 2 and are used for further analysis.

2.3. AGB Map

The reference pan-tropical biomass map of this study is derived from LUCID (Land use, Carbon & Emission Data, http://lucid.wur.nl/datasets (accessed on 15 January 2021). This map is based on two datasets of AGB [48,49]. The first dataset is obtained by combining remote sensing data bases and ground data [48]. The second AGB map is based on multi-sensor satellite data [49]. The independent reference dataset of field measurements and locally calibrated high-resolution biomass maps can provide the estimations of AGB for the tropics (40°N–40°S) with nearly unbiased estimations for mean bias 5 Mg dry mass ha⁻¹ and 1 km resolution [50]. To facilitate comparison with CYGNSS data, the resolution of the AGB map is reduced to 0.05° × 0.05° as displayed in Figure 3.

2.4. CH Map

The CH map derived from ICESat/GLAS and environmental modeling [51] is used for evaluating the proposed algorithm’s sensitivity to the pan-tropical biomass. The map is based on data from the Geoscience Laser Altimeter System (GLAS) aboard ICESat (Ice, Cloud, and land Elevation Satellite) in 2005. It can provide the CH for global scale with a resolution of 10 km, and can be downloaded at https://landscape.jpl.nasa.gov/ (accessed on 15 January 2021). Like the AGB map, the resolution of the CH map is reduced to 0.05° × 0.05° as displayed in Figure 4.

3. Methods

3.1. Observables Derived from CYGNSS

The total GNSS-Reflectometry (GNSS-R) receiver measures the GNSS reflected signals scattered from the Earth’s surface. The results are recorded in the so-called DDMs ($⟨{\left|Y\left(\tau ,f\right)\right|}^2⟩$). which can be modeled as [52,53]:

(1)$⟨{\left|Y\left(\tau ,f\right)\right|}^2⟩=\frac{{\lambda }^2{T}_i^2}{{\left(4\pi \right)}^3}{\displaystyle {\iint }_A\frac{P_t{G}_t{G}_r{\Lambda }^2\left(\tau \right)S^2\left(f\right)}{R_t^2{R}_r^2}\left({\sigma }_0^{coh}+{\sigma }_0^{incoh}\right)dA}$

where $⟨{\left|Y\left(\tau ,f\right)\right|}^2⟩$ corresponds to the GNSS power scattered by the ocean surface as a function of the time delay $\tau$ and the frequency offset $f$. $P_t{G}_t$ represents the GNSS Effective Isotropic Radiated Power (EIRP) and antenna gain of the transmitter. $\lambda$ is the carrier wavelength. $T_i$ is the coherent integration time. ${\Lambda }^2$ and $S^2$ are the components of the Woodward Ambiguity Function (WAF) in delay and delay Doppler frequency, respectively. $R_t$ and $R_r$ are the transmitter-to-surface and surface-to-receiver ranges, respectively. $dA$ denotes the surface element of the scattering area $A$. ${\sigma }_0^{coh}$ and ${\sigma }_0^{incoh}$ symbolizes BRCS, which is related to the glistening zone. When the reflected surface is relatively smooth or has little vegetation density, the coherent component (${\sigma }_0^{coh}$) appears in the scattering power, while the incoherent component (${\sigma }_0^{incoh}$) increases with the increases of surface roughness and vegetation density. Previous works have proposed several algorithms to account for the contributions of the coherent component and incoherent component to the DDMs [53,54,55,56]. These algorithms are based on the assumption that the reflected power decreases quickly away from the specular reflection point. Therefore, the reflected power based on the radar equation is rewritten as:

(2)$⟨{\left|Y_{coh}\left(\tau ,f\right)\right|}^2⟩=\frac{P_t{G}_t{G}_r{\lambda }^2{\Lambda }^2\left(\tau \right)S^2\left(f\right)}{{\left(4\pi \right)}^2{\left(R_t+R_r\right)}^2}{\sigma }_0^{coh}$

where ${\sigma }_0^{coh}$ is equivalent to the traditional Reflectivity (${\Gamma }_{surface}$). Therefore, ${\Gamma }_{surface}$ can be derived from DDM by:

(3)${\Gamma }_{surface}\approx \frac{P_r{\left(R_t+R_r\right)}^2}{P_t{G}_t{G}_r}{\left(\frac{4\pi }{\lambda }\right)}^2$

where $P_r$ is the peak power of DDM. $P_r$ is approximate to SNR in the previous work. Because ${\Gamma }_{surface}$ is affected by surface reflectivity Fresnel reflectivity (${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)$), vegetation attenuation, and surface roughness, (3) has been successfully applied to retrieve soil moisture, surface roughness, vegetation and other surface physical parameters. ${\Gamma }_{surface}$ is given by the following:

(4)${\Gamma }_{surface}={\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)\cdot \gamma \cdot \upsilon$

where ${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)$ represents the surface Fresnel reflectivity in polarization LR (right- hand circular incidence and left-hand circular scattering for CYGNSS). ${\epsilon }_s$ is surface relative complex permittivity, which is related to soil moisture. ${\theta }_i$ is incident angle. $\gamma$ is the vegetation loss term. $\upsilon$ is the attenuation factor due to the surface roughness. ${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)$ is given by:

(5)${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)={\left|\frac{1}{2}\left({\Re }_V-{\Re }_H\right)\right|}^2$

with:

(6)${\Re }_V=\frac{{\epsilon }_s\cos {\theta }_i-\sqrt{{\epsilon }_s-{\sin }^2{\theta }_i}}{{\epsilon }_s\cos {\theta }_i+\sqrt{{\epsilon }_s-{\sin }^2{\theta }_i}}$

(7)${\Re }_H=\frac{\cos {\theta }_i-\sqrt{{\epsilon }_s-{\sin }^2{\theta }_i}}{\cos {\theta }_i+\sqrt{{\epsilon }_s-{\sin }^2{\theta }_i}}$

where ${\Re }_V$ and ${\Re }_H$ are the vertical and horizontal reflection coefficients, respectively. Therefore, to estimate the vegetation loss term ($\gamma$) related to AGB/CH, the ${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)$ and $\upsilon$ have to be calculated first. The ${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)$ is computed by the Dobson model [57]. Under the Dobson model, the information of soil moisture and incident angle is needed. Here, soil moisture is derived from SMAP, while incident angle is derived from CYGNSS level 1 data. Figure 5 shows the dependence of ${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)$ on soil moisture and incident angle.

Figure 5 demonstrates that ${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)$ increases with the increase of soil moisture, while it decreases with the increase of incident angle. In addition, when the incident angle is less than 40°, the change in ${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)$ is very small. Therefore, the CYGNSS observables are used for further analysis with the incident angles less than 40°.

Due to the resolution of this paper of 0.05° × 0.05°, which is approximately 5 km², the surface roughness varies little in this range. Therefore, $\upsilon$ is assumed constant. Thus, ${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)$ is given by:

(8)$\gamma ={\sigma }_0^{coh}/{\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)$

Since the vegetation loss term ($\gamma$) decreases with the increase of AGB/CH, $1/\gamma$ is more suitable for AGB/CH retrievals. Therefore, the correction Reflectivity (${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)/{\sigma }_0^{coh}$) is used for further analysis and ANN input feature. ${\sigma }_0^{coh}$ is the peak point of BRCS from CYGNSS Level 1 data.

The traditional Reflectivity (${\Gamma }_{surface}$) and the correction Reflectivity (${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)/{\sigma }_0^{coh}$) are depicted in Figure 6 and Figure 7, respectively. For the convenience of comparison, both of the observables are normalized. It is noted that, in the regions with lower AGB and CH as shown in Figure 3 and Figure 4, the change of the correction Reflectivity changes less than that of the traditional Reflectivity, especially in Argentina. The reason may be that the traditional Reflectivity includes changes in other physical parameters. The sensitivity of the traditional Reflectivity and correction Reflectivity to the target parameters are summarized in Table 1. The results show that the correlation coefficients are −0.48 and 0.58 compared with AGB, while the correlation coefficients are −0.54 and 0.60 compared with CH for traditional Reflectivity and correction Reflectivity, respectively. Therefore, compared with traditional Reflectivity, correction Reflectivity is more sensitive to AGB/CH.

3.2. Artificial Neural Network (ANN)

In theory, the relationship between the GNSS-R observables and surface geophysical parameters is complex. because Artificial Neural Networks can learn the complex relationship between GNSS-R observables and surface geophysical parameters [54], ANN is applied in this study. Some studies have successfully used ANN for GNSS-R applications [41,58,59,60,61]. In the application of forest biomass, the traditional method uses the traditional Reflectivity and location as input features. Since the AGB is different in different places, the location is also an important information parameter. Here, the correction Reflectivity and location are the input features of the improved method. In addition, by comparing Figure 2 with Figure 3 and Figure 4, both AGB and CH are correlated with soil moisture. Therefore, soil moisture is the last input feature of the improved method.

Before the ANN learning processing, the input features and the corresponding true ABG/CH are first randomly divided into two datasets for training and testing. After testing, it is found that too low test samples have poor accuracy and high efficiency, while high test samples have high accuracy but poor timeliness. When it is increased to 5%, the accuracy is not significantly improved. After comprehensive consideration, 5% is selected as the training sample. To balance the accuracy of the model and efficiency as well as resolution, the data sets are randomly divided into two datasets that contain 5% and 95% for training and testing; the training samples are randomly selected globally at latitudes 40°S to 40°N. To verify the integrity and independence of the algorithm, the testing datasets are used as true values for comparison without participating in training. Figure 8 shows the structure and the difference input features between the traditional method and the improved method of ANN.

In Figure 8, Lat and Lon are latitude and longitude. The ANN structure contains an input layer, hidden layers, and an output layer. After testing, the RMSE and the correlation coefficient are optimal when the number of hidden layers is two. The training datasets are input into input layer. Each layer contains numerous neurons, which are interconnected to each other neuron of the closest layer. Through these interconnections, the weights-array of each layer can be trained by the forward and backward mechanisms. After testing the performance of the algorithm from one neuron to multiple neurons in each hidden layer, it is found that the best performance is six neurons in each hidden layer. Therefore, the ANN structure here contains two hidden layers and six neurons in each layer. Assuming that the training dataset samples of this study are m, and denote as $T=\left\{\left(x_1^1,\dots x_n^1,y^1\right),\dots ,\left(x_1^m,\dots x_n^m,y^m\right)\right\}$, where $x_1^i,\dots x_n^i,i=1,\dots ,m$ represent different input features and $y^i$ denotes the corresponding AGB/CH. If the ith layer has N_i number of neurons, the corresponding weights-bias-array ($Wx+b$) at the ith layer is (N_i+1 * N_i). The activation of N_i+1 neurons at the i+1th layer can be realized through the product of the nonlinear function $\sigma$ and a linear combination of the N_i neurons at the ith layer. The equation of this activation is expressed as follows:

(9)$y=\sigma \left(Wx+b\right)$

where w and b are training parameters; $\sigma$ is the nonlinear function, e.g., the sigmoid function, hyperbolic tangent function, and rectified linear units (RELUs), which are given as:

(10)$\sigma (a)=\left\{\begin{array}{l}\frac{1}{1+e^{-a}}, \mathrm{sigmoid}\\ \frac{e^a-e^{-a}}{e^a+e^{-a}}, \tan \mathrm{gent}\\ \max \left(a,0\right), \mathrm{RELU}\end{array}\right.$

In this study, the tangent is used in the three layers. These functions establish a nonlinear relationship between the weights-bias-array of the connected neurons and transmit a new weights-biases-array to the next neuron. By repeating this process at each neuron until the initial AGB/CH is estimated. This process is defined as forward propagation. The information of the error by updating the training weights and bias of each neuron is also transmitted back to the upper connected neurons. By using a stochastic gradient descent algorithm, the loss function is defined. When increasing the number of iterations until the loss function reaches a threshold minimum value or the maximum iterations are reached, the final ANN training model is obtained. Here, the mean square error function serves as the loss function. Then, the test datasets are used to input the ANN training model to produce the final AGB/CH estimations.

3.3. Strategy of the Validation

Based on the proposed method, the RMSE, the relative RMSE (rRMSE), Bias ratio and correlation coefficient are used as validation parameters. The global AGB/CH results are first shown. Then, the results over some subsets based on different land covers are shown in the discussions. The RMSE is given as:

(11)$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^n{\left(X_{true,i}-X_{modle,i}\right)}^2}}{n}}$

where n is the total samples. $X_{true,i}$ is the corresponding truth value. $X_{modle,i}$ is the retrieval value. rRMSE is given by:

(12)$\mathrm{rRMSE}=\mathrm{RMSE}/\mathrm{mean}$

where mean is the average value of the selected area.

The Bias ratio is equal to ($\left|\mathrm{Bias}\right|/\mathrm{ture}$), where true is the corresponding truth value.

4. Results and Discussions

4.1. AGB Retrievals

To evaluate the performance of the improved method, the RMSEs and the correlation coefficients of AGB retrievals based on the traditional method and the improved method are summarized in Table 2. The AGB retrievals scatter plots based on the two methods are depicted in Figure 9. In order to more intuitively see the retrieval effects and details, Figure 10 provides the pan-tropical AGB map derived from CYGNSS, while Figure 11 presents the absolute errors based on the two methods.

Table 2 shows that the RMSEs of the traditional method and the improved method are 73.38 Mg/ha and 64.84 Mg/ha, while the correlation coefficients are 0.76 and 0.80, respectively. Compared with the traditional method, the RMSE of the improved method is decreased by 11.63%, while the correlation coefficient is increased by 5.26%, respectively. The AGB retrievals scatter plots depicted in Figure 9 show that the AGB retrievals of the traditional method have a saturation at AGB values of about 400 Mg/ha, while the values are 450 Mg/ha for the improved method. Note that, since the training set and the test set are set to 5% and 95% of the total sample, the resolution of the retrieval maps is lower than that of Figure 3 and Figure 4. Although the resolution is reduced, the distributions in Figure 10 are roughly the same as that in Figure 3. Although the absolute errors of the improved method are smaller than that of the traditional method, as shown in Figure 11, there are still some larger values in the improved method. This may due to the small training sets over 450 Mg/ha and the time difference between the reference map (2016) and the used datasets (2019). These areas with highest AGB are mainly located in the north-eastern part of the Amazon basin, the central and western part of the Congo basin, the north part of Indonesia, the western part of Malaysia, Borneo, Papua New Guinea, and the northern part of Myanmar. From the mean SMAP soil moisture from Jan to Jun, 2019, as shown in Figure 2, these areas with larger absolute errors are mainly located in the areas with large soil moisture. Besides the high AGB areas, there are high absolute errors in the top sub-figure of Figure 11, which are located in low AGB areas, such as the eastern part of Peru and the northern part of Chile. Note that, compared with Figure 3 and Figure 4, although AGB and CH vary greatly in the western part of Peru and the northern part of Chile, Figure 6 shows that there is little change in both the CYGNSS correction Reflectivity and traditional Reflectivity. The reason may be the influence of the Andes Mountains. Figure 12 presents the three Landsat maps located in the Andes Mountains, which are located in the red box in the top subfigure of Figure 11. Figure 11 shows the variability of biophysical attributes from vegetation to bare soil. Therefore, a biophysical attributes mutation may lead to the stability of the CYGNSS observables.

However, the analysis of the absolute errors alone cannot reflect the overall performance of the improved method, because the absolute errors are relatively larger where the true values are larger. Therefore, an observable named Bias Ratio is proposed for facilitating the analysis as shown in Figure 13. It is notable that, compared with the traditional method, the improved method shows a better performance in low AGB areas in the range of 0 to 100 Mg/ha. From global land cover maps, as shown in Figure 14 [62] and the soil moisture of Figure 2, the traditional method shows a poor performance from the areas covered by broadleaved evergreen trees compared to other low vegetation areas as well as from high soil moisture areas to low moisture areas. The improved method significantly ameliorates this insufficiency. To see the land cover attributes of the selected subsets more intuitively, the locations of the selected subsets (marked with red rectangles) are also depicted in Figure 14 (Congo Forest (Latitude = [−4, 4]°, Longitude = [9, 28]°), Amazon Forest (Latitude = [−10, 5]°, Longitude = [−75, −54]°), Kalimantan (Latitude = [−4.5, 7.3]°, Longitude = [108, 120]°), Brazil (Latitude = [−20, −13]°, Longitude = [−51, −43]°), Tanzania (Latitude = [−12, −0.9]°, Longitude = [29, 40]°), and Australia (Latitude = [−39, −11]°, Longitude = [112, 155]°)).

4.2. CH Retrievals

To further assess the sensitivity of the proposed algorithms to vegetation, this section presents the results of CH estimations based on the traditional method and the improved method. Likewise, the RMSEs and the correlation coefficients of CH retrievals based on the traditional method and the improved method are summarized in Table 3. The CH retrievals scatter plots based on the two methods are depicted in Figure 15. Figure 16 provides the pan-tropical CH map derived from CYGNSS, while Figure 17 and Figure 18 present the absolute errors and Bias Ratio, respectively, based on the two methods.

Table 2 shows that the RMSEs of the traditional method and the improved method, are 6.83 m and 5.97 m, while the correlation coefficients are 0.79 and 0.83, respectively. Compared with the traditional method, the RMSE of the improved method is decreased by 12.59%, while the correlation coefficient is increased by 5.06%. Similar to AGB retrievals, the performance of the improved method is better than that of the traditional method as shown in Figure 15, Figure 16, Figure 17 and Figure 18, especially in the areas with CH in the range of 0 to 10 m. In addition, compared with the AGB retrievals, CH retrievals points are closer to the 1:1 line and show strong linear correlation. Hence, it could indicate that the performance of both methods to estimate CH is better than estimating AGB.

4.3. Discussions

The results above show a good performance of the improved method on global AGB/CH retrievals. To further investigate the performance of the improved method on some subsets and particularities of some environments, Table 4 and Table 5 present the AGB/CH retrievals performance over Congo Forest, Amazon Forest, Kalimantan, Brazil, Tanzania, and Australia. Since AGB and CH are different in different areas, rRMSEs (as shown in Table 4 and Table 5) are used to assess the retrieval effects. In general, the rRMSEs are smaller in high vegetation areas (Congo, Amazon, Kalimantan), while the rRMSEs are larger in low vegetation areas (Brazil, Tanzania). The reason for the large rRMSEs in low vegetation areas (Brazil, Tanzania) is the small training sample size and small value. Moreover, GNSS reflected signal attenuation is more serious in high vegetation areas, while the GNSS reflected signal is mainly affected by surface parameters (e.g., soil moisture, surface roughness, etc.) in low vegetation areas. Therefore, the performance of the improved method in areas with high vegetation is better than areas with low vegetation. Since Australia has forests, grasslands and deserts, the performance of the improved method over Australia is at an intermediate level.

It is worth noting that Carreno-Luengo et al. [40] also analyzed the AGB retrievals by CYGNSS in different areas, and their results show that the RMSEs reached 30–40 Mg/ha for the Congo and Amazon forests. However, one should be aware that such smaller RMSEs are a bit optimistic because a moving averaging filter has been introduced in their statistics.

5. Conclusions

An improved method for retrieving AGB/CH from CYGNSS based on an observable that is more sensitive to vegetation, soil moisture and location is proposed. Here, the soil moisture data is derived from SMAP. The reference AGB map is derived from LUCID and the CH map is derived from NASA. For better comparisons, all datasets are mapped into a 0.05° × 0.05° grid within 40°N–40°S. To balance the resolution and the stability of the algorithms, after testing, the training set and the testing sets are set to 5% and 95% of the overall samples, respectively.

The sensitivity analysis of correction Reflectivity to target parameters shows that the correction Reflectivity is more sensitive to AGB/CH than the traditional Reflectivity. The CYGNSS AGB/CH estimations based on the observables are compared with the reference AGB/CH data. The results show that, compared with the traditional method, the improved method has a good correlation and low RMSE with reference AGB/CH values. It demonstrates that the use of soil moisture for AGB/CH retrievals is essential. The analysis of the improved method in different areas shows that the performance of the improved method over areas with abundant vegetation is better than areas with low vegetation. Therefore, to obtain higher AGB/CH retrievals accuracy, different areas need different models. The work presented here further strengthens the capability of GNSS-R for global AGB/CH retrievals as well as different land cover areas.

Author Contributions

Conceptualization, F.C. and F.G.; methodology, F.C.; software, F.C.; validation, F.C., F.G., L.L. and Y.N.; formal analysis, F.C. and Y.N.; investigation, F.C.; resources, F.C.; data curation, F.C.; writing—original draft preparation, F.C.; writing—review and editing, F.C., F.G., L.L. and Y.N.; visualization, F.G.; supervision, F.G.; project administration, F.G.; funding acquisition, F.G. All authors have read and agreed to the published version of the manuscript.

Funding

The research was funded by National Key Research and Development Program of China (2017YFB0503402) and the National Natural Science Foundation of China (42074029, 41825009 and 42064002).

Data Availability Statement

The CYGNSS data is available at https://podaac.jpl.nasa.gov (accessed on 5 January 2021). The SMAP data is available at https://nsidc.org/data/smap/smap-data.html (accessed on 7 January 2021). The AGB map can be downloaded at http://lucid.wur.nl/datasets (accessed on 15 January 2021). The CH map can be downloaded at https://landscape.jpl.nasa.gov/ (accessed on 15 January 2021).

Acknowledgments

The authors would like to thank NASA for providing CYGNSS, SMAP data, CH map and Landsat 8 maps as well as LUCID for providing the reference AGB map. We would also like to thank ECMWF for providing land cover map. The numerical calculations in this paper have been done on the supercomputing system in the Supercomputing Center of Wuhan University.

Conflicts of Interest

The author declares no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Mean SMAP soil moisture from January to June, 2019, before interpolation.

Figure 2. Mean SMAP soil moisture from January to June, 2019, after interpolation.

Figure 3. AGB map derived from LUCID.

Figure 4. CH map derived from ICESat/GLAS and environmental modeling.

Figure 5. Relationship between surface Fresnel reflectivity, incidence angle and soil moisture for typical loam soil.

Figure 6. CYGNSS mean traditional Reflectivity from January to June, 2019.

Figure 7. CYGNSS mean correction Reflectivity from January to June, 2019.

Figure 8. The architecture and input features of ANN. (a) traditional method. (b) improved method.

Figure 9. CYGNSS AGBs versus LUCID AGBs based on ANN. (a) traditional method. (b) improved method.

Figure 10. Pan-tropical AGB derived from CYGNSS based on ANN. (a): traditional method. (b): improved method.

Figure 10. Pan-tropical AGB derived from CYGNSS based on ANN. (a): traditional method. (b): improved method.

Figure 11. Absolute errors of AGB retrievals. (a): traditional method. (b): improved method.

Figure 12. The change of surface attributes located in the Andes Mountains from 31 May 2019, to 5 June 2019, of Landsat 8. (a): Point a at the top of Figure 11. (b): Point b at the top of Figure 11. (c): Point c at the top of Figure 11.

Figure 13. Bias Ratio of AGB retrievals. (a) traditional method. (b) improved method.

Figure 14. 300 m annual global land cover time series from 1992 to 2015 [62] and locations of the selected subsets.

Figure 15. CYGNSS CHs versus ICESat/GLAS CHs based on ANN. (a) traditional method. (b) improved method.

Figure 16. Pan-tropical CH derived from CYGNSS based on ANN. (a): traditional method. (b): improved method.

Figure 17. Absolute errors of CH retrievals. (a): traditional method. (b): improved method.

Figure 18. Bias Ratio of CH retrievals. (a): traditional method. (b): improved method.

Table 1

Sensitivity of the CYGNSS observables to the target parameters.

Target Parameters	CYGNSS Observables	Correlation Coefficients
AGB	${\Gamma }_{surface}$	−0.48
	${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)/{\sigma }_0^{coh}$	0.58
CH	${\Gamma }_{surface}$	−0.54
	${\Gamma }_{LR}\left({\epsilon }_s,{\theta }_i\right)/{\sigma }_0^{coh}$	0.60

Table 2

AGB retrievals performance of the traditional method and improved method.

Methods	RMSEs (Mg/ha)	Correlation Coefficients
Traditional method	73.38	0.76
Improved method	64.84	0.80

Table 3

CH retrievals performance of the traditional method and improved method.

Methods	RMSEs (m)	Correlation Coefficients
Traditional method	6.83	0.79
Improved method	5.97	0.83

Table 4

AGB retrievals performance of the improved method over Congo Forest, Amazon Forest, Kalimantan, Brazil moist, Tanzania, and Australia.

Areas	Samples	RMSEs (Mg/ha)	rRMSEs	Correlation Coefficients
Congo	53,879	84.15	0.27	0.77
Amazon	111,592	45.56	0.18	0.79
Kalimantan	22,140	81.95	0.30	0.78
Brazil	20,689	14.00	0.95	0.58
Tanzania	46,828	23.52	1.30	0.46
Australia	10,4036	27.24	0.59	0.82

Table 5

CH retrievals performance of improved method over Congo Forest, Amazon Forest, Kalimantan, Brazil moist, Tanzania, and Australia.

Areas	Samples	RMSEs (m)	rRMSEs	Correlation Coefficients
Congo	53,879	4.86	0.16	0.80
Amazon	111,592	3.88	0.12	0.83
Kalimantan	22,140	6.42	0.23	0.78
Brazil	20,689	2.7	0.43	0.46
Tanzania	46,828	6.14	0.23	0.76
Australia	104,036	2.79	0.74	0.81