Introduction

GNSS-IR is a method for estimating environmental parameters around a geodetic-quality GNSS site. Unlike other reflection techniques, where an antenna is designed to measure reflection signals (Löfgren et al. 2011; Camps et al. 2013) or a geodetic antenna is rotated to improve its ability to measure reflections (Anderson 2000), GNSS-IR uses data collected with (nominally) multipath-suppressing geodetic-quality GNSS antennas in an upright orientation. GNSS-IR has been demonstrated and validated for measuring surface soil moisture (Larson et al. 2008), snow depth (Larson et al. 2009; Nievinski and Larson 2014c, d), permafrost melt (Liu and Larson 2018), tides (Larson et al. 2013, 2017; Löfgren et al. 2014; Roussel et al. 2015), ice-up (Strandberg et al. 2017), firn density (Larson et al. 2015), and vegetation water content (Wei et al. 2015). In addition to these practical demonstrations, Felipe Nievinski developed a simulator that allows a user to test the reflection characteristics for different experimental configurations and surfaces (Nievinski and Larson 2014a, b).

At some level, all GNSS-IR studies are based on the analysis of SNR patterns created by the interference of direct and reflected (or multipath) GNSS signals. There is significant literature on the inherent frequencies in GNSS multipath, and we will not repeat it here (Georgiadou and Kleusberg 1988; Ge et al. 2000; Ray and Cannon 2001; Axelrad et al. 2005; Bilich and Larson 2007). While the multipath frequencies for a planar reflector change as a satellite rises or sets, Axelrad et al. (2005) proposed a simple change of variable (using \(\sin e\) rather than e) that yields one multipath frequency per rising/setting satellite arc. Ignoring the direct signal contribution, SNR data for a single satellite and receiver can be modeled as:

$${\text{SNR}}(e)=A\left( e \right)\sin \left( {\frac{{4\pi {H_R}}}{\lambda }\sin e+\phi } \right),$$
(1)

where e is the GNSS satellite elevation angle with respect to the horizon, \(\lambda\) is the GNSS wavelength, \(\phi\) is a phase constant, HR is the vertical distance between the GNSS antenna phase center and the horizontal reflecting surface, and A(e) represents the amplitude of the SNR data. To be clear, this representation of SNR data is time dependent because e is a function of time. A fuller discussion of the contributions to A(e) can be found in Nievinski and Larson (2014a).

When HR is fixed, surface soil moisture can be derived from the estimated changes in \(\phi\) (Larson et al. 2008; Chew et al. 2016). Using similar assumptions, A can be used to measure vegetation water content (Wei et al. 2015). Changes in A are also important for applications such as sea ice detection (Strandberg et al. 2017). In this short note we ignore A and \(\phi\) and focus on the inherent multipath frequency, \(2{H_{\text{R}}}/\lambda\). By estimating the multipath frequency, we have a simple way to determine HR. This has also been called the reflector height (Larson and Nievinski 2013), although to be clear, it is not a “height” in a geodetic sense. Here, we will mostly use HR so as to limit the confusion between HR and orthometric and ellipsoidal heights. In the next sections, we will discuss our software for GNSS-IR frequency relationships. In particular, we are providing software to do the following:

  • Translate GNSS data files stored in the RINEX format.

  • Map GNSS-IR reflection zones.

  • Calculate the average Nyquist frequency for a GNSS-IR installation.

  • Estimate dominant frequencies (and thus HR) from GNSS SNR data.

Extracting SNR observations needed for GNSS-IR from a RINEX file

Most GNSS networks provide carrier phase and pseudorange data to users in the RINEX format (Gurtner and Estey 2007). These dual-frequency ranging observations are then used with high-precision geodetic or surveying software, with precise ephemerides, to compute daily (or more frequent) Cartesian positions. Unfortunately, the parameters needed for GNSS-IR cannot always be easily extracted from the outputs of high-precision geodetic or surveying software. Here we provide Fortran 77 code that translates the GPS (and more generally GNSS) observations stored in a RINEX file into a format usable for reflections research. In the first (RinexSNR), the elevation and azimuth angles of GPS satellites with respect to the local horizon are computed using the GPS navigation message. The latter is also stored in the RINEX format. The times of the observations (in GPS time, seconds of the day) are also extracted and the SNR observations on the L1, L2, and L5 frequencies are saved. The output format is in columns and thus the data can be easily loaded using other programming languages such as Matlab and python. The current version of the code only reads RINEX version 2.11 format and files with no more than 15 observation types. Signals from non-GPS constellations are ignored.

