計測自動制御学会産業論文集
Vol.10, No.20, 167/174 (2011)

特集「移動体の誘導制御」

A Precise Point Positioning Filtering Method
for Navigation and Orbit Determination of Low Earth Orbiters
Takahiro YOSHIOKA** and Masaaki MURATA*A new precise point positioning (PPP) technique for on-orbit navigation (NAV) and precision orbit determination (POD) of
LEO satellites is proposed. A dual-frequency PPP filter is developed by assuming simple white noise processes for both position
deviation (from the pre-determined reference orbit) and receiver clock. The filter performance is assessed by processing one-day
GPS data collected with a dual-frequency receiver onboard the LEO satellite CHAMP. The results show that dual-frequency
PPP in space is capable of producing NAV solution with 3D positional accuracy of about one meter, and POD solution of about
0.6 meters, both in 95 percentile.
Key Words: PPP, GPS, satellite navigation, EKF, UKF, CHAMP

1. Introduction
A new technique called Precise Point Positioning (PPP)
has been attracting precision geo-science and navigation
community, as it would have a potential of providing highly
accurate point positioning with only a single GPS receiver.
The PPP accuracy, which is almost comparable to that ofkinematic GPS (KGPS), has been reported for terrestrial
and aerial navigation users1)~3). Such a PPP technique
became possible with the availability of precise orbits andclock correction information for the GPS satellites, which
the International GNSS Service (IGS) has begun generating
and disseminating since 1994 as an international civil
cooperation. Compared to DGPS techniques4), PPP has theobvious advantages of simplified system design,operational flexibility (i.e., no specific operational
restrictions), and cost reduction. PPP could therefore be avital alternative to the conventional carrier phase DGPS.
For space PPP, a number of studies on LEO precision
orbit determination (POD) in post-mission mode have been
carried out5)~8), but very few studies have focused on realtime on-orbit navigation application (NAV) for LEO
satellites. In this paper we propose a simple and efficientapproach to LEO NAV as well as LEO POD, which relies
solely on raw un-differenced pseudorange and carrier phasedata. Therefore it differs from the previous investigations(primarily dedicated to LEO POD) in the algorithm: forexample, Bock et al.6) used un-differenced pseudorange and* Professor and Corresponding Author, Aerospace
Engineering, National Defense Academy, Yokosuka,
Kanagawa (E mail to: murata@gnss.co.jp)** Graduate Student, Aerospace Engineering, National
Defense Academy, Yokosuka, Kanagawa(Received October 12, 2010)

epoch-by-epoch carrier phase difference data, and also
Bisnath et al.7) used double-differenced pseudorange andtriple-differenced carrier phase data.
In 2006 and 2007, we developed a basic PPP software at
National Defense Academy (NDA) for terrestrial and aerial
vehicular navigation applications and assessed its
performance indicating the horizontal accuracy at a few
centimeters for static antennas and at better than a few
decimeters for moving antennas9)~10). This PPP software
was intensively modified for space navigation users andevaluated by applying to the 24-hour dataset collected on
CHAMP, German LEO satellite. This paper gives design
philosophy of a dual-frequency PPP filter and analysis
results of positioning performance for both NAV and POD.
2. PPP Filter Design
Implicitly assumed is that the onboard receiver is dualfrequency so that ionospheric-free pseudorange and carrier
phase data are constructed. In that case the observation
equations are described in the form of the linear
combinations for each satellite:P

f12 P1  f 22 P2
f12  f 22

(1)r
   c(dt   r )  c(dt SV   SV
)   PRCf12 C1  f 22 C 2
f12  f 22

(2)

~
r
   c(dt   r )  c(dt SV   SV
)  N   CP

~ cf N  cf 2 N 2
N  1 12f1  f 22

C 2010SICE
TRIA 020/11/1020○(3)

where Pi

for Gaussianity with zero mean, and that it is only anassumption made for technical convenience. The othermodels, such as random walk, first-order Markov model,
etc., could of course be used but the improvement inpositional performance will be marginal.
The state equations are thus described in a simple linear
discrete form as follows:is the measured pseudorange on carrier

wave Li (i=1,2); C i is the measured carrier phase on

