Open-loop GPS signal tracking at low elevation angles from a ground-based observation site

A one-year data set of ground-based GPS signal observations aiming at geometric elevation angles below +2◦ is analyzed. Within the “GLESER” measurement campaign about 2600 validated setting events were recorded by the “OpenGPS” open-loop tracking receiver at an observation site located at 52.3808◦N, 13.0642◦E between January and December 2014. The measurements confirm the feasibility of open-loop signal tracking down to geometric elevation angles of −1◦ to −1.5◦ extending the corresponding closed-loop tracking range by up to 1◦. The study is based on the premise that observations of 5 low-elevation events by a ground-based receiver may serve as test cases for space-based radio occultation measurements, even if the latter proceed at a significantly faster temporal scale. The results support the conclusion that the open-loop Doppler model has negligible influence on the derived carrier frequency profile for strong signal-to-noise density ratios above about 30 dB Hz. At lower signal levels, however, the “OpenGPS” receiver’s dual-channel design, which tracks the same signal using two Doppler models with a 10 Hz offset, uncovers a notable bias. The repeat patterns of the GPS orbit traces in terms of 10 azimuth angle reveal characteristic signatures in both, signal amplitude and Doppler frequency with respect to the topography close to the observation site. On the other hand, vertical refractivity gradients extracted from ECMWF meteorological fields correlate moderately well with observed signal amplitude fluctuations at negative geometric elevation angles emphasizing the information content of low-elevation GPS signals with respect to the atmospheric state within the planetary boundary layer.

G. Beyerle and F. Zus: Low-elevation GPS signal tracking ties of retrieving dual-frequency carrier phase data, when the ray tangent point enters the lower troposphere (Rocken et al., 1997a).More specifically, Rocken et al. (1997a) noticed a significant negative refractivity bias in the lower troposphere at tropical latitudes.At low altitudes the GNSS signals experience multipath beam propagation (see e.g.Gorbunov, 2002b;Hocke et al., 1999).The resulting optical path length differences lead to signal scintillations and these amplitude fluctuations increase the probability of an early loss of tracking lock.To address these issues new signal tracking methods were developed and implemented.Whereas the "fly wheeling" tracking method of JPL's "Blackjack" GPS receivers mounted on CHAMP and GRACE (Ao et al., 2003;Hajj et al., 2004) showed some progress, significant improvements with respect to probing of the planetary boundary layer (PBL), in particular at low latitudes, were obtained with the introduction and implementation of open-loop (O/L) tracking (or raw-sampling) techniques (see e.g.Sokolovskiy, 2001b;Sokolovskiy et al., 2006;Ao et al., 2009;Bonnedal et al., 2010).
The open-loop signal tracking mode successfully resolves the problem of premature loss of signal in the lower troposphere at low latitudes (Sokolovskiy, 2001b;Sokolovskiy et al., 2006).In contrast to closed-loop (C/L) tracking, a receiver operating in open-loop mode partially or completely disregards the tracking loop feedback values from the carrier and code discriminators, but instead steers the corresponding numerically controlled oscillators (NCOs) using a priori parameters.In the following, these O/L parameters, which are usually derived from an atmospheric climatology, are referred to as "O/L model".The time duration, which the receiver operates in open-loop tracking mode, may be controlled by predetermined threshold values in terms of tangent point altitude, elevation angle and/or signal-to-noise ratio (SNR).While these general considerations apply to all open-loop and raw-sampling implementations, the specific realisations vary in detail (see e.g.US patents US6731906 B2 and US6720916 B2).
A key requirement for the adoption of O/L signal tracking in operational GNSS-RO missions is the insensitivity of derived carrier phase paths and code pseudoranges on the particular choice of O/L model.According to Sokolovskiy (2001b) the requirement is met, provided the true atmospheric Doppler profile deviates by not more than half of the sampling frequency from the O/L Doppler model.With a typical sampling rate of f s = 50 Hz this requirement translates into a maximum frequency deviation of 25 Hz.
Recently, Beyerle et al. (2011) claimed, on the basis of GNSS-RO observations recorded by the "IGOR" receiver aboard the TerraSAR-X spacecraft, that the O/L Doppler model may influence the derived refractivity values for low signal amplitudes below about 25 V/V and potentially contribute to the negative refractivity bias.In order to substantiate this hypothesis they proposed to track each GNSS-RO signal with two O/L models, separated in frequency space by a predefined offset.We note, however, that other causes certainly contribute to the observed refractivity bias as well (see e.g.Gorbunov et al., 2015;Sokolovskiy et al., 2010, and references therein).Since IGOR firmware modification aboard TerraSAR is unfeasible, two indirect methods were investigated.First, a measurement campaign was conducted in late 2012 with the GNSS-RO receivers aboard the TerraSAR-X and TanDEM-X satellites recording signals from the same setting GPS satellites, but using different O/L models (Zus et al., 2014b).Second, within the framework of the GLESER (GPS low-elevation setting event recorder) measurement campaign a ground-based experiment was devised and established at an observation site on the "Albert Einstein" science campus in Potsdam, Germany (52.3808 • N, 13.0642 • E) in December 2010.The GLESER campaign targets signals from GPS satellites at low elevations as they set beyond the horizon at elevation angles below +2 • .The measurement hardware, the single-frequency "OpenGPS" receiver, is an in-house development based on C. Kelley's "OpenSourceGPS" concept (Kelley, 2002).
The occurrence of multiple signal paths ("multipath") connecting transmitter device and receiver instrument distinguishes ground-based from space-based GPS measurements.Local multipath, caused by signal reflections in the direct vicinity of the receiving antenna, cannot be avoided in most ground-based observations (see e.g.Parkinson and Spilker, 1996;Hofmann-Wellenhof, 1997); on a spaceborne platform, however, these local reflections may be eliminated to a large extent by careful spacecraft design and suitable antenna placement (cf.Gaylor et al., 2005).Within the framework of space-based radio occultation the term "multipath" assumes an alternative connotation and generally refers to tropospheric signal propagation close to the ray tangent point, thousands of kilometres away from the receiver.Local multipath, i.e. reflections in the vicinity surrounding the antenna, can be described using geometric optics (see e.g.Elósegui et al., 1995;Anderson, 2000;Larson et al., 2008); tropospheric multipath is a diffraction phenomenon and requires wave optical methods to analyse quantitatively (see e.g.Gorbunov, 2002a;Sokolovskiy, 2001b;Jensen et al., 2003).In the following, however, we will argue that for elevation angles below about +2 • wave optical effects may contribute in ground-based observations as well.
Whilst ground-based observations of low-elevation setting events do not allow to derive bending angle profiles for ray tangent points above the receiver altitude (Zuffada et al., 1999;Haase et al., 2014;Healy, 2002;Sokolovskiy et al., 2001), these measurements nevertheless are useful to investigate receiver tracking behaviour under multipath conditions with strongly fluctuating SNRs.In addition, the signal excess phase paths have been shown to be sensitive to the local refractivity field (see e.g.Lowry et al., 2002;Zus et al., 2015).
In the present study we focus on low-elevation setting events with geometric elevation angles decreasing from +2 • at the measurement start to about −1 to −1.5 • at the end of the observation.Here, geometric elevation angle is defined as the line-of-sight elevation angle at the receiver antenna disregarding atmospheric refraction effects.Setting events last for about 10 to 15 min on average; hence their durations are about an order of magnitude longer than typical spacebased radio occultation measurements (see, e.g.Kursinski et al., 1997).Even if a ground-based observation does not lend itself to the derivation of bending angle profiles (Zuffada et al., 1999;Haase et al., 2014;Healy, 2002;Sokolovskiy et al., 2001), from the signal tracking perspective we regard GLESER observations as useful for investigating openloop tracking of signals with strongly fluctuating SNR.The present study is restricted to the observation and analysis of setting events; an extension towards rising events, however, is technically feasible and may be considered for future implementations.
The paper is sectioned as follows.First, closed-loop and open-loop tracking methods are briefly reviewed and the capabilities of GFZ's OpenGPS receiver are illustrated with some example profiles.Second, the measurements conducted during the GLESER campaign are introduced, and the data processing algorithms and analysis methods are discussed.In the main section of this paper the measurement results are discussed and put into perspective using results from multiple phase screen (MPS) simulations.Finally, the OpenGPS hardware and software are described in the Appendix.