To support GNSS-IR, we have provided a separate piece of Fortran code (RinexSNR_GNSS) that extracts GNSS SNR observations from a RINEX file. The main difference between the two programs is that the first uses real-time navigation messages to compute satellite ephemerides and the second uses precise ephemerides. The current constellations supported are GPS, GLONASS, GALILEO, and BEIDOU. Because all GNSS constellations nominally number their satellites the same way, we rename the non-GPS satellites by adding 100 (GLONASS), 200 (GALILEO), or 300 (BEIDOU). Additional instructions are provided in the readme file.

Both Fortran translation programs use the Cartesian station location in the RINEX header to compute the satellite/station elevation angle. If the RINEX file doesn’t have station coordinates, the code stops. It is not necessary that the station location be extremely accurate—but we suggest it be within 50 m of its true position. Most RINEX files created by geodesists and surveyors have much better station coordinates than this, with the exception of some data from the cryosphere. In these cases, the archives often use a single station location for all station files, even though the station is moving rapidly (e.g., hundreds of meters per year). For these users, we allow time-varying receiver coordinates to be read from an external text file. The user should input a Cartesian position at a given epoch and a Cartesian velocity, in meters and meters/year, respectively. This option is not necessary for other GNSS-IR sites.

GNSS-IR reflection zones

The equations for a Fresnel zone near the surface of the Earth are given in the appendix of Larson and Nievinski (2013). The sizes of these elliptical sensing zones are directly sensitive to HR, the satellite elevation angle (e) and the GNSS transmitter frequency (L1, L2, or L5). The orientation of the Fresnel zone with respect to the GNSS antenna depends on the azimuth angle of the satellite. The Fresnel zones gets smaller and closer to the antenna as the elevation angles increase.

We provide two sets of Matlab codes to facilitate mapping the reflection zones. The first, mapview_fresnel_toolbox.m, plots the Fresnel zones for a given GNSS site in a plain “horizontal” map view. Figure 1 shows example L1 Fresnel zones for a GNSS station located in Boulder, Colorado. Two values of HR are used to demonstrate the fundamental relationship to this quantity. The second code, googleEarthFresnel.m, provides similar information, with an output kml file that can be loaded into Google Earth. Figure 2 shows a screen view from Google Earth for the Boulder, Colorado example. Figure 3 shows a screen view for a GNSS site in Alaska. This particular GNSS site is far above sea level, about 68 m, and the input HR value used reflects that.

Fig. 1
figure 1

First Fresnel zones in mapview for GNSS site P041 near Boulder, Colorado. Elevation angles (°) are green (5), cyan (10), blue (15), magenta (20), and red (25). Two values of HR are used: 2 m (top) and 10 m (bottom)

Full size image
Fig. 2
figure 2

Screenshot of first Fresnel zones for GNSS site P041. For clarity, only the reflection zones for elevation angles of 5 and 10° are shown on the Google Earth image. An HR value of 2 m was used

Full size image
Fig. 3
figure 3

Screenshot of first Fresnel zones for GNSS site AC12 and elevation angles of 5°, 7°, and 10° projected on a Google Earth image. An HR value of 68 m was used

Full size image

Both mapview_fresnel_toolbox.m and googleEarthFresnel.m need to know the approximate azimuth of the rising and setting satellites for the GNSS station in question. If the user does not know this information, values can be computed using a separate piece of Matlab code: do_azims.m. The googleEarthFresnel.m code allows the user to either manually set HR or to use mean sea level as the reflecting surface. If desiring mean sea level, the code uses the ellipsoidal GNSS station height and the EGM96 geoid correction (Lemoine et al. 1998). For more information on Fresnel zones in GNSS-IR, the reader is directed to Roussel et al. (2014) and Nievinski et al. (2016).

Frequency extraction from GNSS-IR: theoretical discussion

