計測自動制御学会産業論文集
Vol.6, No.10, 78/84 (2007)

Development and Assessment of GPS Precise Point Positioning Softwarefor Land Vehicular Navigation
Masatoshi HONDA*, Masaaki MURATA* and Yukio MIZUKURA*

GPS precise point positioning (PPP) has attracted precision GPS community as it can provide the accuracy at thecentimeter-decimeter level, which is almost comparable to that obtained by kinematic GPS. In this paper we haveoriginally developed a PPP software for vehicular navigation applications, in that models are pursued so as to beconsistent to those used in the IGS precise ephemeris generation. Assuming pseudorange and carrier phase data from asingle dual-frequency receiver and using IGS precise satellite ephemeris and clock information, the software estimates a
position vector of a moving (and static) antenna, a clock bias of the user receiver, tropospheric delay parameters, andionospheric-free (float) ambiguities in the undifferenced carrier phase data. In particular, all unknowns are assumed to
be random walk processes. Our PPP filter design approach with software has been validated by using a series of GPS
experiment data for static and moving antennas and analyzing data in post-processing mode. The results have shownthat the developed software could yield a 2D accuracy of a few centimeters for static mode and at a decimeter level for
kinematic mode, also reconfirming the reported superior performance of PPP for land vehicular navigation.
Key words: GPS, precise point positioning, IGS, precise orbit and clock, EKF

1. Introduction

We have originally developed a PPP software forvehicular navigation applications as well as static
control survey, which uses both pseudorange andcarrier phase data from a single dual frequency
receiver and IGS products. In this paper we present
the feature of the developed software and its validation
results by processing, in post-mission mode, the GPS
experiment data acquired under static and moving
antenna environments.

Recently, GPS precise point positioning (PPP) has
attracted precision GPS
community as thecentimetr-decimeter accurate positioning results havebeen demonstrated 1),2),3). Such a PPP technique has
become possible with the availability of precise orbits
and clock correction information for the GPS satellites
from IGS (International GNSS service)4). Thisaccuracy level is almost comparable to that obtainablewith carrier phase DGPS, i.e. kinematic GPS3). While
DGPS requires at least two receivers and radio link
between them, and its use is limited generally in therange of about 20 km in baseline length, PPP can work
anywhere with a single GPS receiver only, thus
yielding significant operational flexibility and cost
reduction.
Therefore PPP will be an alternative to theconventional DGPS, and it would find wide areas ofapplication such as remote sensing, monitoring
natural hazards, vehicular navigation, etc. The
application to precision autonomous navigation for
unmanned aerial or land vehicles is especiallyinteresting.

2. PPP
FilterDevelopment

DesignandSoftware

Implicitly assumed is that the receiver is
two-frequency so that ionospheric-free pseudorange
and carrier phase data are available. In that case theobservation equations are written in the form of thelinear combination:PR =f 12 ⋅ PR1 − f 22 ⋅ PR 2
f 12 − f 22

Φ =f 12 ⋅ Φ1 − f 22 ⋅ Φ 2~
N=

f 12 − f 22

= ρ + cdt +d trop+ε PR~= ρ + cdt +d trop +N+ εΦ

cf 12 ⋅ N 1 − cf 22 ⋅ N 2
f 12 − f 22

where PRi is the measured pseudorange on Li ; Φ i is

* Department of Aerospace Engineering, National
Defense Academy, HASHIRIMIZU 1-10-20,
YOKOSUKA,
KANAGAWA-PREF.
239-8686,
JAPAN(Received March 23, 2007)

the measured carrier phase on Li ; ρ is the truegeometric range; c is the speed of light; dt is thereceiver clock bias; dtrop is the tropospheric delay; f i is78the frequency of Li ; Ni is the integer phase ambiguity
on Li ; ε is the measurement noise.
The total troposheric delay is modeled 5) by

