On the Synergy of Airborne GNSS-R and Landsat 8 for Soil Moisture Estimation

Abstract

While the synergy between thermal, optical, and passive microwave observations is well known for the estimation of soil moisture and vegetation parameters, the use of remote sensing sources based on the Global Navigation Satellite Systems (GNSS) remains unexplored. During an airborne campaign performed in August 2014, over an agricultural area in the Duero basin (Spain), an innovative sensor developed by the Universitat Politècnica de Catalunya-Barcelona Tech based on GNSS Reflectometry (GNSS-R) was tested for soil moisture estimation. The objective was to evaluate the combined use of GNSS-R observations with a time-collocated Landsat 8 image for soil moisture retrieval under semi-arid climate conditions. As a ground reference dataset, an intensive field campaign was carried out. The Light Airborne Reflectometer for GNSS-R Observations (LARGO) observations, together with optical, infrared, and thermal bands from Landsat 8, were linked through a semi-empirical model to field soil moisture. Different combinations of vegetation and water indices with LARGO subsets were tested and compared to the in situ measurements. Results showed that the joint use of GNSS-R reflectivity, water/vegetation indices and thermal maps from Landsat 8 not only allows capturing soil moisture spatial gradients under very dry soil conditions, but also holds great promise for accurate soil moisture estimation (correlation coefficients greater than 0.5 were obtained from comparison with in situ data).

1. Introduction

Recent studies have shown that the Global Navigation Satellite Systems Reflectometry (GNSS-R) signal has high sensitivity to the terrain dielectric constant, both from ground-based and airborne platforms [1,2,3,4]. This suggests that GNSS-R techniques could be potentially used to monitor vegetation and soil variables related to the water content. The most recent experiments are based on the attenuation and scattering that vegetation cover produces into the GNSS signal before it impinges on the ground and after it is reflected to the receiving antenna. Due to the novelty of this technique, the performance of airborne reflectometers for vegetation/soil applications needs further research, and remains a challenge of increasing interest, especially with the approval of a new dedicated GNSS-R space-borne mission such as the NASA CYGNSS [5] to be launched in October 2016 and the ESA GEROS-ISS [6], under phase A studies.

The synergy between low-resolution passive microwave observations and optical data at higher spatial resolution is the basis of several proposed multi-resolution soil moisture retrieval approaches for the ESA SMOS mission [7,8]. In this regard, downscaling or fusing techniques are applied to combine L-band radiometric data with optical/infrared imagery from multispectral sensors (e.g., SMOS and MODIS). These disaggregation methods generally combine passive microwave observations with visible and near infrared (VNIR) indices, and little research has been done considering the use of alternative spectral bands, such as the short-wave infrared (SWIR) [9] or indices in the red, green, and blue (RGB) space [10]. Since green vegetation spectra in the SWIR region are dominated by water absorption effects, SWIR-based indices are sensitive to changes in leaf water content. The results of the research of [9] supported the use of MODIS-derived high resolution land surface temperature (LST) and SWIR-based vegetation indices to disaggregate SMOS observations. The use of RGB bands, in turn, was proposed early in the history of satellite remote sensing by Kanemasu [11], but abandoned after several studies which seemed to demonstrate that the combinations of near infrared (NIR) and red were far superior than the combinations of green and red. The VIS/IR spectral response of vegetated or soil areas shows a complex pattern due to soil brightness, environmental effects, shadow, soil color and moisture. Still, by calculating a vegetation index from the red, green, and blue bands of a digital image it is possible to track vegetation status, which may be useful information in VIS/IR and passive microwave data fusion or downscaling algorithms.

The Light Airborne Reflectometer for GNSS-R Observations (LARGO) [12] is an innovative instrument, GNSS-R-based, designed and engineered by the Universitat Politècnica de Catalunya (UPC) for real-time measurements of the ground reflectivity. A new approach to enhance the information provided by LARGO reflectivity by using thermal and optical maps from Landsat 8 is presented. The aim of this work is to show how the combined use of optical and thermal data with the GNSS-R signal improves the sensitivity to soil moisture obtained with each of the sensors separately. A total of two RGB-based indices and five NIR and SWIR-based indices were tested as proxies of vegetation and soil moisture status, together with the land temperature, all of them from Landsat 8. The optical and thermal datasets provide the spatial framework needed to adjust the multi-resolution GNSS-R observations and the resulting soil moisture estimates to a fixed grid of 30 m spatial resolution. Thus, the LARGO observations were converted into reflectivity maps matching the Landsat 8 spatial resolution which correspond to 30 m for VNIR/SWIR bands and native 100 m resampled to 30 m pixel size for TIR bands. This study provides new insight on the fusion of GNSS-R data with optical and infrared datasets, in an effort to obtain spatially coherent and accurate soil moisture estimates.

2. Airborne GNSS-R Observations and Field Campaign