Although a GNSS receiver will track signals at even time intervals, the interval between 2 samplings of sin(e(t)) will be uneven during any given observing window. Furthermore, different satellite tracks will have different sampling intervals, as GNSS satellites that stay low in the sky move more slowly than others that pass higher in the sky. Figure 4 shows how sample data intervals vary for a site in Greenland. We use the Lomb Scargle Periodogram (hereafter LSP), designed to detect periodic signals with unevenly spaced observations, to extract the SNR spectral content (Lomb 1976; Press et al. 1992). To express the SNR spectral frequencies directly in terms of H in meters, we scale the sampling variable sin(e(t)) by a wavelength factor. The code lomb.m computes the normalized periodogram of a sequence of SNR data sampled at X(t) = 2sin(e(t))/λ, during an arc over a given elevation angle range [emin emax].

Fig. 4
figure 4

Interval between two sin(e(t)) sampling variables for station GLS2 with a receiver sampling rate of 15 s. The markers are projected on the satellite ground track for a 2-m antenna height. A marker is plotted every minute because the markers are very close for higher satellite elevation angles

Full size image

A LSP requires the user to provide the highest frequency factor (hifac) and an oversampling factor (ofac). We have recast these inputs into variables that are more intuitive for a GNSS-IR analyst. Instead of hifac and ofac, the user chooses the maximum value of H to be calculated and the desired precision, both in meters. Thus, the Lomb normalized periodogram is calculated at the frequencies H = [0: desired Precision: maxHeight]. These two user inputs are then converted into hifac and ofac for input to lomb.m. More details can be found in the Matlab code get_ofac_hifac.m.

The spectral grid can be as large and fine as you want at the expense of longer computation time. A user may want to compute frequencies up to the pseudo-Nyquist limit (see the next section). For a 1-s GPS sampling rate the pseudo-Nyquist limit is ~ 450 m with N = 1800 samples collected in a 30 min arc. Using a grid precision of 1 cm, the number of frequencies would be Nf = 45,000. The scaling of the lomb.m algorithm, provided in this GPS Toolbox contribution, is of the order of Nf ↔N ~ = 108. For a 5-s GPS sampling rate, this value is reduced by 25 but is still quite large. Instead, for snow accumulation, we have prior knowledge of the frequency range from which to expect a spectral peak. For a GNSS antenna initially set at 2 m, we might select a maximum grid frequency maxHeight of 6 m. Note that faster implementations of the LSP exist if you would like to check all frequencies (Press and Rybicki 1989).

There is no guarantee that the highest peak corresponds to the best frequency HR. And even if you are on the correct peak, the grid spacing is not the precision of the spectral peak. For instance, the spectral resolution scales inversely to the length W of the GNSS type observing window, where W = 2(sin(emax)–sin(emin))/λ in units of inverse meters. It is advisable to have at least one cycle of SNR data in your survey window. Significant frequencies below 1/W need to be analyzed carefully as they could be low-frequency residuals from an imperfect removal of the direct SNR signal. Finally, the peak estimate from one LSP may not be statistically reliable, but by taking the daily average or median over tens to hundreds of tracks one can increase the quality of the HR estimate over a region.

For even sampling, any frequency above the Nyquist frequency will be folded back (aliased) into the lower frequencies. With uneven sampling, the Nyquist-like limit, the limit beyond which no further information from the spectral content of the sampled signal can be extracted, can be much larger than the “Average-Nyquist” frequency computed for the same number of data uniformly sampled during the same time span. This makes intuitive sense because when the spacing varies we collect extra information in the space between two even samples and this can remove the aliasing ambiguity (Press et al. 1992; VanderPlas 2017).

In the context of GNSS SNR data, the samples in the observing window are structured, i.e., they are not random. Any structure in the sampling interval will be reflected in the LSP of the signal. Qualitatively the LSP patterns can be represented as the convolution between the window spectrum and the spectrum of the true signal. To discuss a Nyquist-like limit, we will look at the spectral characteristics of the observing windows. We provide this discussion so that potential GNSS-users will have an understanding of what GNSS time sampling interval to use to appropriately resolve HR from their SNR data. Figure 5d shows L2C SNR observations from the GNSS station GLS2. Operators of the station set its receiver-sampling rate to 15 s and there was no elevation mask. Figure 5e shows three tracks with the same SNR characteristics but with different survey windows in the elevation angle range [5°–20°]. The window spectra are computed up to H of 300 m (To estimate the form of the window power, compute a LSP on a series of unit measurements, but turn off the automatic pre-centering of the data in lomb.m).

