<figure><img src="image_1.png" /></figure>

# Processing of GNSS constellations and ground station networks using the raw observation approach

Sebastian Strasser<fnref n="1" /> ID · Torsten Mayer-Gürr<fnref n="1" /> ID · Norbert Zehentner<fnref n="1" />

Received: 9 July 2018 / Accepted: 2 December 2018 / Published online: 13 December 2018
© The Author(s) 2018

## Abstract

This article describes the raw observation approach as implemented at Graz University of Technology to determine GNSS products like satellite orbits, clocks, and station positions. To assess the performance of the approach, 15 years (2003–2017) of observations from a network of 245 globally distributed IGS stations to the GPS constellation were processed on a daily basis using the IGS14 reference frame and antenna calibrations. The resulting products are evaluated against those determined by IGS analysis centers. Orbit fit quality relative to the IGS combination is comparable to the best-fitting solutions used for evaluation. Starting from early 2017, when the IGS switched to IGS14, the determined orbits fit better to the IGS combination than any other considered solution. Midnight discontinuities show good internal orbit consistency and no noticeable satellite block-dependency. Satellite clocks are comparable to the considered IGS analysis center solutions. Station positions differ from the IGS combination on a similar level to the solutions they were evaluated against. The temporal repeatability of station positions is slightly better than that of the IGS combination. The quality of resulting GNSS products confirms that the raw observation approach is well suited for the task of determining satellite orbits, clocks, and station positions. It provides an alternative to well-established approaches used by IGS analysis centers and simplifies the introduction of additional observables from new and modernized GNSS.

**Keywords** Orbit determination · Satellite clock · Station position · Processing strategy · GPS · IGS

## 1 Introduction

Many global navigation satellite system (GNSS) applications require high-precision GNSS products like satellite orbits and clocks or station positions. The analysis centers (AC) of the International GNSS Service (IGS) (Johnston et al. 2017) routinely generate such products by processing observations from a global GNSS station network. For the longest time, the US Global Positioning System (GPS) and to some degree the Russian Global Navigation Satellite System (GLONASS) have been the primary GNSS in these processing efforts. Many methods and techniques have been developed by the ACs over time with the characteristics of these systems in

mind, e.g., for dual-frequency GPS and/or GLONASS processing (Dach et al. 2009; Loyer et al. 2012). Generally, IGS ACs utilize zero- or double-difference approaches (Weiss et al. 2017) based on the ionosphere-free linear combination (Hauschild 2017b).

With the advent of new GNSS like the Chinese BeiDou Navigation Satellite System (BDS) and the European system Galileo as well as the modernization of GPS and GLONASS in the recent years, these existing processing approaches faced the challenges of incorporating new observables and related parameterizations. In the last few years, a lot of progress has been made within the scope of the IGS Multi-GNSS Experiment (MGEX) (Montenbruck et al. 2014, 2017) with respect to GNSS processing, among other aspects. Some IGS ACs have successfully implemented processing strategies to facilitate multi-GNSS within the MGEX, e.g., the Center for Orbit Determination in Europe (CODE) (Prange et al. 2017), the German Research Centre for Geosciences (GFZ) (Uhlemann et al. 2015), or Wuhan University (Guo et al. 2016). Still, the performance of IGS

<figure><img src="image_2.png" /></figure>

<footnote>Sebastian Strasser
sebastian.strasser@tugraz.at</footnote>

<footnote>Torsten Mayer-Gürr
mayer-guerr@tugraz.at</footnote>

<footnote n="1">¹ Institute of Geodesy, Graz University of Technology, Steyrergasse 30/III, 8010 Graz, Austria</footnote>
multi-GNSS products is not yet competitive to that of GPS and GLONASS products and much work remains to be done.

One of the key tasks concluded by Montenbruck et al. (2017) was for the IGS to establish multi-GNSS processing standards and products for precise point positioning (PPP) using undifferenced and uncombined observations. This reflects the main motivation that led to the development of the GNSS processing approach presented in this article. Direct utilization of raw observations in GNSS processing was first demonstrated by Schönemann et al. (2011) and detailed by Schönemann (2014). Zehentner and Mayer-Gürr (2014, 2016) used a PPP approach based on raw GPS observations to determine kinematic orbits of low-Earth-orbit (LEO) satellites. This raw observation approach was then developed further at Graz University of Technology and generalized to also enable processing of observations from GNSS station networks. The key difference to well-established GNSS processing approaches is the fact that observations are used directly as observed by the receiver. This allows full exploitation of the information contained in each individual observation type and preserves the original measurement accuracy. Observation equations are set up individually for each observable without explicitly forming any linear combinations or differences. This simplifies the inclusion of new observables provided by the new and modernized GNSS.

The aim of this article is to describe the raw observation approach in detail, including all corrections and effects that need to be considered and how to handle the large number of parameters. It should enable others to implement the approach for their GNSS processing. For this reason, the processing strategy used at Graz University of Technology will be presented as well. Aspects of special interest are the handling of the ionospheric influence and the related code bias estimation. Integer ambiguity resolution and the handling of phase biases will also be discussed. In order to evaluate the approach, 15 years of observations from the global IGS station network to the GPS constellation were used to derive satellite orbits, clocks, and station positions. Dual-frequency processing was chosen to be consistent with standard IGS products and to allow for comparisons with IGS ACs.

This article is structured into five major sections. After the introduction, Sect. 2 describes the raw observation approach in detail and illustrates the steps required to estimate a solution. In Sect. 3, the processing strategy implemented at Graz University of Technology is detailed. The resulting GNSS products are presented and evaluated against those of the IGS ACs in Sect. 4. Finally, Sect. 5 summarizes the main points and draws some conclusions.

# 2 Raw observation approach

## 2.1 Observation equations

The observation equations of a standard, single-system GNSS measurement model may be expressed as

$$R_{r,j}^s - \Delta R_{r,j}^s = \rho_r^s + c (\delta_r - \delta^s) + T_r^s + I_{r,j}^s \\ + B_{r,j} + B_j^s + \epsilon_{r,j}^s \quad (1)$$

and

$$L_{r,j}^s - \Delta L_{r,j}^s = \rho_r^s + c (\delta_r - \delta^s) + T_r^s + I_{r,j}^s \\ + \lambda_j (b_{r,j} + b_j^s + N_j) + \epsilon_{r,j}^s \quad (2)$$

for code and phase observations, respectively (Hauschild 2017a). The indices used in this notation relate a term to a receiver r, transmitter s, and signal type j. The frequency is implied in the signal identifier j. Time dependence was omitted for the sake of readability. The terms used in (1) and (2) are: geometrical distance $\rho_r^s$, speed of light c, receiver clock error $\delta_r$, transmitter clock error $\delta^s$, tropospheric delay function $T_r^s$, ionospheric delay function $I_{r,j}^s$, receiver code bias $B_{r,j}$, transmitter code bias $B_j^s$, wave length $\lambda_j$, receiver phase bias $b_{r,j}$, transmitter phase bias $b_j^s$, integer ambiguity $N_j$, and observation noise $\epsilon_{r,j}^s$.