Li ;  is the true geometric range; c is the speed of light;
dt is the receiver clock bias;  r is the relativistic effectr
on receiver clock; dt SV is the transmitter clock bias;  SV is
the relativistic effect on transmitter clock; f i is thefrequency of Li ; N i is the integer phase ambiguity on~
Li ; N is the floating ambiguity; and  i represents the 
rj 1  wr , j , wr , j ~ N (0, qr I 33 )(5)

b j 1  wb , j , wb , j ~ N (0, qb )




~
~
N j 1  N j  wN~ , j , wN~ , j ~ N (0, q N~ I1212 )residual measurement error including random noise. The
measurement error is assumed to be normally distributedwith mean zero and elevation-dependent standard deviation

(6)

(7)

(SD) of k el   0,ksin el 2  a k2

, k  P, C

where I nn denotes an n×n unit matrix. The magnitude of(4)
covariance matrices qr , qb and q N~ for the white noises

where  0 ,k is a nominal SD along the zenith direction, el
is the elevation angle above local horizon, and a k is atuning parameter of constant (set zero throughout thisstudy). Navigation users can use the precise orbits andclocks produced by IGS. Currently the IGS's official
products for precise orbits and clocks are in three
categories: UltraRapid, Rapid, and Final. Among theseproducts, only the UltraRapid orbit and clock can be used inreal-time. In this paper, however, we used the Final productof the 15-minute precise orbits and 30-second clock
corrections, which were computed at the CODE (Center for
Orbit Determination in Europe), one of the IGS analysis
centers. When necessary they were interpolated using
Lagrange's polynomial, that is, the 9th order polynomial fororbits and the 2nd order polynomial for clock corrections.
Also, the correction for the center-of-mass offset to thetransmitting antenna was applied (Fig.1).
The unknown parameters in the PPP filter nominally
wr , wb

and wN~ istuningparameters.

Hereafter,
x ~ N (m, P) means that x is normally distributed with

mean m and covariance matrix P . As the means ofconstructing the PPP filter, we chose the conventional
Extended Kalman Filter (EKF)11)~12) and the Unscented
Kalman Filter (UKF)13) , an alternative to the EKF.
3. Analysis Results
3.1 LEO Data
The CHAMP (CHAllenging Mini-satellite Payload forgeophysical research and application) is a German LEO
satellite designed for geo-science and atmospheric research
and managed by the GeoForschungs Zentrum (GFZ). Its
main scientific mission goals are to provide highly precise
gravity
consist of the position components ( r ) of the receiverantenna, a receiver clock bias ( b ), and ionospheric-free
~(non-integer) ambiguities ( N , max.12) of the un-field,

magnetic

field,

radio

occultation

measurements, and GPS altimetry over a 5-year period.
Fig.1 shows the satellite configuration with the location of
the GPS receiver antennae for POD (and NAV). One-day
dataset collected on November 18, 2006, were focused indifferenced carrier phase data. The position componentsthe present paper. The data included dual-frequency( r ) are the positional deviation from a pre-determinedreference trajectory of the LEO satellite. The float
ambiguities are modeled as a piece-wise constant corrupted
with white noises. All unknowns except for the float
ambiguities are modeled as Gaussian white noise processes
with mean zero and specified variances. The clock model
like this is often called a white noise clock. The white noisemodels are the simplest mathematically and so with lessdemanding for computation. Note that there is no rationale

pseudorange and carrier phase taken by a "BlackJack" 12channel receiver. They were unedited with the samplinginterval of 10 seconds. While the cataloged values of raw,
un-differenced data precision are 10-20 cm in pseudorange
and 0.05 cm in carrier phase, we assigned the values of 50
cm and 1 cm, respectively, as the nominal observation error
standard deviation (see Eq.(4)) for the ionosphere-freecombination in the zenith direction, by taking into account168

conventional terrestrial system (CTS) or equivalently in the
Earth-centered, Earth-fixed frame (ECEF). Here thecoordinate solution is derived by using the pseudorange
data only. We call this set of the temporal coordinate
solution as the pseudo-observations. Satellite dynamics forthe batch filter is expressed in the J2000 inertial system as

the effect of the remaining residual errors. It should be
noted that no attitude information was available for thisresearch. At that time the CHAMP orbited at an altitude of
350 km, with a period of about 90 minutes, an eccentricity
of smaller than 0.001, and an inclination of 89 degrees. The
corresponding satellite visibility is plotted in Fig.2. TheR    R    3  f (t , R, R, p)
R

