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

Authors: Yang Yang; Xiaokui Yue; Jianping Yuan
Journal: International Journal of Aerospace Engineering
 (2015)

---

## Abstract

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.


---

## Body

## 1. 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–5] and reduced-dynamic method (RPOD) [6–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 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–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: Section2 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.


## 2. System Models in GPS Based Orbit Determination

### 2.1. Functional Model of GPS Observations

The GPS code phase and carrier phase observation equations are written as follows [11]:(1)PRrs,f=ρrs+cclr-cls+Irs,f+Trs+Br,PRf-BPRs,f+Mr,PRs+νr,PRs,f,CPrs,f=ρrs+cclr-cls-Irs,f+Trs+Br,CPf-BCPs,f+Nrs,fλf+ϕrsλf+Mr,CPs+νr,CPs,f,where the superscripts s and f indicate the GPS satellite and the signal frequency, respectively, and subscript r 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).c: speed of light: 299 792 458 (m/s).cl: clock error, including relativistic effects (s).I: ionospheric path delay (m).T: tropospheric path delay (m).B: hardware delay bias (m).λ: signal wavelength (m).N: carrier phase integer ambiguity (cycle).ϕ: impact of relative rotation between transmitter and receiver antennas (cycle).M: multipath effects (m).ν: unmodeled code phase and carrier phase errors (m).


#### 2.1.1. 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]:(2)PRr,IFs=f12f12-f22PRrs,1-f22f12-f22PRrs,2,CPr,IFs=f12f12-f22CPrs,1-f22f12-f22CPrs,2.Combining (1) and (2), the IF observations appear as(3)PRr,IFs=ρrs+cclr-cls+Br,PRIF-BPRs,IF+Mr,PRIFs+νr,PRIFs,CPr,IFs=ρrs+cclr-cls+Br,CPIF-BCPs,IF+Nr,IFsλIF+λIFϕr,IFs+Mr,CPIFs+νr,CPIFs.In the IF observation equation, the carrier phase ambiguity,Nr,IFs, does not keep integer nature any more. The receiver clock time offset clr 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.


#### 2.1.2. 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](4)MWrs=1f1-f2f1CPrs,1-f2CPrs,2-1f1+f2f1PRrs,1+f2PRrs,2=λWLN1,rs-N2,rs+Br,MW-BMWs+νr,PRNLs+νr,CPWLs=λWLNWL,rs+Br,MW-BMWs+νr,MWs.


#### 2.1.3. Geometry-Free Combination of GPS Observations

GF observable is formed via the difference between the dual-frequency observations, which is also independent of geometry, receiver clocks, and satellite clocks [11]:(5)PRr,GFs=PRrs,1-PRrs,2=Ir,GFs+Br,PRGF-BPRGFs+Mr,PRGFs+νr,PRGFs,CPr,GFs=CPrs,1-CPrs,2=-Ir,GFs+λGFfNr,GFs+Br,CPGF-BCPGFs+λGFϕr,GFs+Mr,CPGFs+νr,CPGFs.


### 2.2. 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].


### 2.3. 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](6)acentrifugal=-Ω×Ω×rECEF,aCoriolis=-2Ω×vECEF,where Ω is the angular velocity vector and r and v denote the position and velocity of the satellite in the ECEF frame.The motion of the satellite along with the timet is modeled by the conventional orbit dynamics in the Cartesian coordinates:(7)ddtxt=ft,xt,pt=vtat,rt,vt,where the p is the orbital parameter vector, including the solar radiation pressure and air drag coefficients (Cr and Cd). x comprises the position r and v vectors. a is the acceleration vector acting on the LEO satellite. The orbit can be propagated using the integration as follows:(8)xt=xt0+∫t0tfx,p,tdt,where t0 is the initial epoch. Differentiating x(t) with respect to x(t0) results in(9)∂xt∂xt0=I+∫t0t∂f∂xt∂xt∂xt0+∂f∂pt∂pt∂xt0dt.Note that∂p(t)/∂x(t0)=0; hence, the state transition matrix can be obtained by integrating (10) with an initial value equaling the identify matrix:(10)Φ˙t,t0=FΦt,t0,Φt0,t0=I,where (11)Φt,t0=∂xt∂xt0,F=∂f∂xt.Similarly, the sensitivity matrix can be calculated by integrating(12)Ψ˙t=FΨt0+C,Ψt0=0,where Ψ(t)=∂x(t)/∂p(t) denotes the sensitivity matrix and C=∂x(t)/∂p(t) denotes the Jacobian matrix with respect to the considered parameters. Comparatively, the sensitivity matrix is initialised with zero matrix Ψ(t0)=06×(6+2).