The airborne field experiment took place in an agricultural area near the Guareña River, tributary of the Duero River in the center of the Iberian Peninsula (5.36°W; 41.30°N). The main land use is vineyard (100 ha), but some areas of pasture, fallow, irrigated crops, and forest are also present (Figure 1). The airborne campaign was carried out during the 6th August 2014. The flight height was between 500 and 700 m, and lasted from 8 AM to 10:30 AM. During the flight, field measurements of soil moisture (0–6 cm) and surface temperature were taken at 119 ground points (Figure 1). ML3 Theta probes (Delta-T Devices Ltd., Burwell Cambridge, United Kingdom) were used, and after a laboratory test experiment using soil monoliths taken in the area, the manufacturer calibration equation for mineral soils was used. The probes were inserted vertically into the soil and the measurement integrates the soil water content of the first 6 cm of the top layer. The in situ measurement locations were randomly selected along the study area (Figure 1). Four measurements were acquired at each location within an area of 1 m² and geolocated in the field using a differential Global Positioning System instrument. Surface soil temperature was also measured with an infrared thermometer TFA-ScantTemp 380. The in situ measurements were taken from 8 AM to 10:30 AM, at the same time of the flight.

Figure 1. Location of the study area with the main land uses and the sampling layout.

Figure 1. Location of the study area with the main land uses and the sampling layout.

As it is usual under Mediterranean conditions in summer, the in situ observed soil moisture was very low in the top layer, ranging from up to 3% vol. in the vineyard, and up to 14% vol. in the wetter areas such as the pasture field near the river. The field measurements were originally designed to characterize rainfed coverages, such as vineyard; thus, no data over well-watered soils (e.g., irrigated crops) were taken.

The LARGO instrument used in this field campaign is an improved version of a previous LARGO used before in [12], containing an additional switching matrix for calibration purposes. The LARGO instrument is a dual-channel instrument that estimates in real-time the coherent reflectivity ($|\Gamma {({\theta }_{inc})}_{RL}{|}^2$) from transmitted RHCP (Right Hand Circular Polarization) GNSS signals and received LHCP (Left Hand Circular Polarization) reflected GNSS signals using 1 ms of coherent integration time and 1000 incoherent averages. The direct and reflected waveforms are computed by cross-correlating the sampled data with a clean replica of the satellite code, also known as conventional GNSS-R (cGNSS-R). The reflectivity is estimated by dividing the waveform peaks (after noise floor subtraction) of the reflected and direct signals, which are calculated independently. In this sense, the direct signal is used as a reference for estimating the satellite transmitted power. The calibration matrix allows changing which signal (direct or reflected) is sent to each GNSS-R receiver in order to characterize each of the receiving channels and take into account differences among them. The LARGO instrument includes an Inertial Measurement Unit (IMU) to measure the platform movements and to apply antenna pattern corrections in a post-processing stage. It also includes a GPS receiver to geo-locate the platform, compute the position of the specular reflection points, and aid the correlators to apply correlations only against the satellites under observation (not all satellites). The computation of the specular reflection point/area to generate the reflectivity is done by a geometrical approach. For each of the reflectivity measurements, the position of the corresponding satellite on the sky is known using its azimuth and elevation. As GNSS satellites are more than 20,000 km away from the Earth, the direct GNSS signal and the reflected one arrive to the Earth as a plane electromagnetic wave with the same incidence/elevation angle. So, by knowing the height of the platform (GPS data) and the satellite position, the reflection points can be geolocated assuming flat Earth approximation and the atmospheric path as a linear path [13]. The ground-resolution of the LARGO instrument depends on incidence angle and on the surface scattering conditions. In this work, for land observations and low roughness conditions, it is assumed that the main reflected power comes from the first Fresnel zone [14] under specular reflection conditions. Under flat surface conditions, the incidence angle of the GNSS-R data, which is an important parameter for the reflectivity estimation, is the complementary angle of the GNSS satellite elevation angle. The ground-projected first Fresnel zone is initially determined by an ellipse with semi-major and semi-minor axes of 16.5 m and 11.7 m respectively for an incidence angle of 45°, and of 10.1 m and 9.9 m, respectively, for an incidence angle of 10°. The platform movement should also be considered when calculating the Fresnel zone, since, for any given second, there is movement in the along-track direction and the covered area enlarges. Given the platform speed was about 10 m/s, an approximate ground-resolution of 30 × 20 m is obtained, which is comparable to the ground-resolution of LANDSAT images.

When dealing with the terrain topography, the incidence angle must be corrected by the local incidence angle of the reflection area using a Digital Surface Model (DSM). In this case, for all possible reflection points, the conditions for specular scattering were applied using the DSM, and data values were discarded when conditions were not accomplished [15]. Resulting geolocation corrections were of less than 10 m, therefore not affecting the collocation with Landsat (of 30 m resolution).

3. Data and Methods

3.1. LARGO Pre-Processing

A preliminary reflectivity map is shown in Figure 2, where several distinctive regions can be identified, according to the power of the received GNSS-R signal. The higher the reflectivity, the higher the soil moisture, with the blue and red coloring depicting wet and dry conditions, respectively. On the northwestern part, there is a region where reflectivity ranges from −6 to −10 dB, coinciding with an irrigated plot. Nevertheless, on most areas reflectivity ranges from −12 to −18 dB, showing the extremely dry soil conditions observed during the experiment.

Figure 2. (left) Orthophoto of the study area and (right) preliminary reflectivity (dB) map from the UPC LARGO GNSS-R instrument.

Figure 2. (left) Orthophoto of the study area and (right) preliminary reflectivity (dB) map from the UPC LARGO GNSS-R instrument.