Closed-loop and open-loop signal tracking
The GLESER campaign utilises the OpenGPS instrument, a single-frequency 12-channel GPS receiver.Several copies of this device, which is based on C. Kelley's OpenSourceGPS concept (Kelley, 2002), were built at GFZ and used in various ground-based and airborne GPS measurement campaigns (see e.g.Helm, 2008;Helm et al., 2004).In order to provide a self-consistent description of the OpenGPS instrument and its O/L signal tracking implementation, we begin with a brief review of C/L and O/L tracking techniques.
It is well known that inhomogeneities in the tropospheric water vapour field, in particular at low latitudes, can produce multipath propagation of GNSS signals at low elevation angles (see e.g.Gorbunov et al., 2004;Beyerle et al., 2003).Space-based RO observations show that under these conditions SNR values exhibit strong fluctuations, which early GNSS-RO receivers were unable to track properly (Rocken et al., 1997a).To address premature signal loss GNSS-RO receiver tracking algorithms based on closed feedback loops were replaced by "open-loop" techniques (see e.g.Sokolovskiy, 2001b;Sokolovskiy et al., 2006;Ao et al., 2009;Bonnedal et al., 2010;Wang et al., 2016).If a receiver operates in open-loop mode the feedback loop is opened and the NCO producing the replica signal is steered (in part or fully) from model parameters.This model takes into account the expected signal dynamics from both, the transmitter and In closed-loop mode the input signal is correlated with two replica signals, sin and cos, and the latter is phase-shifted by 90 • with respect to the former.The correlation sums are low-pass filtered and the output is examined for phase deviations between input and replica signal.The discriminator adjustments close the feedback loop.Open-loop tracking (bottom panel) dispenses with the phase discriminator and the numerically controlled oscillator (NCO) is solely steered from model values (Doppler model).The observed carrier phase finally is assembled from the NCO phases and the in-phase and quad-phase correlation samples.receiver orbits, clocks biases and drifts as well as the signal propagation characteristics in the lower troposphere.The latter are typically obtained from an atmospheric climatology (Sokolovskiy, 2001b;Bonnedal et al., 2010).
The schematic in Fig. 1 illustrates the two tracking concepts for carrier phase tracking; corresponding considerations apply to code tracking as well.In standard C/L tracking (Fig. 1, top panel) the down-converted input signal ("input") is correlated with two internal replica signals ("sin" and "cos") generated by the NCO (see e.g.Misra and Enge, 2006).The result is low-pass filtered (represented by the box labelled "average"); the carrier phase discriminator determines phase deviations between the observed and modelled replica and provides appropriate adjustments to the loop filter.The sin and cos replica signals, the latter being phase shifted by a quarter cycle, i.e. 90 • with respect to the former, allow us to distinguish phase advances from phase delays between observed and replica signal.The phase discriminator output is digitally filtered ("loop filter") to prevent unstable loop behaviour (see e.g.Lindsey and Chie, 1981;Thomas, 1989).The receiver output samples ("carrier phase output") combine the NCO model phases and the phase residuals from the discriminator to yield the observed carrier phase.
If signal amplitudes drop below certain threshold levels, the corresponding phase residuals start to be dominated by noise and proper alignment between replica and observed signal can no longer be maintained.Figure 2 shows exemplarily the signal-to-noise density ratio C / N 0 (Badke, 2009;Kaplan, 1996;Parkinson and Spilker, 1996) 3).Correspondingly, C / N 0 drops by about 20 dB down to the noise level (Fig. 2) and never recovers during the last part of this setting event.
Open-loop tracking is immune to transient SNR gaps, even if these breaks stretch across extended time periods.The bottom panel of Fig. 1 schematically illustrates the concept.In O/L tracking mode the feedback loop is removed and the NCO is solely controlled by model values.The correlation output values, produced by the cos and sin branches, are denoted by "in-phase" and "quad-phase" correlation samples, respectively (Sokolovskiy, 2001b;Sokolovskiy et al., 2006;  Ao et al., 2009;Bonnedal et al., 2010).In low-SNR conditions the in-phase and quad-phase samples are dominated by noise.However, since no feedback is present in O/L tracking, these noisy samples cannot produce erroneous control input to the NCO and transient signal gaps do not cause loss of tracking lock as illustrated in Figs. 2 and 3. Red and green lines show C / N 0 (Fig. 2) and the NCO carrier frequency (Fig. 3) recorded by the receiver's two open-loop channels, respectively.In the following the two channels are denoted by the letters A and B. Both channels are attached to the same PRN and therefore track the same transmitter.Their corresponding NCO frequencies are derived from the same model, but on channel A's value an additional −5 Hz, and on B an additional +5 Hz, is added, resulting in a 10 Hz interchannel offset.This 10 Hz shift can clearly be identified in Fig. 3 (insert); here, O/L tracking mode starts at an elevation angle of −0.08 • and reaches the nominal 10 Hz shift after a brief initialization phase at −0.13 • elevation.

