Chinese Journal of Aeronautics 24 (2011) 789-796

Contents lists available at ScienceDirect

Chinese Journal of Aeronautics

journal homepage: www.elsevier.com/locate/cja

journal of

aeronautics

Reduced Dynamic Orbit Determination Using Differenced Phase in Adjacent Epochs for Spaceborne Dual-frequency GPS

GU Defenga'b'*, YI Dongyunb

aState Key Laboratory of Astronautic Dynamics, Xi 'an 710043, China bCollege of Sciences, National University of Defense Technology, Changsha 410073, China Received 13 December 2010; revised 14 February 2011; accepted 19 April 2011

Abstract

We study on reduced dynamic orbit determination using differenced phase in adjacent epochs for spaceborne dual-frequency GPS. This method not only overcomes the shortcomings that the epoch-difference kinematic method cannot be used when observation geometry is poor or observations are insufficient, but also avoids solving the ambiguity in the zero-difference reduced dynamic method. As the epoch-difference method is not sensitive to the impact of phase cycle slips, it can lower the difficulty of slip detection in phase observation preprocessing. In the solution strategies, we solve the high-dimensional matrix computation problems by decomposing the long observation arc into a number of short arcs. By gravity recovery and climate experiment (GRACE) satellite orbit determination and compared with GeoForschungsZentrum (GFZ) post science orbit, for epoch-difference reduced dynamic method, the root mean squares (RMSs) of radial, transverse and normal components are 1.92 cm, 3.83 cm and 3.80 cm, and the RMS in three dimensions is 5.76 cm. The solution's accuracy is comparable to the zero-difference reduced dynamic method.

Keywords: dual-frequency GPS; phase difference in adjacent epochs; satellite; reduced dynamic; orbit determination

1. Introduction

Due to the characteristics of all-weather, high precision, and low cost, spaceborne dual-frequency GPS has many successful applications in satellite precise navigation and Earth gravity field recovery [1-2]. With the development of precise orbit determination technology using spaceborne dual-frequency GPS observations, more and more low earth orbit (LEO) satellites have been equipped with dual-frequency GPS receivers for precise navigation, and spaceborne

*Corresponding author. Tel.: +86-731-84573260.

E-mail address: gudefeng05@163.com

Foundation items: National Natural Science Foundation of China (61002033, 60902089); Open Research Fund of State Key Laboratory of Astronautic Dynamics (2011ADL-DW0103)

1000-9361/$ - see front matter © 2011 Elsevier Ltd. All rights reserved. doi: 10.1016/S1000-9361(11)60093-9

dual-frequency GPS has become one of the most important instruments in LEO satellite orbit determination. Dual-frequency GPS has successful applications in precise orbit determination for many satellite experiments, such as ocean topography experiment (TOPEX) [3], challenging mini satellite pay load (CHAMP) [4], gravity recovery and climate experiment (GRACE) [5-6], and gratifying results have been achieved. However, with the improvement in the quality of receiver and the requirement of higher precision orbit in scientific research, there are still many issues worthy of further research in LEO satellite orbit determination using spaceborne dual-frequency GPS.

At present, the LEO satellite orbit determination methods for spaceborne dual-frequency GPS mainly include double-difference and zero-difference method. The ground station observations are needed for double-difference method, which is more complicated than zero-difference method. With the precision of GPS

clock corrections and GPS ephemeris improvement, zero-difference method has been able to achieve the precision comparable to double-difference method [6]. Depending on whether using the orbit dynamic model information or not, the LEO orbit determination methods can be divided into kinematic method and dynamic method. In 1990, Yunck, et al. [7] combined the merits of kinematic method and dynamic method, and proposed reduced dynamic method, which has been widely studied and applied [5,8-9].