the Moon are required, to which two options areavailable in the software, accurate ephemeris by using
JPL ephemeris tape (DE200) and less accurate one byd trop= H dry mdry (e) + H wet mwet (e)+ [Gnorth cos(a) + Geast sin(a)]mgradient (e)Observation Data File

where Hdry, Hwet are the zenith hydrostatic delay andzenith wet delay; Gnorth, Geast are the horizontal delay
gradients in the north and east directions; mdry(e) is thehydrostatic mapping function; mwet(e) is the wet
mapping function; a, e are the azimuth and elevation
angles. The zenith hydrostatic delay can be predictedwith high fidelity by the Saastamoinen model, which
can be described 6) by
H dry =

Database- IGS Precise Ephemeris File and Satellite Clock Correction File
- Reference Coordinate or Reference Trajectory File (Optional)- Earth Orientation Parameters (EOP)- Time Transformation [UTC-TAI]

Developed SoftwareStart

0.0022768 P 0
1 − 0.00266 cos 2 φ − 0.00028 H

where P0 is the barometric pressure in millibars, φ isInitial Values

the latitude and H
is the geodetic height inkillometers. For the mapping functions mdry(e) andmwet(e),we have selected the Niell’s mapping
functions (NMF)7) because of higher accuracy andprecision, and because no surface meteorology isnecessary for its implementation. For the gradientmapping function mgradient(e),mgradient ( e ) =User Options:- Filter- Elevation Cutoff Angle- Mapping Function- Planetary Ephemerides

Read Observation Data (One Epoch)

Assign New Satellite into the Filter State

Time Update1sin( e ) tan( e ) + 0.0032

Edit Data
Process or Reject ?

No

Yes

is implemented 8).
The zenith hydrostatic delay Hdry can be either fixedor estimated, and the other parameters Hwet, Gnorth, andGeast are the solve-for parameters to be estimatedsimultaneously with other unknowns.
Navigation or positioning user can use the precise
orbits and clocks produced by IGS. The IGS’s official
products for precise orbits and clocks are in three
categories: UltraRapid, Rapid, and Final 4). Amongthese products, only the UltraRapid orbit and clock canbe used in real-time. However, in this analysis we haveused the Rapid orbit and clock products. The Rapid
products of precise orbits and clocks are offered at
15-minute and 5-minute intervals, respectively.
Computing satellite position/velocity and clock
corrections at the epochs of data is performed by using
Lagrange’s polynomial interpolation of 9th order fororbits and of 2nd order for clock corrections.
Among phenomena causing site displacements, only
solid earth tide is modeled 9) in the current version, andthe correction for the center-of-mass offset of thetransmitting antenna is applied. For doing these
corrections the planetary ephemerides for the Sun andMeasurement Update

Repeat [1~NSV]Next EpochEnd

Fig.1 Structure and flowchart of the developed
PPP software using a sequential estimation algorithm
analytical orbital theory. Also, relativistic corrections1)
are all taken into account.
The solve-for parameters are position components(3) of the receiver antenna, a receiver clock bias (1),troposphericdelay
parameters(4),
andionospheric-free (float) ambiguities (12) in theundifferenced carrier phase data. The maximum
number of unknown ambiguities is fixed with thenumber of receiver’s channels (12 as usual), and an
ambiguity of a newly coming up satellite can be
automatically incorporated into the filter state. All
unknowns are modeled as random walks with mean
zero and specified variances:
79

The 24-hour data on the day were processed by
the EKF and the position solution was comparedwith the truth. The temporal behaviors of thepositioning errors from the truth are shown in Fig.
2, and the corresponding accuracy statistics istabulated in Table 1. Here, the dE, dN, and dU
denote the positional error components expressed in
the local east-north-up (ENU) frame. From Fig. 2, itis seen that the filter takes about 3 hours or evenlonger for convergence, so that we have excluded
the first 3 hours of errors in calculating statistics inTable 1. After the convergence, the horizontal
components are accurate at a few centimeters; thevertical component is less accurate at a
sub-decimeter level; and the 3-D RMS error is 7 cm.

