PPP-based Swarm kinematic orbit determination

ESA’s Swarm mission offers excellent opportunities to study the ionosphere and to bridge the gap in gravity field recovery between GRACE and GRACE-FO. In order to contribute to these studies, at 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 5 on the Melbourne-Wübbena and ionosphere-free linear combination. The results show that around 95% cycle slips can be repaired and the majority of the cycle slips occur on L2. After the analysis and careful preprocessing of the observations, one 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 ESA in Swarm Level 2 products and SLR measurements. The differences between IfE kinematic orbits and ESA reduced-dynamic orbits are at the 1.5 cm, 1.5 cm and 10 2.5 cm level in the along-track, cross-track and radial directions, respectively. Remaining systematics are characterised by spectral analyses, showing once per revolution period. The external validation with SLR measurements shows RMS errors at the 4 cm level. Finally, fully populated covariance matrices of the kinematic orbits obtained from 30 s, 10 s and 1 s data rate are discussed. It is shown that for data rates larger than 10 s, the correlation should be taken into account. Copyright statement. TEXT 15


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 (Švehla and Rothacher, 2003;Zehentner and Mayer-Gürr, 2015;Jäggi et al., 2016).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 L. Ren and S. Schön: PPP-based Swarm kinematic orbit determination 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 (van den IJssel et al., 2015;Jäggi et al., 2016).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 mea-surements.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.

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.

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 1 -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.

Observation analysis
As a dual-frequency receiver, the Swarm GPS receiver provides code observations C/A (C1C), P 1 (C1P) and P 2 (C2P) and carrier phase observations L 1 (L1C) and L 2 (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 1/ √ sin(Elev).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 1 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 1 code is disturbed by noise of about 0.8 m (1σ ) for satellites at elevations lower than 15   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 ( 2 L(t) = L(t + 1) − 2L(t) + L(t − 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.
The resulting 2 L 1 and 2 L 2 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 1 and L 2 , 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 secondorder differences of the ionosphere-free linear combination 2 L 3 .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 2 L 1 and 2 L 2 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 1 is slightly lower than on L 2 .The standard deviations (1σ ) of 2 L 1 and 2 L 2 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 1 and 1.7 mm for L 2 , according to the laws of error propagation.

Cycle-slip detection and repair
Strong ionospheric scintillations cause not only large noise but also cycle slips in the carrier phase observations.There-fore, 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 where L MW represents the Melbourne-Wübbena linear combination, f 1 and f 2 the carrier frequency, L 1 and L 2 the carrier phase observation in meters, b 1 and b 2 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: where is the operator for epoch differencing, and c 1 = 1 can be identified in the peak close to an integer (difference smaller than 0.1) in the black curve in this example.
Since with the Melbourne-Wübbena linear combination only the difference of cycle slips on L 1 and L 2 can be detected and simultaneous occurring of same magnitude cycle slips on L 1 and L 2 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: where L 3 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 1 and b 2 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 1 and L 2 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 3 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 Together with the equation from Melbourne-Wübbena linear combination, the cycle slip on L 1 and L 2 can be determined: In this example we can calculate the following: b 1 ≈ 0.023 ≈ 0 and b 2 ≈ −1.If the difference of b 1 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 2 .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 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.

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 3 used in cycle-slip detection can also be used to detect the outliers.Considering that the L 3 noise under ionospheric scintillation can reach the 10 cm level, if L 3 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 3 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.

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 variancecovariance matrices of the kinematic orbits obtained from the adjustment are shown and discussed.

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) (Dach et al., 2016;Bock et al., 2009).The observation and error modelling used for the kinematic orbit determination is listed in Table 3.
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 am- biguities 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 obser-vations (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 sinesquared elevation-dependent weight matrix is applied to code observations after fixing instead of a sine elevationdependent weight matrix before fixing.The variances are selected with 6 m (before fixing the converter issue)/0.6m (after fixing the converter issue) for the code and 6 mm for the carrier phase according to the investigation described www.ann-geophys.net/36/1227/2018/Ann.Geophys., 36, 1227-1241, 2018 1236 L. Ren and S. Schön: PPP-based Swarm kinematic orbit determination in Sect.2, taking the error propagation for the ionospherefree linear combination into account and considering the phase noise under ionospheric scintillation.Accordingly, the weight matrix per epoch has a diagonal structure.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.

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 alongtrack, 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.
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 reduceddynamic 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 twiceper-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 −3 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 −1 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

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 be- tween 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.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 cycleslip repairing of around 95 %.It is worth mentioning that almost all the cycle slips occur on L 2 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 onceper-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.

Figure 1 .
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 .Figure 3 .
Figure 2. (a)P 1 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 4
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 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 b 1 − b 2 = 1.

1233Figure 7 .
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 3 between two successive epochs 0.29 = 0.4844 b 1 − 0.3775 b 2 and (b) L 3 between two epochs separated by 100 epochs 0.38 = 0.4844 b 1 − 0.3775 b 2 .

1234L.
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 inTable 2.

Figure 9 .Figure 10 .
Figure 9. Outlier detection and anomalous GPS satellite for Swarm A on DoY 271, 2015, here PRN26 (purple) and PRN32 (green); (a) L 3 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 11 .
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.

Figure 14 .
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. 1239

Figure 15 .
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.

Figure 16 .
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.

Table 1 .
Update of Swarm carrier loop bandwidth

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

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

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.