PPP-based Swarm kinematic orbit determination

The Swarm mission of the European Space Agency (ESA) offers excellent opportunities to study the ionosphere and to provide temporal gravity field information for the gap between the Gravity Recovery and Climate Experiment (GRACE) and its follow-on mission (GRACE-FO). In order to contribute to these studies, at the Institut für Erdmessung (IfE) Hannover, a software based on precise point positioning (PPP) batch least-squares adjustment is developed for kinematic orbit determination. In this paper, the main achievements are presented.

The approach for the detection and repair of cycle slips caused by ionospheric scintillation is introduced, which is based on the Melbourne–Wübbena and ionosphere-free linear combination. The results show that around 95 % of cycle slips can be repaired and the majority of the cycle slips occur on L₂. After the analysis and careful preprocessing of the observations, 1-year kinematic orbits of Swarm satellites from September 2015 to August 2016 are computed with the PPP approach. The kinematic orbits are validated with the reduced-dynamic orbits published by the ESA in the Swarm Level 2 products and SLR measurements. The differences between IfE kinematic orbits and ESA reduced-dynamic orbits are at the 1.5, 1.5 and 2.5 cm level in the along-track, cross-track and radial directions, respectively. Remaining systematics are characterized by spectral analyses, showing once-per-revolution period. The external validation with SLR measurements shows RMSEs at the 4 cm level. Finally, fully populated covariance matrices of the kinematic orbits obtained from the least-squares adjustment with 30, 10 and 1 s data rate are discussed. It is shown that for data rates larger than 10 s, the correlation between satellite positions should be taken into account, for example, for the recovery of gravity field from kinematic orbits.

1 Introduction

The Swarm mission was launched on 22 November 2013 and is the first constellation of satellites of the European Space Agency (ESA) to study the dynamics of the Earth's magnetic field and its interaction with the Earth system (Friis-Christensen et al., 2008). This mission consists of three identical satellites in near-polar low Earth orbits (LEO). The two satellites, Swarm A and C, fly almost side-by-side at an initial altitude of 460 km while the third Swarm satellite B flies in a higher orbit of about 530 km. All the three satellites are equipped with a set of six core instruments: Absolute Scalar Magnetometer (ASM), Vector Field Magnetometer (VFM), Electric Field Instrument (EFI), Star Tracker (STR), Accelerometer (ACC) and GPS Receiver (GPSR). In order to take full advantage of the data information provided by this constellation, Precise Orbit Determination (POD) is necessary, which is obtained with data from the high precision eight-channel dual-frequency GPS receiver. In addition, each satellite has a laser retro-reflector, which makes the independent validation of the GPS-derived orbits with satellite laser ranging (SLR) possible. In recent years, kinematic orbit determination has attracted much attention (Jäggi et al., 2016; Zehentner and Mayer-Gürr, 2015; Švehla and Rothacher, 2003). Compared to the reduced-dynamic approach, the kinematic approach is a purely geometrical approach without using any information on LEO satellite dynamics (e.g. gravity field, air drag, etc.). Therefore, kinematic orbits are preferred for independent gravity field determination.

A first performance evaluation of the RUAG 8-channel GPS receiver was carried out by Zangerl et al. (2014), showing unexpected results and loss of lock due to ionospheric scintillations. Buchert et al. (2015) and Xiong et al. (2016) empirically linked the loss of the GPS signal to ionospheric plasma irregularities. Results from more than 1 year of data using reduced-dynamic and kinematic approaches were reported in van den IJssel et al. (2015) and Jäggi et al. (2016). A comparison with SLR allows an assessment of the orbit accuracy of about 2.5 cm rms for ESA reduced-dynamic orbits and 4 cm for ESA kinematic orbits. The 3-D differences between the reduced-dynamic and kinematic orbits are at the 4–5 cm level, with larger differences close to the geomagnetic poles and along the geomagnetic Equator (van den IJssel et al., 2015). An azimuth and elevation dependent PCV map is generated during flight using the residual approach (Jäggi et al., 2016; van den IJssel et al., 2015). In order to collect more observations and improve the robustness under ionospheric scintillation, the receiver settings have been adjusted several times during the Swarm mission. The impact of the first update on POD performance is analysed in van den IJssel et al. (2016), showing that a lower GPS elevation cut-off angle and wider bandwidths of the tracking loops improve the performance for both orbit types. Gravity field recovery also benefits from the update (Dahle et al., 2017).

In this contribution, the developments for kinematic POD of Swarm satellites at the Institut für Erdmessung (IfE) are reported. This approach is based on the precise point positioning (PPP) technique (Zumberge et al., 1997). This differs from the raw measurements approach by the Institute of Geodesy (IfG) of TU Graz (Zehentner and Mayer-Gürr, 2015), phase-only approach by the Astronomical Institute of University Bern (AIUB) (Jäggi et al., 2016) or the Bayesian weighted least-squares estimator, which is implemented in the GPS High-precision Orbit Determination Software Tools (GHOST) developed at the Deutsches Zentrum für Luft- und Raumfahrt in close cooperation with TU Delft and used for the ESA official orbits (van den IJssel et al., 2015).