Measurements
The GLESER measurement campaign started in December 2010 and since then data acquisition of low-elevation setting events operates almost continuously.Gaps of up to several days occurred and are caused by hardware or software problems, operator errors or other technical reasons.The following discussion is restricted to observations recorded in 2014.The instrument provided on average 8.3 observations per day throughout this year with three data gaps (29 January to 1 February, 29-31 August and 18-22 December 2014).Between 15 July and 6 September the OpenGPS receiver malfunctioned due to an operator error and 437 observations from that time period are removed from the data set, leaving 2581 low-elevation events.The instrument is housed on the top floor of the former water tower on the "Albert Einstein" science campus in Potsdam, Germany (52.3808 • N, 13.0642 • E).An active L1(1575.42MHz) GPS antenna (NovAtel/ANTCOM 2G15A-XTB) is mounted on the tower's observation deck and tilted by about 45 • towards the western horizon to increase recorded signal strength at very low elevation angles.The antenna includes a low-noise amplifier with 33 dB gain and is sensitive to right hand circular polarisation only.A 10 m low-loss signal cable connects the antenna to the OpenGPS receiver.According to the manufacturer's documentation the cable loss is 14 dB per 100 m at 1 GHz.
Figure 4 depicts the panoramic view of the western horizon from the observation deck between azimuth angles of 180 (south) and 360 • (north).The azimuth angles are marked in black at the top side of the photograph.Thick red lines indicate the azimuth angles of those PRNs, which are analysed in the present study.The picture was taken a few metres east of the receiving patch antenna, which is visible to the right of the 290 • azimuth tick mark.
Data acquisition is activated once a satellite enters the observation window which ranges from 200 to 340 • in azimuth (thin broken red lines in Fig. 4) and from +2 to −2 • in elevation.Any additional satellite entering the observation window during this time period is disregarded and not tracked in O/L mode.Raw observables stored on disk include pseudoranges ρ C/L,n , carrier phases ϕ C/L,n and time tags t n for all C/L channels.In addition, raw O/L measurements (inphase and quad-phase correlation sums I n and Q n ; NCO ranges ρ NCO,n and NCO carrier frequencies f NCO,n from both O/L channels) are written to separate files.Here, the subscript n = 1, 2, . . .enumerates the successive TIC events (for an explanation of TIC event, see Appendix A1) recurring every 20 ms.The correlation sums I n and Q n are not sampled at TIC instants but are coherently integrated over a time period of 20 ms preceding the TIC event.The corresponding time offset between I n , Q n and ϕ NCO,n is disregarded in the following analysis.Raw data accumulation rate is about 360 MB day −1 or about 2.5 GB week −1 .
In order to quantify the performance gain of O/L in comparison to C/L detection we compare in Fig. 5 the lowest elevation angles observed in the two tracking modes, O/L off and C/L off .They are defined as the minimum elevation in a given setting event with the smoothed density ratio C / N 0 still exceeding 30 dB Hz.Smoothing is performed by a running mean filter of 1 s width.We find that out of a total of 2581 setting events 2368 and 2366 were recorded by the O/L channel A and B, respectively.In the remaining 8 % of the observations closed-loop tracking stopped too early and O/L was never activated.With active channels A and B, however, in about 86 % of the O/L measurements (2017 out of 2368 and 2069 out of 2366 for channel A and B, respectively) end elevation angles from O/L tracking yielded smaller values compared to the C/L results.Additionally, in 21 % of the measurements (499 out of 2368 and 505 out of 2366 for channel A and B, respectively) O/L tracking extended more than 0.25 • further down than the corresponding C/L observation.

Data processing and analysis
Data processing is performed on an event-by-event basis.Typically, setting events last about 700-900 s; however, measurements extending over a time period of up to 1400 s are occasionally observed.During the initial stage OpenGPS, raw data files are converted to MATLAB ® binary format to facilitate the subsequent processing steps.First, the observed in-phase and quad-phase samples sums are demodulated: (1) This done to allow calculation of the residual phase samples: Here, the four-quadrant arctangent atan2(x, y) denotes the principal value of the angle of the complex number x + i y.
The literature discusses two methods for the determination of the data bits D n , "internal" and "external" demodulation.On the one hand, D n can be extracted from the observations I n and Q n themselves (Sokolovskiy et al., 2009).On the other hand, external demodulation extracts D n from independent observations, such as GFZ's NavBit data base, registered under doi:10.1594/GFZ.ISDC.GNSS/GNSS-GPS-1-NAVBIT (Beyerle et al., 2009).In the following, external demodulation is used since in low-elevation events the modulus of the difference between adjacent the carrier phase residuals, ϕ res,n − ϕ res,n−1 , frequently reaches and exceeds ±90 • .In these cases a clear separation between propagation-induced phase fluctuations and phase changes due to a sign change of D n is difficult to achieve.
Accumulated residual carrier phase samples res,n then follow from with the unwrapping term C n defined by if n > 1 and C n=1 = 0 (Beyerle et al., 2011).In addition, during this first processing stage the receiver's 1 Hz navigation solution, which includes an estimate of the receiver's clock bias, is retrieved as well.
Second, elevation and azimuth angles for each tracked satellite are calculated from broadcast ephemeris data using www.atmos-meas-tech.net/10/15/2017/Atmos.Meas.Tech., 10, 15-34, 2017 the bias-corrected receiver clock time and linearly interpolated from 1 to 50 Hz.From the bit-corrected phase samples res,n the observed carrier frequencies are derived with superscript A and B, indicating the corresponding O/L channel.Figure 6 shows f (A,B) obs,n (Eq.5) corresponding to the event plotted in Fig. 3  the C/L frequency briefly crosses the displayed frequency range again.Below about −1.2 • elevation Earth's limb completely shadows the signal from PRN 7 (see Fig. 2) and from there on O/L and C/L outputs contain noise only.
In addition, with the observed in-phase and quad-phase correlation sums, I n and Q n (Eq.1), the signal-to-noise density ratio C / N 0 , in units of dB Hz, is calculated according to Badke (2009), Kaplan (1996) and Parkinson and Spilker (1996): Here, x is the mean value of x and T c = 0.02 s denotes the coherent integration time.Var Q C/L is the variance of the observed quad-phase correlation sums of the associated mas-ter channel running in C/L mode between 0 and +2 • elevation angle.
If no GPS signal is present and the receiver input signal is pure noise, I and Q are uncorrelated and since Q = 0.Under these conditions the density ratio C / N 0 (dB Hz) decreases to a value of 10 • log 1 T c ≈ 17 dB Hz (8) using Eq.(6) (dashed black line in Fig. 2).

