1 Introduction

GPS code measurements are not only affected by frequency-dependent group delays resulting in the so-called differential code biases (DCB), but also show frequency- and elevation-dependent group delay variations (GDV). These are most pronounced in case of the exceptional satellite SVN49 with meter level differences between horizon and zenith. SVN49 is the only GPS Block IIR-M satellite with the addition of an L5 signal generation unit. The extremely large GDV are caused by satellite internal reflections of the L1 and L2 signals at the auxiliary port used for L5 (Lake and Stansell 2009).

In the course of their SVN49 investigations, Springer and Dilssner (2009) identified further GPS satellites with considerable GDV, especially SVN55, SVN43, and some other Block IIR satellites. Their GDV are almost one order of magnitude smaller than those of SVN49. The findings of Springer and Dilssner (2009) were confirmed by Haines et al. (2010, 2012, 2015) who determined GDV of the ionosphere-free linear combination from GPS measurements collected onboard the low-Earth orbiting GRACE satellite pair and found the largest GDV for some of the Block IIR satellites, larger than those of the GPS Block IIA and IIF satellites.

GDV are also caused by receiving antennas. Murphy et al. (2007) reported GDV of more than 1 m for two antenna types. These results were obtained from laboratory tests. Wübbena et al. (2008) determined elevation-dependent GDV for six receiving antenna types with their robot field calibration device which is mainly intended for the calibration of carrier phase center variations. Absolute GDV showed patterns with maximum variations of 0.5 m between horizon and zenith. Differences of the six antenna types with respect to each other were much smaller and reached up to 0.2 m. With a similar robot-based field calibration procedure, Kersten et al. (2012) determined GDV of some individual receiving antennas, with an emphasis on rover antennas. They found azimuth- and elevation-dependent GDV of up to several decimeters.

There is a severe difficulty and limitation when trying to derive satellite and receiving antenna group delay corrections. Pseudorange biases caused by signal distortions depend on correlator spacing and partly also on the receiver model itself (Hauschild and Montenbruck 2015, 2016). This had also been observed for elevation-dependent GDV in case of the anomalies of GPS SVN49 (Ericson et al. 2010; Hauschild et al. 2012). Hence, no group delay corrections exist which are valid for all receiver types and settings of the tracking channels. Although using observation data from several networks of continuously operating GPS reference stations employing various receiver brands and without documented receiver settings, we were surprised by the fairly good consistency of the GDV calibration results. Only few receivers/stations/antennas did not fit to the majority of the results. The cause for these outliers may be found in the receiver settings or in severe local code multipath effects. The fairly good consistency of the results may partly be due to the requirement of some of the network coordinators to disable multipath mitigation in the GNSS receivers (IGS 2015).

The objective of this paper is to determine GDV correction values for GPS satellite and receiving antennas. This is realized by calibrating GDV with respect to dual-frequency phase observations. This technique was already applied by Wanninger and Beer (2015) for second-generation BeiDou medium Earth (MEO) and inclined geosynchronous orbiting (IGSO) satellites which exhibit GDV of up to 1.5 m. However, the algorithm had to be considerably refined since GDV of GPS satellites are usually smaller by about one order of magnitude.

This paper is organized as follows. In Sect. 2, we present the applied algorithms. Section 3 is dedicated to the determination of satellite GDV corrections and corresponding corrections for various geodetic receiving antenna types. Finally, we apply these corrections to GPS code measurements and demonstrate the successful removal of satellite and receiving antenna GDV for three different applications: single-frequency precise point positioning (PPP; Sect. 4.1), ambiguity fixing based on the Melbourne–Wübbena linear combination (Sect. 4.2), and determination of ionospheric total electron content (TEC; Sect. 4.3). In Sect. 4, we also discuss the effect of these GDV corrections on fractional cycle biases (FCBs) and DCB. The final section summarizes the outcome of this paper.

Throughout the paper we will use RINEX 2 conventions (Gurtner and Estey 2007) to identify GPS signals with a two-character code consisting of an observation code (C for C/A or civil, P for P- or Y-code) and the frequency number (1 or 2). We thus avoid the more complex three-character coding of the RINEX 3 conventions (RINEX WG 2015) which provides no advantage to this paper that deals with traditional GPS signals only.

2 Method

2.1 Code analysis based on MP observable

Certain characteristics of GNSS code observations can be analyzed by forming a special linear combination of single-frequency code and dual-frequency phase observations which usually is called multipath (MP) observable. To our knowledge, this linear combination was first mentioned by Rocken and Meertens (1992). It is widely used to characterize the impact of high-frequency multipath effects on code observations. Its main advantages are that observations of only a single GNSS receiver are required and that code observations on the different frequencies can be analyzed separately.

The MP observable is derived from the code measurement C [m] at frequency i and a linear combination of the carrier phase measurements \(\varPhi \) [m] at frequencies i and j:

$$\begin{aligned} {\text {MP}}_i =C_i +(m_\textit{ij} -1)\cdot \varPhi _i -m_\textit{ij} \cdot \varPhi _j -B \end{aligned}$$
(1)

with

$$\begin{aligned} m_\textit{ij} =\frac{2 \lambda _i^2 }{\lambda _i^2 -\lambda _j^2 } \end{aligned}$$
(2)