Since the GPS data quality is decisive for the obtainable orbit accuracy, we first analyse the in-flight performance of the Swarm onboard GPS receivers in Sect. 2, e.g. the tracking performance and observation noise, especially under the influence of ionospheric scintillation. A mandatory step in the preprocessing is detecting and, if possible, repairing the cycle slips in carrier phase observations in order to reduce the number of ambiguities and therefore strengthen the orbit. Due to the large phase noise caused by ionospheric scintillations, cycle slips in Swarm receivers are difficult to detect at ionospheric-active areas. An approach based on the Melbourne–Wübbena and ionosphere-free linear combination is introduced for this task. In Sect. 3, we introduce the adopted method and models for precise kinematic orbit determination by PPP in a least-squares adjustment. The obtained orbit as well as its fully populated formal covariance matrix are presented. Unlike the CHAllenging Minisatellite Payload (CHAMP) or the Gravity Recovery and Climate Experiment (GRACE), the Swarm GPS receiver provides observations every second since 15 July 2014, which gives us the opportunity to study the temporal correlation of the coordinates with 1 s observation rate. The kinematic orbits are compared with the ESA reduced-dynamic orbits and SLR measurements. Thanks to a more strict data screening strategy and cycle-slip detection, our kinematic orbits show a smaller RMSE compared to ESA kinematic orbits. Section 4 gives a short conclusion of this paper.

Table 1 Update of Swarm carrier loop bandwidth

2 Swarm GPS data qualities

Being a purely geometrical approach without using any dynamic models on LEO satellites, the quality of the kinematic orbits relies completely on the GPS observations as well as on the introduced GPS orbit and clock products, and it is sensitive to measurement errors and the observation geometry (Švehla and Rothacher, 2003). High-quality data is the prerequisite for precise kinematic orbits. Therefore, in this section, we first analyse the Swarm GPS receiver data quality and introduce our approach for cycle-slip detection and repair as well as outlier screening.

2.1 Tracking performance

The Swarm onboard dual-frequency GPS receivers are developed by RUAG Space (Zangerl et al., 2014). They have 8 channels for tracking GPS satellites, which is less than the GRACE and Gravity field and steady-state Ocean Circulation Explorer (GOCE) receivers, with 10 and 12 channels, respectively. In order to collect more observations, the field of view was increased step by step from 10^∘ elevation cut-off angle at the beginning of the mission to 2^∘ cut-off angle since 6 May 2015 for all three Swarm satellites (van den IJssel et al., 2015). Several updates on the bandwidth of the tracking loops were applied to improve the robustness of the tracking in difficult environmental situations. The fast reacquisition of the L₁-carrier tracking was enabled by raising the retry counter from zero to five (ESA, 2015). These updates are shortly summarized in Table 1 (Dahle et al., 2017).

The daily average number of tracked GPS satellites was increased from 7.6 to 7.7 after the first update of receiver bandwidth (van den IJssel et al., 2016). For LEO satellites, the number of satellites in view is usually larger than 8. But due to the limitation of receiver channels, only 8 GPS satellites maximum can be tracked. There are two reasons for epochs within which less than 7 GPS satellites are tracked.

The first reason is the reacquisition time of the Swarm receiver. When a GPS satellite sets, the receiver cannot immediately track another GPS satellite in view and needs around 84 s before the first receiver update (red) and 55 s after the update to acquire another satellite (blue), shown in Fig. 1a, for example, for DoY 250, 2015. Thus, during the reacquisition phase, the number of satellites is less than 8. If some satellites set during short time together, the number of tracked satellites is even less than 7. This shortened reacquisition time is also the main reason for the increased number of tracked GPS satellites.

Another reason for tracking less than 7 GPS satellites is the loss of the lock of signal. There are also two reasons for the loss of the GPS signal, and one is due to the weak signal strength for the GPS satellites at low elevations. However, the loss of lock on GPS satellites at high elevations occurs mainly at polar and equatorial areas under the influence of ionospheric scintillations, which is evidenced by the change of the Total Electron Content (TEC). The occurrence of loss of lock in September 2015 is shown in Fig. 1b, for example.

Although the amount of the loss of lock on GPS satellites at higher elevations was decreased from 37 to 11 after the update, the receiver becomes unstable for the signals from lower elevations. The amount of the loss of lock increased almost fourfold, from 58 to 256. The reacquisition time in this situation was decreased from around 17 to 11 s after the update.

Figure 1 Reduction of reacquisition time for Swarm GPS receivers (a) due to setting and rising satellites on DoY 250, 2015 and (b) due to loss of lock in September 2015, red for Swarm A and blue for Swarm C with different tracking settings.

Figure 2 (a)P₁ code noise with respect to elevation for Swarm A on DoY 100 before fixing the RINEX converter issue and (b) DoY 104, 2016 after fixing the issue.

Figure 3 (a)Δ₂L₁ and (b)Δ₂L₂ noise for Swam A from 18:00 to 24:00 on DoY, 333, 2015, colour coded per PRN (from blue PRN1 to red PRN 32). The time series are intentionally shifted by PRN × 1 cm. The local time is 08:30 for the descending arc and 20:30 for the ascending arc. The right legend shows the geographic latitude.

2.2 Observation analysis

As a dual-frequency receiver, the Swarm GPS receiver provides code observations C∕A (C1C), P₁ (C1P) and P₂ (C2P) and carrier phase observations L₁ (L1C) and L₂ (L2P) in the Swarm Level 1b product in the RINEX 3.0 format. A solid analysis of the observation noise is mandatory for an adequate stochastic model. Typically, models have identical weights or elevation dependent weights, e.g. 1∕sin(Elev) or $\mathrm{1}/\sqrt{\sin \left(\mathrm{Elev}\right)}$. The code accuracy is analysed with the multipath linear combination (Estey and Meertens, 1999), which contains multipath effects, ambiguities, code and phase noise. Since carrier phase noise is smaller than code noise by two or three orders of magnitude, it can be neglected. Ambiguities are subtracted as mean value. Then, only code noise and multipath effects remain.