The LARGO reflectivity retrieved at point scale was converted into a 30 m map following the inverse distance weighted (IDW) method in order to integrate it with Landsat 8 imagery (Figure 3). The IDW method was chosen since it assumes that each measured point has a local influence that diminishes with distance, therefore matching the high point density cloud obtained from LARGO (Figure 2 and Figure 3). To obtain the maps, the LARGO reflectivity was binned between 40° and 60° elevation angles (hereafter LARGO 40-60) and between 60° and 80° (hereafter LARGO 60-80). The elevation angle results to be the complementary angle to the incidence one. The binning is done because, for the same roughness conditions, the surface appears rougher for lower incidence angles than for larger ones [15]. This leads to lower coherent reflectivity values for lower incidence angles (10°–30°) than for larger ones (30°–50°). Hence, for the same soil moisture and roughness conditions, reflectivity from larger incidence angles is higher than for lower ones, and consequently less affected by noise. Furthermore, it should be noted that the Fresnel reflection coefficients at circular polarization are practically flat in the 10°–50° incidence angle region (this is not the case in linear polarization). Consequently, differences in the Fresnel reflection coefficients can be used to estimate the roughness conditions for that particular resolution cell.

Figure 3. Results of the pre-treatment of the datasets: (a) Landsat 8 true-color composite, (b) Landsat 8 temperature, (c) LARGO observations binned between 40° and 60° incidence angles, (d) LARGO observations binned between 60° and 80° incidence angles, (e) reflectivity map from LARGO interpolation at 40°–60° incidence angles, (f) reflectivity map from LARGO interpolation at 60°–80°.

Figure 3. Results of the pre-treatment of the datasets: (a) Landsat 8 true-color composite, (b) Landsat 8 temperature, (c) LARGO observations binned between 40° and 60° incidence angles, (d) LARGO observations binned between 60° and 80° incidence angles, (e) reflectivity map from LARGO interpolation at 40°–60° incidence angles, (f) reflectivity map from LARGO interpolation at 60°–80°.

3.2. In situ Measurements Mapping

A total of 119 in situ measurements of soil moisture were available (each one computed as the average of the four independent measurements taken at each location). In order to make a spatial comparison analysis with the remotely sensed soil moisture estimates, the in situ measurements were interpolated to a 30 m map covering the sampled area (Figure 4). The same interpolation method applied to the LARGO maps was chosen for a fair comparison between datasets.

Figure 4. Field-based soil moisture map interpolated to 30 m.

Figure 4. Field-based soil moisture map interpolated to 30 m.

3.3. LANDSAT 8 Processing and Indices Retrieval

A Landsat 8 scene (202 path-031 row) of 12 August 2014, was selected for this study. The scene was provided at Level 1T, geometrically and terrain corrected. Landsat 8 carries two instruments: the Operational Land Imager (OLI) sensor, including 8 multispectral bands at 30 m, and the Thermal Infrared Sensor (TIRS), providing two thermal bands at 100 m, resampled to the same spatial resolution of the multispectral bands.

Bands blue (450 nm–510 nm), green (530 nm–590 nm), red (640 nm–670 nm), NIR (850 nm–880 nm), SWIR1 (1570 nm–1650 nm), and SWIR2 (2100 nm–2290 nm) from OLI were selected. Ground reflectance was retrieved using atmospheric correction with the module ATCOR of PCI Geomatics 2013 [16], applying the file metadata.

TIRS acquires data in two spectral channels covering 10.60–11.19 μm and 11.50–12.51 μm (bands 10 and 11). The results of vicarious calibrations of TIRS bands [17] suggested that Band 11 data should not be used where absolute calibration is required, due to an out-of-field stray light in the telescope. Thus, only band 10 data was selected to convert from spectral radiance to brightness temperature and further LST.

A wide collection of vegetation water/drought indices have been obtained from different combinations of NIR and shortwave infrared (SWIR) bands and proposed as advanced and direct vegetation moisture indicators [18] from a suite of remote sensors such as Terra/Aqua, Landsat TM-ETM+, SPOT Vegetation, and AVHRR [19,20,21]. The use of SWIR bands as an indicator of vegetation water status builds from the work of [9], which demonstrated the feasibility of SWIR indices to disaggregate passive L-band observations owing the water absorption effects in this band. Likewise, taking advantage of the SWIR bands of Landsat 8, the so-called Normalized Difference Water Indices [22] were calculated as follows:

NDWI1-NIR = (ρ_NIR − ρ_SWIR1)/(ρ_NIR + ρ_SWIR1)

(1)

NDWI2-NIR = (ρ_NIR − ρ_SWIR2)/(ρ_NIR + ρ_SWIR2)

(2)

where ρ_NIR, ρ_SWIR1, and ρ_SWIR2 refer to the reflectivity of bands NIR, SWIR1, and SWIR2 of Landsat 8, respectively. A high value of NDWI1 and NDWI2 is a consequence of a higher reflectance in the NIR band than in the SWIR1 and SWIR2 bands, regions of water absorption. Thus, this indicates sufficient quantities of water in the canopy for photosynthetic activity and, thus, green and healthy vegetation. Conversely, a negative value indicates canopy stress, dry vegetation or bare soil [9]. Whilst less frequent, the normalized ratios with respect to the red band were also calculated (Equations (3) and (4)).

NDWI1-red = (ρ_red − ρ_SWIR1)/(ρ_red + ρ_SWIR1)

(3)

NDWI2-red = (ρ_red − ρ_SWIR2)/(ρ_red + ρ_SWIR2)

(4)