upper (thick) skyline shows the maximum number of GPS
satellites all-in-view from CHAMP’s receiver, indicatingthat approximately 8.2 satellites on average were observed

where t is the Terrestrial Time (TT), μis the gravitational

at an epoch. The lower (thin) line shows the number of GPSＸ(8)parameter of the Earth, p is the unknown constantparameter vector to be estimated, and f denotes the

GPS Antenna

perturbation acceleration vector. For the Earth gravity
computation the NASA's Earth Gravitational Model 1996
(EGM96), truncated up to degree and order 36, was adoptedoin this research. For the atmospheric density, we can use

ＺＹeither exponentially decaying model or the Jacchia-Roberts
1971 model with tabular data consisting of F10 .7 solar flux

Center of Gravity

and 3-hourly planetary geomagnetic index K P . The Earthorientation parameters (EOP) are required as well. The

Fig. 1 CHAMP configuration (GFZ Web)other forces acting on the satellites, such as solid Earth tide,
ocean tide, and earth radiation pressure and relativistic
effects, are marginal for a relatively short arc of LEO
satellites and hence are not taken into account.
The observation equation is thus described asrk , obs  rk ,true   k

where  k ~ N (0, Rk )(9)(k  1,2,...) . The single pointpositioning solution vector rk , obs is calculated in the ECEF
Fig.2 Satellite visibility

frame. The 3×3 positioning error covariance matrix Rksatellites with the elevation angle of less than 10 degrees,
approximately 0.7 satellites on average.

varies with satellite geometry, but we assume they are a

3.2 Reference Orbit Generation
Satellite missions for LEO require real-time or post-flightorbit determination (OD). As usual, post-flight OD is
performed by the use of large-scale orbit determination
program. On the other hand, the proposed PPP filter is
executed assuming a reference orbit with moderate
accuracy. For the case that any reference orbit is notavailable in advance, we have proposed14) to generate itusing a classical least-squares in which the observation data
are an epoch-to-epoch sequence of the single point
positioning solution, i.e. the 3D position coordinates in the

In summary, a reference orbit can be calculated byapplying the following two-step procedures: first, 1) single
point positioning solution (3D coordinates) at each epoch is
obtained by the conventional least-squares in which precise
orbits and clock corrections are used; and secondly, 2) theorbit is dynamically fitted to the 3D coordinates used as thepseudo-observations.
The above OD scenario was evaluated with varying
dynamic models based on the 24-hour CHAMP data, andthe resulting orbit accuracy was investigated by the direct
difference from the GFZ precise ephemeris. As an example,

diagonal  P I 33 where  P = 1 m is adopted.
2169