xi , j +1 = xi , j + wi , jwhere wi , j ~ N (0, Q j (t j +1 − t j )) . The xi , j denotes
the i -th component of the state vector at t j where~~x = (r T b H dry H wet Gnorth Geast N1 L N12 )T .
For parameter estimation, the Extended Kalman
Filter (EKF) was implemented (see Fig.1 for the
structure and flowchart of the developed PPP
software).
In particular, zenith wet delay is a random walkwith a priori mean and standard deviation (STD) of
zero and 10 cm, respectively, driven by white noise of
0.6 cm / h STD. The tropospheric gradients are alsomodeled as random walks driven by white noises of
0.6 cm / h STD. Each of the carrier phase float
ambiguities has a 10 m a priori STD, and is treated asa random walk driven by 0.0006 cm / h white noise.

1.0

(m)

0.5
0.0
-0.5

3. Numerical Results-1.0172800The software validation along with PPP
performance assessment has been conducted, in thatthree datasets, two static (Tsukuba and NDA) and one
kinematic (Kashiwa), were analyzed in post-processing
mode using nominally a full setup of models as
described in the previous section, with elevation cutoff
angle of 10 degrees, unless otherwise noted. In theanalysis the average pressure calculated from themeasured or reported barometric pressure record wereapplied to computing the a priori zenith hydrostaticdelay, so that the estimated zenith wet delay andgradients would reflect the real residual troposphericdelay.

dE
194400dN
216000dU
237600259200

GPS TIME ( sec )Fig. 2 Time history of the position errors (Tsukuba)
0.4

(m)

0.2
0.0

H wet

-0.2
-0.4
172800194400

G northG east216000237600
GPS TIME ( sec )259200Fig. 3 Time history of tropospheric delay parameters
Tsukuba data
The first data to be analyzed are those observed on
22 November 2005 at 30-second interval by the
Geographical Survey Institute (GSI) at Tsukuba,
Ibaragi Pref. This is an IGS site. This dataset wasavailable from the IGS data archives at the CDDIS ofNASA Goddard Space Flight Center, and it covered thewhole one day. No meteorological data have beenreported from GSI. The receiver is an AOA (Allen
Osborne Associates) benchmark ACT, and is driven bya cesium external oscillator. The precise coordinates for
Tsukuba site on that day are also available from IGS
Rapid products, which can be used as the “truth.”

The results for the tropospheric parameter
estimation are plotted in Fig. 3. We observe that thezenith wet delay is significantly time-varying, withmean (STD) of 9 cm (4 cm), and also that the meanvalues (STD) for the gradients Gnorth, and Geast are -0.3
cm (1.2cm), -0.1 cm (1.1cm), respectively. Although nowater vapor radiometer (WVR) data are available with
Tsukuba data, the tropospheric wet delay inferredfrom these estimates should be highly correlated with
WVR data on that day at the site.
The difference from the IGS Rapid very precise clock
corrections for Tsukuba is plotted in Fig. 4. The
80

RMS and the mean (STD) are 0.57 and 0.43 (0.37),
in nanoseconds, indicating the good coincidence
within sub-nanosecond level.

1.0

(m)

0.5

2.0

-0.5

1.0
( nsec )

0.0
dE

-1.0
172800

0.0
-1.0

dN

dU

