International Journal of Aerospace Engineering

Volume 2015 (2015), Article ID 723414, 11 pages

http://dx.doi.org/10.1155/2015/723414

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

^{1}School of Astronautics, Northwestern Polytechnical University (NPU), Xi’an 710072, China^{2}National Key Laboratory of Aerospace Flight Dynamics (AFDL), NPU, Xi’an 710072, China

Received 29 March 2015; Accepted 31 May 2015

Academic Editor: Paolo Tortora

Copyright © 2015 Yang Yang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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

#### 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]: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: : code phase and carrier phase observation (m). : geometric distance between receiver and GPS satellite antenna centres (m). : speed of light: 299 792 458 (m/s). : clock error, including relativistic effects (s). : ionospheric path delay (m). : tropospheric path delay (m). : hardware delay bias (m). : signal wavelength (m). : carrier phase integer ambiguity (cycle). : impact of relative rotation between transmitter and receiver antennas (cycle). : 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]:

Combining (1) and (2), the IF observations appear as

In the IF observation equation, the carrier phase ambiguity, , does not keep integer nature any more. The receiver clock time offset 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]

###### 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]:

##### 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]where is the angular velocity vector and and 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 is the orbital parameter vector, including the solar radiation pressure and air drag coefficients ( and ). comprises the position and vectors. is the acceleration vector acting on the LEO satellite. The orbit can be propagated using the integration as follows:where is the initial epoch. Differentiating with respect to results in

Note that ; 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 integratingwhere denotes the sensitivity matrix and denotes the Jacobian matrix with respect to the considered parameters. Comparatively, the sensitivity matrix is initialised with zero matrix .

###### 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]: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:

#### 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:which consists of LEO satellite position vector and velocity vector in ECEF coordinate frame, orbital parameters , empirical acceleration vector in RIC (Radial, In-track, and Cross-track) coordinate frame, receiver clock offset , and IF carrier phase ambiguity . Only ambiguities of visible GPS satellites at current epoch are to be estimated, so the total dimension of the states is , where 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, and 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:where is the state covariance matrix and is the process noise matrix; subscript is the time epoch.

To use RPOD approach, state transition matrix and process noise covariance matrix can be expressed aswhere , , and can be integrated using the formulations in [25, 28]; 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:where is the Kalman gain, is the design matrix of measurement equation, and is the observation noise matrix. is the observation, while is the calculated one. The formulation of matrix can be written as

#### 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: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 can be estimated using IGS clock products. Then GPS satellite and user receiver HDBs , will be absorbed into ambiguity , 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 ambiguities 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 ; 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 () in pseudorange measurement and CNES-CLS IPCs product is used for satellite clock offsets corrections () 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 inNote 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 depends on SDBS MW combination, which can be formed using (4):Before separating from , a time-smoother is necessary to mitigate and eliminate multipath effects and thermal noises in MW observations as formulated in [14, 18]where subscript indicates the th epoch, and the corresponding variance can be express asThe wide-lane satellite biases 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 -ratio method (see (27)), the following formulation is working simultaneously to test and validate the right integer ambiguity resolution [14, 35]:where is the real ambiguity, is its variance, is the nearest integer of float , and is the second nearest integer of float . Giving the confidence value , for example, , the float ambiguity could be fixed as the nearest integer when is larger than while condition in (27) should be satisfied as well ( 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 (), 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 is the float SDBS narrow-lane ambiguity to be solved. The corresponding variance is given asWith 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 Figure 1.