#### 2.3.1. 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]:(13)ω˙t+1τωt=ut,where ω is the variable, τ is the correlation time, and u(t) is the white Gaussian noise. Hence, the solution of the variable ω can be expressed as follows:(14)ωt+Δt=ωte-Δt/τ,where Δt is the processing interval between two epochs. Define d=e-Δt/τ as the mapping function, and the process noise of the variable can be formulated as an uncorrelated random sequence with zero mean and variance:(15)Eu2=σ21-d2.


## 2.1. Functional Model of GPS Observations

The GPS code phase and carrier phase observation equations are written as follows [11]:(1)PRrs,f=ρrs+cclr-cls+Irs,f+Trs+Br,PRf-BPRs,f+Mr,PRs+νr,PRs,f,CPrs,f=ρrs+cclr-cls-Irs,f+Trs+Br,CPf-BCPs,f+Nrs,fλf+ϕrsλf+Mr,CPs+νr,CPs,f,where the superscripts s and f indicate the GPS satellite and the signal frequency, respectively, and subscript r 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).c: speed of light: 299 792 458 (m/s).cl: clock error, including relativistic effects (s).I: ionospheric path delay (m).T: tropospheric path delay (m).B: hardware delay bias (m).λ: signal wavelength (m).N: carrier phase integer ambiguity (cycle).ϕ: impact of relative rotation between transmitter and receiver antennas (cycle).M: multipath effects (m).ν: unmodeled code phase and carrier phase errors (m).


### 2.1.1. 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]:(2)PRr,IFs=f12f12-f22PRrs,1-f22f12-f22PRrs,2,CPr,IFs=f12f12-f22CPrs,1-f22f12-f22CPrs,2.Combining (1) and (2), the IF observations appear as(3)PRr,IFs=ρrs+cclr-cls+Br,PRIF-BPRs,IF+Mr,PRIFs+νr,PRIFs,CPr,IFs=ρrs+cclr-cls+Br,CPIF-BCPs,IF+Nr,IFsλIF+λIFϕr,IFs+Mr,CPIFs+νr,CPIFs.In the IF observation equation, the carrier phase ambiguity,Nr,IFs, does not keep integer nature any more. The receiver clock time offset clr 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.


### 2.1.2. 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](4)MWrs=1f1-f2f1CPrs,1-f2CPrs,2-1f1+f2f1PRrs,1+f2PRrs,2=λWLN1,rs-N2,rs+Br,MW-BMWs+νr,PRNLs+νr,CPWLs=λWLNWL,rs+Br,MW-BMWs+νr,MWs.


### 2.1.3. Geometry-Free Combination of GPS Observations

GF observable is formed via the difference between the dual-frequency observations, which is also independent of geometry, receiver clocks, and satellite clocks [11]:(5)PRr,GFs=PRrs,1-PRrs,2=Ir,GFs+Br,PRGF-BPRGFs+Mr,PRGFs+νr,PRGFs,CPr,GFs=CPrs,1-CPrs,2=-Ir,GFs+λGFfNr,GFs+Br,CPGF-BCPGFs+λGFϕr,GFs+Mr,CPGFs+νr,CPGFs.