The a priori code and phase corrections $\Delta R_{r,j}^s$ and $\Delta L_{r,j}^s$ comprise all effects that can be adequately modeled and are not contained in any of the right-hand side terms of (1) and (2). Both transmitter and receiver move in inertial space during signal travel time. Space-time curvature caused by Earth's gravitational field also affects the signal. Relativistic effects influence the nominal frequency of transmitter clocks. Hofmann-Wellenhof et al. (2008) provide formulas to correct for those effects. Antenna phase centers have an offset to their reference point, which in turn can have an offset to a satellite's center of mass or a station marker. In addition, the phase center varies depending on the frequency and direction of a signal (Schmid et al. 2005, 2007). Phase observations are also affected by phase windup, which is caused by changes of the mutual orientation of transmitting and receiving antennas (Kouba 2009a). Tidal and loading displacements affect station positions on a sub-daily scale and have to be corrected for as well.

The raw observation approach is a zero-difference method, which means parameters depending only on a transmitter or receiver do not cancel out like in case of the classic double-difference method (Hauschild 2017b). Therefore, clock errors are set up as parameters for each transmitter and receiver for every epoch directly in the combined least squares adjustment. When both transmitter and receiver clocks are estimated at the same time, there is a rank deficiency in the system of equations. To solve this, a zero-mean
constraint is added to the transmitter clocks. This constrains the mean value over those clock parameters to zero at every epoch. Alternatively, one or more very stable clocks (e.g., stations tied to timing labs) can be held fixed as a reference while all others are estimated relative to those (Weiss et al. 2017).

Transmitter and receiver instrumental biases have to be considered as well when using a zero-difference method. Håkansson et al. (2017) give a good overview of code and phase biases in GNSS processing. Code biases cannot be directly estimated in the combined least squares adjustment. There is a rank deficiency since they are fully correlated with clocks and ionospheric total electron content (TEC). Section 3.3 will explain how they can be estimated in an intermediate step during processing. Phase biases prevent integer ambiguity resolution for zero-difference methods. This issue will be touched upon further in Sect. 3.4.

Following the International Earth Rotation and Reference Systems Service (IERS) conventions (Petit and Luzum 2010), the tropospheric delay function can be expressed as

$$T_r^s = m_h(e) D_{zh} + m_w(e) D_{zw} + m_g(e) [G_N \cos a + G_E \sin a] \quad (3)$$

for line-of-sight observations between a station and satellite at elevation e and azimuth a. It consists of a zenith hydrostatic delay Dzh, a zenith wet delay Dzw, and a horizontal delay gradient with north (G_N) and east (G_E) components. They are mapped into line-of-sight with hydrostatic, wet, and gradient mapping functions mh, mw, and mg, respectively. The hydrostatic component accounts for the majority of the tropospheric delay and can be modeled adequately. Modeling the wet component is only possible to some degree. Therefore, a residual wet delay has to be estimated per station. Estimating horizontal delay gradients decreases systematic errors in GNSS solutions when observations at low elevation are used (Bar-Sever et al. 1998).

## 2.2 Ionospheric delay

One of the main characteristics of the raw observation approach is that it does not explicitly use linear combinations of observations. The ionosphere-free combination is commonly used in GNSS processing to eliminate the first-order ionospheric delay (Hauschild 2017b). It is formed separately for code and phase observations. This implies that the ionospheric delay is eliminated independently for code and phase. The combined observable also has an increased standard deviation, e.g., by a factor of 3 for GPS L1/L2 (Misra and Enge 2011). With the raw observation approach, a common slant total electron content (STEC) parameter is estimated for

all code and phase observations between a transmitter and a receiver at an epoch (referred to as observation group hereafter) via least squares adjustment. The ionospheric delay function can be expressed as

$$I_{r,j}^s = I_{r,j}^{s,(1)} + I_{r,j}^{s,(2)} + I_{r,j}^{s,(3)} + I_{r,j}^{s,EPL} + I_{r,j}^{s,dTEC} \quad (4)$$

when using the raw observation approach. All terms in (4) are functions of STEC. Code and phase are differentiated by the signal identifier j, which also implies the frequency. The terms $I_{r,j}^{s,(1)}$, $I_{r,j}^{s,(2)}$, and $I_{r,j}^{s,(3)}$ are the first-, second-, and third-order ionospheric correction functions, respectively. They differ for code and phase measurements, delaying the former and advancing the latter (Fritsche et al. 2005). $I_{r,j}^{s,EPL}$ is the excess path length of a GNSS signal's curved path through the ionosphere that is not described by the first three terms. $I_{r,j}^{s,dTEC}$ is what Hoque and Jakowski (2008) call the range error due to TEC difference at different frequencies. Since a signal's path through the ionosphere is not the same for different frequencies, the respective TEC values also differ slightly for each frequency. Hoque and Jakowski (2008) provide empirical formulas to correct for both $I_{r,j}^{s,EPL}$ and $I_{r,j}^{s,dTEC}$.

GNSS observation equations are nonlinear and have to be solved iteratively. $I_{r,j}^{s,(3)}$, $I_{r,j}^{s,EPL}$, and $I_{r,j}^{s,dTEC}$ are nonlinear functions of STEC, so their partial derivatives in the design matrix require STEC values. Since STEC parameters are estimated when using the raw observation approach, they can be used to consider these corrections during parameter estimation beginning from the second iteration. Therefore, no ionosphere models are required to correct for these effects. The second-order ionospheric delay can amount to a few centimeters in zenith direction, while the other effects are on the millimeter level (Hoque and Jakowski 2008).

## 2.3 Implementation

Figure 1 visualizes a way to estimate a solution using the raw observation approach starting from Eqs. (1) and (2). It is structured into four levels. Level 1 deals with observation equations for one receiver-transmitter pair at a single epoch. Since memory limits make it unfeasible to set up normal equations containing millions of STEC parameters per day (see Table 2), they are eliminated at this level. This does not affect the estimation of other parameters. In the special case of exactly two phase or code observation equations with different frequencies, eliminating the STEC parameter is the same as using the ionosphere-free combination. With the raw observation approach, eliminating the STEC parameter reduces the number of observation equations per observation group by only one, independent of how many code and phase observables are used.
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Flowchart for estimating a solution with an implementation of the raw observation approach</figcaption></figure>

Processing large networks with hundreds of stations at a high sampling rate (e.g., 30 s) leads to hundreds of thousands of clock parameters per day (see Table 2). They can be eliminated on the normal equation level to reduce memory consumption. Clock parameter elimination happens in two steps. At level 2, the receiver clock parameters are eliminated independently for each receiver at every epoch. After accumulating all receiver normal equations for an epoch, the transmitter clock parameters are eliminated for the respective epoch at level 3. To solve the rank deficiency of the clock parameters, a zero-mean constraint is added to the transmitter clocks before they are eliminated. The epoch normal equations are then accumulated at the global level 4. When both satellite orbits and station positions are estimated at the same time, the system can be freely rotated without affecting the internal geometry of the system. To solve this singularity, a no-net rotation constraint must be added to the normal equations (Weiss et al. 2017). It is usually applied to a well-distributed subset of stations, e.g., the IGS core network. In addition, a no-net translation constraint is required when estimating geocenter coordinates. After solving the global normal equation system, all clock and STEC parameters can be reconstructed as if they had been estimated directly.

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Processing strategy at Graz University of Technology. Each Estimate solution step represents the procedure shown in Fig. 1</figcaption></figure>

# 3 Processing strategy

Figure 2 shows a flowchart of the current processing strategy at Graz University of Technology. It represents the procedure for a standard single-system dual-frequency code and phase processing using observations of a global station network at a moderately high sampling rate, e.g., 30 s. The data are processed in daily 24-h intervals. The following subsections will detail the individual processing steps.

## 3.1 Preprocessing and initialization