where the linear factors \((m_\textit{ij} -1)\) and \(m_\textit{ij} \) are calculated from the wavelengths \(\lambda \) [m] of the two frequencies. The linear factors are selected in such a way that ionospheric and tropospheric delays and all geometric contributions (clocks, orbits, and antenna movements) cancel out. Since phase measurements are involved, the MP combination also contains carrier phase ambiguities which change with cycle slips and reappearance of a satellite. The ambiguities cannot be separated from hard- and software-induced delay differences between the different observables. These biases are lumped together in B which is considered to be constant for each continuous ambiguity block. Thus, no absolute MP values are known but only MP variations within such blocks. Differences between the various biases B have to be taken into account when processing several MP sequences in a combined analysis. Besides, a condition has to be applied in order to separate the absolute level of the MP values from the biases. Usually, a zero-mean condition over all MP values is selected (e.g., in Fig. 1). In the following, we prefer the condition to fix the MP value at \(0{^{\circ }}\) nadir angle at the satellite antenna (corresponding to an elevation angle of \(90{^{\circ }}\) at the receiving antenna) to zero.

MP variations mainly reflect high-frequency multipath effects and the tracking noise of the code measurements. They are often quantified by single RMS values for each individual code observable of a receiver. If the elevation-dependence is considered, RMS values are usually calculated for elevation bins of, e.g., \(1{^{\circ }}\) (Fig. 1). In this paper, we are not interested in the high-frequency variations of the MP series but only in their low-frequency variations. The latter contain information on GDV but are contaminated by code multipath effects. We model these low-frequency variations as a function of the nadir angle at the satellite or the elevation angle at the receiving antenna (Fig. 1).

Fig. 1
figure 1

Examples of MP values as a function of the elevation angle (blue dots), elevation-dependent RMS values (dashed green line), and regression model (solid red line)

In contrast to standard MP calculations according to Eq. (1), we slightly modified the algorithm with respect to the effect of phase wind-up (Wu et al. 1993). We corrected the carrier phase observations for the wind-up effect due to satellite rotation which removed a systematic effect from the MP values of up to a few centimeters. Another difference between carrier phase and code observations concerns higher-order ionospheric effects which are ignored in Eq. (1). And since they usually do not exceed 1 cm (Bassiri and Hajj 1993), we did not introduce any corrections for these effects.

When trying to identify an appropriate mathematical model to represent the low-frequency variations of the MP values, we implemented the estimation of coefficients of a series of spherical harmonic functions of maximum degree \(n_{\mathrm{max}}\) and maximum order \(m_{\max } \le n_{\max }\) to describe the dependence on the elevation e and the azimuth \(\alpha \) by

$$\begin{aligned} {\text {MP}}(e,\alpha )&=\sum _{n=1}^{n_{\max } } \sum _{m=0}^n {P_{nm} (\sin e)}\nonumber \\&\quad {\cdot \, (a_{nm} \cos m\alpha +b_{nm} \sin m\alpha )} , \end{aligned}$$
(3)

where \(P_{nm}\) are normalized associated Legendre functions, and \(a_{nm}\) and \(b_{nm}\) are the unknown parameters to be estimated.

Previous GDV studies revealed significant azimuth-dependent variations for GPS satellite antennas (Haines et al. 2004) and azimuth-dependent GDV of receiving antennas could be demonstrated by Kersten et al. (2012). Nevertheless, in this first attempt to determine GDV we preferred to work with a simpler model which only takes nadir- or elevation-dependent GDV into account. This can be achieved by setting \(m_\mathrm{max}\) to zero and thus dropping the modeling of any azimuth-dependence. All results presented in this paper refer to a maximum degree of \(n_\mathrm{max}=5\) and \(m_\mathrm{max}=0\) and are thus based on five coefficients.

2.2 Separation of satellite GDV from receiving antenna GDV

Nadir-dependent GDV of the satellite antenna and elevation-dependent GDV of the receiving antenna sum up to a combined effect. The nadir angle \(\eta \) at the satellite and the elevation angle e at the receiving antenna on the Earth’s surface are directly related by (Schmid and Rothacher 2003)

$$\begin{aligned} \sin \eta =\frac{R}{A}\cdot \cos e \end{aligned}$$
(4)

where R is the Earth’s radius and A the geocentric distance of the satellite which is identical to the semimajor axis of the quasi-circular GPS satellite orbit. All GPS satellites have nearly the same semimajor axis, and thus, the factor R / A in Eq. (4) becomes effectively identical for all GPS satellites.

The elevation ranging from \(0{^{\circ }}\) to \(90{^{\circ }}\) at a receiving station on the Earth’s surface corresponds to a nadir angle ranging from \(0{^{\circ }}\) to almost \(14{^{\circ }}\) at the satellite (Fig. 2). Whereas for high elevation angles an approximately linear relationship exists with respect to the corresponding nadir angles, the lowest 30% of the elevation range are squeezed into the upper 15% of the nadir angle range.

Fig. 2
figure 2

Relation between nadir angles at the GPS satellites and elevation angles at receiving stations on the Earth’s surface

We are able to extract the combined GDV effect from GPS observations. The relationship described in Eq. (4), however, prevents the separation in contributions from satellite and receiving antennas as long as all receiving stations are located on the Earth’s surface. To achieve a separation, we used a set of four receiving antenna models, all of Dorne–Margolin type, to define a reference level for our calibration. Thus, all GDV results for satellite (Sect. 3.1) and receiving antennas (Sect. 3.2) are relative corrections referring to Dorne–Margolin type antennas. The combined satellite and receiver antenna corrections (Sect. 3.3), however, can be considered as independent of the reference antenna type as long as they are used only for observations recorded by receiving stations on the Earth’s surface.

2.3 Carrier phase corrections for satellite and receiving antennas