## 2.1.1. 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]:(2)PRr,IFs=f12f12-f22PRrs,1-f22f12-f22PRrs,2,CPr,IFs=f12f12-f22CPrs,1-f22f12-f22CPrs,2.Combining (1) and (2), the IF observations appear as(3)PRr,IFs=ρrs+cclr-cls+Br,PRIF-BPRs,IF+Mr,PRIFs+νr,PRIFs,CPr,IFs=ρrs+cclr-cls+Br,CPIF-BCPs,IF+Nr,IFsλIF+λIFϕr,IFs+Mr,CPIFs+νr,CPIFs.In the IF observation equation, the carrier phase ambiguity,Nr,IFs, does not keep integer nature any more. The receiver clock time offset clr 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.


## 2.1.2. 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](4)MWrs=1f1-f2f1CPrs,1-f2CPrs,2-1f1+f2f1PRrs,1+f2PRrs,2=λWLN1,rs-N2,rs+Br,MW-BMWs+νr,PRNLs+νr,CPWLs=λWLNWL,rs+Br,MW-BMWs+νr,MWs.


## 2.1.3. Geometry-Free Combination of GPS Observations

GF observable is formed via the difference between the dual-frequency observations, which is also independent of geometry, receiver clocks, and satellite clocks [11]:(5)PRr,GFs=PRrs,1-PRrs,2=Ir,GFs+Br,PRGF-BPRGFs+Mr,PRGFs+νr,PRGFs,CPr,GFs=CPrs,1-CPrs,2=-Ir,GFs+λGFfNr,GFs+Br,CPGF-BCPGFs+λGFϕr,GFs+Mr,CPGFs+νr,CPGFs.


## 2.2. 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].


## 2.3. 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](6)acentrifugal=-Ω×Ω×rECEF,aCoriolis=-2Ω×vECEF,where Ω is the angular velocity vector and r and v denote the position and velocity of the satellite in the ECEF frame.The motion of the satellite along with the timet is modeled by the conventional orbit dynamics in the Cartesian coordinates:(7)ddtxt=ft,xt,pt=vtat,rt,vt,where the p is the orbital parameter vector, including the solar radiation pressure and air drag coefficients (Cr and Cd). x comprises the position r and v vectors. a is the acceleration vector acting on the LEO satellite. The orbit can be propagated using the integration as follows:(8)xt=xt0+∫t0tfx,p,tdt,where t0 is the initial epoch. Differentiating x(t) with respect to x(t0) results in(9)∂xt∂xt0=I+∫t0t∂f∂xt∂xt∂xt0+∂f∂pt∂pt∂xt0dt.Note that∂p(t)/∂x(t0)=0; hence, the state transition matrix can be obtained by integrating (10) with an initial value equaling the identify matrix:(10)Φ˙t,t0=FΦt,t0,Φt0,t0=I,where (11)Φt,t0=∂xt∂xt0,F=∂f∂xt.Similarly, the sensitivity matrix can be calculated by integrating(12)Ψ˙t=FΨt0+C,Ψt0=0,where Ψ(t)=∂x(t)/∂p(t) denotes the sensitivity matrix and C=∂x(t)/∂p(t) denotes the Jacobian matrix with respect to the considered parameters. Comparatively, the sensitivity matrix is initialised with zero matrix Ψ(t0)=06×(6+2).


### 2.3.1. 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]:(13)ω˙t+1τωt=ut,where ω is the variable, τ is the correlation time, and u(t) is the white Gaussian noise. Hence, the solution of the variable ω can be expressed as follows:(14)ωt+Δt=ωte-Δt/τ,where Δt is the processing interval between two epochs. Define d=e-Δt/τ as the mapping function, and the process noise of the variable can be formulated as an uncorrelated random sequence with zero mean and variance:(15)Eu2=σ21-d2.