Due to a software issue in the RINEX converter (ESA, 2016), the code observations in the Level 1b product before 11 April 2016 (DoY 102) suffer from high noise. The P₁ code noise for all the GPS satellites with respect to the elevation for Swarm A on DoY 100, 2016 and DoY 104, 2016 is shown in Fig. 2 as an example. A large code noise of about 1.3 m (1σ) can be seen in the left figure even at high elevations, which is caused by the software issue. After fixing the converter issue, the code noise shows a strong elevation-dependence. The P₁ code is disturbed by noise of about 0.8 m (1σ) for satellites at elevations lower than 15^∘. The noise is smaller than 0.2 m (1σ) for the satellites above 30^∘. A 1∕sin(Elev) curve is used here to approximate the noise behaviour. The standard deviations (1σ) of C∕A code and P₂ code for the satellites above 15^∘ are around 0.33 m and 0.21 m, respectively, after fixing the converter issue. The C∕A code noise measured onboard based on 1 Hz data rate is twice larger than the noise (0.16 m) measured on the ground with double-differences on a zero baseline experiment (i.e. two receivers are connected to one antenna via a signal splitter, but averaged over 10 s), reported in Zangerl et al. (2014). The P₂ code noise with 0.21 m onboard is smaller than the noise in the on-the-ground test with 0.31 m.

The quality of the kinematic orbit determination depends mainly on the available dual-frequency carrier phase observations. To analyse the phase noise, the second-order differences (${\mathrm{\Delta }}_{\mathrm{2}}L\left(t\right)=L(t+\mathrm{1})-\mathrm{2}L\left(t\right)+L(t-\mathrm{1})$) of successive phase observations are used after removing the geometric distances, which are computed with the Swarm Level 2 reduced-dynamic orbits. The receiver-clock offset and drift is removed during differencing. Due to the differencing process involved, this is a qualitative evaluation method rather than a quantitative one.

Figure 4 (a)STEC and (b)its rate for Swam A from 18:00 to 24:00 on DoY, 333, 2015. The right legend shows the geographic latitude.

Figure 5 Carrier phase noise for Δ₂L₁ and Δ₂L₂ at midlatitude areas on DoY 333, 2015, plotted against elevation.

The resulting Δ₂L₁ and Δ₂L₂ phase time series of Swarm A from 18:00 to 24:00 in the GPS time system on DoY 333, 2015 for all the GPS satellites are shown in Fig. 3 with respect to time. In addition, the geographic latitude is given. In order to know whether the phase observations are influenced by ionospheric scintillations, the absolute slant total electron content (STECF), which is defined as integrated electron density along the line of sight from Swarm to GPS satellites and provided in the Swarm Level 2 product, with its rate at Swarm A is shown in Fig. 4. The absolute STEC is derived from the geometry-free linear combination of L₁ and L₂, after corrections for differential code biases of GPS satellite transmitters and Swarm receivers (Swarm TEC Product, 2017).

Figures 3 and 4 highlight the variability of the phase noise with respect to the rate of STEC. Carrier phases on both frequencies are simultaneously disturbed under the influence of ionospheric scintillations for all GPS satellites, even at high elevations. Similar large fluctuations also exist in the second-order differences of the ionosphere-free linear combination Δ₂L₃. The widening of the carrier tracking loops during the updates reduces the impact of ionospheric fluctuations significantly (van den IJssel et al., 2016). The large phase noise, even at decimeter level, can almost always be found at high geographic latitude areas (above 60^∘, see red boxes in Fig. 3), for example, where the STEC is low but varies rapidly. At midlatitude areas (20–60^∘), the ionosphere is less active, and these effects are rare. At equatorial areas, the impact depends on the local time and the activity of the local ionosphere. Exemplarily, from 18:00 to 20:00, the STEC is high but changes slowly, which does not disturb the phase observations much. On the contrary, after 20:00 in ascending arcs, with local time around 20:30 after sunset, the changes in STEC are rapid and the phase observations are strongly perturbed (compare black boxes in Fig. 3). This large noise degrades the accuracy of the kinematic orbit significantly and makes it difficult to detect small cycle slips with standard methods.

The Δ₂L₁ and Δ₂L₂ phase noise at midlatitude areas with respect to elevation is shown in Fig. 5. At a low level of ionospheric activity, the phase noise can reach mm level and the noise on L₁ is slightly lower than on L₂. The standard deviations (1σ) of Δ₂L₁ and Δ₂L₂ are 3.9 and 4.2 mm, respectively. There is no strong dependence on the elevation. Assuming independent observations, the standard deviation (1σ) of the undifferenced phase observations can be assessed to be about 1.6 mm for L₁ and 1.7 mm for L₂, according to the laws of error propagation.

2.3 Cycle-slip detection and repair

Strong ionospheric scintillations cause not only large noise but also cycle slips in the carrier phase observations. Therefore, a proper cycle-slip detection is necessary to set up correctly the ambiguities. Since the typical GPS satellite visibility for LEOs is mostly 40 min, the necessity of estimating additional ambiguities for even shorter segments will further decrease the geometric strength of the positioning.