In order to be able to calibrate GDV with respect to dual-frequency carrier phase observations, we have to refer all these observables to common reference points at the satellite and the receiving antennas. As we aim at obtaining corrections on an accuracy level of a few centimeters, the effects of incorrect antenna reference points on the GDV should be smaller than this accuracy level.

For the receiving antenna, a common reference point of the dual-frequency carrier phase observations is realized by applying absolute antenna corrections to the observations. Such corrections are published by the International GNSS Service (IGS; Dow et al. 2009) for more than 200 antenna types (Schmid 2016). Originating from field calibrations they are considered to be accurate on the subcentimeter level (Wübbena et al. 1997, 2006; Mader 1999).

For the satellite antennas, the situation is different. The IGS determines and publishes carrier phase satellite antenna corrections. Those consist of block-specific x- and y-offsets from pre-launch measurements on the one hand and estimates based on the analysis of GPS observations of the IGS network on the other hand: block-specific nadir-dependent variations and satellite-specific z-offsets. Those corrections referring to the ionosphere-free linear combination (Schmid et al. 2007, 2016) are provided in the ANTEX format (Rothacher and Schmid 2010). This format is not designed to store corrections for ionosphere-free linear combinations but only for the original frequencies. To overcome this limitation, all ANTEX models distributed by the IGS contain identical satellite antenna correction values for L1 and L2. The ionosphere-free linear combination of identical L1 and L2 values results in the same values again. Since we need realistic corrections for the original frequencies, the IGS values are not useful for our application.

A carrier phase z-offset of a satellite antenna being incorrect by \(\varDelta z\) [m] has an impact \(\varDelta v\) [m] on the nadir-dependent variations. \(\varDelta v\) can be estimated from the difference between the effects at the minimum nadir angle (\(\eta _{\min } =0^{\circ })\) and the maximum nadir angle (\(\eta _{\max } \approx 14^{\circ })\) by

$$\begin{aligned} \Delta v=(\cos \eta _{\min } -\cos \eta _{\max } )\cdot \Delta z=0.030 \Delta z. \end{aligned}$$
(5)

To determine GDV with an accuracy of a few centimeters, the z-offset error \(\varDelta z\) should not exceed 0.5 m.

In particular, the ionosphere-free IGS z-offset values for Block IIA satellites of around 2.5 m greatly differ from the geometric distance between the satellite’s center of mass (CoM) and the antenna, which only amounts to about 1.0 m. Therefore, the difference between the ionosphere-free and the frequency-specific z-offset values is most probably bigger than 0.5 m. This is not acceptable for our application. To obtain improved approximate satellite antenna phase center corrections for L1 and L2, we estimated frequency-specific z-offset values, but did not modify the IGS x- and y-offset values and the nadir-dependent phase center variations.

We assume a common geometrical distance \(z_{12}\) between the CoM and the base of the antenna elements to be a first rough estimate for the L1 and L2 phase center z-offsets. Furthermore, we take into account the ionosphere-free z-offset estimates \(z_{0}\) provided by the IGS. Together with the assumption that the arbitrarily chosen \(z_{12}\) is exactly in the middle between the actual z-offset values for L1 and L2, we get the following equations for \(z_{1}\) and \(z_{2}\):

$$\begin{aligned} z_{12}= & {} (z_1 +z_2 )/2 \end{aligned}$$
(6)
$$\begin{aligned} z_0= & {} \frac{f_1^2 }{f_1^2 -f_2^2 }\cdot z_1 -\frac{f_2^2 }{f_1^2 -f_2^2 }\cdot z_2 \end{aligned}$$
(7)

Rearranging these two equations leads to:

$$\begin{aligned} z_1= & {} z_{12} -d \end{aligned}$$
(8)
$$\begin{aligned} z_2= & {} z_{12} +d \end{aligned}$$
(9)

with

$$\begin{aligned} d=\frac{f_1^2 -f_2^2 }{f_1^2 +f_2^2 }\left( {z_{12} -z_0 } \right) \end{aligned}$$
(10)

where \(f_{1},f_{2}\) are the GPS signal frequencies.

To obtain reliable and fairly accurate (uncertainty of a few dm) geometrical distances \(z_{12}\) between the CoM and an assumed common phase center, we evaluated information on the physical dimensions of the GPS satellites from several sources: from the International Laser Ranging Service (ILRS; Pearlman et al. 2002) on the two Block IIA satellites with laser retro-reflectors, documented, e.g., by Bar-Sever et al. (2009), from the ground calibration of a Block IIA satellite antenna (Mader and Czopek 2002), and from the Los Angeles Air Force Base on IIR and IIF satellites (LAAFB 2012a, b).

Eventually, we decided to use a \(z_{12}\) value of 1.0 m for all Block IIA, of 1.5 m for all Block IIR and of 1.6 m for all Block IIF satellites. Based on the IGS corrections \(z_{0}\), we calculated satellite-specific offsets \(z_{1}\) and \(z_{2 }\) using Eq. (8) and (9). These values are summarized in Table 1. Please note that the IGS values of 1.6000 m for some of the Block IIF satellites are preliminary values due to the lack of satellite-specific estimates.

Table 1 z-offsets of the GPS satellite antennas
Fig. 3
figure 3

Selected observation stations for four different types of receiving antennas

The differences between these frequency-specific z-offsets \(z_{1}\), \(z_{2}\) and the IGS values \(z_{0 }\) exceed the meter level for the IIA satellites. They are smaller than 1 m for IIR and negligible for IIF satellites. When we perform the calibration of GDV with respect to dual-frequency carrier phases applying the IGS corrections \(z_{0}\) for both frequencies, the results for Block IIA satellites are significantly worse than those of other satellite blocks. After considering the frequency-specific z-offsets, we obtained results with similar accuracies for all types of satellites.