Simulations
In order to interpret the observed C / N 0 fluctuations we performed a series of MPS simulations (Knepp, 1983;Martin and Flatté, 1988;Grimault, 1998).The propagation of a plane wave through the lower troposphere is modelled by a series of 500 non-equidistant phase screens covering the range from the receiver to a distance of 500 km.On each phase screen the wave suffers a phase delay determined by the inter-screen distance z and the refractivity height profile (Sokolovskiy, 2001a) with height h, scale height h s = 8 km, PBL top height h tp and PBL top transition zone h zn = 50 m.The interaction of the wave with the ground surface is modelled by applying a raised-cosine filter with a 6 dB steepness of 25 m (i.e.within 25 m the filter weight decreases by 6 dB) at zero altitude.
The phase screens extend vertically from −20 to +20 km with a 5 km wide raised-cosine filter applied at the upper and lower boundary to suppress spurious diffraction effects; the receiver altitude is taken to be 50 m.The variation of elevation angle between −2 and +2 • is modelled by tilting the ground surface and its overlying atmosphere correspondingly.
Results from four simulation runs are shown in Fig. 7; it displays the normalised signal amplitude as a function of elevation angle.Signal absorption at the ground surface produces characteristic diffraction patterns for elevation angles below 0 • (red and blue lines).Without ground absorption the diffraction patterns almost disappear and the profiles approximate step functions expected from geometric optics (green and black).The MPS simulations did not produce C / N 0 fluctuations for horizontally oriented PBL tops (parameterised by the parameter h tp in Eq. 9).However, if the top layer tilts towards the receiver, substantial signal deviations at elevation angles above 0 • are observed.Figure 7 illustrates this phenomenon for a PBL top layer ascending Signal absorption at the surface is taken into account (red and blue lines); green and black lines show the result without ground absorption.Two boundary layers are modelled: a horizontal boundary layer top at 2 km altitude (red and black) and a layer top increasing from 1 to 2 km between 30 and about 60 km distance from the receiver (blue and green).For legibility the red, green and blue lines are shifted by an additional +1, +2 and +3 dB offset, respectively.For details see text.
from 1 km at 30 to 2 km at about 60 km distance (green and blue lines); at distances below 30 and above 60 km h tp remains fixed at 1 and 2 km, respectively.The MPS simulation results, plotted in Fig. 7, highlight ground diffraction effects below about 0 • elevation angle (blue and red) and PBL-induced C / N 0 variations above about 0 • (blue and green).The results suggest that these C / N 0 fluctuations are independent from each other and tend to be separated in elevation angle space.Finally, we note that the addition of irregularities on spatial scales characteristic for turbulence to the refractivity profiles did not produce significant C / N 0 changes.

