GPS Based Reduced-Dynamic Orbit Determination for Low Earth Orbiters with Ambiguity Fixing

With the ever-increasing number of satellites in Low Earth Orbit (LEO) for scientific missions, the precise determination of the position and velocity of the satellite is a necessity. GPS (Global Positioning System) based reduced-dynamic orbit determination (RPOD) method is commonly used in the post processing with high precision.This paper presents a sequential RPOD strategy for LEO satellite in the framework of Extended Kalman Filter (EKF). Precise Point Positioning (PPP) technique is used to process the GPS observations, with carrier phase ambiguity resolution using Integer Phase Clocks (IPCs) products. A set of GRACE (Gravity Recovery And Climate Experiment) mission data is used to test and validate the RPOD performance. Results indicate that orbit determination accuracy could be improved by 15% in terms of 3D RMS error in comparison with traditional RPOD method with float ambiguity solutions.


Introduction
Precise orbit information of satellites is the basis for many space missions, such as the earth gravity recovery, atmosphere sounding, and ocean surveillance.From the beginning of 1980s, Global Positioning System (GPS), which can provide continuous, all-day, and high-precision observations, has been widely used in Positioning, Navigation, and Timing (PNT) applications.The concept of spaceborne GPS receivers tracking GPS signals for orbit determination has been proposed since then.During the past 30 years, a large number of satellites have been launched into Low Earth Orbit (LEO) with GPS receivers.GPS based Precise Orbit Determination (POD) has been considered as the fundamental and primary technique for spacecraft missions.
Several GPS based orbit determination strategies have been proposed in the literatures and applied in the real space missions as well.Generally, they are categorised into three groups: kinematic method (KPOD) [1,2], dynamic method (DPOD) [3][4][5] and reduced-dynamic method (RPOD) [6][7][8][9][10].Among them, the reduced-dynamic method is commonly used for centimetre-level of orbit determination precision, which adopts empirical accelerations, impulse velocities, or other random parameters to compensate for the uncertainty in the orbit model so as to achieve very accurate orbit prediction information.On the other hand, GPS system provides two types of ranging observations, namely, pseudorange and carrier phase observations.It is well known that carrier phase observations generally have several millimetres-to centimeter-level of precision (e.g., for geodetic receivers), which are much more precise than pseudorange observations (typically about two orders of magnitude).Hence, the integrity of accurate orbit prediction information together with high-precision GPS measurements in the sequential Kalman filter or batch Least Squares framework will result in good performance of orbit position and velocity estimation.However, carrier phase measurements are contaminated with unknown ambiguous integers [11].If such ambiguities are fixed, the carrier phase measurements become unambiguous pseudoranges but accurate at the level of few millimetres.Usually carrier phase ambiguities are taken as unknown parameters and estimated as float values with other states, for example, orbit positions, velocities, orbit parameters, empirical accelerations, and receiver clock offsets 2 International Journal of Aerospace Engineering in the traditional float PPP (Precise Point Positioning) based RPOD without introducing differencing technique [12,13].Obviously, the loss of integer property of these ambiguities will deteriorate the orbit determination performance.During last decade, many researchers have been focusing on PPP Ambiguity Resolution (PPP-AR).Several strategies have been proposed, which have been tested and validated in terrestrial, marine, and space scenarios with better positioning performance [14][15][16][17][18].This paper will review these approaches and use Integer Phase Clocks (IPCs) based PPP-AR method proposed by Laurichesse et al. [18] for the LEO satellites RPOD demonstration.
The rest of the paper is organised as follows: Section 2 introduces GPS observation functional model and satellite orbit model for propagation.Ionosphere-free combination is used in the PPP processing; Section 3 presents the Extended Kalman Filer (EKF) framework, which combines orbit prediction information with GPS measurement update; PPP-AR strategies are revisited in Section 4, emphasising on IPCs based method.Static data collected from IGS (International GNSS (Global Navigation Satellite System) Services) station are used for testing this method; IPCs based PPP-AR method is applied in RPOD applications in Section 5 and GRACE (Gravity Recovery And Climate Experiment) satellite flight data is used for algorithm testing and validation.Summary of this paper and some conclusions are given in the last section.