Whereas the SWIR-NIR-based indices are related to the vegetation water content, these SWIR-red indices are related to the soil water content, since the reflectance absorption in the red band from the soil water content is negligible, whereas a wet soil exhibits reflectance absorption in the SWIR bands. Thus, it is expected that the wetter the soil, the higher the absorption in SWIR bands, and the lower the reflectance. Accordingly, a positive value of NDWI-red indicates wet soil and, conversely, negative values denote dry soils. Nevertheless, one caveat that should be mentioned on the use of these red-based indices is that they require a prevalence of the bare soil response in comparison with the vegetation canopy response. Otherwise, the absorption both on the red and SWIR bands from well-watered vegetation can lead to the same results (negative values) than a dry soil. Therefore, this index seemed adequate for the vineyard plantation, i.e., since the vines are trellised and separated 2.5 m, their spectral response is dominated by the soil.

NDWI indices based in the SWIR region may be considered as water indices, depicting either the vegetation or the soil water status owing the implicit water absorption in this region. The difference between NDWI-NIR and NDWI-red is that the former refers to vegetation water content, and the latter to soil water content, given the distinctive response of the NIR-band to plants and of the red-band to soil.

Finally, indices in the RGB region were also tested, as they have been widely used to monitor vegetation status. The ratio of the reflectance of green and red bands is sensitive to the ratio between chlorophyll and anthocyanin [11]. Its normalized version is known as the Green-Red Vegetation Index (GRVI) [23]:

GRVI = (ρ_green − ρ_red)/(ρ_green + ρ_red)

(5)

where ρ_green and ρ_red are the reflectance of visible green and red, respectively. The response of GRVI to various ground covers may be simple to interpret because, for dense vegetation, the reflectance in the green band is higher than in the red one, reaching high values.

The Greenness index (Equation (6)) computes the proportion of green reflectance in the whole RGB space. Low index indicates barely vegetated covers. As a well-known and widely used index, the NDVI [24] was also used in the present study as a reference (Equation (7)). NDVI refers to the vegetation greenness and vigor. Green and dense vegetation has large NDVI, whereas it is null or negative for bare soil and water, respectively.

Greenness = ρ_green/(ρ_green + ρ_red + ρ_blue)

(6)

NDVI = (ρ_NIR − ρ_red)/(ρ_NIR + ρ_red)

(7)

3.4. Statistical Analysis

3.4.1. Correlations

Spatial relationships between the LARGO reflectivity, the field-based interpolated soil moisture map, and the Landsat-derived indices were quantified via spatial correlations. Comparisons on a pixel-by-pixel basis were applied, then the correlation coefficient (R of Pearson) associated to the whole image was used as an indicator of the similarity of the spatial patterns between each pair of variables. In addition, the values from the imagery at the ground sampling locations were extracted and compared to the in situ soil moisture values. This analysis would reveal the indices that better correlate with soil moisture, surface temperature and LARGO reflectivity.

3.4.2. Linking Model

Several soil moisture retrieval algorithms have been developed based on the unique relationship between soil moisture, NDVI, and LST for a given region under specific climatic conditions and land surface types [25]. This relationship can be expressed through a regression formula that links soil moisture with vegetation indices and LST. In this work, LARGO reflectivity was also included in the multiple regression, following the approach of [7,26] who showed that adding L-band brightness temperatures to the NDVI/LST space was needed to capture soil moisture gradients at high resolution. Hence, the synergies of Landsat-derived LST, water and vegetation indices, LARGO reflectivity, and in situ soil moisture are expressed through a linear linking model as follows:

Y = a₀ + a₁X₁ + a₂ X₂ + a₃ X₃

(8)

where the dependent variable Y represents soil moisture, and the independent variables X₁, X₂, and X₃ are the LST, the LARGO reflectivity and the corresponding index, respectively. a₀, a₁, a₂, and a₃ are the coefficients of the regression equation. All variables were normalized between their maximum and minimum values.

The in situ sample was randomly divided into the training dataset and the validation dataset. The criteria for assigning points to each dataset were to have a balanced number of sampled points of each vegetation cover (i.e., vineyard, fallow, and pasture) both in the training and validation datasets. The training subsample (one half of the total ground measurements, n = 60) was used to establish the coefficients of the regression in the linking model. Hence, a system of linear equations with an equation per location of in situ measurement was set up by extracting the different variable values from their corresponding image. The system was solved to obtain the regression coefficients. Statistical metrics of the regression analysis (multiple correlation, determination coefficient, coefficients of the polynomial fit at a p-value of 95% and root mean square error, RMSE), were used to evaluate the linking model fitting with the training sample and how well the set of variables can predict the soil moisture.

The other half of the samples (validating subsample, n = 59) were used to validate the soil moisture maps obtained after applying the linking model. The statistics chosen for evaluate the comparison between the estimated soil moisture vs. in situ soil moisture were R, the bias or mean difference (observed-estimated), and the root mean square difference (RMSD). Note, that for validation, RMSD was chosen instead of the RMSE chosen for training since, after applying the model, the estimated soil moisture was compared to the independent validation sample.

4. Results and Discussion

4.1. Landsat Reflectance and Indices for the Main Land Covers