## 2.3.1. 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]:(13)ω˙t+1τωt=ut,where ω is the variable, τ is the correlation time, and u(t) is the white Gaussian noise. Hence, the solution of the variable ω can be expressed as follows:(14)ωt+Δt=ωte-Δt/τ,where Δt is the processing interval between two epochs. Define d=e-Δt/τ as the mapping function, and the process noise of the variable can be formulated as an uncorrelated random sequence with zero mean and variance:(15)Eu2=σ21-d2.


## 3. 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:(16)x=r;v;p;a;c·clr;Awhich consists of LEO satellite position vector r=x;y;z and velocity vector v=vx;vy;vz in ECEF coordinate frame, orbital parameters p=Cr;Cd, empirical acceleration vector a=aR;aI;aC in RIC (Radial, In-track, and Cross-track) coordinate frame, receiver clock offset c·clr, 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+n, where n is the number of visible GPS satellites.


### 3.1. Filter Initialisation

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,Cr=1.3 and Cd=2.3 for GRACE satellites. Other initial values for RPOD demonstration with GRACE data can be found in Table 3.


### 3.2. Time Update

In the time update step of EKF, both the states and their covariance matrix should be predicted:(17)x^k-=fx^k-1+,p^k-1+,0,P^k-=ΦkP^k-1+ΦkT+Q,where P is the state covariance matrix and Q is the process noise matrix; subscript k is the time epoch.To use RPOD approach, state transition matrixΦ and process noise covariance matrix Q can be expressed as(18)Φk=Φr,vΨCr,CdΨa000100000d·I3×300000100000In×n,Q=000000000000σa21-d2·I3×300000σc·clr2τc·clrΔt000000,where Φr,v, ΨCr,Cd, and Ψa can be integrated using the formulations in [25, 28]; d=e-Δt/τ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 first-order Gauss-Markov process, while the receiver clock offset is treated as a stochastic variable with an underlying random walk process model [28].


### 3.3. 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.


### 3.4. 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:(19)Kk=P^k-HkTHkP^k-HkT+R-1,x^k+=x^k-+Kkyk-hx^k-,0,P^k+=I-KkHkP^k-,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)H=∂h∂r∂h∂v∂h∂p∂h∂a∂h∂c·clr∂h∂A=losn×3T0n×30n×20n×31n×10n×nlosn×3T0n×30n×20n×31n×11n×n.


## 3.1. Filter Initialisation

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,Cr=1.3 and Cd=2.3 for GRACE satellites. Other initial values for RPOD demonstration with GRACE data can be found in Table 3.


## 3.2. Time Update

In the time update step of EKF, both the states and their covariance matrix should be predicted:(17)x^k-=fx^k-1+,p^k-1+,0,P^k-=ΦkP^k-1+ΦkT+Q,where P is the state covariance matrix and Q is the process noise matrix; subscript k is the time epoch.To use RPOD approach, state transition matrixΦ and process noise covariance matrix Q can be expressed as(18)Φk=Φr,vΨCr,CdΨa000100000d·I3×300000100000In×n,Q=000000000000σa21-d2·I3×300000σc·clr2τc·clrΔt000000,where Φr,v, ΨCr,Cd, and Ψa can be integrated using the formulations in [25, 28]; d=e-Δt/τ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 first-order Gauss-Markov process, while the receiver clock offset is treated as a stochastic variable with an underlying random walk process model [28].


## 3.3. 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.


## 3.4. 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:(19)Kk=P^k-HkTHkP^k-HkT+R-1,x^k+=x^k-+Kkyk-hx^k-,0,P^k+=I-KkHkP^k-,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)H=∂h∂r∂h∂v∂h∂p∂h∂a∂h∂c·clr∂h∂A=losn×3T0n×30n×20n×31n×10n×nlosn×3T0n×30n×20n×31n×11n×n.