System Models in GPS Based
Orbit Determination where the superscripts  and  indicate the GPS satellite and the signal frequency, respectively, and subscript  indicates the receiver.The terms in the equations are as follows: PR, CP: code phase and carrier phase observation (m).
: geometric distance between receiver and GPS satellite antenna centres (m).

𝑁: carrier phase integer ambiguity (cycle).
: impact of relative rotation between transmitter and receiver antennas (cycle).

Ionosphere-Free Combination of GPS Observations.
According to (1), the GPS observations are contaminated by several errors coming from different sources.In order to obtain an accurate positioning, it is necessary to correct, or at least eliminate, most of them.The first order ionospheric error which accounts for more than 99.9% of the total ionospheric delay can be eliminated via (2) using Ionosphere-Free (IF) combination when the dual-frequency observations are available.The second and higher terms are commonly neglected [19]: Combining ( 1) and ( 2 In the IF observation equation, the carrier phase ambiguity,   ,IF , does not keep integer nature any more.The receiver clock time offset cl  in (3) needs to be determined as well, which is taken as a single white noise estimation parameter and estimated based on the kinematic navigation algorithm [20].
Another two types of observable combinations are commonly used for GPS data quality screening as shown below [21].One is Melbourne-Wübbena (MW) combination and the other is Geometry-Free (GF) combination.

Melbourne-Wübbena Combination of GPS Observations.
The MW combination can eliminate the effects due to the ionosphere, the geometry and the satellite, and receiver clocks, which is given by [22,23]  (5)

GPS Observation Errors Correction.
In addition to ionospheric path delay, other errors in GPS observation are necessary to be corrected, which are including relativistic biases in clocks, phase wind-up effects, antenna phase centre offsets, and multipath errors.Note that in the IF observation (2), the tropospheric path delay is neglected in comparison with the observation equation used for terrestrial applications.This is due to the fact that the LEO satellite is flying above the tropospheric level.Hence, the GPS signal will not be affected by the troposphere.More details on how to correct these errors especially for LEO satellite orbit determination can be found in [24].

Dynamic Model for Orbit Prediction.
For the spacecraft in near-earth orbit, a Newton-Kepler system is traditionally used to describe the orbit for the two-body case.Real orbit modeling, however, should take into account additional gravitational and nongravitational perturbations.In general the accelerations acting on the satellite consist of terms for the geopotential gradient, the 3rd-body gravitational attraction of the sun and moon, the solar radiation pressure, and atmosphere drag on the spacecraft, if no active orbital manoeuvre is performed.The exact formulations for each term can be obtained from classical books, for example [25].
Generally there are two main coordinate systems involved in the GPS-based LEO satellite orbit determination: the earth centred inertial (ECI) system and the earth centred earth fixed (ECEF) system.The former one provides a formulation of the satellite propagation and coordinates of other spatial bodies, for example, the sun and the moon.The solar radiation pressure acting on the spacecraft is also formulated in this frame.On the other hand, the ECEF frame provides modeling of the earth gravity field and atmospheric drag affecting the spacecraft.Aside from that, GPS orbits and its signals are given in this frame.With aim of implementing the sequential filtering algorithm for spacecraft orbit determination, the coordinate systems have to be consistent.In this sense, all the acceleration components are calculated or transformed into the ECEF system in this paper.Additionally, two types of "fictitious" forces are necessary to be added into this rotating reference frame, that is, the centrifugal and Coriolis accelerations, which are written as [20] where Ω is the angular velocity vector and r and k denote the position and velocity of the satellite in the ECEF frame.
The motion of the satellite along with the time  is modeled by the conventional orbit dynamics in the Cartesian coordinates: where the p is the orbital parameter vector, including the solar radiation pressure and air drag coefficients (  and   ).x comprises the position r and k vectors.a is the acceleration vector acting on the LEO satellite.The orbit can be propagated using the integration as follows: where  0 is the initial epoch.Differentiating x() with respect to x( 0 ) results in Note that p()/x( 0 ) = 0; hence, the state transition matrix can be obtained by integrating (10) with an initial value equaling the identify matrix: where Similarly, the sensitivity matrix can be calculated by integrating where Ψ() = [x()/p()] denotes the sensitivity matrix and C = [x()/p()] denotes the Jacobian matrix with respect to the considered parameters.Comparatively, the sensitivity matrix is initialised with zero matrix Ψ( 0 ) = 0 6×(6+2) .
International Journal of Aerospace Engineering

Empirical Accelerations.
Assume that uncertainty of the orbital model could be absorbed into empirical accelerations, which are always modeled as a first-order, stationary, Gauss-Markov process [26,27]: where  is the variable,  is the correlation time, and () is the white Gaussian noise.Hence, the solution of the variable  can be expressed as follows: where Δ is the processing interval between two epochs.Define  =  −Δ/ as the mapping function, and the process noise of the variable can be formulated as an uncorrelated random sequence with zero mean and variance:

Extended Kalman Filter Design
This research uses EKF to determine the LEO satellite positions and velocities epoch by epoch [25].The states vector in the filter is written as follows:  3.

Time Update.
In the time update step of EKF, both the states and their covariance matrix should be predicted: where P is the state covariance matrix and Q is the process noise matrix; subscript  is the time epoch.
To use RPOD approach, state transition matrix Φ and process noise covariance matrix Q can be expressed as where Φ r,k , Ψ   ,  , and Ψ a can be integrated using the formulations in [25,28];  =  −Δ/ a is mapping function of empirical accelerations defined in Section 2.3.1.The process noise matrix is diagonal, with process noises being added on empirical accelerations and the receiver clock offset.As described above, empirical accelerations accord with firstorder Gauss-Markov process, while the receiver clock offset is treated as a stochastic variable with an underlying random walk process model [28].

Data Quality Control.
Traditional GPS data quality control methods by means of MW and GF observation combinations (see (4)∼( 5)) are used for gross errors and cycle slips detection [21].In addition, a receiver clock estimation based GPS data screening approach is also used in this research [29].Within this approach, the a priori orbit positions predicted by orbit models are used to estimate a more precise receiver clock offset compared with conventional single point positioning solution at every epoch.Then the estimated receiver clock offset and its uncertainty can be taken as threshold to screen gross errors and cycle slips.

Measurement Update.
Within the measurement update step of the EKF, the predicted states and covariance matrix are required to be updated using new GPS observations: where K is the Kalman gain, H is the design matrix of measurement equation, and R is the observation noise matrix.y is the observation, while h is the calculated one.The formulation of H matrix can be written as ) . (20)