180000187200
GPS TIME ( sec

194400Fig. 5 Time history of the position errors (NDA)-2.0172800194400216000237600259200

Table 2 Static position accuracy ( NDA )GPS TIME ( sec )Fig. 4 Time history of the clock bias error
Table 1 Static positioning accuracy (Tsukuba)
dE

dN

dU

3D

RMS(m)

0.03

0.02

0.06

0.07

Mean(m)

0.02

-0.01

0.05

STD(m)

0.02

0.01

0.04

dE

dN

dU

3D

RMS(m)

0.03

0.02

0.08

0.09

Mean(m)

0.01

-0.02

-0.06

STD(m)

0.03

0.02

0.05

Kashiwa data
On 7 June 2006, we conducted a GPS experiment fora moving antenna. A van was driven along a rural roadin Kashiwa, Chiba Pref. Two GPS receivers, both
Novatel’s two-frequency OEM4s, participated in theexperiment, one on the vehicle and another at a fixedlocation, to collect the observation data at the rate ofone second. The vehicle’s trajectory and thecorresponding time record of the magnitude of thevehicle velocity vector, which were derived using acommercial kinematic GPS analysis software, areshown in Figs. 6 and 7.

NDA data
On 22 November 2005, the static data were collectedby using a dual- frequency GPS receiver (Novatel’sOEM4) and an antenna with a ground plane for 6hours at one second interval on the roof of the building
at NDA, Yokosuka, Kanagawa Pref. Six hours of NDA
full-rate data were analyzed by the developed software
using the EKF assuming the same nominal setup as in
the previous case. In order to obtain the reference
coordinates or ‘truth’ of the antenna phase center, towhich the PPP solution was compared, this dataset
was submitted to the AUSLIG Online GPS Processing
Service (so-called AUSPOS) operated by the National
Mapping Agency, Australia. The positional errors fromthe truth is plotted in Fig. 5. The position accuracy isgiven in Table 2, where the statistics was calculated
over 5 hours excluding the first one-hour period for thefilter convergence. Again, we observe that: thehorizontal error components, east and north direction,
are accurate at the level of a few centimeters; thevertical error component is less accurate at a
sub-decimeter level; and the 3-D RMS error is 9 cm.
Although not shown here for the sake of economy, themean values (STD) for Hwet, Gnorth and Geast were 1 cm(6cm), -0.03 cm (0.7cm), and -0.3 cm (0.9cm),
respectively.

1500

1000

North ( m )

500
C

A

0

-500

B

-1000
A

B

C

-2000

-1500

-1500
-3500

-3000

-2500

-1000

-500

East ( m )Fig. 6 Trajectory of the vehicle
As guessed from both figures the vehicle repeatedconfigured runs, including frequent go-and-back, andwith the longest straight run of 1.5 km. The total of
4-hour data was collected; in particular one hour atpoint A in static mode; and immediately then the81

vehicle began driving to collect data for 3 hours inkinematic mode on the way of routes B and C. Notethe average vehicle velocity of 5 m/s (for 2 hours) and10 m/s (for 1 hour) over each route. The road is
two-sided, and usually very few cars are running. Theview along the route is mostly open to the sky.
25
B

2.00
1.00
0.00
-1.00
-2.00
266400

C

C

20

(m/s)

dE
5deg

10deg

15deg

(m)

A

origin fixed at the reference site (base station) are
shown in Fig. 9, in which three cases of 5-, 10-, and15-degree cutoff angles are plotted, and the positioning
accuracy for individual cutoff is given in Table 3, where
the statistics was developed over the period of last 2
hours because of the filter convergence.

15
B
10

270000273600277200
GPS TIME (sec)280800

5
A

2.00
1.00
0.00
-1.00
-2.00
2664000270000273600277200280800

GPS TIME (sec)Fig. 7 Time history of the vehicle velocity
The distance between the fixed site and the vehicle
was never more than about 3 km from each other. So,
it was possible to resolve the phase ambiguities in L1
and L2, obtaining the very precise, L1-only, kinematic
solution. This analysis was conducted by using akinematic software called KINGS 10). The resulting
trajectory was used as the ground truth, to which thepost-processing point positioning solutions wereassessed.
The four hours of data were analyzed with
10-degree cutoff elevation angle, except for thedecreasing white noises of 0.06 cm / h STD for the
Hwet and gradients. In addition, another run was
attempted to investigate the effect on the positioning
accuracy of the changing cutoff angles: 5-degrees and
15-degrees. The number of lower-elevation-angle
satellites versus the total number of visible satellites
are shown in Fig. 8.
5〜10 ( deg )10〜15 ( deg )