Fig.3 shows the time history of the orbit difference between
the generated reference orbit (no drag included) and thetruth (i.e. GFZ's precise orbit). The orbit difference fromthe truth is expressed in the local radial-transverse-normal(RTN) frame (see Appendix A). The error in the transverse
direction is the biggest, being the maximum error of about
600 m towards the initial and final epochs. Obviously, thisis due to the lack of atmospheric drag force model in the
OD process. This orbit was used in this research as thereference orbit for the PPP filter.

empirically or through intensive parameter tuning on thebasis of running simulation program of the PPP filter. The
selected values are

qr  qb  106 m2/s and qN~  1010 m2/s for Q, and  1,   2,  3  L for UKF, where L denotes thedimension of the filter state (see Appendix B for details) .
The kinematic CHAMP orbit solution is very sensitive todata editing schemes applied, and therefore online data
screening is a prerequisite for the NAV and POD filters.
The editing strategies implemented in this paper are shownin Table 1. The Melbourne-Wuebbenna (M-W) linear
combination is defined as wide-lane carrier phase minus
narrow-lane pseudorange, i.e.
~
 

1
1
( f11  f 2 2 ) ( f1 P1  f 2 P2 )f1  f 2
f1  f 2(13)c( N1  N 2 )f1  f 2

and this eliminates the effect of the ionosphere, thegeometry, the clocks and the troposphere. Thus, thiscombination can be used to detect cycle slips in the undifferenced carrier phase data. These editing strategies areable to screen the data in the observation domain, with noneed of post-fit residual analysis. Strictly speaking, models
for precise orbits and precise clocks have to be consistentFig.3 Position difference between the reference orbit and
GFZ precise orbit
3.3 PPP Results
The 24 hours of observation data were processed by the
EKF and the UKF, assuming a nominal model setup andwith the mask angle of 10 degrees. The time update portionof the EKF takes the form as follows:x j 1   j xˆ j(10)

Pj 1   j Pˆj Tj  Q

(11)

 
where x  (r T

Table 1 Edit criteria for NAV and POD pre-processing
Criterion
Content
Satellite status*
Satellite health, and SV accuracy(from navigation message)***
Cutoff angle
Reject data with elevation angle lessthan 10°
Radio signal
Reject data with SNR < 6 (RINEX
strength*
standard)
ErroneousReject unless PRmin < PR < PRmax
Pseudorange*
Residuals*
Reject if the moving average of theresiduals weighted with thereciprocal of the estimation error
covariance is larger than thespecified level.
Cycle slip*
Reset the filter if the time differenceof Melboune-Wuebbena linear
combination is larger than thespecified tolerance.
Isolated epoch of
Removed by hand (only for PODdatadata).

~b N T )T is 16-dimensional and qr I 33

Q 0
 0


0
qb
0



0 
q N~ I1212 
0(12)In Eqs. (10) and (11) (and also in Appendix B), the upper
bar and hat denote the a priori and a posteriori estimate
respectively, and the transition matrix  j is the identity
matrix, except for the position and clock bias portions thatare zero for the above baseline white noise model. The* Applied for both NAV and POD data.
** PRN15, decommisioned on 17 Nov. 2006

process noise intensity matrix Q can be set up either

170

Table 2 Positional error of NAV-EKF (Unit:m)parameter dRDTdN
3D rssmedian
0.15
0.13
0.14
0.31
95%
0.87
0.73
0.40
1.13

each other, and it was satisfied in this paper by the use of
15-minute orbits and 30-second clocks both computed by
the CODE. Positioning errors between the EKF navigation
solution and the truth, i.e. precise dynamic orbit of the
GeoForschungsZentrum (GFZ), both expressed in the local
radial-transverse-normal (RTN) frame, were computed inwhich the center-of-mass offset correction against receiver's
phase center was applied by assuming that the spacecraft
attitude is nominally aligned along the RTN directions.

Table 3 Positional error of NAV-UKF (Unit:m)parameter dRDTdN
3D rssmedian
0.15
0.12
0.14
0.30
95%
0.75
0.76
0.38
1.07

Fig. 5 Histogram of 3D rss positioning errors(a) dR

For the purpose of saving paper space, only the plots of the
NAV-EKF results are shown in Fig.4. Fig.4 (a)-(c) plot thepositioning errors, i.e. dR, dT, and dN at each of the RTN
directions, and Fig.4 (d) shows the 3-dimensional root-sumsquares (3D rss). At first it should be noticed that, opposite
to very slow convergence found for the terrestrial or aerial
PPP navigation9)~10), the space PPP filter is converging
much faster, within several number of observation epochs.
This is most likely due to satellite's orbital motion causing
rapid change in the measurement geometry. Second, there
are many spikes in all plots, some exceeding hundred
meters, and they are mostly due to the frequent filter reinitialization when a new satellite is coming up or a cycle
slip occurs. In particular, the procedure when a cycleslip
has occurred goes as follows: calculate first the predictionof the new ambiguity (with use of the M-W linear
combination); and replace the corresponding filter statecomponent with the predicted value and at the same time
initialize the corresponding component of the covariancematrix; and then incorporate that carrier phase data into thefilter.
Third, the POD result does not change so much from theNAV one, but with rather less spikes because of additional
hand editing of separated data and short passes. Fig.5
shows the histogram of positioning errors for NAV case.
Table 2 and 3 tabulate the median and the 95-percentile
point of positional errors in each direction. For NAV
solution, positioning accuracy of less than one meter inboth radial and transverse components and of less than one-(b) dT(c) dN(d) 3D rss
Fig. 4 NAV-EKF positioning errors