3 Determination of GDV

3.1 Satellite antenna GDV

The estimation of satellite antenna GDV requires a global network of observing stations so that each individual satellite contributes observations over the whole nadir angle range from \(0{^{\circ }}\) to almost \(14{^{\circ }}\). Since GPS satellites repeat their ground tracks every (sidereal) day, except for orbit maneuver periods, an extension of the observation period to several weeks or months does not provide additional information. From first test computations, we concluded that the GDV are highly stable in time, and thus, we restricted the data processing to observations of a single GPS week (GPS week 1843: DoY 123-128/2015, May 3–9, 2015).

Fig. 4
figure 4

GDV of 31 GPS satellites based on the observation data of all 43 stations shown in Fig. 3

On the other hand, we generated redundant results by selecting four independent networks, each consisting of 10–12 stations (Fig. 3, see also Table 3). In each network, identical models of receiving antennas are in use. All the four models of receiving antennas are of Dorne–Margolin type and, therefore, expected to produce similar GDV on the receiver side.

The observation data had to meet the following additional requirements: data availability on at least 6 of the 7 days, elevation mask of \(5{^{\circ }}\) or lower, and RMS (MP) values of smaller than 0.5 m for elevation angles between \(10{^{\circ }}\) and \(90{^{\circ }}\). In order to be able to find a sufficiently large number of well-distributed stations, we had to accept mixed networks with respect to

  • antenna domes (various types of radomes or no radome at all),

  • receivers (various manufacturers, various receiver types of the same manufacturer, probably different settings as regards signal tracking),

  • the uniformity of four very similar Ashtech antenna types: ASH700936?_M with “?” being A, B, C, or D.

Whereas all stations deliver C1 and P2 code observations, only some of them provide P1, C2, or C5. Therefore, we restricted our analyses to C1 and P2 and to the ionosphere-free linear combination of C1 and P2.

Table 2 Correction values for GDV of GPS satellites in millimeters

Since we assume that all the four selected antenna models possess similar characteristics with regard to their elevation-dependent GDV, we also produced a combined solution. Thus, the main results of our data processing are satellite-specific GDV for C1 and P2 as shown in Figs. 4, 8 and Table 2. These results do not provide absolute correction values but relative values with respect to the selected reference receiving antenna models.

Many satellites exhibit larger GDV for C1 than for P2 (Fig. 4). In case of C1, they exceed 20 cm for a few satellites. The effect is amplified by forming the ionosphere-free linear combination with values of up to 80 cm for SVN55.

To evaluate the results from the four independent networks with different antenna models, we calculated the differences between the individual and the combined solution based on all 43 stations (Fig. 5). The RMS values of the differences for the ionosphere-free linear combination range from 5.5 to 6.9 cm which indicates that all the four solutions are of similar quality. The corresponding RMS values for C1 and P2 are in the range of only 1.8–2.4 cm and of 1.9–2.3 cm, respectively. With differences on the level of a few centimeters, our assumption of similar GDV for these four antenna types seems to be justified.

Fig. 5
figure 5

Differences between results based on single antenna types and the overall solution based on all 43 stations: satellite-specific GDV for the ionosphere-free linear combination of C1/P2

Haines et al. (2010, 2012, 2015) published results of GPS satellite GDV for the ionosphere-free linear combination based on P1/P2 observations obtained with the receivers of the GRACE satellite pair between 2002 and 2013. These GDV were estimated from observation residuals resulting from the GRACE orbit determination. The antennas onboard the GRACE satellites are choke ring antennas (Haines et al. 2015).

The data sets of Haines et al. (2010, 2012, 2015) and our results have 22 satellites in common. Figure 6a shows our results for the ionosphere-free C1/P2 linear combination from Fig. 4 again but restricted to these 22 satellites. Figure 6b depicts the corresponding P1/P2 results of Haines et al. (2010, 2012, 2015) with each curve constrained to zero in nadir direction. Figure 6c shows the differences after removing an individual bias per satellite. The satellite-specific differences exhibit an overall RMS value of only 6.6 cm. This reflects an excellent agreement considering the different kinds of data sets (ground stations vs. low-Earth orbiters), antennas, time periods, algorithms, multipath conditions, and also code observables (C1/P2 vs. P1/P2).

Fig. 6
figure 6

Comparison of a our GDV results with b those of Haines et al. (2010, 2012, 2015), and c the remaining satellite-specific differences

Figure 7 shows the ionosphere-free GDV for Block IIA, IIR, and IIF satellites. It illustrates that the IIR satellites suffer from larger variations than the IIA or IIF satellites. No distinct differences could be identified between Block IIR-A on the one hand and Block IIR-B/IIR-M on the other.

Fig. 7
figure 7

Ionosphere-free GDV for the different GPS satellite blocks

An exceptional satellite is SVN49 whose signals were observed by 9 of the 43 selected stations in GPS week 1843 although being set unhealthy. The measurements were taken by 6 different receiver models. SVN49 is well known for its large GDV (Springer and Dilssner 2009; Ericson et al. 2010; Hauschild et al. 2012) and has, therefore, not become part of the operational GPS constellation.

Ericson et al. (2010) and Hauschild et al. (2012) demonstrated that the magnitude of the GDV of SVN49 depends on the properties of the tracking channel, especially on the activation of special multipath mitigation techniques. Our results show a GDV range of 1.7 m for C1 and of 0.4 m for P2 (Fig. 8) and, thus, match very well the results for standard early-late correlators. We conclude that probably none of the stations we used for our analyses had multipath mitigation techniques activated. Please note the good agreement of our results with older determinations demonstrating the stability of the satellite GDV with time.