Fig. 5
figure 5

Effect of GNSS type survey window with uneven sampling on the LSP of SNR data. Three windows for L2C data are shown spanning the elevation range [5–20] from the GNSS station GLS2. The latter has a receiver sampling-rate of 15 s. Satellite 27’s southeast track (blue) has nearly uniform samplings. Satellite 30’s northeast track (green) has moderately uneven samplings. Satellite 27’s northeast track (red) has the most uneven samplings found at this site. a LSP of SNR data between 0 and 300 m. b LSP of sampling windows. c Zoom of LSP of SNR data between 0 and 6 m. d SNR data with the same characteristics, as a function of the sampling variable X(t) = 2sin(e(t))/λ in units of inverse meters. e Interval between observations dX as a function of X. Note that in this figure we use X as the variable because it is this sampling variable which expresses the spectral frequencies directly in meters

Full size image

As shown in Fig. 5e, the southeast track for satellite 27 has a nearly uniform sampling. This results in a power spectrum of the window close to regularly narrow spaced spikes with a cadence of 69 m (Fig. 5a). Yet as the frequency H increases, the shape of the peaks changes. The spikes slightly decrease in power and their tail end spreads. As the sampling becomes less uniform the distortion increases. For the most severe case of uneven sampling found at GLS2 (satellite 27 northeast track) the window spectrum becomes noise-like after 250 m. This window structure is reflected in the LSP of the SNR (red traces in Fig. 5a, b).

The distance (cadence) between two spikes stays relatively constant. We will use half of this distance as the pseudo-Nyquist frequency. Any single frequency above this limit will be an imperfect version folded back into the lower frequencies, with a spectrum spread distortion depending on the degree of unevenness in the samplings. For moderately uneven to strong uneven sampled noisy data, the interaction of high-frequency noise bands with the convolution window can be complex with an increased spectral background noise level and eventually create spurious high-leveled peaks. Figure 5c shows that SNR power spectrums for the three tracks in the range [0–6 m] are similar. This is because there is no significant high-frequency noise folded back into this region.

We should clarify that with this type of LSP window, when the power of the second spike is less than ~ half the power of the zero-frequency peak, one should be able to detect significant spectral peaks beyond this pseudo-Nyquist frequency, and the true Nyquist-like limit is much larger. For further qualitative details on the effect of uneven-sampling and understanding the LSP the reader is referred to VanderPlas 2017.

To conclude, the exact structure of the LSP window will vary with each site, track, elevation range, L-band frequency and GPS receiver sampling rate, but the window LSP signature will vary between the 2 extremes mentioned above. This is one reason for the quality of the LSP across sites and tracks, especially in the presence of noise. We can get an approximate order of magnitude pseudo-Nyquist frequency using the average-Nyquist frequency, which is around 30–50 m for GLS2, using a GNSS sampling rate of 15 s at the L2 frequency (Table 1). For a survey window of length W and N observations the average-Nyquist frequency is N/2W.

Table 1 The pseudo-Nyquist and average-Nyquist frequencies for the three survey windows from Fig. 5
Full size table

Matlab code median_avg_nyquist.m is provided to allow a user to calculate the median average-Nyquist frequency. The user provides a receiver sampling interval (in seconds), the station location, GPS frequency, and elevation angle limits. The code simulates rising and setting satellite information for that GNSS site and outputs the median average-Nyquist frequency in meters.

GNSS-IR examples