Discussion and interpretation
We begin the discussion of the GLESER results by comparing the performance of internal versus external data bit demodulation and define the quantity as Here, |x| denotes the modulus of x, the superscript X = A or X = B indicates the O/L channel and "int" or "ext" characterises the demodulation method.The quantity E X n is sensitive to differences between sign changes from one sample to the next, rather than differences between the bit values, since only the former affect the derived carrier frequency.We note that E X n = 1 for randomly chosen D X,int n and D ext n .For 242 observations in January 2014 with signal density ratios C / N 0 ≥ 30 dB Hz the parameters E n (Eq.10) are grouped into an elevation grid with 0.1 • bin size and averaged.The result is shown in Fig. 8.At elevation angles just below 0 • there is good agreement between internal and external data bit retrieval; at lower elevations, however, the failure rate increases significantly.Results for elevation angles below −1.2 • are statistically not significant due to the strongly decreasing number of observations (red line).To eliminate potential errors caused by internal data bit removal, we restrict the following discussion to data processed using external demodulation.
Table 1 lists all 19 PRNs recorded during the 311-day observation period in 2014 and the corresponding number of setting events.The third column gives the mean azimuth angle and its 1σ standard deviation at 0 • elevation.Enhanced changes in azimuth angle of PRN 8 are most likely related to the decommissioning of GPS space vehicle number 38 in October 2014.The last column in Table 1 shows the fractions of profiles in which O/L tracking reached a lower elevation angle than the corresponding C/L result.For example, a value of 70 % indicates that in 30 out of 100 observations the C/L channel tracked to a lower elevation than the corresponding O/L channels.The low O/L enhancement values for PRN 28 might be caused by signal reflections at the water surface of Templiner Lake (see Fig. 4 at about 260 • azimuth angle).As described in the Appendix (Sect.A2), the O/L model is initialised within the elevation angles range between 0 and +2 • .It appears feasible that surface reflections induce C / N 0 fluctuations at these elevation angles (see e.g.Anderson, 2000;Larson et al., 2008), degrading the quality of the O/L model and thereby causing poor O/L performance.The panorama photograph Fig. 4 shows the view to the local horizon within the 140 • wide observation window ranging from 200 to 340 • and centred at 270 • azimuth (west).It exhibits a substantial variation of surface properties in azimuth angle with enhanced building densities of the city of Potsdam in a westerly to north-westerly direction (about 280 to 350 • ), mostly forests at azimuth angles below 270 • and water surfaces from Templiner Lake between about 250 and 260 • .Also indicated in Fig. 4 are the approximate zero-elevation azimuth angles of setting GPS satellites (red marks).Due to the specific geometry of the GPS constellation the orbit traces exhibit a characteristic azimuth angle repeat pattern when viewed from a fixed ground-based location (Axelrad et al., 2005;Choi et al., 2004), i.e. a specific PRN sets at the same azimuth angle with small deviations of less than 1 • (see centre column of Table 1).
In the following the observations are grouped PRN-wise.The azimuth angle repeat pattern implies that these subsets refer to almost the same azimuth angles.For example, settings of PRN 18 invariably take places at lake Templiner Lake, whereas PRN 32 sets across the urbanised area of Pots- dam throughout the year.The non-uniform topography most likely contributes to the observed variability in C / N 0 (see e.g.Fig. 2) and, as will be discussed below, Doppler frequency with respect to azimuth angle.The occurrence time of the setting event, however, shifts by about 4 min day −1 , which corresponds to 24 h per average year (365.25 days).Thus, the statistical analysis performed in the present study essentially averages out any potential diurnal variation.The investigation of these temporal variations is left to future research.
Figure 9 provides a general overview of the observed mean signal density ratios as a function of elevation angle.We note that by expressing C / N 0 in units of dB Hz, Eq. ( 11) constitutes in effect a geometric mean value in terms of the ratios Eq. 6). Figure 9 shows C / N 0 for PRNs 2, 7, 13, 14, 17, 18, 22, 23 and 32.The following analysis is focused on this set of nine PRNs; they were selected according to data availability and azimuth angle coverage.The individual panels in Fig. 9 and the following figures are arranged row-wise according to increasing azimuth angle (see Fig. 4).In Fig. 9  Examining density ratio averages (Eq.11) in lieu of C / N A 0 and C / N B 0 is justified, since the two values agree in the majority of observations.As a matter of principle, however, significant deviations are conceivable, because the density ratio depends on frequency offset f (A,B) between the O/L NCO frequency f (A,B) NCO and the signal's true carrier frequency f 0 , Provided the observed and replica signal are perfectly aligned in the pseudorange domain, the amplitude loss induced by f (A,B) is given by with a coherent integration time of T c .For illustration Fig. 10 shows the normalised C/A code correlation function as a function of code lag and frequency offset in the vicinity of the correlation maximum and the loss function L( f (A,B) ) at zero delay (red line).Thus, in the absence of other factors affecting C / N A 0 and C / N B 0 individually, the density ratios of channel A and B will differ by The statistics of C / N 0 is plotted in Fig. 11  The observed density ratio difference C / N 0 can be analysed quantitatively.For this purpose the observation f (A,B) obs ≈ f 0 is assumed to be a valid approximation for the true carrier frequency f 0 .Using the defining Eq. ( 5) we obtain Figure 12 shows the correlation between the observed density ratio differences C / N 0 and f (A) res (gray data points).Here, the data set is restricted to the subset of typically 210 000 (PRN 2) to 370 000 (PRN 7) samples tracked in O/L mode.The expected result derived from Eq. ( 16) is overlaid in dark blue; dashed lines mark the residual frequency of −5 Hz.
The resulting patterns exhibit a marked dependence on PRN, i.e. on azimuth angle.Whereas the observations of PRN 14 (azimuth angle at about 226 • ) indicate only moderate excursions in terms of residual frequency with most data points clustered at −5 Hz (dashed line), the signals from PRN 32, arriving from an azimuth angle of 318 • , show strong deviations causing residual frequencies exceeding the Nyquist value of f s /2 = 25 Hz.Here the agreement with the theoretical results (blue line, Eq. 16) is evident.A substantial number of observations deviate from the O/L model by more than f s /2 and are therefore affected by aliasing.The light blue curves, which are derived from Eq. ( 16) but shifted by ±50 Hz, show good agreement with the observations and substantiate this interpretation.The aliasing effect is stronger for negative f (A) res since channel A is shifted with respect to the O/L model by an additional −5 Hz.
Correspondingly, aliasing for positive residual frequencies is stronger for channel B since this channel is shifted with respect to the O/L model by +5 Hz as illustrated by Fig. 13, which shows the corresponding result derived O/L channel B data.Again, the dark blue lines (and the aliased curves in light blue) mark the theoretical result derived from Eq. ( 17).
The preceding discussion is based on the assumption 13. Same as Fig. 12, however, showing the corresponding results from O/L channel B.
Figure 14.Same as Fig. 9, however, showing the difference between the two observed frequencies obtained from O/L channel A and B as a function of the mean signal-to-noise density ratio (Eq.11).Mean and 1σ standard deviations, calculated from C / N 0 bins 2.5 dB Hz wide, are marked in green.The fraction of data points exceeding f obs > +40 Hz is indicated as ρ 40 Hz .The result of the statistical analysis excluding this subset still exhibits a positive bias if C / N 0 30 dB Hz (red).
Referring to Fig. 14, however, we deduce that Eq. ( 18) is not strictly fulfilled, in particular for low signal-to-noise density ratios.The panels of Fig. 14 show the observed frequency difference obs as a function of mean density ratio for the selected nine PRNs.Their mean values and 1σ standard deviations, sorted into 2.5 dB Hz bins, are superimposed on the individual data points (gray dots) as green lines.
Table 2. Mean and standard deviations of f obs (centre column) for nine values of the carrier signal-to-noise density ratio C / N 0 between 10 and 50 dB Hz; averaging bin size is 5 dB Hz.The statistics is based on all PRNs shown in Fig. 14 In all panels a small fraction of the data set, denoted by ρ 40 Hz in Fig. 14, populates the frequency band f obs + 40 Hz with a mean value of about +50 Hz and signal levels ranging from 10-20 dB Hz to more than 40 dB Hz.Again, this +50 Hz offset is caused by aliasing.If the true signal frequency f 0 occurs within the frequency range between +20 and +25 Hz, it is correctly tracked by channel A but aliased to the −30 to −25 Hz frequency window by channel B, since channel B's NCO is shifted by −10 Hz with respect to channel A's NCO and therefore its frequency range extends from −30 to +20 Hz.Thus, the frequencies observed by channel A and B differ by +50 Hz.Conversely, signals appearing at frequencies between −20 and −25 Hz are properly recorded by channel B but suffer aliasing in channel A again, producing a +50 Hz offset in f obs .Formally, the observed frequency difference f obs as a function of true frequency f 0 is We note that in both cases the offset is +50 Hz and therefore signals differing from f 0 by integer multiples of 50 Hz cannot be distinguished using dual-channel O/L tracking.
Even though the fraction ρ 40 Hz remains below 4 %, the corresponding samples bias the mean value (green lines in Fig. 14) towards positive frequencies, but they are not the sole cause for the observed bias.If all observations with f obs ≥ +40 Hz are removed from the data set, the corresponding mean values are still biased, albeit significantly less (red lines in Fig. 14).
Table 2 summarises the results shown in Fig. 14.The centre column lists the mean and 1σ frequency differences averaged over all nine PRNs displayed in Fig. 14.The corresponding results for the frequency differences excluding the fraction ρ 40 Hz are given in the last column.
The magnitude the observed frequency bias can be motivated in the following way.In the limit of vanishing signal levels the residual frequencies are dominated by noise and thus Under these circumstances the observed frequency difference between channel A and B becomes lim The mean values observed in Fig. 14 are consistent with this estimate.However, the figure also shows that for C / N 0 30 dB Hz the frequency difference f obs is bias free; the deviations occur solely at signal levels C / N 0 30 dB Hz.We note that this bias is independent from the sampling frequency f s .The numerical values of C / N 0 , however, which characterise the transition zone between biased and bias-free samples, depend on the antenna gain and other receiver-specific parameters.They cannot be compared to observations from space-based GNSS-RO payloads which typically are equipped with higher-gain antennas.
Apart from clusters appearing at +50 Hz, Fig. 14 indicates the presence of frequency offsets at about +25 Hz for PRN 7, 13, 17 and 23.A convincing cause for the presence of these +25 Hz clusters could not be identified; most likely they are related to atmospheric effects on the propagating signal and not receiver induced, since the effect strongly depends on PRN, i.e. azimuth angle.This issue requires further investigations.
The occurrence of strong SNR fluctuations during the last phase of most setting events (see Fig. 2) independent of azimuth angle suggests propagation-induced causes in addition to topographic, i.e. surface interaction effects.Hypothetically we relate these fluctuations to multipath ray propagation within the PBL.
The MPS simulations (Fig. 7) suggest that at negative elevation angles diffraction effects caused by the ground surface dominate the observed C / N 0 fluctuations; at higher elevations atmospheric multipath seems to be more relevant.This hypothesis is tested for the 4.5-month time period from March to mid-July 2014 by correlating the standard deviation of C / N 0 between elevation angles of +1 and +2 • with the mean refractivity gradient dN / dz .The calculation of dN / dz is restricted to the altitude range from 1 to 3 km.Figure 15 shows the results for nine PRNs.
The vertical refractivity profiles N (z) are extracted from European Centre for Medium-Range Weather Forecasts (ECMWF) meteorological fields.Their horizontal resolution is 1 • ×1 • (about 110 km in meridional and 69 km in zonal direction at the receiver location) with 137 height levels ranging from 0 to about 80 km; the averaging interval of 1 to 3 km   E) is selected, which is located about 99.8 km in the north-western direction (314.2 • true bearing).The standard deviation of the carrier signal-to-noise density ratio, σ (C / N 0 ), calculated within the elevation angle range +1 • < < +2 • , is taken as proxy for the signal amplitude fluctuation.Each panel of Fig. 15 includes information on the correlation; the Pearson and Spearman coefficients are quoted in the top and bottom line, respectively (see e.g.Press et al., 1992), the corresponding significance parameters are given in brackets.The values indicate that dN / dz and the standard deviation of C / N 0 are weakly to moderately correlated.With the exception of PRN 17 (top right panel) all calculated correlations are significant on the 5 % level.The (negative) correlations range from −0.17 to −0.40.We note that ECMWF refractivity profiles below 1 km frequently exhibit strong gradients.Their inclusion into the calculation of dN / dz significantly decreases the correlations or even renders them insignificant.The MPS simulations (see Fig. 7) also suggest that the (negative) correlation between σ (C / N 0 ) and dN / dz weakens if elevation angles close to or below the horizon are included.Figure 16 confirms this prediction.In contrast to Fig. 15 in this figure the elevation angle range, used for the calculation of σ (C / N 0 ), is extended downwards to −2 • .With the exception of both correlation coefficients for PRN 7 and the Pearson coefficient for PRN 18, the derived correlations are no longer significant; i.e. based on these results the null hypothesis cannot be rejected on the 5 % significance level.
An investigation of the relationship between signal amplitude and frequency fluctuations, the local atmospheric refractivity field and, potentially, effects arising from surface topography is beyond the scope of this paper and will be addressed with higher resolution meteorological data in future studies.