171half meters in the normal component, with 3D accuracy of
about one meter, can be achieved (all in the sense of 95
percentile). The NAV solution is not biased so much, at thelevel of smaller than a decimeter. For POD, positioning
accuracy of about one-half meters in all components, with
3D accuracy of about 0.6 meters, can be achieved (in 95
percentile). As seen from Table 3, the UKF provides
slightly better positioning results than the EKF, but by
paying the extra computational cost. Indeed, it is wellrecognized that the UKF in its standard form takes much
more computational burden than the conventional EKF.
Fig.6, a microscopic view of Fig.4 focused over the shortperiod (18 min.) of high quality estimation, shows that theposition errors in each direction are more or less oscillating
in a saw-shape manner with a period of about 150 s andwith a peak-to-peak range of 20 cm for dR, some cm for dT
and dN. The oscillatory variation in the radial direction
suggests pitching motion around the center of mass of thespacecraft, although the cause of 150 s periodicity is not
identified. The result evidences superior accuracy of the
PPP filters. Fig.7 compares convergence property of two
filters before and after reset over a period of data gap.
Before reset, both filters are in a steady-state with almost
equal position errors. Immediately after reset, position
errors increase abruptly but decrease gradually as
measurements are processed. Clearly, the UKF is moreaccurate and faster in convergence than the EKF. This isbecause of the second order accuracy of the UKF
unconditionally.
4. Conclusions
A design methodology of a PPP filter for LEO navigationand orbit determination is proposed and its performance isevaluated in post-mission mode using 24 hours of GPS
tracking data collected onboard CHAMP. The potentiality
and effectiveness of the PPP filtering approach aredemonstrated: the 3D rss positioning accuracy (95percentile) of about one meter for NAV and of about 0.6
meters for POD can be achieved. The filter generally shows
fast convergence due primarily to rapid change of thegeometry. However, data editing scheme applied online for
NAV or offline for POD extremely affects the positioning
performance and should be paid special efforts in its
implementation.
Acknowledgements
We thank GeoForschungsZentrum (GFZ), Germany, for
our access to the CHAMP data, the precise science orbit,
and the additional information at the CHAMP Information
System and Data Center (ISDC). In particular, the second
author is grateful to Dr.C.Reigber and Dr.G.Gendt of GFZ
for their cooperation with informative discussions.
References
1) Y.Gao and K.Chen: Performance Analysis of Precise
Point Positioning Using Real-Time Orbit and Clock
Products, Journal of GPS, Vol.3, No.1-2, 95/100 (2004)
2) J.Kouba and P.Heroux: Precise Point Positioning Using
IGS Orbit and Clock Products, GPS Solutions, Vol.5, No.2,
12/28 (2001)
3) O.L.Colombo et al.: Evaluation of Precise, Kinematic
GPS Point Positioning, Proc.of ION GNSS-2004 Meeting,Long Beach, California (2004)
4) B.Hofmann-Wellenhof, H.Lichtenegger, and J.Collins:
GPS-Theory and Practice, 5th revised edition, SpringerVerlag (2001)
5) S.B.Bisnath and R.B.Langley: Precise A Posteriori
Geometric Tracking of Low Earth Orbiters with GPS,
Canadian Aeronautics and Space Journal, Vol.45, No.3,
245/252 (1999)
6) H.Bock,et al.: Efficient Precise Orbit Determination of
LEO Satellites Using GPS, Adv.Space Res., Vol.30, No.2,
295/300 (2002)

Fig.6 Microscopic view of positional errors

Fig.7 Convergence feature before and after reset172

7) S.Bisnath: Precise Orbit Determination of Low Earth
Orbiters with a Single GPS Receiver-Based, Geometric
Strategy, Technical Report No.220, UNB (2004)
8) O.L.Colombo and S.B.Luthcke: Kinematic Point
Positioning of a LEO with Simultaneous Reduced-Dynamic
Orbit Estimation, Proc.of the ION GNSS-2004 Meeting,Long Beach, California (2004)
9) M.Honda, M.Murata, and Y.Mizukura: Development andAssessment of GPS Precise Point Positioning Software forLand Vehicular Navigation, SICE e-Journal, Vol.6, No.10,
78/84 (2007)
10) Y.Shimizu and M.Murata: Flight Evaluation of GPS
Precise Point Positioning Software for Helicopter
Navigation, SICE JCMSI, Vol.1, No.5, 362/367 (2008)
11) A.H.Jazwinski: Stochastic Process and Filtering Theory,
Academic Press, New York (1970)
12) B.D.Tapley, B.E.Schutz, and G.H.Born: Statistical
Orbit Determination, Elsevier Academic Press (2004)
13) E.A.Wan and R.M. van der: The Unscented Kalman
Filter, in Kalman Filtering and Neural Networks (edited by
S.Haykin), John Wiley & Sons (2001)
14) T.Yoshioka: Study on GPS Precise Point Positioning
for Navigation and Orbit Determination of LEO Satellite,
Master Thesis, National Defense Academy (2010) (inJapanese)