Fig. 8
figure 8

GDV of GPS satellite SVN49

3.2 Receiving antenna GDV

In a second analysis step, we applied the corrections contained in Table 2 to the observation data of 132 stations equipped with 13 different antenna models in order to estimate GDV for these types of geodetic antennas. We reused the 43 stations from Fig. 3 and added 6–12 stations per antenna model for the remaining 9 antenna types (Table 3).

Table 3 The 13 antenna types and observation stations selected for the GDV determination
Fig. 9
figure 9

GDV of the 4 receiving antenna types used to determine the satellite corrections

One of the selection criteria was the availability of observations in the elevation range from at least \(5{^{\circ }}\) up to (almost) \(90{^{\circ }}\). Stations with elevation mask angles set to \(10{^{\circ }}\) and stations from polar regions were excluded. No other restrictions were applied with respect to the geographical distribution of the stations equipped with certain antenna types. Therefore, we do not show maps of the station distribution but make available the list of selected stations (Table 3). Also the receiving antenna GDV were derived from observations of GPS week 1843. For each station, we required data availability on at least 6 of the 7 days and RMS (MP) values below 0.5 m in the elevation range from \(10{^{\circ }}\) to \(90{^{\circ }}\). In order to be able to find a sufficiently large number of stations, we had to accept mixed networks with regard to antenna radomes and receiver models again.

Ideally, one would expect corrections of size zero for those four antenna types used for the determination of the satellite GDV corrections. Deviations from zero above the noise level can have various reasons: systematic multipath at the receiving stations, existing differences between the four antenna models, and differences between individual receivers or receiver settings. Figure 9 and Table 4 reveal that, above an elevation angle of \(20{^{\circ }}\), the differences are smaller than 2 cm for C1, 4 cm for P2, and 6 cm for the ionosphere-free linear combination. This confirms that the antennas have similar properties above an elevation of \(20{^{\circ }}\). At lower elevation angles, the agreement gets worse, especially for P2 and, thus, also for the ionosphere-free linear combination. We are not able to identify the major cause of these deviations and cannot exclude any of the possible reasons mentioned above. But the good agreement above \(20{^{\circ }}\) elevation confirms the similarity of the four antenna types with respect to GDV. Hence, it is justified to use this set of antenna models as a common reference for the determination of satellite GDV as described in Sect. 3.1.

Table 4 Correction values for GDV of 13 receiving antenna types in millimeters

The remaining nine antenna types show various levels of GDV (Fig. 10; Table 4). Some models exhibit GDV very similar to the set of four reference antenna types, i.e., deviations of less than a few centimeters in C1 and P2 for elevation angles above \(20{^{\circ }}\): e.g., TRM41249.00 and TRM59800.00. Three antenna types show significant deviations from zero above an elevation of \(20{^{\circ }}\): JAV_RINGANT_G3T, LEIAR25.R3, and LEIAR25.R4. Please note that the two revisions of the Leica AR25 have similar physical dimensions, similar carrier phase center offsets and variations, but still exhibit different GDV, especially in C1 and, thus, also in the ionosphere-free linear combination of C1 and P2.

Fig. 10
figure 10

GDV of all 13 receiving antenna types

3.3 Combined GDV

The correction values presented in Tables 2 and 4 refer to a set of reference antenna types, i.e., the four receiving antenna models used to determine satellite GDV. Thus, they can be considered as relative corrections with respect to the set of reference antenna types. Adding such relative corrections for satellite and receiving antennas yields combined correction data sets which are of course only suitable for specific pairs of satellite and receiving antennas. They are only valid for receiving antennas located on the Earth’s surface. Using 31 relative correction data sets of Table 2 (SVN49 ignored) and 13 relative correction data sets of Table 4, we are able to calculate \(31 \times 13 = 403\) combined correction data sets. They are shown in Fig. 11 for C1, P2, and their ionosphere-free linear combination.

The overall error due to GDV can reach some decimeters in C1 and P2, but almost 1.0 m in their ionosphere-free linear combination (Fig. 11). Since GDV are also of importance for GPS-based time transfer applications, it should be mentioned that these maximum GDV in the ionosphere-free linear combination correspond to 3.3 ns. The combinations of GPS satellite and receiving antennas yielding the largest errors are labeled in Fig. 11.

Fig. 11
figure 11

Combined satellite and receiver antenna GDV corrections for all combinations of the 31 satellites and 13 receiver antenna types (403 data sets per panel)

4 Application of GDV correction models

4.1 Single-frequency ionosphere-free PPP

We performed several tests to demonstrate the effects of the obtained GDV corrections. The first test shows results of single-frequency precise point positioning (PPP) based on the ionosphere-free linear combination \(\varPhi _0\)[m] of code observations C[m] and carrier phase observations \(\varphi \)[cy] (Yunck 1993; Choy 2011):

$$\begin{aligned} \varPhi _0 =\frac{C_i +\lambda _i \cdot \varphi _i }{2} \end{aligned}$$
(11)

where the index i indicates a specific frequency and \(\lambda \)[m] is the associated carrier wavelength. This combination exploits the fact that code and carrier phase are affected by the same amount of first-order ionospheric effects, but with opposite sign (code delay and phase advance) so that the average value forms an ionosphere-free linear combination. In PPP, \(\varPhi _0 \) is treated as a carrier phase measurement with half wavelength ambiguity and a noise level of half the code noise.