Conclusions
For more than a decade the OpenGPS receiver is used at GFZ in several ground-based and airborne measurement campaigns.Owing to its open hardware and software architecture the device can be adapted to address specific signal tracking issues in GNSS radio occultation, reflectometry, scatterometry or related fields.Here, a subset of low-elevation events recorded during the long-term GLESER campaign are introduced and discussed.Between 1 January and 31 December 2014 the instrument recorded 2581 validated setting events at an observation site located at 52.3808 • N, The analysis of open-loop data is performed on demodulated in-phase and quad-phase correlation samples.The present study suggests that navigation message demodulation using external information, e.g.extracted from GFZ's NavBit data base (doi:10.1594/GFZ.ISDC.GNSS/GNSS-GPS-1-NAVBIT), is preferable to internal demodulation.Open-loop signal tracking results are insensitive to the receiver-internal Doppler model for carrier signal-to-noise density ratios C / N 0 30 dB Hz; below this value reconstructed Doppler frequencies gravitate towards the model value and thus potentially constitute a bias source.The present study did not address potential contributions of the ionospheric signal propagation and/or local multipath to the observed C / N 0 fluctuations.These issues need to be addressed in future work preferably using dual-frequency receivers.
MPS simulations suggest that the observed weak to moderate correlations between mean vertical refractivity gradients within the mean PBL and the C / N 0 fluctuations at positive elevation angles are related to multipath induced by tilted layers of refractivity gradients.At lower elevation angles, below 0 • , signal diffraction at the ground surface dominates the observed amplitude variations.