Orbit modeling is done by a standard variational equation method (Montenbruck and Gill 2000). Orbits are integrated as daily 24-h arcs at a 60-s sampling rate. The integrated orbits are fitted to approximate orbits by estimating the initial state vector and a set of solar radiation pressure (SRP) parameters from the Empirical CODE Orbit Model (ECOM) (Arnold et al. 2015). The seven estimated SRP parameters comprise constant ($D_0, Y_0, B_0$), once-per-revolution ($B_{1C}, B_{1S}$), and twice-per-revolution ($D_{2C}, D_{2S}$) accelera-
<table><caption>**Table 1** Force models used for orbit integration</caption><thead><tr><th>Effect</th><th>Model</th><th>Note</th><th>References</th></tr></thead><tbody><tr><td>Earth's gravity field</td><td>GOCO05s</td><td>n<sub>max</sub><sup>static</sup> = 60, n<sub>max</sub><sup>time-variable</sup> = 10</td><td>Mayer-Gürr et al. (2015)</td></tr><tr><td>Astronomical tides</td><td>JPL DE421</td><td>Sun and all planets</td><td>Folkner et al. (2009)</td></tr><tr><td>Solid Earth tides</td><td>IERS 2010</td><td></td><td>Petit and Luzum (2010)</td></tr><tr><td>Ocean tides</td><td>FES2014b</td><td>n<sub>max</sub> = 10</td><td>Carrere et al. (2016)</td></tr><tr><td>Pole tides</td><td>IERS 2010</td><td></td><td>Petit and Luzum (2010)</td></tr><tr><td>Ocean pole tides</td><td>IERS 2010</td><td>n<sub>max</sub> = 10</td><td>Petit and Luzum (2010)</td></tr><tr><td>Atmospheric tides</td><td>AOD1B RL06</td><td>n<sub>max</sub> = 10</td><td>Dobslaw et al. (2017)</td></tr><tr><td>Atmos. and ocean mass variations</td><td>AOD1B RL06</td><td>n<sub>max</sub> = 10</td><td>Dobslaw et al. (2017)</td></tr><tr><td>General relativity</td><td>IERS 2010</td><td></td><td>Petit and Luzum (2010)</td></tr><tr><td>Solar radiation pressure</td><td>Box-wing</td><td></td><td>Rodriguez Solano (2014)</td></tr><tr><td>Earth radiation pressure</td><td>Box-wing</td><td>Based on CERES monthly averages</td><td>Rodriguez Solano (2009)</td></tr><tr><td>Antenna thrust</td><td>IGS model values</td><td>Approximate values per GPS block</td><td>IGS (2011)</td></tr></tbody></table>

$n_{max}$ is the maximum spherical harmonics degree

tions along axes *D*, *Y*, and *B* of a spacecraft-Sun coordinate system. *D* points from the satellite toward the Sun, *Y* points along the solar panel axis, and *B* completes the right-handed system. While a priori box-wing models cover a significant part of the accelerations due to SRP, their accuracy is not sufficient. Therefore, the listed ECOM parameters are additionally estimated to cover residual effects. All force models used for orbit integration are listed in Table 1. Satellite attitude during eclipse seasons is modeled after the *eclipses.f* Fortran routine provided via the IGS (Kouba 2009b; Dilssner 2010). Satellites experiencing an outage according to Notice Advisory to Navstar Users (NANU) messages are discarded for the respective day.

Observation preprocessing is done independently for each receiver. Stations with more than 25% of epochs missing in a day are disabled. Observations below an elevation of 5° are ignored. Initial observation weights depend on the zenith angle *z* and are determined by

$$ \sigma_j(z) = \frac{\sigma_{0,j}}{\cos z} \qquad (5) $$

with $\sigma_{0,code} = 22$ cm and $\sigma_{0,phase} = 1$ mm. Observation groups without code and phase observations on at least two frequencies are discarded. To guarantee overdetermination, epochs with less than five observed satellites are discarded as well. Initial receiver clock errors are estimated from code observations. Epochs where the station position diverges by more than 100 m during initial clock estimation, e.g., due to code cycle slips, are disabled. Continuous tracks are formed from the remaining observations. Cycle-slip detection is performed for each track. It uses two complementary methods to find cycle slips that are insensitive to one of the methods. The first method is based on the Melbourne-Wübbena (MW) combination (Melbourne 1985; Wübbena 1985). It

uses a total variation denoising algorithm (Condat 2013) to find cycle slips in the MW time series. The second method is based on the geometry-free linear combination (Hauschild 2017b). A bidirectional moving-window polynomial prediction is used in this case to detect cycle slips. Each track is split at all epochs with detected cycle slips. Tracks that never exceed an elevation of 15° or contain less than 30 estimable epochs are discarded because a stable ambiguity estimation is not guaranteed for those tracks. Finally, outliers are detected and downweighted for each track based on a robust Huber M-estimator (Huber 1981; Koch 1999).

Since station positions are estimated as constant per day, sub-daily station position changes have to be modeled. Solid Earth, pole, and ocean pole tides are modeled according to the IERS Conventions (Petit and Luzum 2010). FES2014b (Carrere et al. 2016) is used to model ocean tidal loading. Nontidal atmospheric and ocean loading variations are corrected by the model AOD1B RL06 (Dobslaw et al. 2017). Atmospheric tidal loading is also modeled using AOD1B RL06. All models are applied to their respective maximum temporal and spatial resolution. No-net rotation and translation constraints are applied to a set of selected stations, e.g., the IGS14 core network.