Before the correlation analysis, both the spectral response and the indices behavior under different soil and vegetation cover conditions were explored (Figure 5). The reflectance curves (Figure 5a) showed a quite similar behavior for all the main land covers, except for the irrigated areas. The fallow land showed the typical bare soil response, exhibiting the highest reflectance in all bands. This is in agreement with its low water content and the absence of vegetated cover. As shown in the literature, a decrease in water content is reflected as a reflectance increase in the optical wavelength range (0.4–2.5 μm), especially for the red and SWIR spectrum [18]. Many research studies have proved that the whole soil reflectance spectrum diminishes with an increase in soil moisture, with soil reflectance being determined primarily by soil moisture, among other soil characteristics [27,28]. This feature was not easy to identify in this study, due to two main reasons. First, the generally dry conditions of soils during the field experiment (less than 14% vol. for the sampled areas), which were noticeable even for the typically-wet areas of pasture. Due to these conditions, the soil signal was mixed with senescence vegetation, which has similar spectra [29,30]. Second, is the inherently mixed response between soil and plants over the vegetated areas at the Landsat spatial resolution. This effect is especially noticeable in the presence of non-herbaceous plants, e.g., vineyard and forest, which have a sparse distribution. For example, the forested areas are actually scattered trees above a dry soil underneath. Even for the vineyard, which was at the growing season, the spectral response of the plants was mixed with the bare soil one, considering the particular geometry of vineyard plantations (Figure 5c). The vineyard and forested area responses were very similar in all bands (Figure 5a), excepting for those in the SWIR reflectance, where the vineyard registered slightly minor reflectance (higher water content and, thus, more absorption). The strongest relationship between reflectance and water content was observed in the irrigated areas. In this case, no mixed response between soil and plants occurred, owing a more homogeneous growth and a more homogeneous distribution of vegetation, and the covers seemed well-watered according to the very small SWIR reflectance and high NIR reflectance exhibited (Figure 5a).

Figure 5. Landsat 8 bands reflectance (a) and indices (b), both of them averaged for the main land covers in the study area. In (c), several examples of these land covers over a 5-4-3 composite.

Figure 5. Landsat 8 bands reflectance (a) and indices (b), both of them averaged for the main land covers in the study area. In (c), several examples of these land covers over a 5-4-3 composite.

The dominant dry conditions during the experiment are well-captured by all the indices (Figure 5b), each according to its particular behavior. Except for the irrigated areas, vegetation-related indices, such as GVRI, presented negative results, whereas NDVI and Greenness were very low indicating, in all cases, limited vegetation activity.

Regarding the vegetation-water indices, the NDWI1-NIR and NDWI2-NIR reached very low values, close to zero in the first case, in contrast to the watered, canopy-dense irrigated areas. The water-soil related indices, NDWI1-red and NDWI2-red, were both negative, indicating poor water absorption in the SWIR bands and, thus, very dry soil conditions in the case of fallow, vineyard, forest, and pasture. However, as suggested by the rationale of this index, the index does not work properly for the irrigated cover, since it takes the red (i.e., the soil reflectance) as the reference band. For irrigated canopies, such as corn, the soil response is negligible as it remained covered up by the canopy. Thus, the difference between red (null reflectance) and SWIR (small reflectance due to the high water absorption) results were also negative, showing the same trend of the rest of categories. However, these values are actually not soil dryness-induced and hence not comparable to the rest of categories.

4.2. Correlation between Datasets

The relationships between these indices and the soil moisture, both at point scale and for the entire map (Table 1 and Figure 6), showed that indices, such as GVRI, NDWI2-red, and NDWI2-NIR, had good correlations. It can be inferred, thus, that the SWIR2 band results better describe soil moisture content than the SWIR1. Regarding the correlation between indices and temperature, only NDWI1-NIR and NDWI2-NIR showed significant correlation with temperature higher than 0.4 (negative in all cases), both with the Landsat 8 TIR band (LST) and with direct measurements (T). Regarding the LARGO reflectivity, a significant correlation between LARGO 40-60 and soil moisture was found (0.30 < R < 0.40), suggesting a certain potential of using LARGO reflectivities for soil moisture sensing. We hypothesize that the combined use of a water/vegetation index with LARGO reflectivity and LST improves the soil moisture estimation with respect to the use of these variables for independent estimation, as will be tested in the following sections.

Table 1. Correlation coefficients (Pearson) between variables. * Non-significant correlation.

R	Image (30 m)		In Situ (n = 119)
	Soil Moisture	LST (L8)	Soil Moisture	T
GVRI	−0.41	0.06	−0.69	0.07 *
Greenness	−0.05 *	0.14	0.07 *	0.29
NDVI	0.01 *	−0.24	−0.26	−0.27
NDWI1-red	0.20	−0.09	0.37	0.02 *
NDWI2-red	0.58	−0.17	0.59	−0.06 *
NDWI1-NIR	0.25	−0.49	0.06	−0.56
NDWI2-NIR	0.55	−0.41	0.57	−0.48
LARGO 40-60	0.34	−0.22	0.37	−0.28
LARGO 60-80	0.25	−0.20	0.13	−0.05 *

Table 1. Correlation coefficients (Pearson) between variables. * Non-significant correlation.

R	Image (30 m)		In Situ (n = 119)
	Soil Moisture	LST (L8)	Soil Moisture	T
GVRI	−0.41	0.06	−0.69	0.07 *
Greenness	−0.05 *	0.14	0.07 *	0.29
NDVI	0.01 *	−0.24	−0.26	−0.27
NDWI1-red	0.20	−0.09	0.37	0.02 *
NDWI2-red	0.58	−0.17	0.59	−0.06 *
NDWI1-NIR	0.25	−0.49	0.06	−0.56
NDWI2-NIR	0.55	−0.41	0.57	−0.48
LARGO 40-60	0.34	−0.22	0.37	−0.28
LARGO 60-80	0.25	−0.20	0.13	−0.05 *