A classical method for cycle-slip detection is to use the TurboEdit method (Blewitt, 1993) based on the Melbourne–Wübbena linear combination (Melbourne, 1985; Wübbena, 1985). This is represented as

$$\begin{array}{}\text{(1)}& L_{\mathrm{MW}}={\displaystyle \frac{f_{\mathrm{1}}{L}_{\mathrm{1}}-f_{\mathrm{2}}{L}_{\mathrm{2}}}{f_{\mathrm{1}}-f_{\mathrm{2}}}}-{\displaystyle \frac{f_{\mathrm{1}}{P}_{\mathrm{1}}+f_{\mathrm{2}}{P}_{\mathrm{2}}}{f_{\mathrm{1}}+f_{\mathrm{2}}}}={\mathit{\lambda }}_{\mathrm{w}}(b_{\mathrm{1}}-b_{\mathrm{2}}),\end{array}$$

where L_MW represents the Melbourne–Wübbena linear combination, f₁ and f₂ the carrier frequency, L₁ and L₂ the carrier phase observation in meters, b₁ and b₂ the ambiguities and λ_w the wide-lane wavelength.

However, the code observations before DoY 102, 2016 contain large noise due to the RINEX converter issue (ESA, 2016), see Fig. 2. The corresponding standard deviation (1σ) of the Melbourne–Wübbena combination is about 0.8 cycles (wide-lane wavelength). The threshold used in the TurboEdit method for the cycle-slip detection is typically 4σ, which means that cycle slips smaller than 3 cycles cannot be detected, see Fig. 6a.

Thus, a forward and backward moving window averaging (FBMWA) method is applied to reduce the influence of large noise (Cai et al., 2013). In this algorithm, the wide-lane ambiguity is smoothed in both forward and backward directions with a specified size m of the smoothing window in each direction. This differs from the regular TurboEdit algorithm, where only the backward smoothing is performed and the window size continuously grows with the number of epochs. For Swarm satellites, m is set as 50 in order to reduce the noise of wide-lane ambiguity to 0.1 cycles. For example, PRN23 is shown in Fig. 6b, the cycle slip is the difference between the forward moving window average at epoch k and the backward moving window average at epoch k−1, which is represented as follows:

$$\begin{array}{ll}{\displaystyle }{c}_{\mathrm{1}}& {\displaystyle }={\displaystyle \frac{\mathrm{\Delta }{L}_{\mathrm{MW}}}{{\mathit{\lambda }}_{\mathrm{w}}}}={\displaystyle \frac{\mathrm{1}}{m{\mathit{\lambda }}_{\mathrm{w}}}}\left(\sum _{i=k}^{k+m-\mathrm{1}}{L}_{\mathrm{MW}}\right(i)-\sum _{i=k-\mathrm{1}}^{k-m}{L}_{\mathrm{MW}}(i\left)\right)\\ \text{(2)}& {\displaystyle }& {\displaystyle }=\mathrm{\Delta }{b}_{\mathrm{1}}-\mathrm{\Delta }{b}_{\mathrm{2}},\end{array}$$

where Δ is the operator for epoch differencing, and c₁=1 can be identified in the peak close to an integer (difference smaller than 0.1) in the black curve in this example.

Figure 6 Illustration of the proposed cycle-slip detection approach based on Melbourne–Wübbena combination, PRN23 of Swarm A on DoY 333, 2015: (a) TurboEdit method, the cycle slip cannot be detected; (b) FBMWA (forward and backward moving window averaging) method $\mathrm{\Delta }{b}_{\mathrm{1}}-\mathrm{\Delta }{b}_{\mathrm{2}}=\mathrm{1}$.

Since with the Melbourne–Wübbena linear combination only the difference of cycle slips on L₁ and L₂ can be detected and simultaneous occurring of same magnitude cycle slips on L₁ and L₂ are not detectable, another linear combination is still required, in order to determine the cycle slips on both frequencies and repair them. Due to the rapid variations of the ionospheric delays and the large phase noise at ionosphere active areas, the geometry-free linear combination is not suitable for Swarm receivers. Here, the time-differenced ionosphere-free linear combination after removing the difference of geometric distances and receiver-clock errors is thus used so that the rapid and large variations of the ionospheric delay can be avoided. This can be represented by the following:

$$\begin{array}{ll}{\displaystyle }\mathrm{\Delta }{L}_{\mathrm{3}}& {\displaystyle }={\displaystyle \frac{{f_{\mathrm{1}}}^{\mathrm{2}}}{{f_{\mathrm{1}}}^{\mathrm{2}}-{f_{\mathrm{2}}}^{\mathrm{2}}}}\mathrm{\Delta }{L}_{\mathrm{1}}-{\displaystyle \frac{{f_{\mathrm{2}}}^{\mathrm{2}}}{{f_{\mathrm{1}}}^{\mathrm{2}}-{f_{\mathrm{2}}}^{\mathrm{2}}}}\mathrm{\Delta }{L}_{\mathrm{2}}\\ \text{(3)}& {\displaystyle }& {\displaystyle }=\mathrm{\Delta }\mathit{\rho }+\mathrm{\Delta }\mathit{\delta }{t}_{\mathrm{r}}+{\displaystyle \frac{c{f}_{\mathrm{1}}}{{f_{\mathrm{1}}}^{\mathrm{2}}-{f_{\mathrm{2}}}^{\mathrm{2}}}}\mathrm{\Delta }{b}_{\mathrm{1}}-{\displaystyle \frac{c{f}_{\mathrm{2}}}{{f_{\mathrm{1}}}^{\mathrm{2}}-{f_{\mathrm{2}}}^{\mathrm{2}}}}\mathrm{\Delta }{b}_{\mathrm{2}},\end{array}$$