Absolute coordinates on accuracy levels of a few centimeters horizontally and some centimeters vertically can be obtained when applying this processing method to continuous data sets of 24 h. For our test, we considered all 470 daily data sets of DoY 160/2015 (about 5 weeks after the calibration period) provided by the CDDIS archive of GPS reference stations. Only stations equipped with one of the 13 antenna models of Table 4 were considered. Furthermore, we excluded all data sets with less than 23 h of observations and also those with less than 5 GPS satellites per measurement epoch on average. The remaining 317 data sets were processed several times:

  1. (1)

    dual-frequency carrier phase-based PPP to get centimeter accurate reference solutions,

  2. (2+3)

    single-frequency PPP, both with C1 and with P2,

  3. (4+5)

    single-frequency PPP with corrections from Tables 2 and 4 applied.

The nadir- and elevation-dependent corrections mainly affect the vertical component. Figure 12 shows the distribution of vertical coordinate errors for the two frequencies, both without and with corrections applied. Biases of several centimeters are visible for single-frequency PPP results without application of the corrections. The biases completely disappear after their application. Moreover, the standard deviations of the vertical component improve from 5.5 to 3.4 cm for C1 and from 6.6 to 5.4 cm for P2 ignoring biases and outliers. The outliers may be caused by severe code multipath effects at the receiving stations.

Fig. 12
figure 12

Histograms of vertical coordinate errors of 317 single-frequency PPP solutions without and with GDV corrections applied

The horizontal components of the single-frequency PPP results change only slightly when GDV corrections are applied. Nevertheless, Fig. 13 demonstrates that coordinates obtained for stations equipped with certain antenna types can be affected by horizontal biases of several centimeters, although GDV corrections are applied: LEIAT504GG on the second frequency and JAV_RINGANT_G3T on both frequencies. Perhaps, the results could be improved by considering azimuth-dependent GDV corrections.

Fig. 13
figure 13

Horizontal coordinate errors of 317 single-frequency PPP solutions despite GDV corrections applied: LEIAT504GG (red triangles), JAV_RINGANT_G3T (green squares), all other antennas (small dots)

Fig. 14
figure 14

Impact of the combined satellite and receiver antenna GDV corrections for all combinations of the 31 satellites and 13 receiver antenna types on the MW linear combination

4.2 Melbourne–Wübbena linear combination

An important application of GPS code observations for precise positioning is the fixing of widelane ambiguities based on the Melbourne–Wübbena linear combination (MW; e.g., Geng et al. 2010; Loyer et al. 2012). Using code measurements C [m] and carrier phase measurements \(\varphi \)[cycles] at the frequencies i and j, the linear combination MW [widelane cycles] is formed by

$$\begin{aligned} \text {MW}_\textit{ij} =\frac{l_\textit{ij} }{\lambda _i }C_i +\frac{l_\textit{ij} }{\lambda _j }C_j - \varphi _i +\varphi _j \end{aligned}$$
(12)

with

$$\begin{aligned} l_\textit{ij} =\frac{\lambda _j -\lambda _i }{\lambda _i +\lambda _j } \end{aligned}$$
(13)

and \(\lambda _i ,\lambda _j \) [m] being the carrier wavelengths. Here, the linear combination is expressed in units of widelane (WL) cycles to demonstrate effects on ambiguity fixing. The combined effect of satellite and receiver antenna GDV can be estimated using the first two summands of Eq. (12), and the corrections are provided in Tables 2 and 4. The largest combined corrections reach about 0.3 cycles for specific combinations of GPS satellite and receiving antenna (Fig. 14), and thus, they should not be neglected in fast and reliable widelane ambiguity fixing.

An example of the MW linear combination for a complete satellite pass is shown in Fig. 15. The combination of satellite SVN60 with receiving antenna model JAV_RINGANT_G3T was selected to demonstrate the mitigation effect of our corrections in case of large GDV. The unfiltered data set is dominated by code noise and cannot reveal the improvements resulting from the corrections. After low-pass filtering, the error mitigation due to the GDV corrections gets visible. With corrections from Tables 2 and 4 applied, the MW time series is much more stable and reflects a constant widelane signal without any elevation-dependence.

This is also an indication for the successful separation of GDV on the L1 and L2 frequencies. Although the corrections were determined based on the MP linear combination, they show an excellent performance for this completely different MW linear combination.

Fig. 15
figure 15

MW linear combination of dual-frequency C1/P2 observations for a whole satellite pass of SVN60 at IGS station RIO2 (Rio Grande, Argentina; antenna type JAV_RINGANT_G3T, DoY 160/2015): original (dots) and low-pass filtered (solid lines) data

Fig. 16
figure 16

\(\hbox {FCB}_{\mathrm{WL,C1P2}}\) for MW linear combination C1-P2, DoY 160/2015: a own determination (green squares) and values from CNES with CODE \(\hbox {DCB}_{\mathrm{P1C1}}\) applied (blue dots), b difference between the two data sets from (a), c own determination without (green squares) and with (red triangles) carrier phase center and GDV corrections applied, and (d) difference between the two data sets from (c)

Undifferenced ambiguity fixing requires fractional cycle bias (FCB) values for all individual satellites. They are estimated from the observation data of global networks of reference stations. As they are found to be fairly stable for the MW linear combination, a daily update seems to be sufficient to account for their temporal variations. CNES-CLS (Centre National d’Etudes Spatiales—Collecte Localisation Satellites, France) determines and publishes \(\hbox {FCB}_{\mathrm{WL}}\) values under the abbreviation WSB (widelane satellite biases) on its FTP server (Loyer et al. 2012).