10deg

15deg

(m)

266400dN
5deg

273600277200
GPS TIME (sec)280800

dU

(m)

2.00
1.00
0.00
-1.00
-2.00
2664002700005deg
270000

10deg

15deg

273600277200
GPS TIME (sec)280800Fig. 9 Time history of the position errors (Kashiwa)
Table 3 Kinematic position accuracy (Kashiwa)
cutoff

dE

dN

dU

3D

5(deg)

0.05

0.14

0.28

0.32

10(deg)

0.19

0.12

0.65

0.69

15(deg)

0.21

0.12

0.92

0.96

5(deg)

-0.01

0.11

-0.17

10(deg)

-0.10

0.11

0.21

15(deg)

-0.07

0.10

0.36

angle
RMS(m)

Mean(m)

STD(m)

Over 15 ( deg )125(deg)

0.05

0.08

0.23

10(deg)

0.16

0.05

0.61

15(deg)

0.20

0.05

0.85

Visible SV

10
8

The filter convergence seems to be much slow, taking
2-3 hours. Lowering cutoff angles have improved thepositioning accuracy overall but not the convergencecharacteristics which have not changed so much. After
convergence the positioning accuracy seems to be at
the almost equal level for three cases of cutoff angles,
being about a decimeter (RMS) horizontally, and a few
decimeters (RMS) vertically. Note that accuracy
figures clearly depend on the convergence definition.
Since the convergence time is critical not only for

6420266400270000273600277200280800

GPS TIME ( sec )Fig. 8 The number of satellites in view
The position errors in the local ENU frame with the82real-time navigation but also for static survey in ashort session, methods for speeding up theconvergence have to be pursued. Among them theproper setting for the magnitude of Q-matrix, use of alow cutoff angle with NMF, and the accurate initial apriori estimate of the state may be significant. The
higher horizontal accuracy with PPP is especiallypromising for vehicular navigation on land. Moreover,the level of accuracy achieved in this study can beconsidered as the worst and does not always imply thepotential performance of the PPP, because of theunordinary trajectory of the vehicle (see Fig. 6), towhich the PPP was tested.

kinematic and longer data.

4. Concluding Remarks
We have developed and evaluated a basic PPP
software workable in both static and kinematic modes,
and reconfirmed the PPP performance using several
GPS datasets. The main conclusions can besummarized as follows.
(1) Our software can achieve the point positioning
accuracy at a few centimeter level in static mode,
and at ten-odd centimeter level in kinematic mode,
horizontally; and at sub-decimeter level in static
mode and at the level of better than two decimeters
in kinematic mode, vertically. These results clearly
exhibit the validity of our design approach of the
PPP filter and the developed software.
(2) The centimeter-decimeter accuracy of the PPP
results can be achieved by estimating troposhericdelay parameters and by possibly allowing lower
elevation cutoff angle, and after the filter reaches thesteady state. However, the filter generally converges
much slowly, taking several tens minutes to a few
hours depending on antenna’s motion andenvironment, and so strategies for speeding up theconvergence will be a future area of research.

0.2( nsec )

0.1
0.0

-0.1
-0.2
266400270000273600277200
GPS TIME ( sec )280800Fig. 10 Time history of the clock bias estimate
0.4

(m)

0.2

Acknowledgements

0.0
-0.2
-0.4
266400H wet
270000

G north

We wish to thank all those who have provided
information or helped during our research. Among
them: Dr. T. Tsujii of JAXA for computing the Kashiwa
double-differencing solutions by his program KINGS,
Dr. H. Sugawara of NDA for supplying themeteorological data, and Dr. H. Namie of NDA for
providing GPS receivers used in the GPS experiments.