Data availability
The GLESER campaign raw data files have been supplied with the digital object identifier doi:10.running in parallel, the kernel module ogrcvr_mod.ko and the user space application ogrcvr.The module ogrcvr_mod.koexecutes tasks controlling the code and carrier phase tracking loops.To maintain tracking lock the carrier phase and code data have to be read from the GP2021 registers, and the necessary loop adjustments have to be determined and written back to the correlator within a time interval of less than 100-200 µs.The real-time extension layer RTAI guarantees this latency performance even for high disk reading and writing operations or network traffic.ogrcvr handles deferrable tasks, such as determination of satellite positions and velocities from ephemeris data, calculation of the navigation solution and data storage.Communication between real-time module and user space process is accomplished using shared memory and FIFO (first in, first out) buffers.Inspection and modifications of relevant signal tracking parameters, such as loop bandwidths and loop orders, O/L channel frequency and code offsets and sampling frequencies, are performed via a proc-based command line interface.
A third process, ogdspl, can optionally be started to display tracking and positioning information on the terminal screen (Fig. A2).In total, the source code of all kernel and user space modules comprises about 20 000 lines of C code.
The OpenGPS receiver operates in two different observation modes.In "monitor mode" up to 10 of the 12 correlator channels are assigned to PRNs of visible GPS satellites.The O/L channels A and B, in Fig. A2 listed as number 11 and 12, remain unassigned.When a satellite, tracked by channel number k, crosses the +2 • elevation boundary from above, the receiver transitions to "measurement mode", channels A and B are assigned to this satellite's PRN and their code and carrier NCOs are aligned to those of channel k as well.In the following this process is called "cloning" of A and B from the "master" channel k and A/B are referred to as "clone" channels.
In measurement mode code and carrier NCO data, inphase and quad-phase correlation sums from both GP2021 correlation branches ("prompt" and "dither") and position and clock bias results from the real-time navigation solution are stored on the local hard disk.The delay between the prompt and dither branch is fixed at 0.5 chips (about 150 m).Code, carrier phases and correlation sums are written with a temporal resolution of T s = 0.02 ms (corresponding to a sampling frequency of f s = 50 Hz); the navigation solution is provided once per second.We note that the stored correlation sums are coherently integrated over T s , whereas the NCO code and carrier phase are instantaneous values sampled at the corresponding TIC event.
Typically, measurement mode lasts for about 10-15 min and ends when the satellite's elevation angle drops below −2 • .During an initialization phase, for elevation angles above 0 • , the code and carrier phase loop adjustments are collected at each TIC event and stored.Thereafter, when elevation angles are below 0 • , the carrier loop is opened and NCO input values are calculated by linear extrapolation of the stored Doppler adjustments.For code O/L tracking a hybrid method is employed; if C / N 0 ≤ 30 dB Hz, the code NCO values are calculated from a second-order polynomial extrapolation of the stored adjustments.Otherwise, the code loop operates in C/L mode; in addition, the difference between actual adjustments and model-derived values is saved and added to the model.Thus, potential deviations between observations at C / N 0 > 30 dB Hz and the extrapolated loop adjustments are used to update the model.
The OpenGPS receiver software determines geometric elevation angles from GPS almanac information.Since the data post-processing analysis is based on the more precise GPS ephemeris data, the elevation angles of the measurement start and the end of the O/L initialization may deviate from the nominal values of +2 and 0 • , respectively, by some tenths of a degree (see e.g.insert in Fig. 3).
Closer inspection of the O/L carrier NCO frequencies in Fig. 3 (red and green lines; the insert shows a zoomed-in view) reveals residual deviation of the O/L channel from the target Doppler model in addition to the nominal values of ±5 Hz.First, at the start of O/L tracking mode (at about −0.08 • elevation angle) the two clone channels are gradually moved into phase alignment to the O/L model by adjusting the corresponding NCO frequencies (OpenGPS's hardware correlator GP2021 does not support carrier phase adjustments).For channel B (+5 Hz) the duration of this initialization process is shorter, since its initial phase happened to be closer to the model phase.Furthermore, even after completion of the alignment process, small fluctuations on the order of a few 100 mHz are still apparent.These fluctuations are caused by the finite bandwidth of a dedicated feedback loop, which keeps the two O/L channels separated by 10 Hz in Doppler space.(Accessing the GP2021's MULTI_CHANNEL_SELECT registers led to adverse side effects in our tests.) Compared to O/L methods employed by space-based GNSS-RO receivers, which are based on predetermined code and carrier parameters derived from atmospheric climatologies (Sokolovskiy, 2001b;Sokolovskiy et al., 2006;Ao et al., 2009;Bonnedal et al., 2010), the OpenGPS O/L technique creates a different O/L model during the initialization phase (elevation angles between +2 and 0 • ) for each setting event.Thus, diurnal variations in the atmospheric refractivity field may be better accounted for in the O/L model.However, under high-humidity conditions strong refractivity gradients could lead to SNR fluctuations and subsequent (transient) loss of tracking lock already during the O/L initialization phase at positive elevation angles producing a faulty O/L model.In this case it is possible that C/L tracking actually outperforms the O/L results (see Fig. 5).

Figure 1 .
Figure 1.Schematic representation of closed-loop (top panel) and open-loop (bottom) carrier signal tracking.In closed-loop mode the input signal is correlated with two replica signals, sin and cos, and the latter is phase-shifted by 90 • with respect to the former.The correlation sums are low-pass filtered and the output is examined for phase deviations between input and replica signal.The discriminator adjustments close the feedback loop.Open-loop tracking (bottom panel) dispenses with the phase discriminator and the numerically controlled oscillator (NCO) is solely steered from model values (Doppler model).The observed carrier phase finally is assembled from the NCO phases and the in-phase and quad-phase correlation samples.

Figure 3 .
Figure 3. Carrier NCO frequency as a function of geometric satellite elevation angle.For clarity, a constant frequency value of f ref = 1.4071MHz is subtracted.The transition between closed-loop and open-loop tracking is depicted in the insert highlighting the 10 Hz offset between the two O/L channels.Same event as shown in Fig. 2.