Finally, we provide Matlab code to compute GNSS-IR results (sample_gnss_ir.m). The GNSS stations are chosen from Greenland, Antarctica, Alaska, and the western U.S. The GNSS sites have one thing in common: they were deployed without any thought that they could be used for GNSS-IR. Although we refer to GNSS here, these particular sites only tracked GPS satellites. Before discussing the output of this code, we will first describe the various steps within it.

  1. 1.

    Because this is a set of tutorial codes, the user is asked to choose a GPS frequency. The L1 GPS signal generated with the C/A code, or L1C, has the advantage that it is almost always available in RINEX files and it is tracked for all satellites. There can be issues with quality, particularly for some receivers; some of these issues were discussed previously by Larson and Nievinski (2013). The L2 GPS GPS signal currently has two codes. The public code (L2C) is far superior to L2P for GNSS-IR. Unfortunately, it is not always tracked at GNSS sites. As of 2016, L2C is available on 18 GPS satellites. For the sample data files we provided, we strongly recommend using L1 or L2C (when available). None of the example GNSS sites in this tutorial tracked L5 signals, but as with L2C, these new data are excellent for GNSS-IR applications (Tabibi et al. 2015).

  2. 2.

    Various defaults have been set. As discussed in the previous section, two LSP parameters need to be defined. The LSP precision is set to 5 mm and the maximum HR allowed is 8 m. HR values smaller than 0.4 m are not allowed, as GNSS-IR breaks down for small values of HR (Nievinski and Larson 2014c). For tower applications, users must change the maximum allowed value of HR.

  3. 3.

    The code needs to know the rising and setting satellite arcs for each site. Rather than predefine these values, the tutorial codes search for all available satellite data within 45-degree azimuth bins. Once the observations for a particular satellite in a given azimuth bin are found, the SNR data are converted to linear units (from dB-Hz to volts/volts) and a low-order polynomial is removed. The latter represents the direct signal component which is of no interest for GNSS-IR. A LSP is then produced from these “flattened” SNR traces. The user can change the polynomial order by modifying the code directly.

  4. 4.

    After a LSP is computed for a given rising or setting arc, the code must decide whether the peak in the LSP is significant. Here we have used a simple peak/noise ratio test. This is certainly not the only- or the best-way to compute the significance of a peak. The code allows the user to define over which frequencies the noise metric is computed and what ratio between the peak and noise is required. Other quality control metrics that could be used include a simple amplitude minimum value, the number of points (which would depend on the sampling interval), and the elevation angle difference.

  5. 5.

    The tutorial code generates two kinds of output. The significant HR results are printed to a text file. Plots can be automatically generated either as a summary for all azimuths or in separate azimuth bins. This is currently set to 45-degree azimuth bins. However, this can be changed by the user.

  6. 6.

    Figures 6 and 7 give you an overview of the GNSS-IR steps used for GLS2, the sample GNSS site in Greenland. In contrast to the Fresnel zones in the Boulder, Colorado case, which had a void in the north, there is 360-degree azimuthal coverage at GLS2. This is typical for GNSS sites in polar regions. Figure 6 shows SNR data for satellite 30, which rose and set twice at GLS2. This means there are four satellite arcs, one in each geographic quadrant. Figure 7 shows the corresponding LSP for these SNR data. Note that at this site the quality of the oscillations is poorer at higher elevation angles. The user can restrict the analysis accordingly. The reflection signal at GLS2 below 10° is quite strong, which serves as a reminder to station operators that elevation masks hinder GNSS-IR applications.

  7. 7.

    sample_gnss_ir.m can be run with six sample files. Here we will only discuss three of the sites: P038, GLS2, and SG27. Figure 8 shows L2C GNSS-IR results for station P038. This station is located at an airport in New Mexico with little or no terrain relief. We used the “Separate plots by azimuth bin” option and show LSP output for azimuths between 180° and 225°. You can see that the peaks of the LSP are very consistent with a horizontal planar surface. Figure 8 also shows complete GNSS-IR results for a GLS2 data record. In this example, we opted to show all azimuths together and the code calculates the median of the peak values, yielding a HR of 2.95 m.

Fig. 6
figure 6

GPS L2C SNR data for GNSS site GLS2 and satellite 30. Satellite track orientations are given in the titles

Full size image
Fig. 7
figure 7

Lomb Scargle Periodograms for the SNR data shown in Fig. 6

Full size image
Fig. 8
figure 8

Top: Lomb Scargle Periodograms (LSP) for L2C SNR data from GNSS station P038 on January 1, 2018 in an azimuth bin 180–225°; Bottom: LSP for L1C SNR data at GNSS station GLS2 on May 24, 2013 for all azimuths

Full size image