where L₃ represents the ionosphere-free combination of carrier phase (other observation errors are corrected), ρ the geometric distance between GPS satellite and receiver position, and δt_r the receiver-clock error in meters where the coefficients before Δb₁ and Δb₂ are approximately 0.4844 and 0.3775 m, respectively.

The geometric distances are computed with the Medium Precise Orbit (MEO) provided by ESA in the Swarm Level 1b product and removed from the right-hand side of Eq. (3). The prerequisite for this approach is that there is no large jump in the a priori orbit.

After removing the geometric distances term, the differences of receiver-clock error and ambiguity on L₁ and L₂ still remain. If there is no cycle slip in the carrier phase, the ambiguity term should be similar to zero-mean white noise and the difference of receiver-clock error can be computed as the average over all tracked GPS satellites, as it is the same for all receiver channels. If there are cycle slips for one GPS satellite, the deviation between the average and this satellite is much larger than the other normal satellites. So the satellite with cycle slips can be detected and the difference of receiver-clock error should be computed again without this one until no large deviation can be found and removed from the right-hand side of Eq. (3).

After removing the geometric distances and receiver-clock error, only the ambiguity term remains. If there is no cycle slip, the ΔL₃ time series spreads around zero in accordance with the phase noise. When a cycle slip occurs, a large jump occurs in the time series. In order to determine the cycle slip exactly and distinguish it from outliers, this jump should be computed as accurately as possible. Instead of calculating the epoch difference of two successive epochs (Fig. 7a), the epoch difference of two epochs separated by n epochs (Fig. 7b) are calculated and then the mean value of these n-epochs differences is formed so that the influence of the large noise under ionospheric scintillation can be reduced (see Fig. 7b). Considering the large noise of ionosphere-free linear combination under ionospheric scintillations (around 10 cm), we propose to select n=100 in order to reduce the noise to 1 cm level. Taking the example of PRN23, we can calculate the following: $c_{\mathrm{2}}=\mathrm{0.4844}\mathrm{\Delta }{b}_{\mathrm{1}}-\mathrm{0.3775}\mathrm{\Delta }{b}_{\mathrm{2}}$, here c₂≈0.38 m.

Figure 7 Illustration of the proposed cycle-slip detection approach based on an ionosphere-free linear combination from PRN23 of Swarm A on DoY 333, 2015. The following are shown: (a) ΔL₃ between two successive epochs $\mathrm{0.29}=\mathrm{0.4844}\mathrm{\Delta }{b}_{\mathrm{1}}-\mathrm{0.3775}\mathrm{\Delta }{b}_{\mathrm{2}}$ and (b) ΔL₃ between two epochs separated by 100 epochs $\mathrm{0.38}=\mathrm{0.4844}\mathrm{\Delta }{b}_{\mathrm{1}}-\mathrm{0.3775}\mathrm{\Delta }{b}_{\mathrm{2}}$.

Together with the equation from Melbourne–Wübbena linear combination, the cycle slip on L₁ and L₂ can be determined:

$$\begin{array}{}\text{(4)}& \left[\begin{array}{c}\mathrm{\Delta }{b}_{\mathrm{1}}\\ \mathrm{\Delta }{b}_{\mathrm{2}}\end{array}\right]={\left[\begin{array}{cc}\mathrm{1}& -\mathrm{1}\\ \mathrm{0.4844}& -\mathrm{0.3775}\end{array}\right]}^{-\mathrm{1}}⌈\begin{array}{c}{c}_{\mathrm{1}}\\ c_{\mathrm{2}}\end{array}⌉.\end{array}$$

In this example we can calculate the following: $\mathrm{\Delta }{b}_{\mathrm{1}}\approx \mathrm{0.023}\approx \mathrm{0}$ and $\mathrm{\Delta }{b}_{\mathrm{2}}\approx -\mathrm{1}$. If the difference of Δb₁ to its nearest integer is smaller than 0.2, it can be rounded and marked as “repaired”.

The statistics of the repaired cycle slips for Swarm satellites during 1 year from DoY 245, 2015 to DoY 244, 2016 are listed in Table 2. Around 95 % of the cycle slips could be repaired. Issues occur by rapid sequences of cycle slips in short time due to strong ionospheric scintillations. The observations during this period can be excluded as outliers. The majority of the cycle slips occur on L₂. The worldwide distribution of the cycle slips are shown in Fig. 8a. As expected, their occurrence is linked with ionosphere active regions such as polar and equatorial areas. Because of a higher orbit (530 km) where the ionosphere is weaker than at the lower orbit (450 km), the number of cycle slips on Swarm B is less than for Swarm A and C. Due to the similar atmosphere around Swarm A and C, they suffer from the cycle slips almost at the same time. Significantly more cycle slips are found around the geomagnetic South Pole than around the North Pole. The reasons are still under investigation.

Table 2 Statistics of the cycle slips for Swarm A, B and C during 1 year from DoY 245, 2015 to DoY 244, 2016