## 4. 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:(21)PRIFs=ρrs+c·clr+Br,PRIF¯-c·cls+BPRIFs¯+νPRIFs,(22)CPIFs=ρrs+c·clr+Br,CPIF¯-c·cls+BCPIFs¯+NIFsλIF+νCPIFs.The clock offsets are written together with Hardware Delay Biases (HDBs). Antenna phase centre offsets and phase wind-up effects are not explicitly shown in the equation, but they are necessary to be corrected. GPS clock offsets c·cls+BPRIFs can be estimated using IGS clock products. Then GPS satellite and user receiver HDBs Br,CPIF, BCPIFs will be absorbed into ambiguity Nr,IFs, 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.


### 4.2.  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 double-differenced 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].(2)
Double-Differenced Ambiguities (DDA) [16].(3)
Decoupled Satellite Clocks (DSCs) [17].(4)
Integer Phase Clocks (IPCs) [18].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.


#### 4.2.1. IPCs Based PPP-AR Method

The basic concept of this method is to separate pseudorange clocks and carrier phase clocks; the float wide-lane ambiguitiesNWL 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 NNL; 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 (c·cls+Br,PRIF) in pseudorange measurement and CNES-CLS IPCs product is used for satellite clock offsets corrections (c·cls+BCPIFs) 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(23)AIFi,j=1λ11f12-f22f12A1i,jλ1-f22A2i,jλ2=f1f1+f2ANLi,j+f1f2f12-f22AWLi,j=f1f1+f2NNLi,j+BNLi,j+f1f2f12-f22NWLi,j+BWLi,j.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 wide-lane ambiguity resolution NWL depends on SDBS MW combination, which can be formed using (4):(24)MWi,j=λWLNWLi,j-BMWi,j+νMWi,j.Before separating NWL,ri,j from BMWi,j, a time-smoother is necessary to mitigate and eliminate multipath effects and thermal noises in MW observations as formulated in [14, 18](25)MW¯ni,j=MW¯n-1i,j+1nMWni,j-MW¯n-1i,j,where subscript n indicates the nth epoch, and the corresponding variance can be express as(26)σMW¯ni,j2=σMW¯n-1i,j2+1nMWni,j-MW¯n-1i,j2-σMW¯n-1i,j2.The wide-lane satellite biases BMWi,j can be extracted and calculated from CNES-CLS WSBs products. Since the wide-lane 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 R-ratio method (see (27)), the following formulation is working simultaneously to test and validate the right integer ambiguity resolution [14, 35]:(27)qNˇ′qNˇ≥Rthresh,(28)P0=1-∑i=1∞erfci-N-Nˇ2σ-erfci+N-Nˇ2σ,(29)erfcx=2π∫x∞e-t2dt,where N is the real ambiguity, σ is its variance, Nˇ is the nearest integer of float N, and Nˇ′ 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 P0 is larger than 1-α while condition in (27) should be satisfied as well (Rthresh 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 (c·cls+BCPIFs), hence, the estimated IF ambiguity is free of these biases. With known integer wide-lane ambiguity, (23) is derived to solve narrow-lane ambiguities:(30)ANLi,j=1λNLAIFi,jλIF-f22λ2f12-f22NWLi,j,where ANLi,j is the float SDBS narrow-lane ambiguity to be solved. The corresponding variance is given as(31)σANLi,j2=1λNL2σAIFi,j2.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 pseudo-observation 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 Figure1.Figure 1
IPCs based PPP-AR procedure.


#### 4.2.2. 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. Table1 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%.Table 1
IGS station solutions comparison.
PPP-AR strategy
IGS station
Static-PPP
Kinematic-PPP
E (m)
N (m)
U (m)
Fixing Ratio
E (m)
N (m)
U (m)
Fixing Ratio
Rounding
AUCK
0.019
0.024
0.038
97.3%
0.023
0.029
0.059
95.5%
No
0.020
0.024
0.039
0.0%
0.039
0.030
0.059
0.0%
Rounding
SHAO
0.028
0.025
0.085
94.5%
0.036
0.029
0.120
88.9%
No
0.028
0.025
0.086
0.0%
0.040
0.032
0.121
0.0%
Rounding
TLSE
0.009
0.019
0.028
97.4%
0.021
0.037
0.059
93.3%
No
0.013
0.020
0.032
0.0%
0.022
0.036
0.061
0.0%
Rounding
WTZR
0.010
0.025
0.020
97.3%
0.018
0.033
0.047
93.6%
No
0.015
0.026
0.020
0.0%
0.025
0.035
0.049
0.0%Figure2 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.Figure 2
WTZR positioning errors comparison using PPP-AR and float PPP strategies in the kinematic mode.
(a)
(b)
(c)


## 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:(21)PRIFs=ρrs+c·clr+Br,PRIF¯-c·cls+BPRIFs¯+νPRIFs,(22)CPIFs=ρrs+c·clr+Br,CPIF¯-c·cls+BCPIFs¯+NIFsλIF+νCPIFs.The clock offsets are written together with Hardware Delay Biases (HDBs). Antenna phase centre offsets and phase wind-up effects are not explicitly shown in the equation, but they are necessary to be corrected. GPS clock offsets c·cls+BPRIFs can be estimated using IGS clock products. Then GPS satellite and user receiver HDBs Br,CPIF, BCPIFs will be absorbed into ambiguity Nr,IFs, 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.


## 4.2.  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 double-differenced 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].(2)
Double-Differenced Ambiguities (DDA) [16].(3)
Decoupled Satellite Clocks (DSCs) [17].(4)
Integer Phase Clocks (IPCs) [18].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.


### 4.2.1. IPCs Based PPP-AR Method

The basic concept of this method is to separate pseudorange clocks and carrier phase clocks; the float wide-lane ambiguitiesNWL 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 NNL; 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 (c·cls+Br,PRIF) in pseudorange measurement and CNES-CLS IPCs product is used for satellite clock offsets corrections (c·cls+BCPIFs) 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(23)AIFi,j=1λ11f12-f22f12A1i,jλ1-f22A2i,jλ2=f1f1+f2ANLi,j+f1f2f12-f22AWLi,j=f1f1+f2NNLi,j+BNLi,j+f1f2f12-f22NWLi,j+BWLi,j.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 wide-lane ambiguity resolution NWL depends on SDBS MW combination, which can be formed using (4):(24)MWi,j=λWLNWLi,j-BMWi,j+νMWi,j.Before separating NWL,ri,j from BMWi,j, a time-smoother is necessary to mitigate and eliminate multipath effects and thermal noises in MW observations as formulated in [14, 18](25)MW¯ni,j=MW¯n-1i,j+1nMWni,j-MW¯n-1i,j,where subscript n indicates the nth epoch, and the corresponding variance can be express as(26)σMW¯ni,j2=σMW¯n-1i,j2+1nMWni,j-MW¯n-1i,j2-σMW¯n-1i,j2.The wide-lane satellite biases BMWi,j can be extracted and calculated from CNES-CLS WSBs products. Since the wide-lane 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 R-ratio method (see (27)), the following formulation is working simultaneously to test and validate the right integer ambiguity resolution [14, 35]:(27)qNˇ′qNˇ≥Rthresh,(28)P0=1-∑i=1∞erfci-N-Nˇ2σ-erfci+N-Nˇ2σ,(29)erfcx=2π∫x∞e-t2dt,where N is the real ambiguity, σ is its variance, Nˇ is the nearest integer of float N, and Nˇ′ 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 P0 is larger than 1-α while condition in (27) should be satisfied as well (Rthresh 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 (c·cls+BCPIFs), hence, the estimated IF ambiguity is free of these biases. With known integer wide-lane ambiguity, (23) is derived to solve narrow-lane ambiguities:(30)ANLi,j=1λNLAIFi,jλIF-f22λ2f12-f22NWLi,j,where ANLi,j is the float SDBS narrow-lane ambiguity to be solved. The corresponding variance is given as(31)σANLi,j2=1λNL2σAIFi,j2.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 pseudo-observation 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 Figure1.Figure 1
IPCs based PPP-AR procedure.


### 4.2.2. 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. Table1 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%.Table 1
IGS station solutions comparison.
PPP-AR strategy
IGS station
Static-PPP
Kinematic-PPP
E (m)
N (m)
U (m)
Fixing Ratio
E (m)
N (m)
U (m)
Fixing Ratio
Rounding
AUCK
0.019
0.024
0.038
97.3%
0.023
0.029
0.059
95.5%
No
0.020
0.024
0.039
0.0%
0.039
0.030
0.059
0.0%
Rounding
SHAO
0.028
0.025
0.085
94.5%
0.036
0.029
0.120
88.9%
No
0.028
0.025
0.086
0.0%
0.040
0.032
0.121
0.0%
Rounding
TLSE
0.009
0.019
0.028
97.4%
0.021
0.037
0.059
93.3%
No
0.013
0.020
0.032
0.0%
0.022
0.036
0.061
0.0%
Rounding
WTZR
0.010
0.025
0.020
97.3%
0.018
0.033
0.047
93.6%
No
0.015
0.026
0.020
0.0%
0.025
0.035
0.049
0.0%Figure2 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.Figure 2
WTZR positioning errors comparison using PPP-AR and float PPP strategies in the kinematic mode.
(a)
(b)
(c)


## 4.2.1. IPCs Based PPP-AR Method

The basic concept of this method is to separate pseudorange clocks and carrier phase clocks; the float wide-lane ambiguitiesNWL 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 NNL; 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 (c·cls+Br,PRIF) in pseudorange measurement and CNES-CLS IPCs product is used for satellite clock offsets corrections (c·cls+BCPIFs) 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(23)AIFi,j=1λ11f12-f22f12A1i,jλ1-f22A2i,jλ2=f1f1+f2ANLi,j+f1f2f12-f22AWLi,j=f1f1+f2NNLi,j+BNLi,j+f1f2f12-f22NWLi,j+BWLi,j.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 wide-lane ambiguity resolution NWL depends on SDBS MW combination, which can be formed using (4):(24)MWi,j=λWLNWLi,j-BMWi,j+νMWi,j.Before separating NWL,ri,j from BMWi,j, a time-smoother is necessary to mitigate and eliminate multipath effects and thermal noises in MW observations as formulated in [14, 18](25)MW¯ni,j=MW¯n-1i,j+1nMWni,j-MW¯n-1i,j,where subscript n indicates the nth epoch, and the corresponding variance can be express as(26)σMW¯ni,j2=σMW¯n-1i,j2+1nMWni,j-MW¯n-1i,j2-σMW¯n-1i,j2.The wide-lane satellite biases BMWi,j can be extracted and calculated from CNES-CLS WSBs products. Since the wide-lane 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 R-ratio method (see (27)), the following formulation is working simultaneously to test and validate the right integer ambiguity resolution [14, 35]:(27)qNˇ′qNˇ≥Rthresh,(28)P0=1-∑i=1∞erfci-N-Nˇ2σ-erfci+N-Nˇ2σ,(29)erfcx=2π∫x∞e-t2dt,where N is the real ambiguity, σ is its variance, Nˇ is the nearest integer of float N, and Nˇ′ 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 P0 is larger than 1-α while condition in (27) should be satisfied as well (Rthresh 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 (c·cls+BCPIFs), hence, the estimated IF ambiguity is free of these biases. With known integer wide-lane ambiguity, (23) is derived to solve narrow-lane ambiguities:(30)ANLi,j=1λNLAIFi,jλIF-f22λ2f12-f22NWLi,j,where ANLi,j is the float SDBS narrow-lane ambiguity to be solved. The corresponding variance is given as(31)σANLi,j2=1λNL2σAIFi,j2.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 pseudo-observation 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 Figure1.Figure 1
IPCs based PPP-AR procedure.


## 4.2.2. 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. Table1 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%.Table 1
IGS station solutions comparison.
PPP-AR strategy
IGS station
Static-PPP
Kinematic-PPP
E (m)
N (m)
U (m)
Fixing Ratio
E (m)
N (m)
U (m)
Fixing Ratio
Rounding
AUCK
0.019
0.024
0.038
97.3%
0.023
0.029
0.059
95.5%
No
0.020
0.024
0.039
0.0%
0.039
0.030
0.059
0.0%
Rounding
SHAO
0.028
0.025
0.085
94.5%
0.036
0.029
0.120
88.9%
No
0.028
0.025
0.086
0.0%
0.040
0.032
0.121
0.0%
Rounding
TLSE
0.009
0.019
0.028
97.4%
0.021
0.037
0.059
93.3%
No
0.013
0.020
0.032
0.0%
0.022
0.036
0.061
0.0%
Rounding
WTZR
0.010
0.025
0.020
97.3%
0.018
0.033
0.047
93.6%
No
0.015
0.026
0.020
0.0%
0.025
0.035
0.049
0.0%Figure2 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.Figure 2
WTZR positioning errors comparison using PPP-AR and float PPP strategies in the kinematic mode.
(a)
(b)
(c)


## 5. 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 (GeoForschungsZentrum, 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.Table 2
Orbit models used in this paper.
Satellite model
Cannonball model (constant surface and mass) [25]
Earth gravity field
GGM03S 75 × 75 [39]
Earth tides
k
20 Solid tides [40]
Atmospheric density
Harris-Priester [25]
Solar radiation pressure
Assuming the surface normal is always aligned with the sunIncluding eclipses and distance from the sun [25]
Planetary ephemerides
Low precision model [25]
Coordinates transformations
IERS1996/IAU1980 transformations [41]Table 3
RPOD EKF parameters settings.
Initial values
A priori variance
Correlation time
Process noise
Observation noise
a
R [nm/s2]
100
σ
r [m]
30
τ
a
(
R
,
I
,
C
) [s]
500
σ
a
R [nm/s2]
200
σ
P
R
IF [cm]
80
a
I [nm/s2]
100
σ
v [m/s]
0.1
τ
c
⋅
c
l
r [s]
100
σ
a
I [nm/s2]
600
σ
C
P
IF [cm]
0.8
a
C [nm/s2]
100
σ
C
r
0.2
—
—
σ
a
C [nm/s2]
400
—
—
C
r
1.3
σ
C
d
0.4
—
—
σ
c
⋅
c
l
r [m]
500
—
—
C
d
2.3
—
—
—
—
—
—
—
—RPOD solutions in RIC coordinates with and without phase ambiguity fixing using 10 hours of GRACE-B satellite data are shown in Figure3. 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 3
GRACE-B data RPOD errors comparison (10 hours on February 20, 2012).Figure4 (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.Figure 4
GRACE-B data RPOD carrier phase ambiguity resolution analysis (10 hours on February 20, 2012).The empirical accelerations estimated in RPOD are plotted in RIC coordinate system in Figure5. All the three components are several hundred nm/s2, 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.Figure 5
GRACE-B satellite empirical accelerations in RIC coordinate system (10 hours on February 20, 2012).The IF observation residuals are plotted in Figure6. The RMS of IF pseudorange observable residuals is 0.545 m, while that of IF carrier phase observable residuals is 8 mm.Figure 6
GRACE-B IF observations postfit residuals (10 hours on February 20, 2012).


## 6. Concluding Remarks

This paper presents PPP ambiguity resolution in reduced-dynamic 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.



---

*Source: 723414-2015-07-07.xml*