Tropospheric delay is modeled using Vienna Mapping Functions 1 (VMF1) (Böhm et al. 2006). This model provides mapping coefficients and zenith delays for the hydrostatic (ZHD) and wet (ZWD) components on a global 2.5° × 2° grid with a sampling rate of 6 h. Zenith delay values are given at grid height. Mapping coefficients for the stations are computed from the grid via bilinear interpolation. ZHDs for each station are computed following Böhm (2017, personal communication). ZHDs at the surrounding grid points of a station are converted to pressure using the inverse model of Saastamoinen (1972). Next, these pressure values are extrapolated to the station height using
$$p = p_B \exp \left( -\frac{g}{R'_d T_v} dh \right) \qquad (6)$$

with pressure at the grid point $p_B$, mean gravity $g = 9.80665 \text{ m/s}^2$, specific gas constant $R'_d$, virtual temperature $T_v$, and height difference $dh$. $R'_d = R_g/M_d$ using the universal gas constant $R_g = 8.3143 \text{ J/K/mol}$ and the molar mass of dry air $M_d = 0.028965 \text{ kg/mol}$. The virtual temperature can be computed via

$$T_v = T \left( 1 + 0.378 \frac{e}{p_B} \right) \qquad (7)$$

using temperature $T$ and water vapor pressure $e$ from the empirical troposphere model GPT2w (Böhm et al. 2015). Bilinear interpolation is then used to compute the pressure at a station from the extrapolated grid values. The interpolated pressure at the station is converted back to ZHD using the model of Saastamoinen (1972). Since residual ZWDs are estimated anyway (see Sect. 2.1), ZWDs for each station are computed at grid height using bilinear interpolation.

Table 2 lists the parameters that are set up for each day in case of dual-frequency GPS processing using code and phase observations. Satellite initial states and SRP parameters are set up again to estimate updates when orbits are fitted to observations. In addition, a velocity change (referred to as pseudo-stochastic pulse) is estimated at the center of each 24-h orbit arc. These parameters are used to cope with residual orbit modeling deficiencies related to nonconservative forces (Hugentobler and Montenbruck 2017). The pulses are constrained to 0.1 µm/s along each axis. Transmitter and receiver clock errors are set up every 30 s if the number of observations is sufficient at the respective epoch. GPS block IIA satellites are disabled for 30 min after they exit Earth's

shadow. Their attitude is uncertain during this time due to a post-shadow recovery maneuver. Therefore, no observation equations are set up for these satellites during this period. Station positions are estimated as constant per day. A slant TEC value is set up for each observation group that is used in processing. Tropospheric zenith wet delays are estimated via degree 1 splines with hourly nodes for each station. Tropospheric delay gradients are set up using two nodes at the start and end of a day per station. Estimated Earth orientation parameters (EOP) comprise length of day (LOD) as well as constant and trend parts for polar motion.

## 3.2 Processing of a network

As shown in Fig. 2, the actual processing starts with a core network, which is a well-distributed subset of 40–50 stations from the full network. This is sufficient to determine all parameters to a point where integer ambiguity resolution is possible. The steps to estimate a solution have been described in Sect. 2 and are visualized in Fig. 1. Solutions are estimated in an iterative process using variance component estimation (VCE) (Koch 1999) to determine observation weighting. This implies that outliers are automatically downweighted. Iteration stops once a convergence threshold is reached. Receiver TEC biases are then estimated in an intermediate step (see Sect. 3.3). This is done to get realistic STEC values that can be used for higher-order and additional ionospheric corrections. Transmitter TEC biases can also be estimated in this step if no a priori transmitter code biases were used. The next step is integer ambiguity resolution, which is described in Sect. 3.4. Once ambiguities are resolved, transmitter and receiver TEC biases are reestimated before iteratively estimating a solution again to update all parameters.

After processing the core network, all additional stations are processed individually. Transmitter-dependent parameters are held fixed during this part of the processing. Since transmitter phase biases have already been determined at this point, ambiguities can be resolved independently for each station. After this step, all ambiguities in the network have been resolved. Finally, the full network is processed together and all parameters, including transmitter-dependent ones, are estimated again. After some iterations, the TEC biases are updated one more time in an intermediate step before continuing with the final iterations. On a desktop computer with a quad-core CPU at 3.6 GHz and 16 GB of RAM, one iteration using the full network takes approximately 11 min with the settings described in this section. The initial iterations including only the core network take around 3 min each.

## 3.3 TEC and code bias estimation

Simplifying Eq. (1) to the terms relevant to code bias estimation leads to

<table><caption>**Table 2** Estimated parameters for daily processing of the GPS constellation and IGS14 network using two frequencies</caption><thead><tr><th>Parameter type</th><th>Count</th><th>Daily avg.</th></tr></thead><tbody><tr><td>Satellite initial state</td><td>6 · n<sub>satellite</sub></td><td>180</td></tr><tr><td>Satellite clock error</td><td>n<sub>epoch</sub> · n<sub>satellite</sub></td><td>86,250</td></tr><tr><td>Solar radiation pressure</td><td>7 · n<sub>satellite</sub></td><td>210</td></tr><tr><td>Pseudo-stochastic pulse</td><td>3 · n<sub>satellite</sub></td><td>90</td></tr><tr><td>Ambiguities + phase b.</td><td>2 · n<sub>track</sub></td><td>22,000</td></tr><tr><td>Code biases (C1P-C1C)</td><td>n<sub>satellite</sub></td><td>30</td></tr><tr><td>Slant TEC</td><td>n<sub>observation</sub>/4</td><td>5,008,000</td></tr><tr><td>Tropo. zenith wet delay</td><td>25 · n<sub>station</sub></td><td>4425</td></tr><tr><td>Tropo. delay gradient</td><td>4 · n<sub>station</sub></td><td>708</td></tr><tr><td>Station position</td><td>3 · n<sub>station</sub></td><td>531</td></tr><tr><td>Station clock error</td><td>n<sub>epoch</sub> · n<sub>station</sub></td><td>500,000</td></tr><tr><td>Earth orientation</td><td>n<sub>x<sub>p</sub>,y<sub>p</sub></sub> + n<sub>LOD</sub></td><td>5</td></tr></tbody></table>

Daily averages represent period 2003–2017
$$R_{r,j}^{s} = \rho_{r}^{s} + c \Delta \delta_{r}^{s} + \frac{40.3}{f_{j}^{2}} \text{STEC}_{r}^{s} + B_{r,j} + B_{j}^{s} \quad (8)$$

where j would be replaced by, for example, C1W, C2W, or C1C in case of dual-frequency GPS processing. Since GNSS observations are relative, there is a rank deficiency when receiver and transmitter code biases are estimated together. This is solved by applying a zero-mean constraint to all transmitter code biases, same as for clock errors. Code biases are not accessible in an absolute sense. Therefore, only differential code biases (DCB) can be estimated. By convention, IGS clock corrections are consistent with respect to C1W and C2W observables. The corresponding DCB can be defined as

$$B_{\text{C1W},\text{C2W}} = B_{\text{C2W}} - B_{\text{C1W}}. \quad (9)$$

If the ionosphere-free linear combination of C1W and C2W is used, then $B_{\text{C1W},\text{C2W}}$ cancels out (Dach et al. 2015). This means a standard dual-frequency PPP user does not have to apply this DCB correction when using IGS products. The pseudo-absolute code biases $B_{r,j}$ and $B_{j}^{s}$ in (8) can be expressed in terms of DCBs, e.g.,

$$B_{\text{C1W}} = \frac{f_2^2}{f_1^2 - f_2^2} B_{\text{C1W},\text{C2W}} \quad (10)$$

$$B_{\text{C2W}} = \frac{f_1^2}{f_1^2 - f_2^2} B_{\text{C1W},\text{C2W}} \quad (11)$$

$$B_{\text{C1C}} = \frac{f_2^2}{f_1^2 - f_2^2} B_{\text{C1W},\text{C2W}} + B_{\text{C1W},\text{C1C}}. \quad (12)$$

When STEC values are estimated based on C1W and C2W observables, it is not possible to separate STEC and $B_{\text{C1W},\text{C2W}}$ in the least squares adjustment. This means the estimated STEC parameters are biased, which can lead to incorrect consideration of higher-order ionospheric corrections; see Sect. 2.2. This issue is solved in intermediate steps during processing, called *Estimate TEC biases* in Fig. 2. Based on the observation equation

$$\text{STEC}_{r}^{s} = m(z) \text{VTEC}_{r} + B_{r,\text{TEC}} + B_{\text{TEC}}^{s} \quad (13)$$

a least squares adjustment is solved using all estimated, biased STEC values as observations. They are mapped to a common vertical TEC (VTEC) value per receiver and epoch based on a single-layer model mapping function (Schaer 1999)

$$m(z) = \frac{1}{\cos z'} \quad \text{with} \quad \sin z' = \frac{R}{R+H} \sin z \quad (14)$$

using zenith angle $z$, radius $R = 6371$ km, and ionosphere height $H = 350$ km. Additionally, TEC biases are estimated

per receiver ($B_{r,\text{TEC}}$) and transmitter ($B_{\text{TEC}}^{s}$) as constant per day. A zero-mean constraint is applied to all transmitter TEC biases to solve the rank deficiency. Estimated TEC biases can be converted to DCBs via

$$B_{\text{C1W},\text{C2W}} = 40.3 \frac{f_1^2 - f_2^2}{f_1^2 f_2^2} B_{\text{TEC}} \quad (15)$$

where $B_{\text{TEC}}$ is expressed in TEC units (TECU) and $B_{\text{C1W},\text{C2W}}$ in meters.

If receivers track further observables in addition to C1W and C2W, e.g., C1C, the related $B_{\text{C1W},\text{C1C}}$ transmitter DCBs can be directly estimated in the main least squares adjustment as additional parameters based on (8) in combination with (12). Cross-correlation receivers provide a linear combination of code observables. This linear combination has to be considered when setting up DCBs in the respective observation equations. Once all DCBs are determined, they can be converted to pseudo-absolute code biases for the individual signals using, for example, (10), (11), and (12) in case of dual-frequency GPS processing. Raw observations can then be directly corrected in further processing steps using these biases.

## 3.4 Ambiguities and phase biases

In case of zero-difference approaches, receiver and transmitter phase biases ($b_{r,j}$, $b_{j}^{s}$) prohibit direct access to integer ambiguities $N_j$ (Hauschild 2017a). This means each track has associated ambiguities

$$n_{r,j}^{s} = b_{r,j} + b_{j}^{s} + N_{j} \quad (16)$$

where j would be replaced by, for example, L1 and L2 in case of dual-frequency GPS processing. Going through the list of tracks, a phase bias parameter is set up per frequency for each transmitter and receiver the first time it appears on a track. This means no integer ambiguity parameters are set up for these tracks, and they are only used to estimate the corresponding phase biases. For all further tracks, integer ambiguity parameters are set up in addition to the phase biases according to (16). Since there is a rank deficiency when both transmitter and receiver phase biases are estimated, a zero-mean constraint is applied to all transmitter phase biases, same as for code biases.

Initially, all ambiguity parameters are estimated as float values in the *Estimate solution* steps (see Sect. 2). Integer ambiguity resolution is then performed based on the least squares ambiguity decorrelation adjustment (LAMBDA) method (Teunissen 1995), specifically the modified algorithm (MLAMBDA) by Chang et al. (2005). First the covariance matrix of the integer ambiguity parameters is computed by eliminating all but those parameters from the
full normal equation matrix and inverting it. Then, a Z-transformation is performed as described by Chang et al. (2005) to decorrelate the ambiguity parameters without losing their integer nature.

The search process follows MLAMBDA and uses integer minimization of the weighted sum of squared residuals. As shown in Table 2, the number of ambiguities is usually in the tens of thousands for a global network. It is computationally infeasible to search a hyper-ellipsoid with dimensions that large. Instead, a blocked search algorithm is performed by moving a window with a length of, e.g., 200 parameters over the decorrelated ambiguities, starting from the most accurate. In each step, the window is moved by half of its length and the overlapping parts are compared to each other. If all fixed ambiguities in the overlap agree, the algorithm continues. Otherwise, both windows are combined and the search is repeated using the combined window, again comparing with the overlapping part of the preceding window. Once the algorithm finishes, all ambiguity parameters are fixed to integer values. In contrast to MLAMBDA, it is not guaranteed that the resulting solution is optimal in the sense of minimal variance with given covariance. This trade-off is necessary to cope with the large number of ambiguities.

The ambiguities are then transformed back to the original parameter space following Chang et al. (2005). Integer ambiguities $N_j$ are now considered resolved and are reduced from the observations as if they were known corrections. Only transmitter and receiver phase biases remain and are estimated in subsequent processing iterations, as shown in Fig. 2. Observations associated with ambiguities that were potentially resolved incorrectly will be downweighted by VCE in further processing steps.

# 4 Results and discussion

The raw observation approach was tested and evaluated by processing real data for a period of 15 years, from the beginning of 2003 to the end of 2017. Observations from the IGS14 station network to the GPS constellation were used. A standard dual-frequency processing using GPS code and phase observations on L1 and L2 was chosen to be consistent with IGS products, which were used for evaluation. The full time period was consistently processed using the same settings and the latest available models (see Sect. 3), e.g., IGS14 antenna calibrations. IGS analysis center (AC) and combined products used for evaluation come from the second IGS reprocessing campaign (repro2) until the end of 2013, and from final products afterward. This was chosen to allow for evaluation against time series that contain relatively recent models and processing techniques. IGS repro2 products are based on the IGb08 reference frame and antenna calibrations. IGS ACs transitioned from IGb08 to IGS14 on

<table><caption>**Table 3** IGS analysis centers used for evaluation</caption><thead><tr><td>Label</td><td>Institution</td></tr></thead><tbody><tr><td>COD</td><td>Center for Orbit Determination in Europe</td></tr><tr><td>EMR</td><td>Natural Resources Canada</td></tr><tr><td>ESA</td><td>European Space Operations Centre</td></tr><tr><td>GFZ</td><td>GeoForschungsZentrum</td></tr><tr><td>GRG</td><td>Groupe de Recherche de Géodésie Spatiale</td></tr><tr><td>JPL</td><td>Jet Propulsion Laboratory</td></tr><tr><td>MIT</td><td>Massachusetts Institute of Technology</td></tr><tr><td>SIO</td><td>Scripps Institution of Oceanography</td></tr></tbody></table>

January 29, 2017, for final products. Evaluation focuses on station positions and GPS satellite orbits and clocks. Products processed in the scope of this article are labeled as TUG (Graz University of Technology). IGS ACs used for evaluation are listed in Table 3. IGS combined products are used as reference in the comparisons since they are assumed to be the best approximation of truth.

## 4.1 Satellite orbits

Figure 3 shows daily GPS orbit RMS values of the TUG and IGS AC solutions relative to the IGS combination. Only orbits that were computed by all institutions were used in this comparison. Therefore, satellites experiencing outages according to NANU messages were not considered in the overall RMS for that day, even though some institutions might have processed them anyway. Considering all available orbits from the IGS combination as 100%, TUG covered 96.4% of those orbits. The spread for IGS ACs is 95.2–99.8%, with a mean of 97.7%. Reference frame differences between an individual solution and the IGS combination were corrected by estimating daily sets of Helmert parameters in the inertial frame. They comprise constant and trend components for translation, rotation, and scale to remove any systematic effects, e.g., due to differing Earth orientation parameters.

Overall, the individual solutions differ from the IGS combination by a median RMS of 9.1–14.8 mm, with the TUG solution on the lower spectrum with 9.5 mm. During the repro2 period, especially after 2004, the relative RMS differences between the solutions are stable over time. This was to be expected due to the consistent processing during that period. Most solutions show periodic signals close to annual or GPS draconitic (351.6 days) periods. IGS AC and combined final products were used from the beginning of 2014 onward, which is clearly visible as jumps in the RMS time series. Since the TUG solution was processed consistently over the full time series, its jump at this point most likely originates in the reference IGS combined orbits. The weights of individual AC solutions in the IGS combination
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Daily GPS orbit RMS relative to IGS combination. Limited to orbits computed by all institutions. Reference frame differences corrected (Helmert). 91-day median-filtered for clarity</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Orbit discontinuity RMS at day boundaries per GPS block for TUG solution. Reference frame differences between days corrected (Helmert). 91-day median-filtered for clarity</figcaption></figure>

differ between the repro2 and final periods. This also explains jumps in the time series of the ACs themselves. Further jumps are visible shortly after in April 2014, when ESA changed their solar radiation pressure modeling (Weiss et al. 2017). Again, this affected the relative weighting of the IGS AC solutions and therefore the IGS combined orbits. From 2014 onward, the TUG solution fits increasingly better to the IGS combination. This is most likely related to modeling and processing technique updates in the AC solutions, which reflect into the IGS combination. After the switch to IGS14 on January 29, 2017, TUG fits best to the combined solution, even though it is not part of the combination.

Figure 4 shows the orbit discontinuity RMS at day boundaries. This is a measure for the internal consistency of a solution. It expresses how well consecutive daily arcs agree at midnight epochs. Reference frame differences between consecutive days were corrected by estimating a set of seven Helmert parameters. The discontinuity RMS is displayed per GPS block to see if there are any block-dependent systematic effects, e.g., related to solar radiation pressure. In general, all blocks show an RMS of around 13 mm over the full time series. A small annual/draconitic signal is visible for

most blocks. Draconitic signals are also present in nearly all standard IGS products. Griffiths and Ray (2013) show that other than mismodeling of orbit dynamics and aliasing of near-sidereal local station multipath effects, sub-daily EOP tide errors could also be the source of these draconitic signals. Large peaks for block IIA come from satellite SVN 29, which ceased operations in October 2007. Other than that, there is no major block-dependency noticeable. Evaluation against IGS AC solutions was not possible for this measure since IGS orbit products do not contain overlapping midnight epochs.

## 4.2 Satellite clocks

Daily GPS clock RMS values of the TUG and IGS AC solutions relative to the IGS combination are shown in Fig. 5. The clocks were synchronized, meaning only epochs with clocks available from all institutions were used in the comparison. This implies a 5-min clock sampling rate up to the end of 2014, since not all IGS ACs provided 30-s clocks before that period and for repro2. COD and SIO are not included in the comparison because no clocks were provided for repro2 or in general, respectively. Solving the rank deficiency for clocks results in a system-wide absolute clock shift when comparing to other solutions. This shift was corrected relative to the reference IGS combination for each epoch during RMS computation.

The median RMS of clocks processed in the scope of this article is 40.5 mm. IGS AC solutions range from 27.8 to 70 mm median RMS. It is generally stable over the full period and does not show large outlier periods as some other solutions do. Generally, comparing clocks can be tricky. They are affected, for example, by the transmitter attitude model used in the processing. Inconsistent attitude modeling during eclipse seasons can result in clock discrepancies between solutions. This can also lead to the rejection of some solutions during the IGS combination process (Weiss et al. 2017).
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Daily GPS clock RMS relative to IGS combination. Limited to clocks computed by all institutions. System-wide absolute clock shifts corrected at each epoch. 91-day median-filtered for clarity</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6** Daily station position RMS relative to IGS combination. Note that the station count differs since all IGS14 stations processed by the respective institution were used. Reference frame differences corrected (Helmert). 91-day median-filtered for clarity</figcaption></figure>

## 4.3 Station positions

Figure 6 shows daily station position RMS values of the TUG and IGS AC solutions relative to the IGS combination. Note that the station count differs for each solution since all IGS14 stations processed by the respective institution were used. Limiting to only stations processed by all institutions at the respective day resulted in a very small subset of stations that would not have been representative. Differences in the reference frames between an individual solution and the IGS combination were corrected by estimating a set of seven Helmert parameters per day. Since nontidal atmospheric and ocean mass variations were considered during processing of the TUG solution, the same model (AOD1B RL06) was used to correct the IGS solution during RMS computation between these two solutions.

Overall, the solution processed in the scope of this article is on the average level of IGS AC solutions with a median RMS of 2.7 mm relative to the IGS combination. AC solutions range between 1.8 and 3.9 mm. During the repro2 period, all solutions are again relatively stable over time. The TUG solution shows a slight increase toward the end of the repro2 period. Further investigations are necessary to determine

whether this is possibly related to ionospheric effects due to the Sun's 11-year solar cycle, which had its latest peak around 2013–2014. With the switch from IGb08 to IGS14 by the IGS ACs on January 29, 2017, TUG's RMS drops from 3 to 2 mm. Since the TUG solution was processed using IGS14 for the full time series, this jump originates solely from the better consistency between all solutions from this point onward.

The temporal repeatability of station positions over the full time series was investigated in order to validate the station positions independent of a reference solution. Overall RMS values for TUG, IGS ACs, and the IGS combination are listed in Table 4. They were computed using all IGS14 stations included in the respective solution with at least 30 days of data over the 15 years. All solutions used for evaluation were corrected for nontidal atmospheric and ocean mass variations using the model AOD1B RL06. This was done to reduce temporal variability and to be consistent with TUG. Each station time series was split into intervals using discontinuities determined for the ITRF2014 reference frame. Additional break points were added at the switch from repro2 to final on January 1, 2014, and from IGb08 to IGS14 on January 29, 2017, for all but the TUG solution. Mean, trend, and annual signals
<table><caption>**Table 4** Overall temporal position repeatability for IGS14 stations included in the respective solution</caption><thead><tr><th rowspan="2">Solution</th><th rowspan="2">Stations</th><th colspan="4">Overall temporal RMS (mm)</th><th rowspan="2">% Outliers removed</th></tr><tr><th>Total</th><th>North</th><th>East</th><th>Up</th></tr></thead><tbody><tr><td>COD</td><td>232</td><td>2.7</td><td>1.7</td><td>1.6</td><td>4.0</td><td>3.9</td></tr><tr><td>EMR</td><td>188</td><td>6.4</td><td>6.3</td><td>5.0</td><td>7.5</td><td>1.6</td></tr><tr><td>ESA</td><td>216</td><td>3.4</td><td>2.0</td><td>2.0</td><td>5.1</td><td>4.7</td></tr><tr><td>GFZ</td><td>206</td><td>3.4</td><td>2.0</td><td>1.9</td><td>5.2</td><td>4.0</td></tr><tr><td>GRG</td><td>185</td><td>8.3</td><td>8.4</td><td>5.4</td><td>10.3</td><td>2.6</td></tr><tr><td>JPL</td><td>227</td><td>6.1</td><td>5.9</td><td>4.6</td><td>7.4</td><td>1.5</td></tr><tr><td>MIT</td><td>242</td><td>5.9</td><td>5.9</td><td>4.0</td><td>7.4</td><td>2.0</td></tr><tr><td>SIO</td><td>238</td><td>8.9</td><td>10.2</td><td>4.7</td><td>10.5</td><td>2.3</td></tr><tr><td>IGS</td><td>243</td><td>2.9</td><td>1.6</td><td>1.6</td><td>4.5</td><td>5.3</td></tr><tr><td>TUG</td><td>245</td><td>2.8</td><td>1.5</td><td>1.5</td><td>4.3</td><td>4.7</td></tr></tbody></table>

Mean, trend, and annual signal reduced considering discontinuities. Outliers removed based on robust 3σ-level

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Temporal position repeatability of IGS14 stations for TUG (a) and IGS (b) solutions. Mean, trend, and annual signal reduced considering discontinuities. Outliers removed based on robust 3σ-level</figcaption></figure>

were reduced from each interval. Outliers were then removed for each station based on a robust 3σ-level over the full time series.

The overall RMS for TUG is 2.8 mm, which is slightly lower than the IGS combined solution. This is noteworthy since the IGS combination is a weighted mean of the IGS AC solutions, which is expected to result in increased stability. Only COD has a lower total RMS of 2.7 mm. It has to be noted that 3-day solutions from COD were used in this article, while all others are 1-day solutions. This might explain the lower temporal variability for COD. The percentage of removed outliers is relatively high for TUG with 4.7%, although it is similar to other solutions with a comparable total RMS. It is noticeable that there seem to be two groups: one with total RMS of around 3 mm and about 4.5% outliers removed, and one with total RMS of around 6–8 mm and roughly 2% outliers. Figure 7 shows the RMS values and percentage of removed outliers per station for TUG and the IGS combination. Some stations with a very high percentage of outliers are clearly visible in both solutions.

# 5 Summary and conclusions

The raw observation approach was described in detail, including the handling of clocks and instrumental biases as well as tropospheric and ionospheric delays. A possible implementation of the approach with all necessary steps to estimate a solution was shown. Furthermore, the processing strategy applied at Graz University of Technology was detailed. It uses a multi-step procedure in which initially only a core network is processed to determine all parameters and enable ambiguity resolution. Additional stations can then be processed individually. After resolving all ambiguities, a final full network processing updates all parameters. The approach includes sophisticated handling of different error sources and a realistic observation weighting scheme based on VCE. The description of the raw observation approach and the applied processing strategy should enable others to implement this approach.

Based on the presented results, it is concluded that the raw observation approach is well suited for determining GNSS products like satellite orbits, clocks, and station positions. The quality of resulting products from dual-frequency GPS
processing using code and phase observations is comparable to or better than that of well-established approaches used by other institutions. The main advantage of the raw observation approach is the fact that observations are used directly as observed by the receiver. This preserves the original measurement accuracy and allows full exploitation of the information contained in each individual observation type. New observables (e.g., GPS L5) can be included in a straightforward way because observation equations are set up individually and no linear combinations are formed. This becomes especially important in view of new signals and systems that are being introduced in GNSS.

Since it has been confirmed that the raw observation approach is competitive for a standard dual-frequency GPS processing, the next step is to implement processing capabilities for further GNSS in order to take full advantage of the approach. Challenges of multi-frequency and multi-GNSS processing, like the presence of inter-frequency and inter-system biases, will have to be dealt with. It is expected that the flexibility of the raw observation approach will simplify the parametrization of those and similar time-variable biases. The software implementation is already generalized in a way to facilitate multi-GNSS processing.

**Acknowledgements** Open access funding provided by Graz University of Technology. The presented work was funded by the Austrian Research Promotion Agency (FFG) in the frame of the Austrian Space Applications Programme (ASAP) Phase 12 (Project No.: 854024). We thank the IGS and its associated institutions for providing GNSS products and observation data. The feedback from three anonymous reviewers was much appreciated.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

## References

Arnold D, Meindl M, Beutler G, Dach R, Schaer S, Lutz S, Prange L, Sośnica K, Mervart L, Jäggi A (2015) CODE's new solar radiation pressure model for GNSS orbit determination. J Geod 89(8):775-791. https://doi.org/10.1007/s00190-015-0814-4

Bar-Sever YE, Kroger PM, Borjesson JA (1998) Estimating horizontal gradients of tropospheric path delay with a single GPS receiver. J Geophys Res Solid Earth 103(B3):5019-5035. https://doi.org/10.1029/97JB03534

Böhm J, Werl B, Schuh H (2006) Troposphere mapping functions for GPS and very long baseline interferometry from European Centre for Medium-Range Weather Forecasts operational analysis data. J Geophys Res Solid Earth. https://doi.org/10.1029/2005JB003629

Böhm J, Möller G, Schindelegger M, Pain G, Weber R (2015) Development of an improved empirical model for slant delays in the troposphere (GPT2w). GPS Solut 19(3):433-441. https://doi.org/10.1007/s10291-014-0403-7

Carrere L, Lyard F, Cancet M, Guillot A, Picot N (2016) FES 2014, a new tidal model - validation results and perspectives for improvements. Presentation at ESA living planet symposium 2016, Prague, 9-13 May 2016

Chang XW, Yang X, Zhou T (2005) MLAMBDA: a modified LAMBDA method for integer least-squares estimation. J Geod 79(9):552-565. https://doi.org/10.1007/s00190-005-0004-x

Condat L (2013) A direct algorithm for 1-D total variation denoising. IEEE Signal Process Lett 20(11):1054-1057. https://doi.org/10.1109/LSP.2013.2278339

Dach R, Brockmann E, Schaer S, Beutler G, Meindl M, Prange L, Bock H, Jäggi A, Ostini L (2009) GNSS processing at CODE: status report. J Geod 83(3-4):353-365. https://doi.org/10.1007/s00190-008-0281-2

Dach R, Lutz S, Walser P, Fridez P (eds) (2015) Bernese GNSS software version 5.2. University of Bern, Bern Open Publishing, Bern. https://doi.org/10.7892/boris.72297

Dilssner F (2010) GPS IIF-1 satellite-antenna phase center and attitude modeling. Inside GNSS 5(6):59-64

Dobslaw H, Bergmann-Wolf I, Dill R, Poropat L, Thomas M, Dahle C, Esselborn S, König R, Flechtner F (2017) A new high-resolution model of non-tidal atmosphere and ocean mass variability for de-aliasing of satellite gravity observations: AOD1B RL06. Geophys J Int 211(1):263-269. https://doi.org/10.1093/gji/ggx302

Folkner WM, Williams JG, Boggs DH (2009) The planetary and lunar ephemeris DE 421. The interplanetary network progress report, vol 42(178), pp 1-34. http://ipnpr.jpl.nasa.gov/progress_report/42-178/178C.pdf. Accessed 19 Dec 2017

Fritsche M, Dietrich R, Knöfel C, Rülke A, Vey S, Rothacher M, Steigenberger P (2005) Impact of higher-order ionospheric terms on GPS estimates. Geophys Res Lett 32(23):L23-311. https://doi.org/10.1029/2005GL024342

Griffiths J, Ray JR (2013) Sub-daily alias and draconitic errors in the IGS orbits. GPS Solut 17(3):413-422. https://doi.org/10.1007/s10291-012-0289-1

Guo J, Xu X, Zhao Q, Liu J (2016) Precise orbit determination for quad-constellation satellites at Wuhan University: strategy, result validation, and comparison. J Geod 90(2):143-159. https://doi.org/10.1007/s00190-015-0862-9

Håkansson M, Jensen ABO, Horemuz M, Hedling G (2017) Review of code and phase biases in multi-GNSS positioning. GPS Solut 21(3):849-860. https://doi.org/10.1007/s10291-016-0572-7

Hauschild A (2017a) Basic observation equations. In: Teunissen PJ, Montenbruck O (eds) Springer handbook of global navigation satellite systems, Springer handbooks. Springer, Cham, pp 561-582. https://doi.org/10.1007/978-3-319-42928-1_19

Hauschild A (2017b) Combinations of observations. In: Teunissen PJ, Montenbruck O (eds) Springer handbook of global navigation satellite systems, Springer handbooks. Springer, Cham, pp 583-604. https://doi.org/10.1007/978-3-319-42928-1_20

Hofmann-Wellenhof B, Lichtenegger H, Wasle E (2008) GNSS - Global navigation satellite systems. Springer, Vienna. https://doi.org/10.1007/978-3-211-73017-1

Hoque MM, Jakowski N (2008) Estimate of higher order ionospheric errors in GNSS positioning. Radio Sci. https://doi.org/10.1029/2007RS003817

Huber PJ (1981) Robust statistics. Wiley series in probability and statistics. Wiley, Hoboken. https://doi.org/10.1002/0471725250

Hugentobler U, Montenbruck O (2017) Satellite orbits and attitude. In: Teunissen PJ, Montenbruck O (eds) Springer handbook of global navigation satellite systems, Springer handbooks. Springer, Cham, pp 59-90. https://doi.org/10.1007/978-3-319-42928-1_3

IGS (2011) Calculated and estimated GPS transmit power levels. http://acc.igs.org/orbits/thrust-power.txt. Accessed 13 Feb 2018

Johnston G, Riddell A, Hausler G (2017) The International GNSS Service. In: Teunissen PJ, Montenbruck O (eds) Springer hand-
book of global navigation satellite systems, Springer handbooks.
Springer, Cham, pp 967–982. https://doi.org/10.1007/978-3-319-42928-1_33

Koch KR (1999) Parameter estimation and hypothesis testing in linear models. Springer, Berlin. https://doi.org/10.1007/978-3-662-03976-2

Kouba J (2009a) A guide to using International GNSS Service (IGS) products. http://kb.igs.org/hc/en-us/article attachments/203088448/UsingIGSProductsVer21_cor.pdf

Kouba J (2009b) A simplified yaw-attitude model for eclipsing GPS satellites. GPS Solut 13(1):1–12. https://doi.org/10.1007/s10291-008-0092-1

Loyer S, Perosanz F, Mercier F, Capdeville H, Marty JC (2012) Zero-difference GPS ambiguity resolution at CNES/CLS IGS analysis center. J Geod 86(11):991–1003. https://doi.org/10.1007/s00190-012-0559-2

Mayer-Gürr T, Pail R, Gruber T, Fecher T, Rexer M, Schuh W-D, Kusche J, Brockmann J-M, Rieser D, Zehentner N, Kvas A, Klinger B, Baur O, Höck E, Krauss S, Jäggi A (2015) The combined satellite gravity field model GOCO05s. Presentation at EGU general assembly 2015, Vienna, 12–17 April 2015. https://doi.org/10.13140/RG.2.1.4688.6807

Melbourne WG (1985) The case for ranging in GPS-based geodetic systems. In: Proceedings of the first international symposium on precise positioning with the Global Positioning System. National Geodetic Survey (NGS), Rockville, pp 373–386

Misra P, Enge P (2011) Global Positioning System: signals, measurements, and performance, 2nd edn. Ganga-Jamuna Press, Lincoln

Montenbruck O, Gill E (2000) Satellite orbits. Springer, Berlin. https://doi.org/10.1007/978-3-642-58351-3

Montenbruck O, Steigenberger P, Khachikyan R, Weber G, Langley R, Mervart L, Hugentobler U (2014) IGS-MGEX: preparing the ground for multi-constellation GNSS science. Inside GNSS 9(1):42–49

Montenbruck O, Steigenberger P, Prange L, Deng Z, Zhao Q, Perosanz F, Romero I, Noll C, Stürze A, Weber G, Schmid R, MacLeod K, Schaer S (2017) The multi-GNSS experiment (MGEX) of the International GNSS Service (IGS) – Achievements, prospects and challenges. Adv Space Res 59(7):1671–1697. https://doi.org/10.1016/j.asr.2017.01.011

Petit G, Luzum B (eds) (2010) IERS conventions. IERS Technical Note 36. Verlag des Bundesamts für Kartographie und Geodäsie, Frankfurt am Main. http://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn36.html. Accessed 19 Dec 2017

Prange L, Orliac E, Dach R, Arnold D, Beutler G, Schaer S, Jäggi A (2017) CODE’s five-system orbit and clock solutionthe challenges of multi-GNSS data analysis. J Geod 91(4):345–360. https://doi.org/10.1007/s00190-016-0968-8

Rodriguez Solano CJ (2009) Impact of albedo modelling on GPS orbits. Master’s thesis, Technische Universität München, München. https://mediatum.ub.tum.de/doc/1368717/file.pdf. Accessed 19 Dec 2017

Rodriguez Solano CJ (2014) Impact of non-conservative force modeling on GNSS satellite orbits and global solutions. Ph.D thesis. Technische Universität München, München. https://mediatum.ub.tum.de/doc/1188612/file.pdf. Accessed 19 Dec 2017

Saastamoinen J (1972) Atmospheric correction for the troposphere and stratosphere in radio ranging satellites. In: Henriksen SW, Armando M, Chovitz BH (eds) The use of artificial satellites for geodesy. American Geophysical Union, Washington, D.C., pp 247–251. https://doi.org/10.1029/GM015p0247

Schaer S (1999) Mapping and predicting the earth’s ionosphere using the Global Positioning System. Ph.D thesis. University of Bern. https://www.researchgate.net/publication/252260542_Mapping_and_Predicting_the_Earth's_Ionosphere_Using_the_GlobalPOSITIONING_SYSTEM. Accessed 7 May 2018

Schmid R, Rothacher M, Thaller D, Steigenberger P (2005) Absolute phase center corrections of satellite and receiver antennas. GPS Solut 9(4):283–293. https://doi.org/10.1007/s10291-005-0134-x

Schmid R, Steigenberger P, Gendt G, Ge M, Rothacher M (2007) Generation of a consistent absolute phase-center correction model for GPS receiver and satellite antennas. J Geod 81(12):781–798. https://doi.org/10.1007/s00190-007-0148-y

Schönemann E (2014) Analysis of GNSS raw observations in PPP solutions. Schriftenreihe der Fachrichtung Geodäsie 42. http://tuprints.ubl.tu-darmstadt.de/3843/. Accessed 27 Nov 2018

Schönemann E, Becker M, Springer T (2011) A new approach for GNSS analysis in a multi-GNSS and multi-signal environment. J Geod Sci 1(3):204–214. https://doi.org/10.2478/v10156-010-0023-2

Teunissen PJG (1995) The least-squares ambiguity decorrelation adjustment: a method for fast GPS integer ambiguity estimation. J Geod 70(1–2):65–82. https://doi.org/10.1007/BF00863419

Uhlemann M, Gendt G, Ramatschi M, Deng Z (2015) GFZ global multi-GNSS network and data processing results. In: Rizos C, Willis P (eds) IAG 150 years, International Association of Geodesy symposia, vol 143, pp 673–679. Springer, Cham. https://doi.org/10.1007/1345_2015_120

Weiss JP, Steigenberger P, Springer T (2017) Orbit and clock product generation. In: Teunissen PJ, Montenbruck O (eds) Springer handbook of global navigation satellite systems, Springer handbooks, pp 983–1010. Springer, Cham. https://doi.org/10.1007/978-3-319-42928-1_34

Wübbena G (1985) Software developments of geodetic positioning with GPS using TI 4100 code and carrier measurements. In: Proceedings of the first international symposium on precise positioning with the Global Positioning System. National Geodetic Survey (NGS), Rockville, pp 403–412

Zehentner N, Mayer-Gürr T (2014) New approach to estimate time variable gravity fields from high-low satellite tracking data. In: Marti U (ed) Gravity, geoid and height systems. International Association of Geodesy symposia, vol 141, pp 111–116. Springer, Cham. https://doi.org/10.1007/978-3-319-10837-7_14

Zehentner N, Mayer-Gürr T (2016) Precise orbit determination based on raw GPS measurements. J Geod 90(3):275–286. https://doi.org/10.1007/s00190-015-0872-7