Phase observation is always used in high-precision GPS applications, but there is the impact of cycle slips and ambiguity in phase observation processing. In the high dynamic spaceborne environment, the probability of phase cycle slip occurrence is much larger than that on the ground, and high-speed movement of LEO satellites makes it more difficult to detect the slip [10]. At present, the cycle slip detection methods mainly include high order difference method, Melbourne-Wubbena combination method, ionosphere-free combination method, geometry-free combination method, and so on [11-13]. In the ground static measurement, the receiver will track a GPS continuously for a long time, so the signal changes slowly, which makes cycle slip detection easier. But in the high dynamic spaceborne environment, the high-speed LEO satellite flying makes the time of receiver's continuous tracking of a fixed GPS satellite shorter (ground measurements may be hours of continuous tracking, while for the CHAMP satellite, the average tracking time of one GPS satellite is only 30 min), the signal changes rapidly, and phase cycle slip detection becomes difficult. Therefore, small cycle slip detection has become a focus in spaceborne GPS application. Phase slip smaller than 5 cycles is hard to detect using ground static measurement method especially [14].

In zero-difference phase observation processing, if one phase cycle slip occurs and the corresponding epoch has not been detected exactly, then observation processing in the whole tracking arcs will be affected. Numerical experiments show that the artificial increase or decrease of a cycle slip by one will result in decimeter impact in single point positioning [13]. But for epoch-difference phase observations, the impact of an undetected cycle slip is limited to one or two epochs around the occurring cycle slips, and can be easily avoided.

Some scholars have studied "kinematic" orbit determination with the phase epoch-difference method. Using this method, simulation shows that decimeter-level precision of position components is possible [15]. Epoch-difference accuracy of the TOPEX satellites' orbit positions is about 2.5 cm in radial component and 1 cm in other two components [16]. Orbit position accuracy of the CHAMP and GRACE satellites is only about 30-40 cm [13,17-18]. The phase epoch-difference method will drop some observations which are not from the same satellite in two adjacent epochs, so observations are not fully utilized. On the other hand, the

kinematic method is particularly sensitive to erroneous measurements, unfavorable viewing geometry and data outages, which sometimes restrict its value in practice. In this paper, reduced dynamic orbit determination method using differenced phase in adjacent epochs is presented, which makes use of the known physical models of the spacecraft motion to constrain the resulting position estimates [19-20]. This allows an averaging of measurements from different epochs and the satellite trajectory can even be propagated across data gaps.

The difference introduces correlation between phase observations in adjacent epochs, which makes the weight matrix of observations and clock error design matrix no longer a diagonal matrix, so the long arc epoch-difference method faces the problem of high-dimensional matrix computation. In order to avoid high-dimensional matrix computation, this long arc observation data will be divided into several short arcs, and epoch-difference is taken inside each short arc, but not between two adjacent arc junctions. Then at the adjacent epoch between two short arc junctions, phase observations are not correlated, weight matrix of phase epoch-difference observation and design matrix of clock error become block diagonal matrices. This solution strategy effectively solves the high-dimensional matrix computing problems, and makes the phase epoch-difference method be successfully applied to satellite dynamic orbit determination.

2. Orbit Determination Model Using Differenced

Phase Observations

2.1. Observation equation

In order to eliminate the first order ionosphere delay, dual-frequency ionosphere-free combination observations are always adopted. For the pseudo code and carrier phase observations, the ionosphere-free combination yields

PF (ti) = —jf P (ti) - -¡f P (ti) =

f J 2 Jl J 2

p1 (ti) + c(8tr(t,) -St1 (ti)) + sp ,IF(ti)

Lf^) = -¡j L (ti) - -¡j L (ti) =

J1 J2 J1 J2

p (ti) + c(Str (ti) - St1 (ti)) + 4fAiF + eiiF (ti)

where subscript "IF" denotes ionosphere-free combination, subscripts "1" and "2" denote different frequencies, superscript "1 " denotes the 1th GPS satellite, P is code observations, L phase observations, f the carrier frequency, p geometric distance from the LEO satellite to the GPS satellite, Str receiver clock error, St GPS satellite clock error, c speed of light, and A,IFAIF the ambiguity of phase ionosphere-free combination,

£p,if and sL IF contain thermal measurement noise, mul-tipath, and all other unmodeled errors. The geometric range