Figure 9 shows sample LSP results for GNSS site SG27 in Barrow, Alaska. Only the southeast quadrant is shown, as this region is unobstructed and planar (Liu and Larson 2018). The L1C LSP results produce a median HR value of 3.60 m. Figure 9 also shows LSP results for L2P. Instead of a single significant peak, there appear to be strong peaks at both the L1 peak and at ~ 4.6 m. The second peak is located at 3.60 multiplied by the ratio of the L2 and L1 wavelengths (0.244 and 0.19 m, respectively). This second peak is to be expected given that geodetic receivers must cross correlate to extract the L2P data. While the secondary peak is straightforward to observe and exclude, we caution that the peaks will be much more difficult to separate at smaller values of HR. We have also provided options in sample_gnss_ir.m so that you can evaluate your own GNSS site and/or modify default assumptions.

Fig. 9
figure 9

Lomb Scargle Periodograms for SNR data from GNSS station SG27 on January 1, 2018 for the southeast quadrant. Top: L1C data; bottom: L2P data

Full size image

Our goal for distributing this software is to make it easier for you to visually understand what reflected GNSS signals look like in SNR data. You cannot immediately use the codes in an operational sense, but they can easily be modified for that use. The main change you will need is to convert sample_gnss_ir.m into a function. Most typically at this stage all plots would be turned off, and this new function would be called with, e.g., station name, year, day of year, frequencies, elevation angles limits, desired precision, maximum HR, and azimuth ranges. One reason we think it is useful to start out evaluating SNR data visually is so that you can see which azimuth and elevation angles are generating usable reflection data. The reflection zone mapping software we have provided earlier gives you a way to validate these azimuth and elevation angle choices. While it is certainly possible to automate the azimuth and elevation angle choices, as we did for PBO H2O (Larson 2016), you will be able to develop a better automation scheme if you start with the raw SNR data. Visually inspecting GNSS-IR periodograms will also encourage you to think about how to decide which LSP retrievals are significant and which are not. The Fortran and MATLAB code for this software, and the sample data, can be accessed via the GPS Toolbox website at http://www.ngs.noaa.gov/gps-toolbox.

Final remarks

We hope that this software will provide guidance for new users of GNSS-IR. The technique is particularly straightforward to use in the cryosphere (Shean et al. 2017; Siegfried et al. 2017). With excellent azimuthal coverage at the poles and large planar surfaces, several hundred GNSS-IR reflection retrievals can easily be made per day, yielding an extremely robust daily average. GNSS receivers/antennas on towers can also provide accurate measurements of the surface over a large spatial region. If the GNSS antenna pole is set in ice, GNSS-IR and traditional GNSS vertical measurements can be used to simultaneously constrain the density of the firn layer and snow accumulation (Larson et al. 2015). Liu and Larson (2018) recently demonstrated that GNSS-IR can also be used to measure the deformation of the transition zone in permafrost regions.

The GNSS-IR technique is increasingly being used as a tide gauge (Larson et al. 2017). The motivation for doing so is primarily its simultaneous ability to measure changes in the water surface and the antenna phase center in a terrestrial reference frame (Santamaría-Gómez and Watson 2017). For the tide gauge application, additional corrections are needed and codes for those corrections are not provided in this GPS Toolbox contribution. First, a correction is needed if the water height changes significantly during a rising or setting satellite arc (Larson et al. 2013), i.e., the so-called \({\dot {H}_{\text{R}}}\) term. Second, a refraction correction must be made. The reader is directed to Williams and Nievinski (2017) for more discussion of refraction in GNSS-IR and how to correct for it. Because the tide gauge application requires subdaily measurements, significant efforts have been made to improve the resolution of a single HR value. We direct the reader to Strandberg et al. (2016), Reinking (2016), and Wang et al. (2018) for additional information on these efforts.

As a final note, we encourage GNSS-IR enthusiasts to take advantage of GNSS-IR simulators such as Nievinski and Larson (2014d) and Roussel et al. (2014) and the mapping/Nyquist tools provided here to properly design new GNSS-IR sites. At a minimum, one should:

  1. 1.

    Track all frequencies.

  2. 2.

    Track all codes (i.e. L2C, L5).

  3. 3.

    Track all constellations.

  4. 4.

    Remove elevation angle masks.

If there are significant cost issues related to telemetry, one can choose a receiver sampling interval that both limits those costs and allows reflection monitoring. We also encourage those station operators that need to place their GNSS antennas on buildings to consider placing them near the edge of the roof so as to enable reflection science.