(Schmid et al., 2007)(van den IJssel et al., 2016)(Wu et al., 1993)(IS-GPS-200H, 2013)(Hofmann-Wellenhof et al., 2008)

Table 3 Summary of the measurement and error models used for Swarm kinematic orbit determination.

Figure 8 Global distribution of cycle slips on Swarm A (red), B (green) and C (blue) (a) from DoY 245, 2015 to DoY 244, 2016 and (b) from September 2015. The black lines are the geomagnetic Equator and geomagnetic parallel of ±20^∘. The results are given in Table 2.

In order to see the impact of update of bandwidth of the tracking loop, the cycle slips in September 2015 are additionally shown in Fig. 8b. It can be seen that the number of cycle slips are significantly reduced after the update in May 2015 on Swarm C. Because the ionosphere has been very weak since May 2016, the tracking is very stable. There were just a few cycle slips, so the impact of the update in June and August 2016 cannot be observed.

2.4 Outlier detection

As a purely geometric approach, the kinematic orbit is sensitive to the GPS measurement errors. As a consequence, it is important to screen the outliers. The epoch-differenced ionosphere-free linear combination ΔL₃ used in cycle-slip detection can also be used to detect the outliers. Considering that the L₃ noise under ionospheric scintillation can reach the 10 cm level, if ΔL₃ is larger than 20 cm, an outlier is detected and should be eliminated. Some systematic errors can also be detected with it. For cases where the orbit of a given GPS satellite is not accurate enough, the time series of ΔL₃ of the corresponding GPS satellite will not be zero-mean but will show a strong trend. Thus, to avoid a systematic error in the LEO orbit and carrier phase residuals, the corresponding satellites have to be removed accordingly. This can be seen in Fig. 9 for PRN 26 (purple) and PRN 32 (green). Such satellites can be found in the IGS weekly summary (e.g. ftp://ftp.igs.org/pub/product/1864/igs18647.sum.Z, last access: 17 September 2018). After removing these satellites, the quality of the kinematic orbits is significantly improved. The same effects can also be found on Swarm C at the same time.

3 Kinematic orbit determination

In this section, the approach for kinematic orbit determination at IfE will be introduced first. Then 1-year kinematic Swarm orbits are validated with the reduced-dynamic orbits and the SLR residuals. The fully populated variance–covariance matrices of the kinematic orbits obtained from the adjustment are shown and discussed.

3.1 Observation modelling

A MATLAB-based PPP software was developed at IfE Hannover using the least-squares adjustment method. The ionosphere-free linear combinations of P code and phase are computed and introduced as observables in the least-squares adjustment model to eliminate the first-order ionospheric effect. The GPS orbits and satellite clock errors are computed based on the CODE final GPS orbits and 5 s clock products provided by the Center for Orbit Determination in Europe (CODE) (Bock et al., 2009; Dach et al., 2016). The observation and error modelling used for the kinematic orbit determination is listed in Table 3.

Figure 9 Outlier detection and anomalous GPS satellite for Swarm A on DoY 271, 2015, here PRN26 (purple) and PRN32 (green); (a) ΔL₃ between epochs separated by 100 s, (b) phase residuals and (c) deviations of the kinematic orbit with respect to ESA reduced-dynamic orbits with (red) or without (blue) the anomalous GPS satellites.

Figure 10 Position residuals for Swarm A (red), B (green) and C (blue) with respect to ESA reduced-dynamic orbits for DoY 333, 2015. The time series are offset by 10 cm for Swarm A and C. The legend indicates the rms (r) and mean value (m) in cm.

The 1 Hz data is processed over 30 h, with 3 h overlap in the previous and following day. In this way, the removal of some GPS observations can be avoided due to the short arc at the day boundary. This increases the stability of the ambiguities at the day boundary. Because the observation time saved in the RINEX file is synchronized to the receiver internal clock, which has a large clock drift, a time offset is added to synchronize the internal clock with GPS time everyday at the day boundary. This causes a clock jump in the code observations and a clock jump as well as a phase jump in the carrier phase observations. In order to get a continuous arc at the day boundary, these phase jumps need to be fixed. The clock jump is first determined from the average of each epoch-differenced code observations (removing the a priori geometric distances) at the day boundary. Then, the phase jump is the average of each epoch-differenced phase observations (removing the a priori geometric distances) minus the clock jump.

The phase noise is weakly elevation-dependent (see Fig. 5) and an identical weight matrix is thus applied to phase observations. The code noise has a stronger elevation-dependent behaviour after fixing the converter issue. Therefore, a sine-squared elevation-dependent weight matrix is applied to code observations after fixing instead of a sine elevation-dependent weight matrix before fixing. The variances are selected with 6 m (before fixing the converter issue)/0.6 m (after fixing the converter issue) for the code and 6 mm for the carrier phase according to the investigation described in Sect. 2, taking the error propagation for the ionosphere-free linear combination into account and considering the phase noise under ionospheric scintillation. Accordingly, the weight matrix per epoch has a diagonal structure.

Figure 11 Differences of kinematic orbits with respect to ESA reduced-dynamic orbits of Swarm C of different institutes from DoY 1 to 7, 2016. The time series are intentionally shifted by 15, 5, −5 and −15 cm for AIUB (red), IfG (green), ESA (blue) and IfE (black), respectively. The legend gives the rms in cm.

Table 4 rms and mean of position residuals of IfE and ESA kinematic orbits with respect to ESA reduced-dynamic orbits for Swarm A, B and C during 1 year (DoY 245, 2015 to DoY 244, 2016). The values for ESA kinematic orbits are given in brackets.

The three unknown position coordinates x, y, z and receiver-clock error δt_r are estimated epoch-wise, while the ambiguities b are constant over one consecutive GPS satellite-to-satellite tracking pass. A batch least-squares adjustment with pre-elimination and back-substitution is used for the estimation of the unknown parameters.

3.2 Kinematic orbit results

One-year orbits from 2 September 2015 to 31 August 2016 are computed. Here only exemplary but typically results are shown. The position residuals for Swarm A, B and C on DoY 333, 2015 with respect to the reduced-dynamic orbits provided by ESA in the Swarm Level 2 product are shown in Fig. 10. For a better illustration, the time series for Swarm A and C are shifted by ±10 cm.

The RMSEs in along-track, cross-track and radial directions are around 1.5, 1 and 2 cm, respectively. The noise caused by the ionospheric scintillations shown in Fig. 3 degrades the position when passing at equatorial and polar regions. It is mainly absorbed in the radial component. In some ionosphere-active areas, the deviations are beyond 10 cm or even 20 cm. Comparing Swarm A and C, a similar signature can be seen due to the similar GPS constellation and atmospheric conditions. Besides the large noise, all three orbits are affected by systematic errors in the along and radial components, predominantly at timescales in the orbital period. These periodic variations may either be caused by modelling deficiencies in the PPP solution or the reduced-dynamic orbits that are used as reference trajectory.

The RMSEs and average of position residuals for Swarm A, B and C of IfE and ESA kinematic orbits from DoY 245, 2015 to DoY 244, 2016 are listed in Table 4. The days with orbit maneuvers and instrument problems and epochs with GDOP values larger than 5 were excluded. The differences between the reduced-dynamic and kinematic orbits are at the 1.5, 1.5 and 2.5 cm level in the along-track, cross-track and radial directions, respectively. Since Swarm B flies in a higher orbit, the influence of ionosphere is weaker than for Swarm A and C. A better kinematic orbit quality can be observed for Swarm B, especially in the radial direction. It is also worth mentioning that the average offset between the kinematic and reduced-dynamic orbit in the cross-track direction is larger than in the along-track and radial directions.

Figure 12 Spectral analysis of Swarm C orbits from and (a) AIUB, (b) ESA,(c) IfG and (d) IfE in March 2016. Slope 0 represents white noise and slope −1 flicker noise.

Figure 13 Daily RMSEs of the ionosphere-free carrier phase residuals and the daily average STEC derived from Swarm C over all GPS satellites for 1 year of data from DoY 245, 2015 to DoY 244, 2016.

Figure 14 SLR residuals over 1 year from DoY 245, 2015 to 244, 2016, for (a) Swarm A, (b) B and (c) C. The residuals for IfE are given in blue and ESA kinematic orbits in red.

Figure 15 Correlation matrix for 30 h time span in X-component obtained with a (a) 30 s, (b) 10 s and (c) 1 s data rate. The identical epochs (30 s interval) are extracted from the respective correlation matrix. Panel (d) shows the mean correlation over 2 h.

The computed kinematic orbits of Swarm C are also compared with the kinematic orbits from ESA (van den IJssel et al., 2015), IfG of TU Graz (Zehentner and Mayer-Gürr, 2015) and AIUB (Jäggi et al., 2016). Again, ESA reduced-dynamic orbits are the reference trajectory. The position residuals of all institutes from DoY 01 to DoY 07, 2016 are shown in Fig. 11. The rms values are given in the legend. Our results are in good agreement with the kinematic orbits derived from the other institutes.

In order to analyse the remaining deviations, a spectral analysis is performed for Swarm C for March 2016 in Fig. 12. The once-per-revolution periodic signatures in the along-track and radial direction are well identified in all orbits except that of IfG, with 93.5 min for IfE and AIUB, 80.3 min for ESA. Besides the once per revolution, a twice-per- revolution periodic signature can be observed in the three orbits probably due to the impact of ionospheric scintillations at polar areas, with 46.2 min for IfE and ESA, 48 min for AIUB. Interestingly, the signatures of the cross track are different; IfE orbits are showing no significant periodicity. A half-day periodicity can be observed in AIUB and ESA orbits. The remaining high-frequency noise in the radial component is clearly visible in IfE and ESA orbits, which is linked to the observation noise during the passage of ionosphere-active areas. For the along- and cross-track component, IfE is at the AIUB level for frequencies larger than 10⁻³ Hz. High-frequency white noise contribution dominates in IfE and AIUB orbits until 250 s and ESA and Graz orbits until 100 s. Then the flicker noise persists in all the components. Between once per revolution and twice per revolution, some along-track and radical components in IfE persist periodically. The cross-track component of IfE is showing a f⁻¹ behaviour.

An additional possibility to assess the orbit quality is to analyse the residuals of GPS observations (van den IJssel et al., 2015). The daily RMSEs of the ionosphere-free carrier phase residuals for Swarm A, B and C for 1 year of data from DoY 245, 2015 to 244, 2016 are shown in Fig. 13, for example. The average of the daily RMSEs are at the 4 mm level. As expected, the daily RMSEs of Swarm B are slightly smaller than for Swarm A and C. Swarm A and C show similar results. The rms values show a decreasing trend with time, which could be related to the decreasing activity of ionosphere. Comparing Swarm A and C in September 2015, a significant improvement due to the update of the bandwidth on Swarm C can be observed.

The comparison with reduced-dynamic orbits and the residuals of GPS observations is an internal validation of the orbits. Independent satellite laser ranging (SLR) measurements provide another external validation of the orbits. The GPS-derived kinematic orbits are first interpolated with a polynomial of degree 5 to the observation time of SLR measurements. Then, the computed ranges between the ground stations of tracking networks and the GPS-derived kinematic orbits are compared with the observed laser ranges. The coordinates of the ground stations are given in ITRF2014 (Altamimi et al., 2016) with consideration to the solid Earth tide and ocean loading (FES2004 Lyard et al., 2006). The positions of the laser reflector in the satellite reference frame are given in Astrium GmbH (2013). The tropospheric delay and relativistic effects are also corrected. Reflector-dependent range corrections are not applied. The SLR residuals for Swarm A, B and C from DoY 245, 2015 to DoY 244, 2016 are shown in Fig. 14. The SLR residuals larger than 30 cm are excluded as outliers. The RMSEs of SLR residuals are at a 3–4 cm level for our orbits and ESA kinematic orbits.

Figure 16 Correlation matrix for 30 h time span in Z-component obtained with a  (a) 30 s, (b) 10 s and (c) 1 s data rate. The identical epochs (30 s interval) are extracted from the respective correlation matrix. Panel (d) shows the mean correlation over 2 h.

3.3 Covariance Information

The least-squares estimation provides not only the coordinates but also the variance–covariance information, which is important for gravity field recovery (Jäggi et al., 2011). Usually, only the mathematical correlations between the coordinates at the same epoch are considered. But the carrier phase ambiguities also correlate coordinates over multiple epochs. The fully populated formal variance–covariance matrices for 30 h PPP are computed with 30, 10 and 1 s GPS observation data rates, respectively. For undifferenced 30 and 10 s GPS observations, the covariance matrix can be directly computed from the inverse of the normal equation. However, due to computational limitation, the covariance of 1 s is computed epoch-wise by the sequential inversion of block submatrices. For a better illustration of the covariance matrix, it is converted into a correlation matrix.

Figure 15 shows the corresponding correlation matrix of the Earth-fixed X-component for all 30 h in the 30 s intervals. The X-coordinates of 30 s rate are strongly correlated at the polar areas (red spots) and the correlation decreases with time (Fig. 15a). The temporal mathematical correlation between two successive 30 s epochs is around 0.6. The correlation time varies between 0.5 and 2.5 h and shows a 24 h variation, which is caused by the period of the GPS constellation. Due to the increased number of observations, the correlation for 10 s observations is much lower than for 30 s. The correlation between two successive epochs is only around 0.4. But it shows a pattern similar to 30 s observations (Fig. 15b). As expected, with increased observations frequency, the correlation for the 1 s data rate is further decreased with a correlation between two successive epochs lower than 0.1 (Fig. 15c). Thus, the increase of data implies an increasing number of coordinates that are less correlated with the ambiguities, since no filter model is used in our PPP approach. The average correlation of X-components with times for the different rates is shown in Fig. 15d. The X-components are still correlated for 30 s data rate after 2 h. For the 10 s data rate, the average correlation is around 0.4 after 10 s and decreases smoothly to 0.1 after 2 h. There is almost no mathematical correlation for 1 s data rate, and the correlation drops immediately to 0.05 after 1 s. The correlation structure of the Earth-fixed Y-component is similar to the X-component.

The correlation structures of the Z-component for 30, 10 and 1 s data rates are shown in Fig. 16. The correlation between two successive epochs is similar to the X-component, around 0.4, 0.2 and below 0.1 for 30, 10 and 1 s data rates, respectively. However, the correlation times are much shorter than for the X-component, shown in Fig. 16. After an hour, the correlations can be neglected for 30 and 10 s data rates.

4 Conclusions

In this contribution, we present the approach for the determination of the IfE Swarm kinematic orbit with PPP. Since the GPS data quality is a prerequisite for the high quality of kinematic orbit, a thorough analysis was first performed. Although only 8 channels are available, the requirements for an accurate positioning are fulfilled. If the reacquisition time can be shortened, the tracking performance can be further improved. Carrier phase observations at both frequencies are strongly disturbed by ionospheric scintillations, where noise larger than 0.2 m degrades the detection of cycle slips. The approach is proposed, which allows for a successful cycle-slip repairing of around 95 %. It is worth mentioning that almost all the cycle slips occur on L₂ due to ionospheric scintillations.

Precise kinematic orbits at IfE are generated with the PPP method. The differences between our kinematic orbits and the ESA reduced-dynamic orbits are at 1.5, 1.5 and 2.5 cm level in along-track, cross-track and radial directions, respectively. The daily rms values of ionosphere-free carrier phase residuals are at 4 mm level, showing strong internal consistency. The spectral analysis of the orbits highlights the once-per-revolution period and higher harmonics as well as the dominant noise processes. A comparison with orbits from other institutes shows very good agreement with the rms level. A comparison with SLR underlines the accuracy of the kinematic orbits of 3–4 cm. Due to the ambiguities in the carrier phase, the coordinates are correlated. The mathematical temporal correlation between two successive epochs is larger than 0.6 for the 30 s data rate. The correlation decreases with increased sampling frequency. For a 1 s data rate, the correlation between two successive epochs is below 0.1.