Figure 4 .Figure 5 .
Figure 4. View towards west from the observation platform.The picture was taken about 2-3 m behind the receiver's GPS antenna, which is tilted by about 45 • towards west (the antenna and its mounting pole is visible at about 290 • azimuth angle).The metal structure at about 220 • is part of a lightning protection system and does not block the antenna's field of view, which ranges from 200 to 340 • .Throughout the observation period a satellite with a given PRN sets beyond the horizon at almost the same azimuth angle.The angles of those nine PRNs discussed in the present study are marked in red.

Figure 6 .
Figure6.Observed carrier frequencies f obs , reconstructed from NCO and residual frequencies (Eq.5), using external data bit demodulation as a function of elevation angle (red and green).In addition, the closed-loop result is plotted in blue.For clarity, a constant offset f ref = 1.4071MHz is subtracted.Same event as shown in Figs.2 and 3.

Figure 7 .
Figure 7. Normalised signal amplitudes as a function of elevation angle derived from several MPS simulations.Refractivities on the individual phase screens are calculated from an exponential profile furnished with a planetary boundary layer in the lower troposphere.Signal absorption at the surface is taken into account (red and blue lines); green and black lines show the result without ground absorption.Two boundary layers are modelled: a horizontal boundary layer top at 2 km altitude (red and black) and a layer top increasing from 1 to 2 km between 30 and about 60 km distance from the receiver (blue and green).For legibility the red, green and blue lines are shifted by an additional +1, +2 and +3 dB offset, respectively.For details see text.

Figure 9 .
Figure9.Carrier signal-to-noise density ratio, averaged over the two O/L channels, as a function of elevation angle.The panels show observations from nine PRNs sorted row-wise according to increasing azimuth angle (see Fig.4).Mean and 1σ standard deviations, calculated from 0.2 • elevation bins, are marked in dark blue.Once the signal is completely blocked by the horizon, C / N 0 drops to about 17 dB Hz (see Eq. 8).

Figure 10 .
Figure 10.Normalised C/A code correlation function as a function of delay and Doppler frequency offset.At zero delay the code correlation function corresponds to the frequency response (solid red line).It closely matches the modulus of the normalised sinc function (dotted red line, shifted by +2.5 chips for clarity).
Figure 11.Same as Fig. 9, however, showing C / N 0 , the difference of the two O/L density ratios, as a function of elevation angle.

Figure 12 .
Figure 12.Same as Fig. 9, however, showing the residual frequency in O/L channel A as a function of C / N 0 , the difference of the two O/L density ratios (gray points).The theoretically expected result (Eq.16) is shown in dark blue; light blue lines are the theoretical result shifted by ±50 Hz highlighting the occurrence of frequency aliasing.The horizontal dashed lines mark the −5 Hz residual frequency.

Figure 16 .
Figure 16.Same as Fig. 15 but the correlation analysis now includes all observations at elevations between −2 and +2 • .With the exception of PRN 7 (and the Pearson's coefficient for PRN 18) the derived correlations are no longer significant.Note the change of scale with respect to Fig. 15.
13.0642 • E. The OpenGPS receiver tracks signals from setting GPS satellites simultaneously in both, closed-loop and open-loop mode down to geometric elevation angles of −1 to −1.5 • .These low-elevation events are characterized by fluctuations of about 10-20 dB Hz in signal-to-noise density ratio and about 10-20 Hz in carrier frequency.Tracking the same event with one closed-loop and two open-loop channels in parallel allows for direct intercomparison of open-loop versus closed-loop performance.Whilst open-loop tracking allows us to follow strongly fluctuating signals to very low elevation angles, in about 14 % of the observations closed-loop tracking outperformed the open-loop channels, since fluctuations in the early phase of the setting event between +2 and 0 • elevation angle prevented proper initialization of the openloop model.

Figure A2 .
Figure A2.Screenshot of a terminal window running the OpenGPS user space application ogdspl.In the lower half characteristic information on the 12 tracking channels is shown.Note that ogdspl maps azimuth angles to the interval [−180, +180 • ] with −90 and +90 • referring to west and east, respectively.At the time of the measurement on 16 October 2015 channel 1 tracked PRN 28 at an elevation angle of +1.64 • (column 6 entitled "elev").Its azimuth value (−98.15• ) is within the selected observation window ([−160, −20 •]) and the O/L channels A and B, here listed as channel 11 and 12, are cloned from channel 1. In-phase and quad-phase samples as well as parameters describing the alignment between master and clone channels are listed at the very bottom of the screen.The top few lines list position and clock solutions from to the real-time navigation solution.
Figure8.Failure rate of internal navigation bit retrieval in terms of averaged performance parameter E (Eq.10) as a function of elevation angle for signal density ratios C / N 0 ≥ 30 dB Hz (blue, left axis).The corresponding number of observations per elevation bin is plotted as well (red, right axis).The statistics is based on 242 observations recorded in January 2014.

Table 1 .
Number of low-elevation events and azimuth angle at zero elevation for all 19 PRNs recorded during 311 days in 2014.Azimuth angle is given as mean and 1σ standard deviation.The last column lists the percentages of observation in which O/L tracking reached a lower elevation angle than the corresponding C/L detection.The first number refers to O/L channel A, the second to channel B.
. The third column lists the corresponding result, neglecting frequency deviations larger than 40 Hz.
Figure15.Standard deviation of C / N 0 at elevation angles between +1 and +2 • versus mean refractivity gradient for nine PRNs extracted from ECMWF (March to mid-July 2014).For PRN 13, 14 and 23 one data point exceeds the axis limit of 4.2 dB Hz; its respective mean refractivity gradient is marked by an arrow.(Of course, these observations are included in the statistical analysis.)In the upper left corner of each panel correlation coefficients are given (top: Pearson's coefficient; bottom: Spearman's coefficient).
5880/GFZ.2016.1.1.002and are available from GFZ's data archive at http://dx.doi.org/10.5880/GFZ.2016.1.1.002.The data are supplemented with documentation describing the measurement data files and an archive containing the OpenGPS receiver software used during the measurement campaign.The OpenGPS receiver software is free software and available at http://dx.doi.org/10.5880/GFZ.2016.1.1.002.You can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License or (at your discretion) www.atmos-meas-tech.net/10/15/2017/Atmos.Meas.Tech., 10, 15-34, 2017