G east273600277200
GPS TIME ( sec )280800Fig. 11 Time history of tropospheric delay parameters
For reference the estimates for clock bias andtropospheric delay parameters are shown in Figs. 10 to
11. The estimated Hwet and gradients are, on average,
about 20 cm and a few millimeters, respectively. On
the contrary to the static case, reliable estimation ofthese nuisance parameters in kinematic mode is verymuch difficult, and so estimating only the zenith wet
delay while navigating is of special concern. Several
runs changing estimation strategy and cutoff angle
were tried, and it turned out that the best result interms of positioning quality and convergence property
was with estimating Hwet and gradients and with
5-degree cutoff, following Hwet -only estimation with
10-degree cutoff. This may not be the whole story, but
we need further evaluation using different sets of[ References ]1) J.Kouba and P.Heroux: Precise Point Positioning
Using IGS Orbit and Clock Products, GPS Solutions,
5-2, 12/28 (2001)
2) Y.Gao and K.Chen: Performance Analysis of Precise
Point Positioning Using Real-Time Orbit and Clock
Products, Journal of GPS, 3-1/2, 95/100 (2004)
3) O.L.Colombo et al.: Evaluation of Precise, Kinematic
GPS Point Positioning, Proc. of ION GNSS-2004
Meeting, Long Beach, California (2004)
4) J.Kouba: A Guide to Using International GPS
Service (IGS) Products (Taken from IGS Website)5) Y.E.Bar-Sever et al.: Estimating horizontal
gradients of tropospheric path delay with a single
83GPS receiver, JGR, 103-B3, 5019/5035 (1998)
6) J.L.Davis et al.: Geodesy by radio interferometry:
Effects of atmospheric modeling errors on estimatesof baseline length, Radio Sciences, 20-6, 1593/1607(1985)
7) A.E.Niell: Global mapping functions for theatmosphere delay at radio wavelengths, JGR,
101-82, 3227/3246 (1996)
8) G.Chen and T.A.Herring: Effects of AtmosphericAzimuthal Asymmetry on the Analysis of Space
Geodetic Data, JGR, 102-B9, 20489/20502 (1997)
9) D.D.McCarthy (ed.): IERS Conventions (1996),
IERS Technical Note 21 (1996)
10) T.Tsujii et al.: Development of kinematic GPS
software KINGS and its application to theobservation of crustal deformations in the Izu-island
district, J. Geod. Soc. Jpn,, 43-2, 91/105 (1997) (inJapanese)

American Astronautical Society.
Yukio Mizukura (Non Member)
He joined the National Defense
Academy of Japan (NDA) in 1967.
He received the B.S. degree from the
Kantou Gakuin University in 1970.
Since 1982, he is a Research
Associate in the Department of Aerospace Engineering
at NDA. His research interests include man machine
system and human pilot of aerospace vehicles.

Masatoshi Honda (Non Member)He received the B.S. degree from the
Nihon University in 2001 and the
M.E. from the Department of
Aerospace Engineering at the
National Defense Academy of Japan
in 2007. He joined Maritime Self-Defense Force in
2001. Since April 2007, he has been with theTechnical Research and Development Institute of the
Ministry of Defense. His main research interests
include applying estimation theory to the modelingand analysis of aerospace vehicles.
Masaaki Murata (Member)He received the B.E. and M.E.
degrees from Kyoto University in1967 and 1969, respectively, andPh.D. from The University ofTexas at Austin in the Department
of Aerospace Engineering in 1982. In 1969, he joined
the National Aerospace Laboratory (now JAXA) andworked on satellite geodesy and GPS applications.
Since 2003, he is a Professor in the Aerospace
Engineering Department at the National Defense
Academy of Japan. His research interests are in thearea of applied estimation, satellite geodesy, precision
orbit determination, and aerospace navigation.
Currently he is developing GPS precise point
positioning methods for aerial and space users. He is amember of the Geodetic Society of Japan and the84