Figure 6. Relationship between the indices the best correlated with soil moisture and soil moisture.

Figure 6. Relationship between the indices the best correlated with soil moisture and soil moisture.

Other combinations of bands may be considered. Particularly, the blue band is considered by Zhang et al. [18] as a the band least sensitive to vegetation moisture variation compared to the other bands, hence the water-sensitive band (SWIR and red) and the less sensitive band (blue) may be used to maximize moisture variation. Furthermore, the combination SWIR-green has been used to separate built-up features from water bodies [31] and SWIR-red to delineate marsh areas [29], among others. In this line, other combinations of SWIR-visible bands were tested with the current dataset. The combination SWIR-blue and SWIR-green, although plausible, exhibited worse results than those for the SWIR-NIR and SWIR-red combinations (R < 0.50).

The reflectivity from LARGO generally showed low correlations with vegetation indices, probably due to the weak relationship between reflectivity and spectral-based indices. The spectral indices indicate plant and soil conditions (e.g., vigor, photosynthetic activity, water content) that seemed difficult to be tracked through the reflectivity along the entire area for all types of vegetation. However, for the category where the indices took extreme values, such as irrigated crops, some patterns were found. LARGO 60-80, more than LARGO 40-60, showed noteworthy correlations, with R~−0.60 for the vegetation indices NDVI, NDWI1-NIR, and NDWI2-NIR (for which the indices were positive and very high) and positive R~0.60 for the soil-water indices NDWI1-red and NDWI2-red (for which the indices were negative and very small). In other words, index values denoting dense vegetation (i.e., either positive and high NDVI and NDWI-NIR, or negative and small NDWI-red) reached high inverse correlation with the reflectivity, showing an influence of the vegetation on the signal: the denser and healthier the canopy, the smaller the reflectivity. The vegetation type and geometry also showed influence on the reflectivity, not through the index values, but considering the plant size. Indeed, the average LARGO reflectivity for woody-high plants (e.g., trees and vines) resulted ~25% smaller than for the small, mixed with soil, fallow, and pasture covers. In summary, the reflectivity decreases as the canopy height (i.e., vineyard and trees) and density (i.e., irrigated) increases, as it was also found in [32].

4.3. Linking Model Calculation and Validation

The indices GVRI, NDWI2-red, NDVWI1-NIR, and NDWI2-NIR exhibit the best correlations with both soil moisture (Figure 6) and temperature and, to a lesser extent, with LARGO reflectivity. Thus, a linking model is built using each of these indices together with thermal and reflectivity data, and was adjusted to the soil moisture observations of the training subsample (Table 2). A good fit was found for all of them except for NDWI1-NIR with LARGO 60-80, a combination that seemed inadequate for the soil moisture estimation. LARGO 40-60 fit performs better than the LARGO 60-80 fit. This occurs because it is less sensitive to the presence of vegetation (as mentioned before) and more sensitive to the soil water content. This was an expected result as, for larger incidence angles, the surface appears to be smoother than for lower incidence angles and, consequently, coherent reflectivity increases [15]. When comparing the results of the linking model (Table 2) to the individual correlations with soil moisture (Table 1), it can be inferred that the integration of reflectivity, temperature and index drastically improved the results. To further analyze this point, an additional analysis of correlation was made (Table 3) in which, applying the training dataset, it can be seen that the R values improved by adding the reflectivity measurements into the linking model. Indeed, the correlation improves for all cases, even switching the non-significant correlation for the NDWI1-NIR into significant. However, as stated before, LARGO 40-60 inclusion in the model performs better than the LARGO 60-80 one. In other words, the synergy between GNSS-R and conventional multispectral VIS/IR sensors allows for soil moisture estimation at high spatial resolution scales, which opens new insights in this promising technique. However, these results must be carefully analyzed because of the large RMSD of the estimation, close to 2% in a few cases. Taking into account the dry soil conditions during the experiment, this error surpasses in most cases the actual soil water content. In fact, the dry conditions during the field campaign hampered the soil moisture retrieval using only GNSS-R, as seen in Table 1. This is in line with the results obtained by Valencia et al. [33] with the Passive Advanced Unit (PAU) GNSS-R instrument in this area who found that, under very dry conditions (soil moisture below 8% vol.), there is a low reflectivity sensitivity to soil moisture. This could explain the low correlation obtained for reflectivity maps and soil moisture in this study, where soil moisture is below 14% vol.

The weight of each coefficient in Equation (8) was tested, resulting in a higher coefficient a₃ in all cases. The higher weight associated with the index indicates the large influence of the vegetation/soil proxy, expressed as a spectral index, in the soil moisture estimation.

Table 2. Statistical results of the linking model training for the two LARGO datasets and the chosen indices. Correlation coefficient (R), determination coefficient (R²), and root mean square error (RMSE). * Non-significant correlation.

Model Training In Situ Soil Moisture (n = 60)	LARGO 40-60			LARGO 60-80
	R	R2	RMSE (%)	R	R2	RMSE (%)
GVRI	0.69	0.48	1.60	0.69	0.48	1.60
NDWI2-red	0.64	0.42	1.69	0.52	0.27	1.89
NDWI1-NIR	0.51	0.27	1.90	0.12 *	0.01	2.20
NDWI2-NIR	0.59	0.34	1.79	0.55	0.30	1.85

