Multi-constellation GNSS interferometric reflectometry with mass-market sensors as a solution for soil moisture monitoring

. Per capita arable land is decreasing due to rapidly increasing population, and fresh water is becoming scarce and 10 more expensive. Therefore, farmers should continue to use technology and innovative solutions to improve efficiency, save input costs, and optimise environmental resources (such as water). In the case study presented in this manuscript, the GNSS-IR technique was used to monitor soil moisture during 66 days, from December 3, 2018, to February 6, 2019, in the installations of the Cajamar Centre of Experiences, Paiporta, Valencia, Spain. Two main objectives were pursued. The first was the extension of the technique to a multi-constellation solution using GPS, GLONASS, and GALILEO satellites, and the second 15 was to test whether mass-market sensors could be used for this technique. Both objectives were achieved. At the same time the GNSS observations were made, soil samples taken at 5 cm depth were used for soil moisture determination to establish a reference dataset. Based on a comparison with that reference data set, all GNSS solutions, including the three constellations and the two sensors (geodetic and mass-market), were highly correlated, with a correlation coefficient between 0.7 and 0.85. GLONASS, and GALILEO constellations solutions and the use and comparison of a geodetic and mass-market antenna solutions.


Introduction 20
Soil moisture is a fundamental component of the hydrological cycle, and a key observable variable for optimising agricultural irrigation management. Additionally, soil moisture monitoring has been one of the main goals of the remote sensing satellite missions Soil Moisture and Ocean Salinity (SMOS), (Kerr et al., 2001), Soil Moisture Active Passive (SMAP), (Chan et al. 2016), and Sentinel-1, (Mattia et al., 2018). SMOS is used to derive global maps of soil moisture every three days at a spatial resolution of about 50 km, SMAP every two-three days with a spatial resolution of about 40 km (gridded to 36 km since the 25 radiometer is the only instrument on board that works), and one Sentinel-1 satellite 12 days (two Sentinel-1 satellites are in orbit which decreases the revisit time) with a spatial resolution of about 1 km.
To obtain information about soil moisture at a very local scale and continuously, Global Navigation Satellite System (GNSS) reflectometry began to be tested as a possible solution (Masters et al., 2002;Zavorotny et al., 2003;Katzberg et al., 2005). This was possible because GNSS satellites transmit in the L-band (microwave frequency), so the GNSS signal reflected by 30 2 nearby surfaces and recorded by the antenna contains information about the environment surrounding the antenna (scale of about 1000 m 2 ). In particular, the ground-reflected global positioning system signal measured by a geodetic-quality GNSS system can be used to infer temporal changes in near-surface soil moisture. This technique, known as GNSS-interferometric reflectometry (GNSS-IR), analyses changes in the interference pattern of the direct and reflected signals, (Fig. 1), which are recorded in signal-to-noise ratio (SNR) data, as interferograms. Thus, GNSS-IR can be considered as another remote sensing 35 technique for monitoring soil moisture in a local scale and continuously, independent from climatological conditions (the technique is valid in raining and foggy conditions) and illumination (day or night). Temporal fluctuations in the phase of the interferogram are indicative of changes in near-surface (depth of about 5-7 cm) volumetric soil moisture content, (Larson et al., 2008a(Larson et al., , 2008b. Commercially available geodetic-quality GNSS receivers and antennas can be used for GNSS-IR. The method has been tested 40 with the Global Positioning System (GPS) satellite constellation, and it has been shown to provide consistent measurements of upper surface soil moisture content, (Larson et al., 2008a(Larson et al., , 2008b(Larson et al., , 2010Larson and Nievinski, 2013;Chew et al., 2014Chew et al., , 2015Chew et al., , 2016, Vey et al., 2015Wan et al., 2015;Chen et al., 2016;Zhang et al., 2017).
With the use of the GPS constellation, the GPS-IR reflection footprint is far from homogeneous, Fig. 2, and some tracks cannot be included in the process and analysis (Vey et al., 2015;Chew et al., 2016). Therefore, GPS-IR needs to evolve to Global 45 Navigation Satellite System reflectometry, GNSS-IR, where multi-constellation observation provides the solution. The integration of new navigation satellite constellations will produce a more homogeneous footprint around the antenna (Fig. 2). Roussel et al. (2016) introduced the GLONASS Russian constellation to retrieve soil moisture over bare soil, but there are no references in the literature for the European GALILEO or Chinese BEIDU constellations. Roesler and Larson (2018) provided a software tool for generating map GNSS-IR reflection zones that support GPS, GLONASS, GALILEO, and BEIDU 50 constellations.
Therefore, the first novelty of this research was to extend, compare and combine the GPS-IR methodology to a multiconstellation scenario (GPS, GLONASS, and GALILEO; BEIDU is not introduced in this research because the antennas used in the experiment are not able to decode BEIDU signals), which will produce a much larger sample set of observations around the antenna than is obtained with only the GPS constellation, as shown in Fig. 2. 55 Additionally, geodetic-quality GNSS receivers and antennas are an expensive solution. If we keep in mind that the final market will be the agricultural market, a technique developed using those devices will never be introduced into the sector. Thus, the (main) second novelty of this research was the introduction of mass-market GNSS sensors as the basis for the technique. If the use of these mass-market devices can be confirmed, it will be possible to use them (one or several at the same time to add redundancies) at a very low cost. 60 3 2 Materials and methods

Location of the experiment
The experiment was conducted in the installations of the Cajamar Centre of Experiences, located in Paiporta, Valencia, Spain The centre began its activities in 1994. Some of the research topics carried out by the centre are the valorisation of agricultural by-products and the use of microorganisms in food, pharmaceuticals, and aesthetics using the latest biotechnology resources; the design of new containers and bio-functional formats for the marketing of healthy foods with high added value; improvement in irrigation automation, biological control management, and agronomic management in organic production; and the introduction of alternative value crops and new varieties that guarantee the sustainability of agricultural sector. 70

Instrumental and observations
A geodetic GNSS receiver (Trimble R10 GNSS receiver, from the Department of Cartographic Engineering Geodesy and Photogrammetry of the Universitat Politècnica de València) and a mass-market receiver (Navilock GNSS receiver based on a u-blox 8 UBX-M8030-KT chipset with a built-in antenna) connected to a Raspberry Pi 3 as a control device and for storing the observations, were used to obtain multi-constellation SNR observables (GPS, GLONASS and GALILEO). Five seconds 75 sample rate observations were obtained simultaneously for both sensors (Fig. 3).
The radio-signal structure of GPS, GLONASS and GALILEO systems are similar. Different carrier signals in the L-band are broadcast, L1 and L2 corresponds with the two main frequencies of the signal emitted from the GPS satellites and E1 and E5 with the two main frequencies of the signal emitted from the GALILEO satellites. In contrast to GPS and GALILEO, GLONASS satellites transmit carrier signals at different frequencies from a basic L frequency, GLONASS L1 frequencies are: 80 where fo = 1602.0 MHz, and ∆! "# = 0.5625 678, and k is the carrier number assigned to the specific GLONASS satellite (Hoffmann et al., 2008). Thus, the frequency for each satellite should be computed and included in the GLONASS file. 85 The frequencies used in the experiment were L1, for the GPS and GLONASS satellite constellations and E1 for the GALILEO constellation. This choice was forced because the mass-market device could not track the L2 or E5 satellite signals. However, Vey et al. (2011) showed that the soil moisture root mean square difference between L2C and L1 was only 0.03 m 3 /m 3 . L2C corresponds to the Civil L2 signal of the block satellites IIR-M and IIF of the GPS constellation, available only since 2005 when the first block IIR-M was launched. This signal is designed specifically to meet commercial needs, which increases 90 robustness of the signal, improve resistance to interference, and improve accuracy (Leick et al., 2015).

4
The GNSS-IR footprint for a single rising or setting satellite is an elongated ellipse in the direction of the satellite track (Fresnel ellipse or zones; Larson et al., 2010;Wan et al., 2015;Vey et al., 2015;Roesler and Larson, 2018). As the satellite rises and the elevation angle increases, the Fresnel zone becomes smaller and closer to the GNSS antenna. Data with elevation angles higher than 30 degrees should be discarded from the SNR series because they contain no significant oscillations and cannot be 95 retrieved reliably. Data with elevation angles lower than 5 degrees should also be discarded in order to avoid strong multipath effects from trees, artificial surfaces, and structures surrounding the antenna. A GNSS satellite takes about one hour to rise from an elevation angle of 5 degrees to an angle of 30 degrees.
The geodetic GNSS receiver store the observations (including SNR data) in the commonly used RINEX files, so the elevation and azimuth of a satellite for an epoch should be computed from the observation RINEX file and the navigation RINEX file, 100 (Hofmann-Wellenhof et al., 2008).
The mass-market receiver uses NMEA GSV sentences to provide integer numbers for elevation, azimuth and signal-to-noise ratio (SNR) directly. NMEA is an acronym for the National Marine Electronics Association. GNSS NMEA is a standard data format supported by all manufacturers to output measurement data from a sensor in a pre-defined format in ASCII. In the case of GNSS, it output position, velocity, time and satellite related data (for the constellations that the antenna can decode). There 105 are quite a few NMEA messages or sentences, specifically, GSV sentences provide integer numbers for elevation, azimuth and signal-to-noise ratio.
The results were compared with soil moisture measurements based on soil samples taken at a depth of 5 cm and weighed before and after being dried (gravimetric method) in a laboratory (Fig. 4). These measurements were considered the reference dataset. The soil samples were taken one per day except weekends and the location, in comparison with the antenna position, 110 can be seen in Fig. 2.
In total, 66 days of measurements, from December 3, 2018, to February 6, 2019, were observed, processed, and analysed. The height of the antennas from the ground was 1.80 m for the geodetic GNSS device and 1.84 m for the mass-market device.
Precipitation data were added in the final plot results. These data were obtained from a meteorological station located in the Cajamar Experiences Centre (100 meters from the GNSS antennas). 115

Theoretical background
The theoretical background is based on the procedure developed by Larson et al., (2010)  2) A low-order polynomial (second degree) is fit to the Slineal in order to eliminate the direct satellite signal, so that the reflected signal is isolated: 9 :;<=>: ?=@:=AB=C , (Wan et al., 2015;Chew et al., 2016).
3) A Lomb-Scargle periodogram (Lomb, 1976;Press et al., 1992;Roesler and Larson, 2018), is then computed from 125 9 :;<=>: ?=@:=AB=C , and the track goes to the next step only if there is a clear signal that reflects a primary wave. Tracks with multiple peaks or low maximum average power (less than four times the background noise) are not included in the next step. If the Lomb-Scargle periodogram is computed using the sine elevation angle as the input X axis, the result converts the frequency into antenna height in the output X axis. Only tracks with computed antenna height consistent with the measured antenna height (less than 0.1 meters difference), go to the next step. 130 4) The selected tracks are modelled using the expression below: The equation means that 9 :;<=>: ?=@:=AB=C can be modelled in terms of the amplitude A and phase offset f of a primary wave.

135
l is the GNSS wavelength (L1 for GPS and GLONASS and E1 for GALILEO), e is the satellite elevation, and h is the antenna height, which is assumed to be a constant due to the low signal penetration on the ground (Chew et al., 2014;Roussel et al., 2016;Zhang et al., 2017). The least squares algorithm (Strang and Borre, 1997;Leick et al., 2015) is used to estimate A and f. The file containing the NMEA observations from the mass-market antenna was used to generate three different input files for the processing process, one for each satellite constellation. However, due to the integer nature of the SNR, elevation, and azimuth observation numbers, an extra processing step was included for the mass-market observation files. This step used the navigation files from the International GNSS Service (IGS) repository (http:/www.igs.org) to compute float numbers for 170 elevation and azimuth values of the observed satellites.
The rest of the processing followed the steps defined in the previous section. Only full GNSS tracks data covering more than 30 minutes and cover more than 10 degrees of elevation in its trajectory were considered in our study.

Results
The geodetic antenna SNR data in volts for satellite GPS number 23 are shown in Fig. 5a, the SNR data with the direct signal 175 removed are shown in Fig. 5b, the Lomb-Scargle periodogram for the SNR reflected signal is shown in Fig. 5c, and the SNR reflected signal with the adjusted wave (Step 4 in the previous section) is shown in Fig. 5d. Fig. 6 portrays the same concepts for the same satellite but using the mass-market antenna observations. Fig. 7 and 8 portray the same concepts for the GLONASS satellite number 5, and Fig. 9 and 10 display these for the GALILEO satellite number 21.
The SNR values from the geodetic antenna and the mass-market antenna for the GPS constellation are similar, as suggested in 180 Li and Geng (2019), because the u-blox chipset uses an active, right-handed, circularly polarized antenna with uniform antenna gain. However, the SNR values for GLONASS and GALILEO present a systematic bias of about 3-5 db-Hz between the geodetic and mass-market antennas (Fig. 7a and 8a and Fig. 9a and 10a).
A linear relationship between reference data and every GNSS constellation and antenna was computed using the methodology proposed by Zhang et al. (2017), the results can be seen in Table 1  The numerical values for Fig. 11, 12, 13 and 14 are listed in Table 2, where MAE is the mean absolute error, RMSE is the root mean square error, mean and Std. are the mean and the standard deviation respectively between the GNSS antennas and the reference values. The Pearson correlation coefficient can be used to summarize the strength of the linear relationship between two data samples. Spearman correlation can be used to summarize if two variables are related by a nonlinear relationship, such 200 that the relationship is stronger or weaker across the distribution of the variables.

Discussion
Based in the results summarized in Table 2 Our RMS results using the a priori slope values of 65.1° are comparable with those obtained by Zhang et al. (2017), who 225 processed six months of continuous observations and obtained a mean standard deviation value of 0.036 m 3 /m 3 , and those of Vey et al. (2015), who processed 6 years of observations and obtained a standard deviation value of 0.06 m 3 /m 3 .
The SNR bias between the geodetic and mass-market antenna for GLONASS and GALILEO constellations (Fig.7b and 8b and Fig. 9b and 10b) has no effects in the final phase offset variations for the adjusted wave.

According to
Step 3 of Section 2.3, the 70% of the GPS tracks recorded by the geodetic antenna were considered valid for 230 processing, as were 73% for GALILEO, and 74% for GLONASS. This percentage is reduced to around a 10% if we consider the tracks recorded by the mass-market antenna. Nonetheless, one of the main important problems in this research is related with the selection of the correct tracks to be processed and adjusted using Step 4 of Section 2.3. Based on the mentioned criteria (tracks with multiple peaks or low maximum average power and computed reflector height consistent with the measured antenna height), some tracks that should not be processed are finally processed (around 8% of all tracks irrespective the 235 constellation). These wrongly processed tracks introduce outliers in the computed VGNSS, which are eliminated in the daily final mean VGNSS computation because they produce a high RMS in the daily computations using all satellites. One way to accomplish this task could be to use good figures, such as those from Fig. 5c Fig. 5d, to produce a valid set of training images and use machine learning tools (image recognition) to decide automatically whether a new track can be considered as a good track (so it can be processed) or not. This idea is currently under development. 240 In situ observations are needed to solve Eq. 3 (VResidual parameter). However, if there are no reference values, this constant cannot be included, and the results will present an offset in comparison with the real values. A possible solution would be the estimation of the parameter based on the soil type (URL 1); though, that requires having a long enough time series to make the assumption that, at some point during the time series, soil moisture was low enough to hit the residual value. However, the results can be used in a relative way, that is, can be used to infer VWC variations from one day to another. This relative 245 comparison can be performed only if the observations are continuous. If there is an interruption in the raw data (because the antenna is turned off) of more than two or three hours, the previous reference is lost and the relative comparisons should start again (from the moment the antenna is turned on again). In situ observations are also needed if we want to adjust the linear relationship between the computed phase offset and the soil moisture, as is developed in Zhang et al. (2017); however, in case the linear relationship is positive, a value of 65.1° can also be used to obtain acceptable results. 250

Conclusions
The case study presented in this research is focused on the GNSS SNR data acquisition and processing using the GNSS-IR technique to monitor soil moisture. The main objectives of this research were the use, comparison and combination of GPS, Eliminado: if there are no reference values, they can be estimated based on the soil type (URL 1). Though, that requires having a long 9 GLONASS, and GALILEO constellations solutions and the use and comparison of a geodetic and mass-market antenna solutions. 260 Independent GPS, GLONASS, and GALILEO solutions were generated to demonstrate that the technique can be extended to a multi-constellation solution. This is necessary because a single constellation solution presents a reflection footprint that is far from homogeneous around the antenna and because 30-35% of the observed satellite tracks of the geodetic antenna are not valid for processing (40-45% if the mass-market antenna is considered).
The use of a mass-market GNSS antenna was confirmed to be a viable tool for GNSS-IR, with the caution of using the IGS 265 navigation files to transform the observed integer numbers obtained in the NMEA messages for the elevation and azimuth of the satellites into floating numbers. With the use of mass-market sensors, it will become possible to design scenarios with several GNSS stations generating redundant observations. Therefore, maps of soil moisture variations by specific and selective areas of soil, cultivation, and/or management can be generated, instead of obtaining only an average value for the entire observation area. 270 GNSS-IR is still a technique with numerous technological challenges in order to becoming a competitive solution with respect to current observation techniques, but it has great potential with regard to continuity of observation (can be implemented in a real or quasi-real time scenario), precision, and measurement acquisition cost if mass-market antennas are used.
Acknowledgements. The authors want to thank the Cajamar Centre of Experiences staff for all their help and collaboration 275 during the field observation. The authors would also like to thank the comments from the editor and the two anonymous reviewers who have helped to improve the original manuscript.

Data availability
GNSS raw observations used to conduct this study are available upon request from the corresponding author (Angel Martin)