The \(\hbox {FCB}_{\mathrm{WL}}\) values provided by CNES-CLS refer to the code observations on P1 and P2. Thus, we call them \(\hbox {FCB}_{\mathrm{WL,P1P2}}\) [WL cy]. As we use C1 instead of P1, we have to apply the P1-C1 differential code biases \(\hbox {DCB}_{\mathrm{P1C1}}\) [ns] in order to obtain \(\hbox {FCB}_{\mathrm{WL,C1P2}}\) [WL cy]:

$$\begin{aligned} {\mathrm{FCB}_\mathrm{WL,C1P2} }= & {} {\mathrm{FCB}_\mathrm{WL,P1P2} -\frac{f_1^2 -f_1 f_2 }{f_1 +f_2 }\cdot \mathrm{DCB}_\mathrm{P1C1} } \nonumber \\&{\approx \mathrm{FCB}_\mathrm{WL,P1P2} -0.196\cdot \mathrm{DCB}_{P1C1} } \end{aligned}$$
(14)

Monthly \(\hbox {DCB}_{\mathrm{P1C1}}\) values are determined and published by CODE (Center for Orbit Determination in Europe, Bern, Switzerland; see Schaer 2014). Making use of these monthly values, we calculated \(\hbox {FCB}_{\mathrm{WL,C1P2}}\) for DoY 160/2015.

Apart from that, we determined \(\hbox {FCB}_{\mathrm{WL,C1P2}}\) values from the observations of the 317 IGS reference stations on DoY 160/2015 that were already used in Sect. 4.1. First, MW FCB values were determined for each station individually. Afterward, they were combined taking into account widelane ambiguities and station-specific signal delays. Thus, we obtained \(\hbox {FCB}_{\mathrm{WL,C1P2}}\) values for each GPS satellite valid for DoY 160/2015.

Zero-mean conditions were applied to make the two \(\hbox {FCB}_{\mathrm{WL,C1P2 }}\) data sets comparable. The comparison (Fig. 16a, b) reveals an excellent agreement between the two sets with an RMS error of only 0.03 widelane cycles. Hence, we can be confident that our algorithm and data processing procedure produce reliable and accurate results.

We repeated the data processing with carrier phase center corrections (Sect. 2.3) and GDV corrections (Sects. 3.1, 3.2) applied. The results differ very much from those of the first processing (Fig. 16c, d). The RMS value of the differences reaches 0.26 widelane cycles. The differences for the Block IIF satellites are very similar which corresponds quite well to the earlier finding that the GDV of the Block IIF satellites are more consistent than those of the Block IIR satellites (cf. Fig. 7).

In a subsequent step, we tested the GDV corrections and corresponding \(\hbox {FCB}_{\mathrm{WL,C1P2}}\) with regard to their effect on MW ambiguity fixing. We reused the observation data sets of 317 stations of DoY 160/2015, subdivided them into observation periods of 15 min each, and computed fractional MW ambiguities. For the complete data set and more than 200,000 estimated MW ambiguities, we see no effect of the GDV corrections on the overall statistics of the fractional ambiguities. We suspect that other error sources, e.g., code multipath and noise, or receiver-specific FCB, dominate over potential improvements by GDV corrections.

A positive effect of GDV corrections, however, is found when selecting subsets of the MW ambiguity estimates. When we restrict the statistical analysis to satellite and receiving antenna pairs with large GDV and to specific receiving antennas types, and, thus, to specific receiver types, we obtain smaller fractional MW ambiguities. Figure 17 shows results of those 17 satellite and receiving antenna combinations which involve LEIAR25.R3 or LEIAR25.R4 and have GDV of larger than 0.2 m. In this example, the percentage of MW ambiguity estimates with deviations of less than 0.15 cy from the closest integer increases from 84.8 to 87.5 when applying the GDV corrections.

Fig. 17
figure 17

Distribution of fractional MW ambiguities without and with application of GDV corrections and corresponding \(\hbox {FCB}_{\mathrm{WL,C1P2}}\) (receiving antenna types: LEIAR25.R3, LEIAR25.R4, only for combined GDV larger 0.2 m, observation periods of 15 min, sample size: 4,852)

The main conclusions with respect to PPP ambiguity fixing are:

  • the GDV corrections are able to slightly improve ambiguity fixing based on the MW linear combination.

  • FCB values significantly change when applying carrier phase center and GDV corrections, and, thus, these corrections should be considered in the FCB determination to maintain consistency.

4.3 TEC determination

Dual-frequency GNSS observations are a valuable source for ionospheric total electron content (TEC) determination. The TEC can be derived from dual-frequency code, dual-frequency carrier phase or even from single-frequency code and carrier phase observations. All these observations are affected by hard-/software-induced signal delays both at the satellite and at the receiver. Furthermore, carrier phase observations are biased due to their ambiguities.

The TEC [TECU] (\(1\hbox { TECU }= 10^{16}\hbox { el}/\hbox {m}^{2})\) can be computed from dual-frequency code observation C [m] at frequencies i and j by (Klobuchar 1996)

$$\begin{aligned} \mathrm{TEC}_{ij} =\frac{1}{40.3\times 10^{16}}\cdot \frac{f_i^2 \cdot f_j^2 }{f_j^2 -f_i^2 }\cdot \left( {C_i -C_j } \right) +\mathrm{Bias} \end{aligned}$$
(15)

and f being the signal frequency in [Hz].