PPP Ambiguity Resolution
4.1.Traditional Float PPP.Carrier phase integer ambiguity resolution was applied in Real-Time Kinematic (RTK) positioning at early stage, since differencing technology could eliminate and mitigate most of the observation errors and biases so that the integer property of the ambiguity could be preserved.In the traditional PPP algorithm, however, it is hard to achieve double-differenced positioning accuracy.
Here the GPS IF observations for LEO satellite orbit determination in (3) are revisited: The clock offsets are written together with Hardware Delay Biases (HDBs).Antenna phase centre offsets and phase windup effects are not explicitly shown in the equation, but they are necessary to be corrected.GPS clock offsets  ⋅ cl  +   PR IF can be estimated using IGS clock products.Then GPS satellite and user receiver HDBs  ,CP IF ,   CP IF will be absorbed into ambiguity   ,IF , which will not hold its integer property any more.Therefore, it is impossible to fix ambiguity parameters as integers in the traditional PPP algorithm.Instead, they are estimated as float values.Hence, the positioning performance will be deteriorated in the traditional float PPP approach.

PPP Ambiguity Resolution Strategies.
Since 2006, researchers started fixing integer ambiguity by correcting initial phase biases (namely, HDBs) in zero-difference PPP simulations.Wang and Gao found that phase biases of single-difference between receivers are unstable, while that of Single-Difference Between GPS Satellites (SDBS) are very stable in several days.Hence, they concluded that PPP ambiguities could be fixed as integers with GPS satellite phase biases well corrected [30,31].Ge et al. [14] also presented the stability of SDBS wide-lane Uncalibrated Phase Delays (UPDs).Based on the UPD corrections, PANDA software was used to process one-day data collected from ground stations; 30% positioning improvements could be achieved in the East component with 80% ambiguity fixing ratio.Geng et al. [15] improved Ge's method by introducing the Fractional Cycle Biases (FCBs) and reducing the broadcasting load of the corrections.This method was applied in marine surveying with long baseline; centimetre-level of accuracy could be achieved with ambiguity fixing using PPP algorithm [32].Another PPP integer ambiguity resolution method was proposed by Bertiger et al. [16].Unlike Ge's UPD and Geng's FCB corrections, the undifferenced float ambiguities estimated in the network solution are directly broadcasting to users.Therefore, users could form doubledifferenced ambiguities to clean all the biases.Additionally, satellite biases are corrected by assimilating them into the clock parameters.Collins et al. [17] introduced Decoupled Satellite Clocks (DSCs) based PPP-AR method, in which pseudorange biases and carrier phase biases are treated differently.Laurichesse et al. [18] proposed Integer Phase Clocks (IPCs) based PPP-AR method.Only GPS satellite clock and biases in IF carrier phase are taken into account in comparison with Collins' method.Both Laurichesse and Bertiger applied their methods into LEO satellite precise orbit determination applications.According to the formulation to separate the hardware biases, these strategies are named as follows: (1) Single Difference Between GPS Satellites (SDBS) [14,15].
In summary, with receiving extra information from network solutions, standalone receivers are able to implement integer ambiguity resolution in PPP algorithm.Since 2009, IGS CNES-CLS (Centre National d'Etudes Spatiales, Collecte Localisation Satellites) started to release GPS satellite IPCs products together with Wide-lane Satellite Biases (WSBs) [33,34].This research will focus on IPCs based PPP-AR algorithm in LEO satellite RPOD demonstration.

IPCs Based PPP-AR Method.
The basic concept of this method is to separate pseudorange clocks and carrier phase clocks; the float wide-lane ambiguities  WL are to be fixed as integers after WSBs corrections and GPS satellite clock offsets are corrected by IPCs product so as to recover the integer narrow-lane ambiguity  NL ; then IF ambiguity could be solved and states of position can be corrected and updated with fixed carrier phase ambiguity.
The IF observations given in ( 21) are used for PPP processing.The IGS precise clock product is used for satellite clock offsets corrections ( ⋅ cl  +  ,PR IF ) in pseudorange measurement and CNES-CLS IPCs product is used for satellite clock offsets corrections ( ⋅ cl  +   CP IF ) in carrier phase measurement.Since the wavelength of IF observation is only 0.6 cm, it is usually separated as wide-lane (86.7 cm) and narrow-lane (10.7 cm) ambiguity as formulated in Note that in order to avoid estimation of receiver clock offsets, which contains corresponding uncalibrated hardware biases at receiver end, SDBS is formed.
Wide-Lane Ambiguity Resolution and Validation.The widelane ambiguity resolution  WL depends on SDBS MW combination, which can be formed using (4): Before separating  , WL, from  , MW , a time-smoother is necessary to mitigate and eliminate multipath effects and thermal noises in MW observations as formulated in [14,18] MW where subscript  indicates the th epoch, and the corresponding variance can be express as ) . ( The wide-lane satellite biases

𝑖,𝑗
MW can be extracted and calculated from CNES-CLS WSBs products.Since the widelane wavelength is around 86.2 cm, a simple rounding method is used to estimate the single-differenced wide-lane ambiguity by (24) after correction.Validation procedure is required to be done before the above ambiguity candidate solutions are fixed as integers.Except for the traditional ratio method (see (27)), the following formulation is working simultaneously to test and validate the right integer ambiguity resolution [14,35]: where N is the real ambiguity,  is its variance, Ň is the nearest integer of float N, and Ň is the second nearest integer of float N. Giving the confidence value , for example,  = 0.01%, the float ambiguity could be fixed as the nearest integer when  0 is larger than 1 −  while condition in (27) should be satisfied as well ( thresh is commonly set with a value of 3).
Narrow-Lane Ambiguity Resolution.IF ambiguities are estimated by means of traditional float PPP algorithm.Since the satellite integer phase clock produced by CENS-CLS contains the IF satellite biases ( ⋅ cl  +   CP IF ), hence, the estimated IF ambiguity is free of these biases.With known integer wide-lane ambiguity, ( 23) is derived to solve narrow-lane ambiguities: where

𝑖,𝑗
NL is the float SDBS narrow-lane ambiguity to be solved.The corresponding variance is given as With solved float narrow-lane ambiguity and its variance, the rounding method is also used to find the integer solution and formulations in ( 27) and ( 28) are used for validation.
After fixing the narrow-lane ambiguity, the IF ambiguity is recalculated together with the integer wide-lane ambiguity.This precise information can be considered as new pseudoobservation and used for positions update in a new run of Kalman filter.The flow diagram of the IPCs based PPP-AR procedure is shown in Figure 1.

IGS Data Tests.
Firstly, the IPCs based PPP-AR method is tested with IGS station observations.Data collected from four IGS stations (AUCK, SHAO, TLSE, and WTZZ) on February 20, 2012, are used for both static and kinematic PPP testing.Table 1 presents positioning errors in ENU (East, North, and Up) coordinate system and ambiguity fixing ratios.Several centimetres of positioning accuracy could be achieved using the static PPP algorithm except for SHAO station with worse positioning performance in the Up direction.However, PPP ambiguity resolution does not improve the positioning accuracy significantly.One decimetre of positioning accuracy could be achieved using kinematic PPP algorithm, and PPP-AR can improve positioning accuracy obviously.In comparison with traditional float PPP solutions, positioning errors in terms of 3D RMS of four stations have been reduced.In particular, positioning error in East direction of AUCK is reduced by 40% and 3D RMS error is reduced by over 9%.
Figure 2 depicts positioning errors of WTZR in ENU coordinates using kinematic PPP algorithm.In the figure, blue plots indicate the positioning errors using traditional float PPP approach, while green plots indicate that using IPCs based PPP-AR approach.It is noted that PPP-AR approach improves positioning accuracy obviously and shortens the convergence time as well.

GRACE Data Analysis
Spaceborne GPS data from the GRACE mission is used for testing RPOD performance with ambiguity resolution, which   has twin LEO satellites flying in a group for space missions of earth gravity field modeling and climate experiments [36].This data can be accessed from JPL's (Jet Propulsion Laboratory, USA) Physical Oceanography Distributed Active Archive Centre (PODAAC) or from GFZ's (Geo-ForschungsZentrum, Germany) Information System and Data Centre (ISDC) [37].GPS flight data products are in the format of GPS1B while GPS Navigation (GNV) products are also released, which contain reduced-dynamic precise orbit determination solutions with a time-varying accuracy of a few centimetres in position [38].In this research, the GNV solutions are taken as the real trajectories of GRACE satellites.
To demonstrate the RPOD with PPP ambiguity resolution described above, both L1 and L2 code phase and carrier phase data from GRACE-B onboard GPS receivers are used, which covers 10 hours on February 20, 2012, at a sampling interval of 10 s.GPS satellite clocks with consistent satellite orbits are provided weekly by the CNES-CLS IGS analysis centre [33], in which the clock product consists of WSBs and phase clock with IF satellite hardware biases.The dynamic models/parameters used for orbit propagation and corresponding state transition matrix and sensitivity matrix are summarised in Table 2. Various error sources in the observation model equations (see (3)) should be carefully corrected, for example, the carrier phase windup effects, the antenna phase centre offsets, and some other systematic noises [24].EKF parameter settings are given in Table 3.
RPOD solutions in RIC coordinates with and without phase ambiguity fixing using 10 hours of GRACE-B satellite International Journal of Aerospace Engineering  Low precision model [25] Coordinates transformations IERS1996/IAU1980 transformations [41] data are shown in Figure 3.In the figure, red plots indicate the traditional float PPP based RPOD results, while blue plots indicate PPP-AR based RPOD results.Except for some loss of accuracy in the radial component, orbit determination performance in other two components using PPP-AR approach have been improved considerably in comparison with float PPP approach, with 0.062 m, 0.025 m, and 0.014 m in terms of RMS error.And 3D RMS error of 0.068 m using PPP-AR strategy is reduced by 15% in comparison with float PPP results of 0.080 m. Figure 4 (up) presents PPP ambiguity resolution results using the specified GRACE-B satellite data.3435 epochs have been fixed with integer ambiguity in the total 3600 epochs (first 100 epochs of the day are excluded for analysis due to the convergence), with a fixing ratio of 95.4%.In Figure 4   (bottom), red stars indicate that the epochs PPP-AR strategy fails in those epochs.
The empirical accelerations estimated in RPOD are plotted in RIC coordinate system in Figure 5.All the three components are several hundred nm/s 2 , higher than the values stated in [10,28,29].The radial acceleration component is affected by atmospheric drag, which leads to orbital decay.The in-track component mainly reflects a mismodeling of atmospheric drag.The cross-track component might be related to albedo effects, mismodeled solar radiation pressure, or unmodeled tidal perturbations.Hence, the orbital models used in this paper should be refined in the future work.
The IF observation residuals are plotted in Figure 6.The RMS of IF pseudorange observable residuals is 0.545 m, while that of IF carrier phase observable residuals is 8 mm.

Concluding Remarks
This paper presents PPP ambiguity resolution in reduceddynamic orbit determination using GPS observations.Several centimetres of accuracy has been achieved.Carrier phase ambiguity resolution could significantly improve orbit determination performance by 15% in terms of 3D RMS in the GRACE-B satellite scenario.In the RPOD algorithm, the empirical accelerations are in accordance with the orbital model precision.Hence, it is important to determine the value of empirical accelerations.On the other hand, the estimated empirical accelerations can help detect mismodeling part in the orbital models.Since GNSS real-time PPP-AR method has attracted much more attention these days, it is potential to apply this technique into real-time precise orbit determination in the future work.

Figure 2 :
Figure 2: WTZR positioning errors comparison using PPP-AR and float PPP strategies in the kinematic mode.
Single point positioning algorithm is used for determining the positions and velocities of LEO satellite and receiver clock offsets at the first epoch.Carrier phase ambiguities are initialised by carrier phase measurement minus pseudorange measurement.As for the orbit parameters, namely, solar radiation pressure and atmospheric drag coefficients, their initial values are given empirically.For instance,   = 1.3 and   = 2.3 for GRACE satellites.Other initial values for RPOD demonstration with GRACE data can be found in Table which consists of LEO satellite position vector r = [; ; ] and velocity vector k = [V  ; V  ; V  ] in ECEF coordinate frame, orbital parameters p = [  ;   ], empirical acceleration vector a = [  ;   ;   ] in RIC (Radial, In-track, and Cross-track) coordinate frame, receiver clock offset ⋅cl  , and IF carrier phase ambiguity A. Only ambiguities of visible GPS satellites at current epoch are to be estimated, so the total dimension of the states is 12 + , where  is the number of visible GPS satellites.3.1.Filter Initialisation.

Table 2 :
Orbit models used in this paper.