Fig. 8 RTN frame

Appendix B: UKF Algorithm
Given discrete-time nonlinear state and observation
equationsYk  G ( X k )   k(B2)

 k ~ N (0, R ) , the UKF algorithm can be written in thefollowing pseudo-code form13).
Initialize withXˆ 0  EX 0 

(B3)



Pˆ0  E ( X 0  Xˆ 0 )( X 0  Xˆ 0 )T



(B4)

For k  1,,  ,(A1)

k 1   Xˆ k 1
Here, Rˆ , Tˆ and N̂ are aligned with the radial, transverse(along-track), and normal (cross-track) directions,
respectively. This set of three triads defines the righthanded rotating coordinate system called the radialtransverse-normal (RTN) frame (Fig.8). The coordinatetransformation matrix from the reference frame to the RTN
frame is given by
T
RTN
C REF
 Rˆ Tˆ Nˆ
(A2)
which is used to convert the position errors in the reference
Cartesian frame into those in the RTN frame.



(B1)

where u k is a known input, wk ~ N (0, Q ) andAppendix A: RTN Frame
Given the position and velocity vectors R and V of thesatellite at an epoch in a reference Cartesian coordinate
system, we can compute the orthogonal unit vectors as
follows:


 
R
R
V
Rˆ   , Nˆ    , Tˆ  Nˆ  Rˆ|R|
| R V |X k 1  F ( X k , uk )  wk

Xˆ k 1   Pˆk 1

Xˆ k 1   Pˆk 1 
(B5)Time update:k / k 1  F k 1 , uk 1 

(B6)

2L

X k  Wi ( m )i , k / k 1



(B7)

i 0

Pk  Wi ( c ) i , k / k 1  X k i , k / k 1  X k   Q
2L

T

i 0

(B8)

Yk / k 1  Gk / k 1 
173

(B9)

Masaaki MURATA (Member)
2L

Yk  Wi ( m )Yi , k / k 1

(B10)

Received the B.E. and M.E. degrees
from Kyoto University in 1967 and1969, respectively, and the Ph.D. degree
in aero-space engineering from the
University of Texas at Austin in 1982.
In 1969, he joined the National
Aerospace Laboratory (now JAXA) andworked on satellite geodesy and GPS
applications. From April 2003 to March
2010, he was Professor of Aerospace
Engineering at the National Defense
Academy (NDA). Currently, he isTechnical Adviser at GNSS
Technologies, Inc. He is a member of
the Geodetic Society of Japan and theAmerican Astronautical Society.

i 0Measurement update:P~y k ~y k  Wi ( c ) Yi , k / k 1  Yk Yi , k / k 1  Yk   R
2L

T

i 0

(B11)

P~xk ~y k  Wi ( c ) i , k / k 1  X k Yi , k / k 1  Yk 
2L

T

i 0

(B12)

K k  P~xk ~y k P~yk 1~y k(B13)Xˆ k  X k  K k Yk  Yk 

(B14)

Pˆk  Pk  K k P~y k ~y k K kT

(B15)

where  L   ,    2 ( L   )  L is a scaling

parameter, L  dimension of the state, Q  process noise
covariance, R  measurement noise covariance, andWi  weights given by

W0( m )   /( L   )W0( c )   /( L   )  (1   2   )
Wi

(m)

 Wi

(c)

(B16)

 1 /2( L   ), i  1,,2 L

The constants  ,  ,   are adjustable scaling parameters
and their values are usually set in the following way:

104    1 ,   0 or 3  L , and   2 (optimal for
Gaussian distributions).
Takahiro YOSHIOKA
Received the B.E. degree fromthe Tohoku University in 2005,
and the M.E. degree from NDA in
2010. He has been with the Japan
Grand Self-Defense Force since
March 2005. His field of interest
is applied estimation and GPS
navigation.

174