Whenever code observations are used for TEC determination, elevation-dependent GDV affect the results. The combined effect of satellite and receiver antenna GDV can be estimated using the first two summands of Eq. (15) and the corrections provided in Tables 2 and 4. The largest combined corrections reach about 4 TECU for specific combinations of GPS satellite and receiving antenna (Fig. 18), and thus, their application is strongly recommended for precise TEC determination based on code observations.

The vertical axis of Fig. 18 is labeled with two different scales. Delay differences between C1 and P2 can be interpreted as an effect of ionospheric refraction and, thus, converted to TEC, or they are considered as hard-/software-induced delays and named DCB. In practice, differences between dual-frequency code observations must be corrected for satellite and receiver DCB to enable the estimation of unbiased TEC values.

Fig. 18
figure 18

Impact of the combined satellite and receiver antenna GDV corrections for all combinations of the 31 satellites and 13 receiver antenna types on the delay difference between C1 and P2

Fig. 19
figure 19

Differences between dual-frequency code- and phase-based TEC determinations for a whole satellite pass of SVN55 at IGS station KRGG (Kerguelen Islands; LEIAR25.R4, DoY 160/2015): original (dots) and low-pass filtered (solid lines) data

Fig. 20
figure 20

\(\hbox {DCB}_{\mathrm{C1P2}}\): a own determination for DoY 160/2015 (green squares) and values from CODE for June 2015 (blue dots), b difference between the two data sets from (a), c own determination without (green squares) and with (red triangles) GDV corrections applied, and d difference between the two data sets from (c)

An example of TEC errors from dual-frequency code observations is shown in Fig. 19 for a complete satellite pass. The combination of satellite SVN55 with receiving antenna model LEIAR25.R4 was selected to demonstrate the mitigation effect of our corrections in case of large GDV. Due to the lack of reliable and accurate reference TEC, we compare our estimates from dual-frequency code observations with the ones obtained from dual-frequency phase observations. The latter are less noisy and less affected by multipath but contain an ambiguity-related bias. Hence, Fig. 19 shows differences between these two TEC determinations with an emphasis on the elevation-dependent variations over the whole satellite pass. We do not intend to illustrate any effects of satellite DCB, but only of the GDV.

The unfiltered differences are dominated by code noise and only partly reveal the improvements resulting from the GDV corrections. After low-pass filtering, the mitigation effect due to the corrections gets visible. With corrections from Tables 2 and 4 applied, the TEC time series does no longer show the anomaly which is visible in the uncorrected data set for elevations above \(60{^{\circ }}\).

Finally, we determined a set of \(\hbox {DCB}_{\mathrm{C1P2}}\) values following the approach of Montenbruck et al. (2014). This procedure makes use of IGS global ionosphere maps (GIM) to correct for the impact of ionospheric refraction on code observations. \(\hbox {DCB}_{\mathrm{C1P2}}\) values for individual satellite-receiver combinations can thus be obtained from averaged ionosphere-corrected C1-P2 differences over whole satellite passes. Using the observations on DoY 160/2015 of 317 globally distributed stations equipped with one of the antenna models of Table 4, we computed \(\hbox {DCB}_{\mathrm{C1P2}}\) values for each combination of satellite and station. In a second step, we computed a common set of satellite biases by separating satellite-specific from receiver-specific \(\hbox {DCB}_{\mathrm{C1P2}}\) by introducing a zero-mean condition for the average of all the satellite biases (Fig. 20a).

To evaluate our results, we also considered the satellite \(\hbox {DCB}_\mathrm{P1P2}\) and \(\hbox {DCB}_{\mathrm{C1P1}}\) values for June 2015 published by CODE (Schaer 2014). The difference between these two data sets from CODE yields \(\hbox {DCB}_{\mathrm{C1P2}}\) values (Fig. 19a). The differences between the results provided by CODE and our estimates have a scatter of 0.10 ns RMS which is rather small (Fig. 20b).

Then, we applied the GDV corrections from Tables 2 and 4 and repeated the determination. \(\hbox {DCB}_{\mathrm{C1P2}}\) values derived from uncorrected and corrected code observations are compared in Fig. 20c, d. Their differences have a scatter of 0.33 ns RMS, and the largest difference amounts to 0.7 ns. This demonstrates that DCB values significantly change when applying the GDV corrections.

The main conclusions with respect to the TEC determination are:

  • GPS code measurements should be corrected for GDV to improve the accuracy of the TEC determination from code observations.

  • \(\hbox {DCB}_{\mathrm{C1P2}}\) values significantly change when applying GDV corrections, and, thus, these corrections should be considered in the DCB determination to maintain consistency.

5 Conclusions and outlook

The combined GDV of the present GPS satellites and a set of 13 receiving antenna types reach up to 0.4 m in C1 and P2, and up to 1.0 m in their ionosphere-free linear combination (corresponding to 3.3 ns), up to 0.4 widelane cycles in the Melbourne–Wübbena linear combination, and up to 1.5 ns (corresponding to 4 TECU) in the difference between C1 and P2. They can be calibrated with respect to dual-frequency carrier phase observations. Data of a large number of globally distributed observation sites are needed to cover the whole nadir angle range of the satellites and to mitigate local code multipath effects. The presented corrections should be used with care since receiver tracking channels with activated multipath mitigation techniques may experience different GDV.

The GDV corrections for GPS satellite and receiving antennas obtained in this paper refer to a set of reference antennas and are, thus, relative by their nature. One way to achieve absolute corrections of receiving antennas are robot calibrations as demonstrated by Wübbena et al. (2008) and Kersten et al. (2012). It is desirable to combine both techniques to achieve sets of absolute corrections for GPS receiving and satellite antennas.