P (ti ) =\r (ti ) - rG (ti -T )|

is simply given as the distance between the antenna phase center position of the GPS receiver, r(t), and the GPS satellite, r£. (t-tj), at the moment of signal reception and transmission, respectively. In Eq. (2), tj is the signal path delay which can be obtained by iterative calculation. The center for orbit determination in Europe (CODE) currently provides the final GPS ephemerides and clock error products [21], ephemeris interval of 15 min, clock error interval of 30 s and 5 s [22].

For continuous tracking arc, we make phase difference in adjacent epochs as follows:

hLLF(ti) = LF(ti+1) - LF(t,) =

P (ti+0 -P (ti) + c(5tr(t,.+,) -8tr(ti)) -

c(5tj (ti+,) + 5tj (ti)) +

Through epoch difference, AIF Aj , ambiguity of phase ionosphere-free combination, is eliminated. Even if a phase cycle slip occurs and the exact epoch has not been detected, the impact of undetected cycle slips is limited to one or two epochs when the cycle slips occur, and can be easily avoided through 3 a edition method of the observation minus computation (O-C) residuals in orbit determination. So phase epoch-difference method can lower the difficulty of phase observation preprocessing in zero-difference method. However, the difference introduces correlation between phase observations in adjacent epochs. In order to estimate orbit parameters successfully, we should pay additional attention to this problem.

2.2. Observation correction

(1) Relativity correction

The relativity correction formula is

2 J ua 2

APei =----e sin E = — (xx + yy + zz) (4)

where ju is the Earth gravitational constant, a semimajor axis of the satellite orbit, e orbit eccentricity, E orbit eccentric anomaly, and (x, y, z) orbit position.

(2) Antenna offset correction for GPS satellite

The antenna center offsets, Ax, Ay and Az, are defined in the GPS satellite body coordinate system. The axis vectors of the GPS satellite body coordinate system, ex, ey and ez, can be calculated in the Earth inertial coordinate system from Eq. (5)

e z = —

e y =\

ex = ey X ez

sion laboratory (JPL) solar system ephemeris. Then, the antenna offset correction vector for GPS satellite is

Ae + Ayfiy + Azez.

(3) Antenna offset correction for LEO satellite The range correction caused by the LEO satellite antenna offset is

APnt = :

rG - r

■ M -Ar

where Arbody is LEO satellite antenna offset vector in satellite body coordinate system, and M attitude rotation matrix from satellite body fixed coordinate to the Earth inertial coordinate.

2.3. Reduced dynamic orbit determination

The major difference between the kinematic method and the dynamic method is the fact that the individual spacecraft positions at each measurement epoch are replaced by the spacecraft trajectory model. We make use of known physical models of the spacecraft motion to constrain the resulting position estimates, and introduce three experienced acceleration components to absorb all other unmodeled perturbation. As the atmospheric density and solar activity are difficult to model accurately, atmospheric drag coefficient and solar radiation pressure coefficient are also estimated. Therefore, parameters to be estimated are: receiver clock error parameter at each epoch, cSt,-; six-dimensional initial satellite position and velocity vector, -«o = [r0T v0T]T; a solar pressure coefficient, CR; an atmospheric drag coefficient, CD, the piecewise linear empirical acceleration [24]

«, (t, ) = to +( j +1T-1'

[aR aT aN],+

t, - to - T

[aR aT a

nJ j+1

in consecutive time intervals (t0+jz, t0 + (j + 1)r) for compensating deficiencies in the applied dynamical models, and an independent set of empirical acceleration parameters (aR, aT, aN) is estimated for entire data arc, j = 0, 1,—, na.

When grouping the estimation parameters in the nT -dimensional vector of GPS receiver clock offsets

T = [Str(to) StT(fi)

Str(tnT -l)]T

where rs is the position of Sun provided by jet propul-

and the «^-dimensional vector concerning the satellite trajectory modeling, also referred to as the dynamic estimation parameters,

X = [CR CD a0 a1 L an„ ]T

The vector with all estimated parameters can be expressed as

y = [ XT Str(t0) Str(t1) l 8tr(tnr -1)]T

Observation equation can be expressed by h(y),

then, around initial values y0, the weighted least squares update value of y is

Ay = (HTWH)1 HTW(z -h(y0)) (7)

where z is observation vector, W — Q- the inverse of

observation weight matrix, and H=dh(y0)/dy0 the design matrix containing the linearized partial derivatives of the modeled measurements with respect to the estimation parameters. In accordance with the partitioned formulation, the design matrix H is split up into a part containing the modeled linearized measurement par-tials with respect to dynamic estimation parameters, HX, and the receiver clock, HT. The GPS observation model is linearized around initial values of the epoch clock offset X0 and dynamic estimation parameters T0:

T = To + AT,

X = X +AX

Here, the initial points T0 and X0 are given by the dynamical smooth filtering of discrete positions through single point positioning. The least squares solution of AT and AX is given by the following formula,

d (To, Xo) d(To, Xo)

d(T), Xo) and the overall design matrix.

W(z - h(To, Xo))

d(T), Xo)

= [ Ht HX ]

is constructed from the two partitioned ones. The least squares estimation is rewritten as

ht-WHx

HxWHt Hx WHx

which can be further reduced to

Ntt N TX "AT " nT

N XT N XX _ AX nX _

hTw ( z - h(To, X o)) hX W ( z - h(To, Xo))

It is now possible to firstly resolve the dynamic estimation parameter updates:

AX = (N XX — NXT NTT NTX ) (nX — NXTNTT nT)

which are subsequently back-substituted to obtain the updates for the clock offsets,

AT = NTlT(nT - NTX AX)

After having obtained the updates for the initial estimates, the newly obtained values are now used as initial values for a second run. Multiple Gauss-Newton

iterations of this kind are required to cope with the non-linearity of the reduced dynamic estimation problem, and convergence is typically achieved within 3 to 4 iterations.

As the clock offsets need to be estimated at each epoch, the dimension of matrix NTT associated with the clock offsets is very large. If the sampling interval of observation is 10 s, the number of clock offsets to be estimated will reach 8 640 for one day, and then the dimension of matrix Njt is 8 640 * 8 640. It is hard to get direct inversion of the full matrix with dimension 8 640 * 8 640. In the zero-difference dynamic orbit determination method, because observations at each epoch are independent, and matrix NTT is a diagonal matrix, the high-dimensional matrix inversion problem does not really exist[1,10]. But the epoch-difference method introduces correlation between phase observations in adjacent epoch. If zero-difference phase observation Lp is time independent and the variance is cr^ IF , then the weight matrix for epoch-difference phase observation is

Qz = cov

" A4 (ti ) '

ALf (t2 ) ALj (ts )

ALf (t -i)_

2 -1 -1 2 -1 -1 2

O O —1

—1 2

Through the phase difference in adjacent epochs, the observation weight matrix Qz is no longer a diagonal

matrix, and matrix

Ntt — HT Q. HT

is also no longer

a diagonal matrix. In order to avoid direct inversion of high-dimensional matrix in Eqs. (11)-(12), the long arc observation data will be divided into several short arcs, and epoch-difference is only taken inside a short arc, but not between two adjacent arc junctions. Then at the adjacent epoch between two short arc junctions, phase observations are not correlated, and observation co-variance matrix Qz and clock offsets design matrix HT become block diagonal matrices (see pig. 1). purther-more, matrix NTT becomes block diagonal matrix. Due to the partitioned formulation of the problem, Eqs. (11)-(12) can be solved more efficiently than solved by direct inversion of the full matrix.

On the one hand, in order to achieve better dimension reduction effect and improve the computational efficiency, a short arc length of each selection should not be too long. On the other hand, we make no epoch-difference at the link epoch between two adjacent short arcs, and this strategy causes observations to be not fully utilized, so the length of each short arc selection should not be too short. The length of each short arc selected in this paper is 1 h (about 360 epochs).

Fig. 1 Block partitioned matrix of matrices HT and Qz.

3. Dynamic Orbit Model and Parameter Selection

Satellite dynamic orbit equation in the geocentric inertial coordinate system can be generally expressed as

r =—r + a (t, r, r, P)

where r and r denote the satellite velocity and acceleration, respectively, P is the vector of dynamic orbit parameter, and a the sum of all perturbation acceleration, except for the two-body center acceleration of gravity.

a = aNS + aNB + aTD + aD + aSR + aRL + aRTN (15)

where aNS is the acceleration of the Earth non-spherical perturbation, aNB the 3rd body gravitational perturbation acceleration including the Sun and the Moon, aTD the acceleration of tidal perturbation mainly consisting of solid tide and the ocean tide, aD the atmospheric drag perturbation acceleration, aSR the acceleration of solar pressure perturbation, aRL the acceleration of relativity perturbation, and aRTN the empirical acceleration.

aRTN = aR eR + aT eT + aNeN

where eR, eT and eN are unit vectors, denoting the satellite radial, along-track, and cross-track flight directions, respectively. The empirical accelerations are considered to be piecewise linear in pre-defined sub-intervals. The entire data arc is divided into na intervals of equal duration t. Intervals of 900 s duration have been selected in this paper.

The models and parameter information used for processing reduced dynamic orbits are in Table 1.

Table 1 Dynamic orbit model and parameters

Item Description

Static gravity field GGM02C 150x150

Solid Earth tide IERS96, 4x4

Polar tide IERS96

Ocean tide CSR4.0

The 3rd body gravity Sun and Moon

Solar radiation pressure Ball model, conical Earth shadow, CR is

estimated

Jacchia 71density model (National

Oceanic and Atmospheric Administra-

Atmospheric drag tion (NOAA) solar flux (daily) and

geomagnetic activity (3 hourly)), CD is

estimated

Relativity Schwarzs child

Precession IAU1976

Nutation IAU1980 + EOPC correction

Earth orientation E0PC04

Solar ephemerides JPL DE405

Orbit determination process is typically conducted in single day (24 h) data batches, where the measurements are processed in 10 s steps. The selected coordinate system is J2000 inertial reference frame, and the selected time system is terrestrial dynamic time (TDT). The Adams-Cowell multi-step integration method is used for orbit integration.

4. Numerical Examples

In reduced dynamic orbit determination with zero-difference phase observations, if a phase cycle slip occurs and the exact epoch has not been detected, the phase observation processing in the whole tracking arcs will be affected. But this impact can be easily avoided by epoch-difference method. Fig. 2 gives a numerical example to show the impact of a cycle slip in dynamic orbit determination of LEO satellite. At the epoch t = 2.5 h, a size of 1 cycle slip is added to the phase observations. Compared to the situation without

Fig. 2 Cycle slip effect on reduced dynamic orbit determination.

cycle slips, the orbit determination precision of phase zero-difference method in the entire arc is obviously deteriorated, and the maximum impact reaches 8 cm near the epoch where the cycle slip occurs. But such big impact cannot be found in phase epoch-difference method, and the orbit determination precision of phase epoch-difference method is just the same as that with phase zero-difference method without cycle slip (see Fig. 2).

Over a period of 31 days from January 1 to 31, 2006, GRACE A dual-frequency GPS observations are processed, and root of mean square (RMS) value of the O-C residuals obtained from reduced dynamic orbit determination using differenced phase is only 6.7 mm (see Fig. 3), which reflects the consistency of the applied models with the GPS observations. The differences between the orbit obtained from reduced dynamic orbit determination using differenced phase and the GeoForschungsZentrum (GFZ) post science orbit are computed in the radial, along-track and cross-track directions at the discrete epoch (see Fig. 4).

For epoch-difference method, the short arc length should not be too long otherwise the dimension of matrix NTT will be too large, meanwhile it cannot be too short otherwise the orbit determination precision will be deteriorated. The reasonable length of each selected short arc is 1 h (see Table 2).

Fig. 4

Comparison between GFZ science orbit and reduced dynamic orbit determination result using differenced phase for GRACE A on January 2, 2006. The RMSs of R, T and N components are 2.1 cm, 3.8 cm and 2.9 cm.

Table 2 Orbit determination with different short arc length for GRACE A on January 2, 2006

Short arc length ATT dimension Orbit precision/cm

300 s 30x30 7.89

900 s 90x90 6.43

0.5 h 180x180 5.60

1 h 360x360 5.30

2 h 720x720 5.24

3 h 1 080x1 080 5.20

Fig. 3 Phase O-C residual of reduced dynamic orbit determination using differenced phase for GRACE A on January 2, 2006 (RMS=6.7 mm).

Over the 31-day period, compared with GFZ post science orbit, for phase epoch-difference reduced dynamic orbit determination method, the RMSs of radial (R), transverse (T) and normal (AO position components are 1.92 cm, 3.83 cm and 3.80 cm, and the RMS in three dimensions is 5.76 cm (see Fig. 5). The solution's accuracy is comparable to the zero-difference reduced dy-

Fig. 5 Comparison between GFZ science orbit and reduced dynamic orbit determination results using differenced phase for GRACE A in January 2006 (Average RMS = 5.76 cm).

Fig. 6 Comparison between GFZ science orbit and zero-difference reduced dynamic orbit determination results for GRACE A in January 2006 (Average RMS = 5.53 cm).

namic method, for which, the RMSs in R, T and N position components are 1.96 cm, 3.41 cm and 3.85 cm, and the RMS in three dimensions is 5.53 cm (see Fig. 6). The orbit determination accuracy of CHAMP and GRACE satellites with phase epoch-difference kinematic orbit determination method is about 3040 cm [13,18]. Compared with kinematic orbit determination method, reduced dynamic orbit determination method can significantly improve orbit accuracy.

5. Conclusions

(1) The proposed reduced dynamic orbit determination method using differenced phase in adjacent epochs for spaceborne dual-frequency GPS is a new attempt. Compared with phase zero-difference method, phase epoch-difference method avoids solving the ambiguity, and is not sensitive to the impact of phase cycle slips, which can lower the difficulty of phase observation preprocessing.

(2) In the solution strategies, the high-dimensional matrix computation problems are solved by decomposing the long observation arc into a number of short arcs, making epoch-difference inside each short arc, but not between two adjacent arc junctions. The reasonable length of each selected short arc is 1 h.

(3) According to the proposed phase epoch-difference reduced dynamic orbit determination method, the GRACE satellite orbit determination results show that the orbit position accuracy of R, T and N position components is 1.92 cm, 3.83 cm and 3.80 cm, and the RMS in three dimensions is 5.76 cm. The solution's accuracy is comparable to the zero-difference reduced dynamic method, for which, the RMSs of R, T and N position components are 1.96 cm, 3.41 cm and 3.85 cm, and the RMS in three dimensions is 5.53 cm.

References

[1] Kroes R. Precise relative positioning of formation flying spacecraft using GPS. PhD thesis, Delft University

of Technology, 2006.

[2] Schwintzer P, Reigber C. The contribution of GPS flight receivers to global gravity field recovery. Journal of Global Positioning Systems 2002; 1(1): 61-63.

[3] Kang Z, Schwintzer P, Reigber C, et al. Precise orbit determination for TOPEX/Poseidon using GPS-SST data. Advances in Space Research 1995; 16(12): 59-62.

[4] van den Ljssel J, Visser P, Patino Rodriguez E. Champ precise orbit determination using GPS data. Advances in Space Research 2003; 31(8): 1889-1895.

[5] Kang Z, Tapley B, Bettadpur S, et al. Precise orbit determination for the GRACE mission using only GPS data. Journal of Geodesy 2006; 80(6): 322-331.

[6] Jäggi A, Hugentobler U, Bock H, et al. Precise orbit determination for GRACE using undifferenced or doubly differenced GPS data. Advances in Space Research 2007; 39(10): 1612-1619.

[7] Yunck T P, Wu S C, Wu J T, et al. Precise tracking of remote sensing satellites with the global positioning system. IEEE Transactions on Geoscience and Remote Sensing 1990; 28(1): 108-116.

[8] Svehla D, Rothacher M. Kinematic and reduced-dynamic precise orbit determination of low Earth orbiters. Advances in Geosciences 2003; 1: 47-56.

[9] Montenbruck O, van Helleputte T, Kroes R, et al. Reduced dynamic orbit determination using GPS code and carrier measurements. Aerospace Science and Technology 2005; 9(3): 261-271.

[10] Gu D F. The spatial states measurement and estimation of distributed InSAR satellite system. PhD thesis, National University of Defense Technology, 2009. [in Chinese]

[11] Liu J Y. GPS satellite navigation/positioning theories and methods. Beijing: Science Press, 2007; 319-329. [in Chinese]

[12] Kong Q L, Ou J K, Chai Y J. Detection and repair of gross errors and cycle slips in LEO-based GPS data on zero level. Journal of Geodesy and Geodynamics 2005; 25(4): 105-109. [in Chinese]

[13] Zheng Z Y. Study and software implementation of GPS data pre-processing and onboard GPS kinematic orbit determination. PhD thesis, Shanghai Astronomical Observatory, Chinese Academy of Sciences, 2004. [in Chinese]

[14] Bae T S, Kwon J H, Grejner-Brzezinska D A. Data

screening and quality analysis for kinematic orbit determination of CHAMP satellite. ION Technical Meeting. 2002; 771-777.

[15] Bisnath S B, Langley R B. Precise a posteriori geometric tracking of low Earth orbiters with GPS[J]. Canadian Aeronautics and Space Journal 1999; 45(3): 245-252.

[16] Bock H, Hugentobler U, Springer T, et al. Efficient precise orbit determination of LEO satellites using GPS. Advances in Space Research 2002; 30(2): 295- 300.

[17] Wu J F, Du P, Wang L, et al. Methods and progress on kinematic orbit determination of LEO based on GPS. Progress in Astronomy 2006; 24(2): 113-128. [in Chinese]

[18] Chen J P, Wang J X. Kinematic precise orbit determination of low Earth orbiter based on epoch-difference strategy. Journal of Geodesy and Geodynamics 2007; 27(4): 57-61. [in Chinese]

[19] Peng D J, Wu B. Zero-difference and single-difference precise orbit determination for LEO using GPS. Chinese Science Bulletin 2007; 52(6): 715-719. [in Chinese]

[20] Gu D F, Zhu S B, Yi D Y. CDGPS relative position determination of distributed SAR satellite formation based on dynamic orbit model. Chinese Journal of Space Science 2009; 29(5): 515-521. [in Chinese]

[21] Dach R, Brockmann E, Schaer S, et al. GNSS procesing at CODE: status report. Journal of Geodesy 2009; 83

(3-4): 353-366.

[22] Bock H, Dach R, Jäggi A, et al. High-rate GPS clock corrections from CODE: support of 1 Hz applications. Journal of Geodesy 2009; 83(11): 1083-1094.

[23] Xu G. GPS: theory, algorithms and applications. Heidelberg: Springer Verlag, 2003; 97.

[24] Jäggi A, Hugentobler U, Beutler G. Pseudo-stochastic orbit modeling technique for low Earth orbiters. Journal of Geodesy 2006; 80(1): 47-60.

Biographies:

GU Defeng Born in 1980, he received B.S. and Ph.D. degrees from National University of Defense Technology in 2003 and 2009 respectively, and then became a teacher there. His main research interests are modern statistical method and measurement data processing. E-mail: gudefeng05@163.com

YI Dongyun Born in 1965, he received B.S. degree from Nankai University and Ph.D. degree from National University of Defense Technology in 2003. Since then he has been a professor at the Department of Mathematics and System Science, National University of Defense Technology. His main research interests are dynamic system analysis and measurement data processing. E-mail: dongyunyi@nudt.edu.cn