Table 2. Statistical results of the linking model training for the two LARGO datasets and the chosen indices. Correlation coefficient (R), determination coefficient (R²), and root mean square error (RMSE). * Non-significant correlation.

Model Training In Situ Soil Moisture (n = 60)	LARGO 40-60			LARGO 60-80
	R	R2	RMSE (%)	R	R2	RMSE (%)
GVRI	0.69	0.48	1.60	0.69	0.48	1.60
NDWI2-red	0.64	0.42	1.69	0.52	0.27	1.89
NDWI1-NIR	0.51	0.27	1.90	0.12 *	0.01	2.20
NDWI2-NIR	0.59	0.34	1.79	0.55	0.30	1.85

Table 3. Results of correlation when adding consecutive parameters in the linking model (only the vegetation index; index and LST; index, LST, and reflectivity). * Non-significant correlation.

	R with the Training Sample (n = 60)
Model Options	Only Index	Index+LST	Index+LST+ LARGO 40-60 (LARGO 60-80)
GVRI	−0.66	0.67	0.69 (0.69)
NDWI2-red	0.52	0.52	0.64 (0.52)
NDWI1-NIR	−0.10 *	0.10 *	0.51 (0.11 *)
NDWI2-NIR	0.43	0.52	0.59 (0.55)

Table 3. Results of correlation when adding consecutive parameters in the linking model (only the vegetation index; index and LST; index, LST, and reflectivity). * Non-significant correlation.

	R with the Training Sample (n = 60)
Model Options	Only Index	Index+LST	Index+LST+ LARGO 40-60 (LARGO 60-80)
GVRI	−0.66	0.67	0.69 (0.69)
NDWI2-red	0.52	0.52	0.64 (0.52)
NDWI1-NIR	−0.10 *	0.10 *	0.51 (0.11 *)
NDWI2-NIR	0.43	0.52	0.59 (0.55)

Once the best proxies for the linking model were determined, and the regression coefficients were retrieved using the training samples, maps of soil moisture at 30 m resolution were produced (Figure 7). These maps were validated at point scale through comparison with the validating subsample in situ measurements (Table 4 and Figure 8) and through comparison with the field-based interpolated soil moisture map described in Section 3.2 (Figure 4) using spatial correlation (Table 5). When comparing the remotely sensed maps in Figure 7 with the field-based soil moisture map in Figure 4, slightly higher soil moisture content is noticeable in the remotely-sensed maps. This is quantitatively confirmed by the negative bias resulting when comparing the field-sampled to the remotely-sensed soil moisture (Table 4 and Table 5) for all the linking model alternatives.

Considering the soil moisture spatial distribution, the wetter areas near the river are detected in all cases. The driest areas (e.g., buildings, roads), in turn, are also detected except for the case of using NDWI1-NIR, where the expected patterns of soil moisture seemed to be switched (dry to wet and vice-versa). Indeed, for this alternative R was non-significant, and the error and bias were the highest of all datasets. Additionally, some artifacts were detected in the building area (northwestern part in Figure 7b–d) for the maps retrieved with GVRI and NDWI1-NIR, where this area was confusedly represented. In view of these results and the corresponding soil maps, more accurate soil moisture estimates can be obtained using the NDWI2-NIR index (SWIR2 band) in the linking model (Figure 7g,h). This is in agreement with the results of the correlation analysis on Table 1. Still, from the validation results at point and map scales (Table 4 and Table 5), there are no substantial differences in the statistical scores obtained with the different proposed model alternatives, except for the above-mentioned NDWI1-NIR. This result confirms the hypothesis that the SWIR1 band is not adequate for soil moisture estimation. Regarding the reflectivity integration in the linking model, the LARGO 40-60 revealed slightly better results for the training dataset (Table 2), but the two selected ranges of incidence angles result in comparable results for the validating dataset and field-based interpolated soil moisture map.

Figure 7. (a) Ground observations distribution. (b to h) Estimated soil moisture (% vol.) maps following the seven alternatives considered.

Figure 7. (a) Ground observations distribution. (b to h) Estimated soil moisture (% vol.) maps following the seven alternatives considered.

Table 4. Statistics of estimated soil moisture resulting from the linking model vs. in situ soil moisture measurements: correlation coefficient (R), mean difference (bias), and root mean square difference (RMSD). * No-significant correlation.

Model Validation In Situ Soil Moisture (Point Scale)	LARGO 40-60			LARGO 60-80
	R	bias (%)	RMSD (%)	R	bias (%)	RMSD (%)
GVRI	0.67	−0.205	2.619	0.70	0.010	2.462
NDWI2-red	0.67	−0.748	2.675	0.67	−0.753	2.650
NDWI1-NIR	0.16 *	−0.353	3.284	---	---	---
NDWI2-NIR	0.66	−0.801	2.894	0.68	−0.771	2.733

Table 4. Statistics of estimated soil moisture resulting from the linking model vs. in situ soil moisture measurements: correlation coefficient (R), mean difference (bias), and root mean square difference (RMSD). * No-significant correlation.

Model Validation In Situ Soil Moisture (Point Scale)	LARGO 40-60			LARGO 60-80
	R	bias (%)	RMSD (%)	R	bias (%)	RMSD (%)
GVRI	0.67	−0.205	2.619	0.70	0.010	2.462
NDWI2-red	0.67	−0.748	2.675	0.67	−0.753	2.650
NDWI1-NIR	0.16 *	−0.353	3.284	---	---	---
NDWI2-NIR	0.66	−0.801	2.894	0.68	−0.771	2.733

Figure 8. Relationship between the estimated soil moisture resulting from the linking model and the in situ soil moisture measurements using GVRI, NDWI2-red, and NDWI2-NIR.

Figure 8. Relationship between the estimated soil moisture resulting from the linking model and the in situ soil moisture measurements using GVRI, NDWI2-red, and NDWI2-NIR.

Table 5. Statistics of estimated soil moisture resulting from the linking model vs. the in situ soil moisture map: correlation coefficient (R), mean difference (bias), and root mean square difference (RMSD). * Non-significant correlation.

Model Validation In situ Soil Moisture (Map)	LARGO 40-60			LARGO 60-80
	R	bias (%)	RMSD (%)	R	bias (%)	RMSD (%)
GVRI	0.37	−0.046	1.632	0.28	0.107	1.844
NDWI2-red	0.57	−0.571	1.394	0.57	−0.315	1.267
NDWI1-NIR	0.23 *	−0.475	1.696	---	---	---
NDWI2-NIR	0.45	−0.287	1.520	0.35	−0.150	1.839

Table 5. Statistics of estimated soil moisture resulting from the linking model vs. the in situ soil moisture map: correlation coefficient (R), mean difference (bias), and root mean square difference (RMSD). * Non-significant correlation.

Model Validation In situ Soil Moisture (Map)	LARGO 40-60			LARGO 60-80
	R	bias (%)	RMSD (%)	R	bias (%)	RMSD (%)
GVRI	0.37	−0.046	1.632	0.28	0.107	1.844
NDWI2-red	0.57	−0.571	1.394	0.57	−0.315	1.267
NDWI1-NIR	0.23 *	−0.475	1.696	---	---	---
NDWI2-NIR	0.45	−0.287	1.520	0.35	−0.150	1.839

The bias in Table 4 and Table 5 is very large, particularly if considering the field point measurements, which led to an unacceptable RMSD. However, the correlation coefficients are encouraging (R > 0.66), slightly lower for the pixel-by-pixel comparison at the interpolated map (Table 5) than for the sampled locations (Table 4), where the soil moisture was directly observed. Again, it should be pointed out that we are considering very specific conditions to depict the relationships between vegetation and water soil status. Also, small variations in soil moisture led to a high standard deviation and RMSD owing the small dynamic range of this variable under those conditions. To illustrate this, the model validation was repeated removing just the four pasture points of the validation subsample. These four points had the highest soil moisture (~13% vol.) as compared to the remaining 55 points (~3% vol.), resulting in a reduction of the RMSD of 40% on average.

The validation exercise compares the estimated soil moisture to the in situ measurement under the assumption that both data correspond to soil moisture within the same soil layer. However, the synoptic model considers LST (skin temperature), L-band reflectivity (5 cm), and a vegetation index (with vegetation status supposedly responding to root-zone water status). Thus, the resulting soil moisture can be considered as an integrated value of different soil depth responses. We hypothesize the mismatch in sensing depth between LST, L-band reflectivity, and in situ probes has a limited impact in the obtained results. However, given the important weight of the vegetation index in the model, the estimated soil moisture could be representative of a deeper soil layer than that sensed by the in situ soil moisture. This mismatch could partly explain the negative bias (overestimation) and the high values of errors resulting for the comparisons. Further research should be directed in determining the sensing depth of the synergistic soil moisture estimates.

5. Conclusions

The objective of this work was to test the capability of an airborne GNSS-R sensor (LARGO) to estimate and characterize both vegetation cover and soil moisture status. In view of the expected development of new GNSS-R-based missions, the potential interactions between soil moisture estimates and observations at high spatial resolution from airborne or satellite sensors allow for the development of synergistic approaches that can be later transferred to sensors on-board satellite platforms. New insight on the synergy between airborne GNSS-R reflectivity and satellite optical, infrared and thermal data is presented, using a multiple regression model that links the soil moisture to land surface temperature, vegetation/water indices and GNSS-R reflectivity altogether. Taken these datasets separately, and considering the particular weather and vegetation conditions of the experiment, LARGO reflectivity did not show distinctive sensitivity to the main vegetation covers present in the area. However, a noticeable relationship to soil moisture was detected. This relationship was strongly reinforced if the reflectivity was merged to the surface temperature and some vegetation/water indices, all of them retrieved from Landsat 8 bands. Indeed, encouraging correlations were found (R > 0.60) when applying a multiple regression model that linked the soil moisture to the temperature, the vegetation index, and the LARGO reflectivity binned by a certain range of incidence angles. These correlations were much higher than taking each variable separately. The remotely-sensed soil moisture maps obtained with the proposed linking model were highly correlated with the ground measurements, even though a high bias and RMSD were observed. Future campaigns, under diverse environmental conditions, should be performed to further our understanding of the limitations and applicability of the proposed fusion scheme. The methodology developed has the advantage of being applicable at a variety of spatial scales, from several meters, using airborne platforms, to kilometers using space-borne observations. This versatility allowed exploring the optimal spatial resolution to transform the point scale GNSS-R observations into reflectivity maps suited for soil moisture information. In this case, the 30 m from Landsat results were adequate, but further analyses are needed to explore higher spatial resolution sensors.