Enhancing PPP with GLONASS Integration
Enhancing PPP with GLONASS Integration
Ian Martin
June, 2013
Abstract
Precise Point Positioning (PPP) provides GNSS navigation using a stand-alone re-
ceiver with no base station. As a technique PPP suffers from long convergence times
and quality degradation during periods of poor satellite visibility or geometry. Many
applications require reliable real time centimetre level positioning with worldwide
coverage, and a short initialisation time. To achieve these goals, this thesis considers
the use of GLONASS in conjunction with GPS in kinematic PPP. This increases
the number of satellites visible to the receiver, improving the geometry of the visible
satellite constellation.
To assess the impact of using GLONASS with PPP, it was necessary to build a real
time mode PPP program. pppncl was constructed using a combination of Fortran
and Python to be capable of processing GNSS observations with precise satellite
ephemeris data in the standardised RINEX and SP3 formats respectively. pppncl
was validated in GPS mode using both static sites and kinematic datasets. In GPS
only mode, one sigma accuracy of 6.4 mm and 13 mm in the horizontal and vertical
respectively for 24 h static positioning was seen. Kinematic horizontal and vertical
accuracies of 21 mm and 33 mm were demonstrated.
pppncl was extended to assess the impact of using GLONASS observations in addi-
tion to GPS in static and kinematic PPP. Using ESA and Veripos Apex G2 satel-
lite orbit and clock products, the average time until 10 cm 1D static accuracy was
achieved, over a range of globally distributed sites, was seen to reduce by up to
47%. Kinematic positioning was tested for different modes of transport using real
world datasets. GPS/GLONASS PPP reduced the convergence time to decimetre
accuracy by up to a factor of three. Positioning was seen to be more robust in com-
parison to GPS only PPP, primarily due to cycle slips not being present on both
satellite systems on the occasions when they occurred, and the reduced impact of
undetected outliers.
i
Acknowledgements
This research was funded by the Engineering and Physical Sciences Research Coun-
cil, and Veripos/Subsea 7.
I would like to thank my supervisors Professor Phil Moore and Dr Nigel Penna for
their help and guidance, as well as their feedback reading this thesis.
In addition, thanks are due to Martin Robertson for his invaluable help during the
data collection, and Sam Webb and Matthew Goode for their assistance in fieldwork.
Thanks to Dr Stuart Edwards for his flexibility during the final stages of the thesis
write up, in parallel with another research project.
Thanks also to Pieter Toor of Veripos for providing Veripos orbit and clock products,
offshore GNSS datasets, as well as an industrial perspective. My thanks to Dr Tim
Springer of ESA for the data and information he provided. The data and products
from ESA, CODE, the IGS and the NERC BIGF network has been invaluable in
completing this research.
Thanks to Captain Neale Smiles, Albemarle Barracks for arranging access to the
airfield, and to members of the Lynceus survey vessel for supporting data collection
in the River Tyne.
Finally, thanks to Clare for the endless support and encouragement, without whom
none of this would have been possible.
ii
Contents
1 Introduction 1
1.1 Research background . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Research motivations and objectives . . . . . . . . . . . . . . . . . . . 2
1.3 Research methodology . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Thesis layout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
iii
CONTENTS
4 PPP method 49
4.1 Extended Kalman filter . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.2 Kinematic PPP measurement models . . . . . . . . . . . . . . . . . . 54
4.3 The state vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.4 Calculating the expected observations . . . . . . . . . . . . . . . . . . 57
4.5 Design matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.6 Observation stochastic modelling . . . . . . . . . . . . . . . . . . . . 59
4.7 Parameter stochastic modelling . . . . . . . . . . . . . . . . . . . . . 61
4.8 Quality control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.8.1 Outlier detection / blunder detection . . . . . . . . . . . . . . 64
4.8.2 Cycle slips . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.9 Feasibility of PPP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.9.1 Dilution of precision and the effect of observation error budget
on estimated position . . . . . . . . . . . . . . . . . . . . . . . 69
4.9.2 Expected positional accuracy . . . . . . . . . . . . . . . . . . 71
4.10 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5 Software development 74
5.1 Creating a PPP navigation program . . . . . . . . . . . . . . . . . . . 74
5.2 Adding GLONASS processing . . . . . . . . . . . . . . . . . . . . . . 75
5.3 Source code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.4 Input formats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.5 Configuration file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.6 Output files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.7 Tools for setup and running . . . . . . . . . . . . . . . . . . . . . . . 81
5.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
6 PPP validation 83
6.1 Static positioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
6.1.1 Estimated tropospheric delay . . . . . . . . . . . . . . . . . . 91
6.1.2 Coordinate convergence . . . . . . . . . . . . . . . . . . . . . 93
6.1.3 Year time series . . . . . . . . . . . . . . . . . . . . . . . . . . 98
6.2 Kinematic positioning . . . . . . . . . . . . . . . . . . . . . . . . . . 102
6.2.1 Common processing elements . . . . . . . . . . . . . . . . . . 105
iv
CONTENTS
8 Conclusion 163
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
8.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
8.2.1 pppncl software . . . . . . . . . . . . . . . . . . . . . . . . . . 163
8.2.2 Static GPS PPP validation . . . . . . . . . . . . . . . . . . . 164
8.2.3 Kinematic GPS PPP validation . . . . . . . . . . . . . . . . . 164
8.2.4 GLONASS carrier phase noise . . . . . . . . . . . . . . . . . . 165
8.2.5 GPS/GLONASS static PPP repeatability . . . . . . . . . . . 165
8.2.6 GPS/GLONASS kinematic PPP . . . . . . . . . . . . . . . . . 165
8.2.7 Software architecture . . . . . . . . . . . . . . . . . . . . . . . 167
8.3 Suggestions for further work . . . . . . . . . . . . . . . . . . . . . . . 168
References 185
v
List of Figures
vi
LIST OF FIGURES
7.1 Ionosphere free carrier phase residuals using ESA products . . . . . . 117
7.2 Comparison of GPS and GLONASS LC carrier phase residual RMS
using 2 degree elevation bins with ESA final navigation products. . . 118
7.3 Distribution of the 87 sites used in the static GPS/GLONASS PPP
repeatability test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
7.4 Mean RMS of static PPP estimated North coordinates with respect
to IGS reference coordinates. . . . . . . . . . . . . . . . . . . . . . . . 121
7.5 Mean RMS of static PPP estimated East coordinates with respect to
IGS reference coordinates. . . . . . . . . . . . . . . . . . . . . . . . . 122
7.6 Mean RMS of static PPP estimated Height coordinates with respect
to IGS reference coordinates. . . . . . . . . . . . . . . . . . . . . . . . 122
7.7 RMS of the static PPP estimated North coordinates with respect to
IGS reference coordinates. . . . . . . . . . . . . . . . . . . . . . . . . 125
7.8 RMS of the static PPP estimated East coordinates with respect to
IGS reference coordinates. . . . . . . . . . . . . . . . . . . . . . . . . 125
7.9 RMS of the static PPP estimated Height coordinates with respect to
IGS reference coordinates. . . . . . . . . . . . . . . . . . . . . . . . . 126
7.10 Overall model test statistics for GPS and GPS/GLONASS PPP with
code phase step bias . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
7.11 3D RMS coordinate error for GPS and GPS/GLONASS PPP with
code phase step bias . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
7.12 Overall model test statistics for GPS and GPS/GLONASS PPP with
a carrier ramp bias . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
7.13 3D RMS coordinate error for GPS and GPS/GLONASS PPP with a
carrier ramp bias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
7.14 Overall model test statistics for GPS and GPS/GLONASS PPP with
a code phase ramp bias . . . . . . . . . . . . . . . . . . . . . . . . . . 131
7.15 3D RMS coordinate error for GPS and GPS/GLONASS PPP with a
code phase ramp bias . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
vii
LIST OF FIGURES
7.16 Sections of pickup truck route processed successfully using Track . . . 135
7.17 Distance of the pickup truck from the nearest base station . . . . . . 136
7.18 Pickup truck test Track-estimated reference coordinate uncertainty . 136
7.19 Difference between GPS PPP and reference coordinates during pickup
test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
7.20 Difference between GPS/GLONASS PPP and reference coordinates
during pickup test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
7.21 The Melbourne-Wübbena combination of code and carrier phase for
GPS satellite 2 following a cycle slip . . . . . . . . . . . . . . . . . . 139
7.22 Route of the commercial survey vessel between Svalbard and North
Norway . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
7.23 Distance of survey vessel from GPS reference station NYA2 . . . . . . 142
7.24 Formal errors of the Track estimated reference coordinates for the
Norway survey vessel . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
7.25 Difference between the GPS PPP coordinates and the reference co-
ordinates, for the offshore survey vessel. . . . . . . . . . . . . . . . . . 145
7.26 Difference between the GPS/GLONASS PPP coordinates and the
reference coordinates, for the offshore survey vessel. . . . . . . . . . . 145
7.27 Difference between the GPS PPP coordinates and the reference co-
ordinates, for the offshore survey vessel when less than 110 km from
the base station. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
7.28 Difference between the GPS/GLONASS PPP coordinates and the
reference coordinates, for the offshore survey vessel when less than
110 km from the base station. . . . . . . . . . . . . . . . . . . . . . . 146
7.29 Difference between GPS and GPS/GLONASS PPP height solutions . 148
7.30 GPS/GLONASS PPP and Veripos Ultra real-time Height solutions
compared . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
7.31 Ground track of the survey vessel Lynceus . . . . . . . . . . . . . . . 150
7.32 Distance of the survey vessel Lynceus from the base station NSLG . . 151
7.33 Tyne survey vessel Track reference coordinate formal errors . . . . . . 152
7.34 Difference between GPS PPP and reference coordinates for the Tyne
survey vessel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
7.35 Difference between GPS/GLONASS PPP and reference coordinates
for the Tyne survey vessel . . . . . . . . . . . . . . . . . . . . . . . . 154
7.36 Large scale difference between GPS PPP and reference coordinates
for the Tyne survey vessel . . . . . . . . . . . . . . . . . . . . . . . . 155
viii
LIST OF FIGURES
ix
List of Tables
6.1 Static PPP error compared to IGS weekly solution for 282 stations
on 10 October 2010 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6.2 The daily bias and standard deviation for pppncl estimated ZTD
relative to the IGS ZTD product . . . . . . . . . . . . . . . . . . . . 93
6.3 Mean, minimum and maximum magnitude of position error . . . . . . 96
6.4 Average 2D and 3D pppncl position error of IGS station UNBJ for 2008100
6.5 Mean horizontal and 3D position error of GAPS and pppncl compared
to IGS weekly solutions for 2008 . . . . . . . . . . . . . . . . . . . . . 102
6.6 Airfield reference solution formal errors from Track . . . . . . . . . . 108
6.7 Mean, standard deviation and RMS of difference between the pppncl
and reference coordinates at Albemarle Airfield . . . . . . . . . . . . 111
6.8 Time taken until kinematic pppncl coordinates are within a given
difference from the reference coordinates, Albemarle Airfield . . . . . 111
x
LIST OF TABLES
B.1 List of sites and their coordinates used in static PPP validation . . . 232
xi
Chapter 1 Introduction
Current Global Navigation Satellite Systems (GNSS) methods used to estimate po-
sition provide varying levels of accuracy and coverage. There is a demand for reliable
real-time centimetre level positioning with worldwide coverage. Relative positioning
methods achieve high accuracy with quick start up times using a stream of data
transmitted from a nearby base station, or generated from a regional network of
GNSS base stations. However, these methods are limited in operational area due to
a decrease in accuracy with increasing range from the base station or regional net-
work. In remote locations without accessible base stations, such as offshore, Precise
Point Positioning (PPP) provides an alternative (Russell-Cargill, 2010). In place
of the data from a reference base station or regional network, PPP uses satellite
position and clock data calculated using a global network of GNSS receivers that is
considerably more accurate than ephemeris data broadcast by the satellites them-
selves. The remaining error sources are then mitigated where possible by linear
combinations of observables, physical and empirical models and estimated parame-
ters (Zumberge et al., 1997). Although able to provide global coverage, PPP suffers
from a convergence period where the attainable positional accuracy progressively
improves.
Bisnath and Gao (2008) suggested that integrating the Global Positioning Sys-
tem (GPS) with other navigation systems such as the Russian Globalnaya Navigat-
sionnaya Sputnikovaya Sistema (GLONASS), European Galileo or Chinese BeiDou
could provide more observations and improve the accuracy of positioning and con-
vergence time of the PPP method. This theory was based on having an increased
number of visible satellites, and improving the geometry for determining the posi-
tion. Of these alternative navigation systems, currently only GLONASS contains
enough satellites for functional operation, with 23 out of the planned 24 satellites
already in operation. As of late 2011, both Galileo and the global segment of BeiDou
are at initial testing stage.
1
Chapter 1. Introduction
GPS/GLONASS was by Cai and Gao (2007); this research showed no significant im-
pact on PPP results with the addition of two to three GLONASS satellites that were
visible at any one time in 2007. Since then the number of operational GLONASS
satellites has increased from 13 to 23. The importance of high rate clock correc-
tions was shown by Hesselbarth and Wanninger (2008) in 1 Hz PPP positioning.
They showed a 1.5 to 2.5 times reduction in time to decimetre accuracy of GP-
S/GLONASS PPP compared to GPS only PPP when using GPS and GLONASS
clock corrections tabulated at 30 s intervals. No such corrections were available, so
they used phase interpolation of the low frequency clocks that were available at the
time. Recently high rate GLONASS clock products have become available as part of
an integrated GNSS navigation product, this allows for combined GPS/GLONASS
kinematic PPP processing.
However, the increased number of GLONASS satellites now in operation, and the
new orbit and clock products available, provide the opportunity for further investi-
gation into the improvement of PPP performance, particularly for kinematic posi-
tioning, using GLONASS satellite observations.
In both scientific research and industrial activity in remote locations, there is a need
for determining accurate positions within a short period of time. Positioning perfor-
mance can be quantified by the following metrics: accuracy, precision, initialisation
time, reliability (in terms of correctly reported accuracy), availability and continuity
of positioning solution (Bisnath and Gao, 2008). These metrics are discussed further
in Chapter 2.
2
Chapter 1. Introduction
Whilst PPP solves the availability problem for high accuracy GNSS positioning in
remote locations, the initialisation time and reconvergence time following loss of
signal tracking are major drawbacks with this technique.
Additional satellites offer the prospect of improved geometry and therefore the po-
tential to reduce both the initialisation time and the reconvergence time. There is
also the possibility of improved accuracy dependent on the noise present in the addi-
tional observations. Availability can be improved in locations with poor sky view due
to the additional satellites present, and also frequency diversity of multiple GNSS
operating on different frequencies, providing resilience against interference in a given
frequency band. The addition of extra satellites by using GLONASS also brings in-
creased operational diversity through the use of two independent systems. There
is also the potential for reliability to be improved by the increase in observations
leading to an increased ability to identify erroneous observations.
The recent release of the European Space Agency (ESA) combined GPS/GLONASS
ephemeris products (Springer, 2010) is the first publicly available opportunity to
assess GPS/GLONASS kinematic PPP. Additionally the almost complete network
of GLONASS satellites (23 out of 24 planned) provides increased additional satellite
numbers compared to previous studies.
The objective of this thesis is to quantify the benefits of using the GLONASS satellite
network in addition to the GPS satellite network. Specifically, this research considers
the improvement brought by these additional satellites with respect to the following
metrics: positional accuracy, time to a given positional accuracy and reconvergence
time following an outage or cycle slip. The choice of these metrics is discussed in
Chapter 2.
Throughout this work the approach has been to complete the above comparison
3
Chapter 1. Introduction
for a range of real world datasets from different transport modes in order to draw
representative conclusions.
The initial phase involved writing a program which estimates the position of the
moving receiver in real-time mode using the PPP method. Although there are several
existing programs, the source code was not available for modification. Therefore a
new program, pppncl, was written to be used as a platform for the research, using
routines from Track, a double difference post processing kinematic GPS program
(Chen, 1998), as a starting point to build a PPP program. The software had to be
capable of:
• Real-time mode processing of GPS datasets giving results in line with or better
than existing published work for static and kinematic modes.
• Producing outputs such that the convergence time and position accuracy could
be compared for different options.
Once developed, pppncl output was compared to published results from other aca-
demic programs as well as the commercially available Veripos Ultra. Both static
and kinematic validations were completed using a range of globally distributed sites
and datasets.
The relative magnitudes of GPS and GLONASS observation noise was assessed. The
impact of GLONASS on accuracy and convergence time was determined using static
and kinematic datasets from several sites. Real world data from global locations was
used to verify that the results are applicable on a global scale. Additionally, datasets
provided by Veripos allowed testing of the software using data from an environment
in which PPP is commercially deployed.
4
Chapter 1. Introduction
Chapter 3 considers the possible error sources the PPP method is subject to, and
describes the relevant mitigation used in each case to enable accurate positioning.
The chapter starts with discussion of satellite ephemeris and clocks, and then follows
tropospheric and ionospheric effects. Next, the relativistic effects experienced by
satellite clocks and the importance of satellite orientation, antenna phase centres and
phase windup are considered. This is followed by an explanation of the deformation
effects of the Earth body tides, ocean tide loading, atmospheric pressure loading,
differential code biases and multipath errors. Finally, the residual range error budget
is calculated.
Chapter 4 goes into more depth on the PPP method, describing the theory and
fundamentals as implemented in pppncl. The chapter begins with an overview of
the Kalman filter, the state vector and the design matrix, moving on to observation
and parameter stochastic modelling and finally, the approach taken to observation
quality control.
Chapter 5 describes the pppncl software created during this investigation. This gives
a high level overview of the software design and in particular the considerations
needed for including GLONASS in a PPP program. A detailed description of the
user interface and input and output formats is given in Appendix A.
Chapter 6 describes the validation of the PPP software, pppncl, created. This
includes static validation across a range of globally distributed sites. Next, validation
of kinematic datasets is described and results compared to published work.
Chapter 7 looks at the effect of incorporating GLONASS observation data into the
PPP method. The GLONASS observation noise is evaluated, and the orbit and clock
products described. Next, the main research objective, understanding the impact
of using GLONASS combined with GPS satellite observation data on positioning
quality, is investigated. The chapter steps through this comparison for a range of
datasets collected from different transport modes.
Finally, Chapter 8 presents a summary of the conclusions drawn from this research
5
Chapter 1. Introduction
6
Chapter 2 Types of kinematic GNSS positioning and
navigation
• Accuracy,
• Integrity,
• Continuity,
• Availability.
Accuracy is the difference between the estimated position and the true position,
specified with a certain confidence.
The integrity of a system is the level of trust that can be placed in the output of
that system. In some applications such as aircraft positioning, integrity includes
the ability of the system to alert the user when the accuracy is worse than defined
operational limits (Ochieng et al., 2003).
Continuity is the probability that accuracy and integrity will be maintained for a
certain period into the future.
With this framework in mind, this chapter examines the principles of GNSS posi-
tioning, and how the error budget defines the performance of a positioning method.
The distinct classes of GNSS positioning techniques used to reduce the error budget
to achieve improved performance are discussed. Additionally the advantages and
disadvantages of each technique are considered, along with further metrics that may
be useful in understanding the limitations and strengths of each system.
The final part of this chapter briefly discusses the current and planned GNSS. In
particular the differences between GPS and GLONASS are reviewed, followed by the
expected effect of combining GPS and GLONASS in PPP. The chapter concludes
with an examination of the particular metrics relevant to characterising the impact
of combined GPS/GLONASS PPP.
7
Chapter 2. Types of kinematic GNSS positioning and navigation
A position can be determined by measuring the range from the rover receiver point to
a set of GNSS satellites in orbit around the Earth, using the signals transmitted from
each satellite. A brief introduction to GNSS positioning is given here, for further
background information to GNSS the reader is referred to Hofmann-Wellenhof et al.
(2008) and Leick (2004).
The satellites contain clocks which are used to control the generation of one or
more coded signals which are then transmitted on one or more radio frequencies. A
receiver that knows the transmitted code and the time of generation can measure
the reception time of that code against its local clock. This allows the receiver to
determine the propagation time of the signal from the satellite to the receiver. The
observed propagation time can be converted to a range measurement by multiplying
by the speed of light in a vacuum.
This measured range is not the geometric range from satellite to receiver, but a
measure containing the receiver satellite range in addition to other biases. For this
reason the measured ranges are known as pseudoranges. If the satellite positions
and clock offsets are known then each pseudorange observation describes the surface
of the sphere upon which the receiver may lie. With three or more observations the
location of the receiver in 3D space can be determined. In practice four or more
observations are needed in order to also account for the receiver clock offset.
The biases in the pseudoranges are often considered as errors, as failing to account
for them leads to a reduction in positioning accuracy. There are several classes
of error including: errors in the values of the satellite position and clock offset;
physical effects that delay or advance the signal compared to a propagation model
based on a straight line path moving at the speed of light in a vacuum; and receiver
measurement errors. There are a further set of errors introduced by the definition
of the coordinate system, both in the accuracy of its definition, and due to Earth
deformation effects that are taken into account in the definition of the reference
frame. Finally, multiple reflected signals interfering at the receiver introduce a
source of noise which is highly dependent on the local environment.
The combination of all these effects gives the error budget for the measurements
made by a receiver. This error budget can be related to the accuracy of the estimated
position as discussed in Section 4.9.1. To improve the accuracy of positioning, the
errors contained within observations must be minimised. There are several methods
8
Chapter 2. Types of kinematic GNSS positioning and navigation
of GNSS positioning that aim to correct, or accept the limitations of ignoring, these
errors. The key concepts of some common GNSS positioning methods are described
in the following sections.
2.2 SPP
Single Point Positioning (SPP) is the simple default form of GNSS positioning. This
form only requires a GNSS receiver, no extra infrastructure is required. Satellite
navigation information (orbits and clocks) is broadcast by the satellites and updated
every 4 hours. The accuracy of this broadcast is low due to limitations of the format
and the widely spaced update time. The quality of the orbits is generally the limiting
factor in positioning so most other error sources are ignored. An ionospheric model
is also broadcast which helps single frequency receivers mitigate approximately 50%
of the range bias due to the ionosphere.
The User Equivalent Range Error (UERE) is the combination of the Signal In Space
(SIS) User Range Error (URE) (an estimate comprising orbit data, satellite clock,
ionospheric and tropospheric delay), errors introduced by the user equipment, as well
as local environment errors. Typical values for the UERE are given in Table 2.1.
The total column is formed from the root sum of squares of the bias and random
quantities. The UERE is an estimate of the standard deviation of the measurement
error.
The receivers measure the code observations giving a pseudorange to the satellite.
With four or more satellites in view, the 3D position and receiver clock offset may be
solved for. This method of GNSS positioning is most commonly used in consumer
applications such as in-car navigation and smart phones.
9
Chapter 2. Types of kinematic GNSS positioning and navigation
Relative positioning makes use of the fact that many of the error sources that bias
the receiver to satellite range measurements are spatially correlated. Using one
or more reference stations at known positions, the errors affecting measurements
at the rover can be corrected for. All relative positioning methods need at least
one simultaneously observing reference station, and for real-time applications, a
data link. This imposes limitations on operating areas, for example offshore use
is not practical as there are no nearby reference stations. Three common forms
of relative positioning are described here: Differential GNSS (DGNSS), Real-Time
Kinematic (RTK) and Network RTK.
2.3.1 DGNSS
The methods described so far use the code phase pseudorange as the observable
in the position estimation calculation. Even high end geodetic receivers can only
measure this with at best decimetre accuracy (Groves, 2008). The code phase is
also susceptible to large errors caused by multipath due to the frequency of the code
chipping rate (Hofmann-Wellenhof et al., 2008). It is also possible to measure the
phase of the carrier wave onto which the code signal is modulated. The carrier phase
10
Chapter 2. Types of kinematic GNSS positioning and navigation
can be measured by receivers with an root mean square (RMS) error of 0.5◦ to 2◦
depending on the signal to noise ratio which corresponds to millimetre accuracy in
range (Braasch and Van Dierendonck, 1999).
The bias b may be split into two parts, b = b0 + N where N is the integer part,
commonly known as the carrier phase ambiguity or simply ambiguity, and b0 is the
fractional part. There are two classes of positioning based on how b is determined.
If b is estimated directly then it is known as ambiguity float, and the integer and
fractional part are estimated as one. If instead the integer and fractional part can
be separated, and the integer part N estimated, the technique is referred to as
ambiguity fixed. The following sections describe the different techniques that make
use of the carrier phase.
2.4.1 RTK
RTK positioning involves a base station transmitting its code and carrier phase
observations to the rover. Receiver specific errors are removed by differencing ob-
servations between two satellites forming a single difference observation. Satellite
specific errors are removed by differencing the previously described single differences
from the rover and the base station; this is known as a double difference observa-
tion. Over short baselines the double difference removes all errors sufficiently that
the correct set of carrier phase integer ambiguities can be found. At longer baselines
the errors seen at the rover and base station, particularly due to the troposphere
and ionosphere, are no longer sufficiently similar for the integer ambiguities to be
resolved. RTK has been shown to work effectively up to baseline lengths of approx-
imately 10 km to 20 km.
11
Chapter 2. Types of kinematic GNSS positioning and navigation
Network RTK (NRTK) extends the concept of RTK to a system with multiple base
stations. This allows the spatial variability in the errors to be accounted for based on
the rover’s location relative to the set of base stations. There are two common forms
of NRTK: Virtual Reference Station (VRS) and Master Auxiliary Concept (MAC).
Both approaches attempt to reduce the amount of information required to be sent
to the rover.
In VRS, the observations are combined in the network processing stage to produce
a new set of observations at a virtual reference station at, or very close to the
receiver position. This virtual set of observations contains the expected errors at
the rover based on the current set of reference station measurements. Standard RTK
is then performed using the VRS as the base station. This approach ensures the
distance between the rover and the reference station will be short enough for integer
ambiguity resolution.
In MAC, the complete set of observations from one master reference station and the
atmospheric, orbit and clock errors estimated by the network for the surrounding
auxiliary stations, are transmitted to the rover (Euler et al., 2002). Thus the spatial
distribution of the errors may be interpolated by the rover relative to the master
station, or the observations at the auxiliary stations may be reconstructed allowing
for a multi-baseline solution.
2.5 PPP
Zumberge et al. (1997) introduced and then Kouba and Héroux (2001) developed the
technique of PPP using undifferenced code and carrier phase observations from dual
frequency receivers. Rather than using differences between receivers and satellites
to reduce the error budget, PPP uses external correction products and models for
the error sources.
The starting point for PPP is orbit and clock products accurate to within a few
centimetres. These can be produced from a global network of receivers; this then
enables positioning of a receiver anywhere in the world. Using ionosphere free linear
12
Chapter 2. Types of kinematic GNSS positioning and navigation
combinations of the code and phase observables, along with models for predominant
error sources, Zumberge et al. (1997) showed it was possible to reduce the error bud-
get sufficiently to enable centimetre accurate static positioning. A full examination
of the error sources and applicable models is given in Chapter 3.
In a review of the status of PPP, Bisnath and Gao (2008) found little in the way
of biases in position estimates, so there is little difference between the precision and
the accuracy. North and East one sigma level position errors were a few centimetres
in static and decimetre level in kinematic mode (Bisnath, 2004; Dixon, 2006; Gao
et al., 2005; Muellerschoen et al., 2001).
When using PPP for position estimation there are some key metrics that can be
used to evaluate its performance, these include accuracy, integrity and continuity,
as introduced at the start of this chapter. There is an additional metric that needs
to be evaluated in the case of PPP, known as convergence time. This is the time
for the estimated position to reach a given accuracy. Various authors have chosen
different metrics for convergence time, so that inter-comparison of results is not
straightforward. Metrics proposed have included time to decimetre-level (Bisnath
and Gao, 2008) or centimetre level accuracy (Martín et al., 2011), and accuracy
after a given time period (Colombo et al., 2004; Zumberge et al., 1997).
Due to the unknown fractional part of the ambiguity at the satellite, the carrier
phase ambiguities are estimated as floating. Initially the position accuracy is en-
tirely dependent on the code phase observations. This results in initialisation times
of approximately 30 min for decimetre level accuracy and significantly longer for cen-
timetre level accuracy, reflecting the time taken for the ambiguities to be sufficiently
determined (Bisnath and Gao, 2008).
Gao and Shen (2002) introduced an alternative observation model that uses com-
bined code and carrier phase observations, in addition to the ionosphere free carrier
phase combination. In this model two float ambiguities are estimated for each satel-
lite. Out of 36 one hour static datasets from the Canadian Active Control Station
PRDS, three of the datasets that did not converge below the metre level with the
standard method described by Héroux et al. (2004), did converge to decimetre ac-
curacy using the alternative observation model.
Colombo et al. (2004) evaluated the performance of kinematic PPP using 30 s clock
products. They noted 3D RMS position errors of 4.6 cm to 9.3 cm with typically a
30 min to 40 min or longer convergence period. They proposed using a “stop and
go” survey method where the receiver is held fixed at known coordinates in order to
13
Chapter 2. Types of kinematic GNSS positioning and navigation
In kinematic positioning there is a need for high rate clock products, ideally at the
positioning estimation rate, so that the short term variability in the satellite clocks
is captured and corrected for (Colombo et al., 2004). Hesselbarth and Wanninger
(2008) examined the stability of the satellite clocks for use with PPP and found that
there was minimal change in the kinematic coordinate error RMS when using 30 s
tabulated clock corrections as opposed to 10 s tabulated corrections.
There are a range of groups, including the International GNSS Service (IGS) analysis
centres and commercial providers, that produce the required precise orbit and clock
products at a range of latencies from near real-time to two weeks delayed. The use
of additional data in processing means the delayed products are more accurate.
Investigations have focussed on reducing the error budget, this includes more accu-
rate orbits and clocks, as well as improved error source models. In order to reduce
the error introduced by the atmosphere, a tropospheric delay parameter is usually
estimated and linear combinations of observables are used to reduced the effect of
the ionosphere (Kouba, 2009a). Using precision atmospheric products could reduce
the number of unknown parameters or remove the need for the linear combina-
tions that increase observation noise, and show the potential to improve positioning
performance (Dodd, 2007; Keshin et al., 2006). Whilst steps are being made in
this direction, for example the European ionosphere forecast (Belehaki et al., 2006),
producing sufficiently accurate range corrections (better than a few centimetres)
remains a challenge (Bisnath and Gao, 2008).
The convergence period present in PPP affects the achievable accuracy at a given
time since the start of navigation. Therefore for a given accuracy, the length of the
convergence period to that accuracy directly affects the availability. Any event which
requires the carrier phase biases to be re-estimated, such as temporary obstruction
of the sky, high receiver dynamics or interference causing loss of lock, introduces
a period of reconvergence. This reconvergence is an interruption in the continuity
and reducing the length of this period would improve the overall performance of the
system.
In kinematic PPP the interruption of a satellite signal is both much more likely than
the static case and is more problematic (Grinter and Roberts, 2011). Therefore there
has been significant work on identifying and correcting cycle slips. For real-time
applications methods involving curve fitting across discontinuities of carrier-phase
observations (Beutler et al., 1984) or of linear combinations of observations (Bisnath,
14
Chapter 2. Types of kinematic GNSS positioning and navigation
Multipath, cycle slips, malfunctioning satellites and erroneous navigation data are
all mechanisms through which the measurement accuracy could exceed the expected
error budget. Including such outliers negatively affects the accuracy, hence detecting
these outliers is required to achieve both accuracy and integrity. Bisnath and Gao
(2008) highlight the current poor state of PPP integrity checking as compared to
relative GNSS methods. There is no external quality reference as in Network RTK
where positions can be computed to several of the network base stations.
Apart from filter position covariance estimates, quantitative quality measures of the
obtained results are limited. The estimation filter covariance estimates tend to be
optimistic by comparison with the true error in positional accuracy. This is largely
due to the time correlated nature of many of the error sources which is ignored
in the filter. Correctly determining the residual observation noise after modelling
allows the generation of more realistic uncertainties for the estimated coordinates.
Therefore the estimation of the residual error budget, discussed in Chapter 3, is
important to achieving integrity.
For real-time PPP the reliability of transport for the correction stream is critical to
the availability of the method. For real-time use the IGS real-time working group
15
Chapter 2. Types of kinematic GNSS positioning and navigation
are distributing open access real-time orbit and clocks products in NTRIP format
over the internet (Weber et al., 2007). Commercial bodies provide real-time clocks
distributed by satellite links to ensure global reach (Rocken et al., 2011).
After initialisation, PPP achieves good availability when used in a continuously un-
obstructed open sky environment. However in an environment where the sky view
is regularly obstructed, for example under trees or bridges, the necessary reinitiali-
sation results in a reconvergence period similar to the initialisation period. Equally
the reduced satellite number and resulting poor geometry in, for example, urban lo-
cations, can significantly increase the convergence time. Despite these shortcomings,
PPP has been used in many areas where centimetre static or decimetre accurate kine-
matic positioning is required, particularly in areas where there is no local network
of reference receivers. Some static geodetic applications include crustal deforma-
tion monitoring (Zumberge et al., 1997) and glacial isostatic rebound measurements
(Thomas et al., 2011); scientific kinematic applications include positioning low earth
orbiters (Bisnath, 2004); commercial applications include precision farming (Dixon,
2006), marine sensor positioning, sea floor mapping, marine construction (Arroyo-
Suarez et al., 2005; Rocken et al., 2011) and airborne mapping (Gao et al., 2005).
In summary, PPP has a wide range of uses, but convergence time and reconvergence
time is still a major issue. Two complementary approaches have been investigated:
fixing the ambiguities to integer values is discussed in the following section, and
integrating GPS with other navigation systems such as GLONASS, covered in Sec-
tion 2.8.
The carrier phase ambiguities that must be estimated during PPP positioning can
be split into three parts. These are an integer part, and two Fractional Cycle Biases
(FCBs), one associated with the transmitting satellite and the other associated with
the receiver, b = bj + bi + N .
If both fractional biases can be determined, then it would be possible to use the same
search techniques as employed in relative carrier phase positioning to determine the
integer part. Fixing the integer part to the correct value is desirable as it removes
the need to estimate the carrier phase ambiguity. The carrier phase can then be
used as an unambiguous measurement in the same way as the code phase, but with
a measurement accuracy at the millimetre rather than metre level.
16
Chapter 2. Types of kinematic GNSS positioning and navigation
Collins (2008), Ge et al. (2008) and Laurichesse et al. (2009) showed that with addi-
tional network processing it is possible to separate out the satellite FCBs which can
then be used by other receivers, not part of the original network, to fix undifferenced
carrier phase ambiguities to integer values.
Ge et al. (2008) decomposed the ambiguity into widelane and narrowlane ambigu-
ities. By taking differences between satellites, they remove the receiver dependant
FCB. Then using average values for the Melbourne-Wübbena linear combination of
code and carrier phase (Melbourne, 1985; Wübbena, 1985) from complete satellite
passes for a set of receivers in a global or regional network, the widelane FCBs can
be determined for each satellite. Widelane FCBs are very stable for days or weeks
at a time (Gabor, 1999; Gao and Wang, 2008) so daily estimated values can be used.
In a similar manner the narrow lane FCBs were determined by averaging the frac-
tional part of estimated narrow lane ambiguities, formed from the combination of
the ionosphere-free carrier phase linear combination, and the previously determined
widelane ambiguities.
As the narrow lane FCBs are not stable with time, producing a set of average
values every 15 minutes was proposed (Ge et al., 2008). Fixing ambiguities in PPP
then follows a similar process to the network fixing step, first, the single difference
widelane biases are fixed, having corrected for the satellite widelane FCBs using
the values previously determined in the network solution. Then the narrow lane
ambiguities and their standard deviations are calculated based on the float values
from the position estimation process and the fixed widelane biases. The narrow lane
ambiguities can now be fixed to integers if the standard deviation is low enough
based on a statistical test that the probability of the narrow lane ambiguity be
equal to its nearest integer (Dong and Bock, 1989). Alternatively, the LAMBDA
decorellation method can be used together with a ratio test between the two best
fit sets of narrow lane ambiguities (Teunissen, 1995).
Laurichesse et al. (2009) proposed an approach that determines the widelane FCBs
in the same way but does not estimate narrow lane FCBs. Instead, during the
network estimation stage, the narrow lane biases are directly fixed to integers and
the narrow lane FCB is absorbed into the estimated satellite clock product. Collins
(2008) demonstrated a method that decouples the carrier phase and code phase
receiver clock estimate by estimating one term for each and then fixing wide and
narrow lane ambiguities as in Laurichesse et al. (2009). Both these methods result
in the narrow lane FCB being absorbed in the satellite clock estimate and hence the
generated clock products have been named Integer Recovery Clocks (IRCs) (Geng
17
Chapter 2. Types of kinematic GNSS positioning and navigation
et al., 2010b).
Fixing the PPP narrow lane ambiguities at a single receiver is then simpler than the
method proposed by Ge et al. (2008), as the narrow lane satellite FCB correction is
already contained in the clock products.
In a comparison of the two approaches Geng et al. (2010b) showed that in static
24 hour positioning both methods performed with equivalent accuracy. The mean
RMS position error (compared to the IGS weekly coordinates) of 350 sites for 2008
reduced in the East, North and Up components from 3.4 mm, 2.2 mm and 6.2 mm,
to 2.0 mm, 2.1 mm and 5.9 mm, respectively.
Geng et al. (2010c) showed that ambiguity fixed PPP can improve the accuracy of
static positioning after one hour, reducing the 3D RMS from 5 cm to 1.6 cm. This is
due to the step change in accuracy at the point of fixing ambiguities. However, the
time taken until ambiguities are fixed can be lengthy; for example Laurichesse et al.
(2009) state this to be around 30 min for static positioning, but around 1.5 h for
kinematic positioning with a static initialisation preferred. Such latencies arise from
the time taken to sufficiently accurately determine the narrow lane biases as float
values to enable them to be reliably fixed to integers, as required by PPP ambiguity
resolution.
Li et al. (2011) showed that interpolating the zero difference atmospheric biases
found at base stations in a regional network during ambiguity fixing could be used
to provide instantaneous PPP ambiguity resolution at a rover within the network.
The average baseline used in their test was 60 km and the validity of this approach
at larger distances (such as offshore on a survey vessel) is unlikely due to the decor-
relation of atmospheric biases with increasing distance.
Geng et al. (2010a) showed the possibility for rapid reconvergence following cycle
slips or loss of lock. First the the ionosphere delay is estimated once the ambigui-
ties have been fixed using linear combinations of the carrier phase. Then, following
loss of lock instead of using the noisy Melbourne-Wübbena combination to fix the
widelane ambiguity, the predicted ionosphere delay is used to correct the widelane
combination. It is then possible to estimate the wide lane ambiguity directly by con-
straining the position using the ionosphere free code phase observation. Using this
technique Geng et al. (2010a) showed an improvement in the number of ambiguity
fixed epochs during a moving vehicle trial from 7.7% to 93.6%.
CNES produce a GPS orbit and clock product using the IRC method suitable for
18
Chapter 2. Types of kinematic GNSS positioning and navigation
Table 2.2 summarises the accuracy and initialisation time of the GNSS kinematic
methods discussed previously. As shown, PPP is the only option for global stand-
alone centimetre positioning. However the method suffers from a long initialisation
time and reinitialisation is required after loss of lock, which is particularly a prob-
lem for real-time users. Whilst ambiguity fixed PPP offers improved accuracy and
the potential for very short reconvergence times, any technique that could improve
the float solution convergence time would not only benefit the user position, but
also decrease the time to ambiguity fixing, by enabling the faster initial accurate
determination of the narrow lane biases.
There are currently four proposed or existing Global Navigation Satellite Systems
(GNSS): the American GPS, Russian GLONASS, European Galileo and the Chinese
BeiDou. These are all independent systems aiming to provide global coverage for
positioning, navigation and time transfer. In this section the differences between
these systems are described and the considerations that must be made in combining
multiple GNSS are discussed.
Galileo has just reached the In-Orbit Validation (IOV) phase with two satellites
launched in mid October 2011. These satellites are designed to validate both the
space and ground based components of Galileo; further satellite launches are planned
to reach Initial Operational Capability (IOC) by mid-decade, and Full Operational
Capability (FOC) by 2020 (ESA, 2011).
19
Chapter 2. Types of kinematic GNSS positioning and navigation
The first satellite of the BeiDou system was launched in April 2007, since then further
launches have consisted of geostationary and high Earth orbit satellites centred over
China. Further satellite launches to medium Earth orbits that will provide global
coverage are planned to be completed by 2020 (Beidou, 2011). Galileo and BeiDou
do not currently provide a global service and will not be discussed in further detail
in this thesis.
The first GPS satellite was launched in 1989 with IOC reached in 1993 when 24
satellites were operational. FOC was achieved in 1995 with 24 satellites once the
initial Block I satellites had all been replaced. GPS satellites are arranged in medium
Earth orbits in one of 6 orbital planes.
GLONASS initially reached FOC in 1996 with 24 satellites, however due to a lack of
funding, the number of available satellites declined to six to eight satellites in 2001
(Hofmann-Wellenhof et al., 2008). Since then the system has been maintained to
reach 23 satellites.
GLONASS GPS
Constellation
Number of orbital planes 3 6
Semi-major axis 25 510 km 26 580 km
Orbital height 19 130 km 20 200 km
Orbital period 11 hr 15.8 min 11 hr 58 min
Inclination 64.8◦ 55◦
Distinguishing satellites FDMA CDMA
Signal Characteristics
Carrier frequencies 1602 + k ∗ 0.5626 MHz 1575.42 MHz
1246 + k ∗ 0.4375 MHz 1227.60 MHz
Code frequencies ST code: 0.511 MHz C/A code: 1.023 MHz
VT code: 5.11 MHz P code: 10.23 MHz
Reference System
Reference frame PZ-90 WGS-84
System time GLONASS time GPS time
The key characteristics of GPS and GLONASS systems are detailed in Table 2.3.
The implications of these characteristics on positioning, and the considerations that
must be made in combining both systems in one positioning solution are discussed
below.
20
Chapter 2. Types of kinematic GNSS positioning and navigation
All GPS satellites use the same set of carrier frequencies and make use of the Code
Division Multiple Access (CDMA) method to identify each satellite. This means
that each satellite generates a different set of code transmissions with which the
receiver can identify each satellite. GLONASS uses Frequency Division Multiple
Access (FDMA) to identify the signal from each satellite. This means that all the
satellites transmit the same pseudo-random noise code signal, but each on a different
frequency. As only half the satellites are ever visible from any point on the Earth,
the satellites at opposite points of an orbital plane share the same frequency (ICD
GLONASS, 2008).
The two frequencies transmitted by each GLONASS satellite are based around two
base frequencies. The nominal values of the L1 and L2 carrier frequencies are defined
as follows:
f 1k = f 10 + k∆f 1 (2.1)
f 2k = f 20 + k∆f 2 (2.2)
where k is the channel number of the corresponding satellite. The constants are
f 10 = 1602 MHz, ∆f 1 = 562.5 kHz, f 20 = 1246 MHz and ∆f 2 = 437.5 kHz.
The channel number of each satellite is provided in the almanac which is broadcast
as part of the GLONASS navigation message. All satellites launched after 2005
use values of k in the range of −7 to +6. The ratio of the two carrier frequencies
f 1/f 2 = 9/7 is constant for all k (Leick, 2004, pg. 87) and is very similar to the
ratio for GPS of 154/120.
In a CDMA system, any frequency dependent bias introduced by the receiver hard-
ware will be the same for a given signal from all satellites. As this bias is constant
across all satellites the effect is inseparable from a bias in the receiver clock. In
FDMA systems, such as GLONASS, this bias would be different for each satellite.
In an analysis of the effect of this bias on the ionosphere free code phase combination
(Equation 4.18), Reussner and Wanninger (2011) conclude that “the effect is mainly
frequency dependent but seems to be receiver individual and cannot be modelled
with a simple linear function”. This effectively reduces the accuracy with which the
GLONASS code phase can be observed with biases of up to 6 m measured.
21
Chapter 2. Types of kinematic GNSS positioning and navigation
similar for both frequencies and for receivers of the same type. The biases appear
stable over time and Wanninger (2011) found no temperature dependence when
comparing bias estimates produced 6 months apart.
These additional biases complicate the use of GLONASS satellites for precise po-
sitioning. If they are ignored then this increases the error budget of the observed
pseudoranges. For float PPP any constant bias can be grouped with the ambiguity
and considered a single term so the carrier phase biases will not affect the positioning
solution.
Due to the wider band of radio spectrum used by the GLONASS FDMA signals,
GLONASS is more resistant against narrow-band interference compared to GPS
(Hofmann-Wellenhof et al., 2008).
GLONASS satellites have an orbital inclination of 64.8◦ , compared to the GPS or-
bital inclination of 55◦ ; this provides better coverage at higher latitudes. Specifically
in Arctic and Antarctic regions GLONASS satellites will reach higher elevation an-
gles.
When specifying the location of a satellite, it must be done with reference to some
underlying coordinate system. Similarly for the satellite clock offset, it must be
specified with respect to a time system. Both GPS and GLONASS satellites broad-
cast satellite orbit and clock information in their own reference systems and time
scales.
The GPS terrestrial reference system is WGS-84 (IS-GPS-200E, 2010). GPS time
was set to match UTC in 1980 but is not corrected to match the rotation of the
Earth by the addition of leap seconds. GPS time therefore is at a constant offset of
−19 s compared with International Atomic Time (TAI). The GLONASS reference
system is known as PZ-90 and the time system is GLONASS time. GLONASS time
is closely related to UTC but has a constant offset of three hours. Apart from this
offset, GLONASS time is kept within 1 millisecond of UTC.
If using GPS and GLONASS satellite coordinates and clock offsets, as transmitted
in broadcast navigation, a transformation must be made so that both sets of data
are in the same time system and referenced to the same time scale. There have been
many attempts to derive transformation parameters between WGS-84 and PZ-90,
in particular the IGEX campaign produced several solutions that were consistent at
the decimetre level (Boucher and Altamimi, 2001). The accuracy with which such a
transformation is able to capture the potentially time varying differences has a direct
22
Chapter 2. Types of kinematic GNSS positioning and navigation
impact on the accuracy of the used satellite coordinates and clock offsets. This must
be taken into consideration when estimating the pseudorange error budget if using
broadcast orbits and clocks.
In September 2007 the PZ-90 reference frame was adjusted to align more closely
with the International Terrestrial Reference Frame 2000 (ITRF 2000) (Revnivykh,
2007). WGS-84 is regularly adjusted to maintain close alignment with the ITRF
and therefore differences in terrestrial coordinates between the two systems are now
consistent to within 2 cm (Hegarty and Chatre, 2008).
The GPS orbital period results in the ground track of a particular GPS satellite
repeating once every sidereal day. This causes a resonance between the satellites’
orbital periods with variations in the Earth’s gravity field. This resonance can create
orbital perturbations of up to 4 m (Ferreira and de Moraes, 2009). To counteract
this, the satellites must be manoeuvred regularly to keep them close to their nominal
23
Chapter 2. Types of kinematic GNSS positioning and navigation
orbit. The ground track of the GLONASS satellites repeats every 8 days, avoiding
the resonance experienced by the GPS satellites.
In summary, compared to GPS, GLONASS has additional biases on the code phase
and carrier phase due to the FDMA method used. There is a receiver inter-system
bias composed of receiver hardware delays and the GPS/GLONASS time system
offset. The GPS and GLONASS terrestrial reference systems are different, but
using precise orbit and clock products for both systems that are produced in the
same reference system eliminates this difference.
Bisnath and Gao (2008) suggested that integrating GPS with other navigation sys-
tems such as GLONASS, or in the future, Galileo or BeiDou, could provide more
observations and improve the accuracy of positioning and convergence time of the
PPP method. This conjecture was based on having an increased number of visible
satellites, and hence improving the geometry for determining the position. Píriz
et al. (2009) undertook a limited test on the accuracy of short-occupation static
PPP, by processing 1 hour of data from 20 IGS stations collected on 14 June 2009
using the magicGNSS software (Píriz et al., 2008), and found that the horizontal
accuracies of static batch least squares 1 hour PPP solutions improved from around
10 cm for GPS-only, to around 5 cm for GPS/GLONASS. However, a similar analysis
by Martín et al. (2011) of 1 hour batch static sessions with the same software, using
one day of data from 1 January 2010 for eight IGS stations, showed no improvement
for GPS/GLONASS over GPS-only.
24
Chapter 2. Types of kinematic GNSS positioning and navigation
Pertinent to the inconclusive outcomes from these previous studies is that dur-
ing 2010, nine new GLONASS satellites started transmitting data (one was with-
drawn) to result in an active satellite constellation of 23, and high-rate (30 second)
GLONASS satellite clock values were made freely available by ESA (Springer, 2010).
These developments provide motivation for further investigation of the impact of
adding GLONASS observations to GPS kinematic PPP, with emphasis on reduc-
tions in kinematic convergence time and solution robustness when supplementing
GPS with GLONASS tracking.
In order to assess the impact of combining GLONASS with GPS on the performance
of PPP it was necessary to select relevant performance metrics. The following
metrics were chosen with reasons for their selection also provided.
The convergence time is the time for the positional accuracy to reach a minimum
required level from the start of positioning. The convergence time was expected
to improve due to the increased number and improved geometry of the combined
satellite constellation.
Sturza and Brown (1990) showed that the combination of GPS and GLONASS im-
25
Chapter 2. Types of kinematic GNSS positioning and navigation
The metrics described here form the framework for the tests used in the analysis of
the performance of combined GPS with GLONASS, presented in Chapter 7.
2.10 Summary
This chapter has outlined the principles of several GNSS methods, and the different
approaches to minimising the error sources in order to achieve the required level
of accuracy, depending on application and environment. This shows the value of
PPP, the only method able to provide global centimetre accurate positioning. An
overview of the research into PPP has been given, as well as examining the recent
developments including ambiguity fixed PPP. Whilst this method offers improved
accuracy, PPP still suffers from long convergence times.
Having identified the importance of minimising the error budget in PPP, a detailed
description of the error sources, relevant models and mitigation strategies is given
in Chapter 3.
26
Chapter 3 Error sources and mitigation in PPP
The details of the position and clock offset of the satellites are central to all GNSS
positioning. For GPS, the satellites transmit Keplerian orbit elements and the clock
offset, rate and drift parameters (IS-GPS-200E, 2010). These are calculated by the
GPS ground segment and uploaded to the satellites periodically. The accuracy and
precision of the data sent from satellites is limited by the satellite storage, downlink
bandwidth and ground control processing strategy. For accurate positioning a better
quality of satellite information is required.
Using a global network of receivers, accurate satellite orbits and clock offsets can be
calculated and distributed to users via the internet or a communications satellite.
This information is available commercially or from the IGS analysis centres which
produce orbit and clock products at a range of latencies and corresponding accuracies
(Dow et al., 2009).
Real-time products are created by predicting the satellite orbits and clock offsets
forwards for a short period. This reduces the accuracy achievable as the prediction
is an extrapolation; this is shown in Table 3.1.
In this thesis, commercially available products were used, the real-time Veripos Apex
GPS-only and GPS/GLONASS (Rocken et al., 2011).
Table 3.1 shows the accuracy and tabulation interval of the various IGS orbit and
clock products, as well as the ESA Final and Veripos Apex products. Both the IGS
Ultra Rapid and Veripos Apex G2 products can be used for real time positioning.
The Ultra Rapid orbits and clocks are produced with a latency of 3 h but include
6 h of predicted data. The Veripos Apex G2 product is delivered to users via a
27
Chapter 3. Error sources and mitigation in PPP
continuous satellite link, rather than using batch processing as is done for the IGS
Ultra Rapid product, they are continuously updated and transmitted to the user
with a 2 s latency.
The remaining residual error in the orbit and clock products is given in the equivalent
range error column in Table 3.1. ESA and Veripos do not publish accuracies for their
products. In comparison between the ESA Final orbit and clock products, and the
IGS Final product, the RMS difference in orbit and clocks was between 1 cm to 2 cm
RMS (Kouba, 2009a).
Satellite positions from the IGS analysis centres are produced at 15 min intervals,
and must therefore be interpolated to determine the satellite position at transmission
time (Schenewerk, 2003). By comparing satellite positions interpolated from 15 min
intervals to those produced at a 5 min interval, Yousif and El-Rabbany (2007) showed
the standard deviation of the error introduced by a 9 point Lagrange interpolation
was below the mm level.
Satellite clocks, however, contain high frequency noise that is not captured by such
a long interpolation interval. Therefore the clock product should ideally be tabu-
lated at the same interval as the required positioning interval. When performing
high-rate positioning (1 Hz or higher) this is not possible due to the computational
requirements of generating such high rate clock estimates in the network solution.
Bock et al. (2009) showed the deterioration in 3D coordinate RMS for 1 Hz static
positioning was only 2% using 5 s clock products, but up to 30% using 30 s clock
products. Therefore, to avoid degradation of potential accuracy, the highest rate
clock products up to the sampling rate of the positioning should be used.
The IGS produces a combined GLONASS orbit product using the solutions provided
by the Center for Orbit Determination in Europe (CODE) and the ESA. However,
there is not an IGS combined GLONASS clock product. GMV produce GLONASS
orbit and clock products but they are only available for use within their online posi-
tioning service, magicGNSS (Píriz et al., 2009). Currently, free of charge GLONASS
clock products can only be obtained from ESA (Springer, 2009). CODE orbits and
clocks are used for kinematic GPS PPP validation due to their high clock tabulation
rate (5 s). ESA orbits provide an integrated GPS/GLONASS product.
28
Provider Latency Clock Clock Orbit Orbit Equivalent range GLONASS
RMS (ns) interval (s) RMS (cm) interval (s) error (cm) clock product
IGS Ultra Rapid predicted 3 900 5 900 90 No
IGS Rapid 17 h to 41 h 0.15 300 2.5 900 5 No
IGS Final 12 d to 18 d 0.075 30 2.5 900 3 No
29
ESA Final 12 d to 18 d - 30 - 900 - Yes
CODE Final 12 d - 5 - 900 - No
Veripos Apex G2 2s - 30 - 30 - Yes
As discussed in Chapter 2, any bias in the receiver clock introduces a common range
error on all received signals. The clock in most receivers is usually a quartz crystal
oscillator due to their low cost. These are affected by changes in temperature and
do not have the long term stability of the atomic clocks used in GNSS satellites,
so often drift compared to the satellite network time scale. Most GPS receivers
maintain their internal clock to within 1 ms of GPS time, either by introducing
receiver clock jumps or steering the internal clock (Hofmann-Wellenhof et al., 2008).
The clock offset can be removed by differencing between satellites or, as is the
standard approach in PPP, estimated as a stochastic parameter along with the
receiver position (Kouba, 2009a).
In a GNSS receiver the internal components responsible for GPS and GLONASS
signal measurement can produce different delays for each network. Therefore a
system time difference parameter is needed to account for this (Cai and Gao, 2007),
and for GPS/GLONASS PPP, two unknown receiver clock parameters are estimated.
3.3 Troposphere
As electromagnetic signals pass through the troposphere they travel slower than they
would in a vacuum. This is due to the refractive index of the neutral atmosphere, and
manifests itself as a range bias compared to the assumption of the signals travelling
at the speed of light in a vacuum. Over the range of frequencies that GNSS satellites
transmit on, this delay is a frequency independent effect. In terms of the refractive
index n the tropospheric path delay is defined as:
Z
δtrop = (n(s) − 1)ds (3.1)
The delay can be split into two parts: the hydrostatic part that follows the laws of
ideal gases, often known as the “dry” delay; and a more variable part that is harder
30
Chapter 3. Error sources and mitigation in PPP
to model. The latter is known as the non-hydrostatic or “wet” delay and is related
to water vapour present in the troposphere (Hopfield, 1969). Due to the increase
in the distance that signals at low elevation angles travel through the atmosphere,
the magnitude of this effect is larger at lower elevation angles. The dry part is
responsible for a zenith delay of approximately 230 cm at sea level. The wet delay
is responsible for up to 40 cm in the zenith direction (Leick, 2004).
It has been shown that the tropospheric delay may be approximated for a given
angle by a zenith delay and a corresponding mapping function such that, for a given
elevation angle E, the tropospheric path delay can be written as:
where Zwet and Zdry are the wet and dry zenith delay, and mwet and mdry are the
wet and dry mapping functions (Hopfield, 1969).
The zenith delay can be modelled as follows using the method of Saastamoinen
(1973) as given by Davis et al. (1985).
0.002277p
Zdry = (3.3)
1 − 0.00266 cos 2φ − 0.0028 × 10−3 h
1255 0.002277e
Zwet = + 0.05 (3.4)
T 1 − 0.00266 cos 2φ − 0.0028 × 10−3 h
given the atmospheric pressure p, partial pressure of water vapour e, both in mil-
libars, temperature T in Kelvin, latitude φ and height h in km.
There are many mapping functions to project the zenith delay to a delay at a given
elevation angle. A commonly used example is the Niell Mapping Function (NMF)
which uses a continued fraction with tabulated seasonal and latitude dependent co-
efficients (Niell, 1996). More recently the Vienna Mapping Functions 1 (VMF1) were
produced by ray tracing through numerical weather models such as the European
31
Chapter 3. Error sources and mitigation in PPP
Centre for Medium-Range Weather Forecasts (ECMWF) global weather model. The
required coefficients of VMF1 are provided for IGS, International VLBI Service (IVS)
and International DORIS Service (IDS) stations, as well as on a global grid (2.5 x 2.0
degrees) with a latency of less than 24 hours. The Global Mapping Function (GMF)
is a spherical harmonic fit to seasonal average VMF1 parameters. It requires only
station coordinates and day of year as inputs, making it suitable for use in real-time
global PPP (Böhm et al., 2006a). Compared to the NMF, the GMF significantly
reduces regional height biases and annual errors. Both the wet and dry mapping
functions take the following form:
a
1+
b
1+
m(E) = 1+c (3.5)
a
sin E +
b
sin E +
sin E + c
where E is the satellite elevation angle. The dry and wet mapping functions are
formed from corresponding sets of parameters for a, b and c designated with indices
h and w. The coefficients ah and aw for any site coordinates and day of year (doy)
can be determined from a mean value a0 and an annual amplitude A of a sinusoidal
function using Equation 3.6.
!
doy − 28
a = a0 + A cos 2π (3.6)
365
where the mean value and seasonal value are determined from the tabulated spatial
spherical harmonic coefficients Anm and Bnm up to degree and order 9 as given in:
n
9 X
X
Pnm (sin φ) [Anm cos(mλ) + Bnm sin(mλ)] (3.7)
n=0 m=0
where Pnm are the Legendre polynomials, φ and λ, the site latitude and longitude.
Parameters b and c in Equation 3.5 were estimated by least squares fit to ray traced
values produced from monthly mean atmosphere profiles for 2001 (Böhm et al.,
2006b):
32
Chapter 3. Error sources and mitigation in PPP
bh = 0.0029 (3.8)
bw = 0.00146 (3.9)
" ! ! #
doy − 28 c11
ch = c0 + cos 2π + Ψ + 1 + c10 (1 − cos φ) (3.10)
365 2
cw = 0.04391 (3.11)
Where Ψ specifies the Northern or Southern Hemisphere, the parameters are given
in Table 3.2.
In this thesis the Saastamoinen wet dry zenith troposphere delay model is used with
input meteorological data generated using the GPT model. This accounts for the
dry delay to sub-millimetre RMS (Mendes and Langley, 1998). Most of the residual
error due to the troposphere is caused by the wet delay. Based on comparisons with
radiosondes and water vapour radiometers, Rocken (2005) suggest the wet zenith
delay can be estimated as a parameter using kinematic GPS to approximately 13 mm
RMS. The residual range error contribution to the error budget is then a combination
both of errors in the estimation of the wet trophospheric delay and in the mapping
function.
A correction to the wet delay is estimated as this is highly variable and hard to
model. The zenith values are mapped to a satellite elevation angle using the GMF.
The effect of the troposphere delays is increased at low elevations so any mismod-
elling in the zenith delay will be increased at low elevations. Therefore, as is done in
33
Chapter 3. Error sources and mitigation in PPP
3.4 Ionosphere
The ionosphere consists of charged particles that affect radio signals in a frequency
dependent way. The Total Electron Content (TEC) is equal to the number of free
electrons in a column of unit area along which the signal travels between satellite
and receiver. The free electrons delay the pseudoranges, known as group delay,
and advance the carrier phases, known as phase advance, by an equal amount.
Coronal mass ejections and extreme ultraviolet solar radiation are the main causes
of the ionisation, however the relationship between TEC and solar flux is irregular,
sometimes showing very poor correlation (Doherty et al., 2000). The spatial and
temporal variability of the TEC rules out the option of eliminating the ionospheric
delay through modelling.
The ionospheric range delay for the code signal transmitted on a carrier frequency
f is given to first order by Leick (2004):
40.30
If,P = TEC (3.13)
f2
leading to the following relationships between code (P ) and phase (Φ) observations
and frequency:
with fL1 the L1 carrier frequency, and fL2 the L2 carrier frequency, as given in
Table 2.3.
As the ionospheric range delay is proportional to the TEC, signals at lower elevation
are affected more than signals at higher elevations, due to the slant angle through
the atmosphere. The several metre-level effect this has on measured pseudoranges
requires elimination to achieve centimetre or decimetre accurate positioning in PPP.
34
Chapter 3. Error sources and mitigation in PPP
formed to eliminate the first order ionospheric effects for both the code:
2 2
fL1 PL1 − fL2 PL2
PLC = 2 2
(3.16)
fL1 − fL2
2 2
fL1 ΦL1 − fL2 ΦL2
ΦLC = 2 2
(3.17)
fL1 − fL2
2 2
The normalisation by fL1 − fL2 is performed so that PLC and ΦLC are also measure-
ments of the satellite to receiver range.
Higher order ionospheric effects are generally ignored in kinematic PPP. Although
the maximum effect on the “ionosphere free” phase combination is of the order
40 mm, the majority of the second order effect is absorbed by the satellite clock
corrections (Petrie et al., 2011). Solar activity rises and falls in cycles, with a
recent low activity period from 2006–2010, as shown in Figure 3.1. The period of
data analysed in this thesis (2006–2011) covers low to medium solar activity; this
corresponds to reduced ionospheric activity and therefore a considerably lower effect
on GNSS signals than the worst case. In the near future the impact of higher order
ionospheric effects on PPP will have a greater impact.
200
Average sunspot number
150
100
50
0
6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
199 199 199 199 200 200 200 200 200 200 200 200 200 200 201 201
Date
Figure 3.1: Solar activity cycle showing a minimum during the pe-
riod 2006-2011. Data from [Link]
[Link]
35
Chapter 3. Error sources and mitigation in PPP
In this thesis the first order ionospheric range error is removed using the ionosphere
free measurement combination described above (Equations 3.16 and 3.17). Higher
order terms are ignored due to their relatively small effect on the range error.
The Sagnac effect occurs due to the time dilation experienced by a clock rotating in
a non inertial frame (Ashby, 2003). The magnitude of the effect on a satellite clock
due to orbital eccentricity, ∆tr , is given by:
rrs · ṙrs
∆tr = (3.18)
c2
where rrs is the receiver to satellite position vector, ṙrs the receiver to satellite velocity
vector, and c the speed of light. This is applicable to the satellite clocks, and as
the effect is corrected for in the generation of orbit and clock products, observations
must be corrected during PPP processing to remain consistent with those products
(Kouba, 2009a).
The Earth’s gravitational field causes a space-time curvature of the satellite signal,
this curvature must be accounted for as the range is modelled as a Euclidean range.
This correction may be represented in the following form (Hofmann-Wellenhof et al.,
2008, p. 123):
2µ ρj + ρi + ρji
δ rel = 2 log j (3.19)
c ρ + ρi − ρji
where µ is the Earth’s gravitational constant. The geocentric distances of satellite j
and receiver i are denoted as ρj and ρi , and ρji is the distance between satellite and
receiver.
is a coordinate system with its origin at the satellite centre of mass that rotates
36
Chapter 3. Error sources and mitigation in PPP
with the satellite. The orientation of GPS and GLONASS satellites is maintained
so that the antenna is pointing at the centre of the Earth; k is the unit vector in
this direction. The solar panels are maintained orientated towards the Sun, so given
the satellite to Sun unit vector ê, then j, defined along the solar panel axis is given
by j = k × ê. The unit vector i, located in the Sun-satellite-Earth plane, completes
the right-handed coordinate system. This orientation is also known as the “nominal
yaw attitude” (Bar-Sever, 1996).
The effective reception point, the phase centre, of an antenna varies depending on
the angle of incidence and frequency of the received signal. This also applies to the
phase centre of transmission. As the phase centre is an electromagnetic rather than
physical property of an antenna, and may be inside or outside the antenna, a mark
is put on the outside of the antenna to allow the location of the phase centre to
be measured relative to this point. This mark is known as the Antenna Reference
Point (ARP). Several methods exist for calibration — relative, robot and chamber
calibrations (Görres et al., 2006). In order to maintain consistency with the orbit
and clock products used, the same set of antenna descriptions should be used in
processing as are used in the generation of the products.
For each of the two GPS carrier frequencies, the description of the position of the
phase centre can be split into two parts: a constant per frequency Phase Centre
Offset (PCO), and an azimuth and elevation dependent Phase Centre Variation
(PCV). The IGS maintains a consistent set of calibrations tabulated at 5◦ intervals
for receivers and 1◦ for satellites; by interpolating these values corrections can be
applied to GNSS observations at any elevation angle.
For the offsets provided by the IGS the following sign conventions are used:
where x̄pc is the mean phase centre position, xARP is the location of the antenna
reference point, and xpco is the phase center offset vector (given in a topocentric
left-handed system: north, east and up component). Ignoring all other propagation
37
Chapter 3. Error sources and mitigation in PPP
effects, the observed distance ρobs is then the sum of the geometric range ρgeom
and the phase centre variation PCV(elevation, azimuth). Here the azimuth counts
clockwise from the North towards the East.
Given an offset xpco in the satellite body frame, the conversion from satellite centre
of mass xcom calculated using precise orbits, to antenna transmission point in the
Earth Centred Earth Fixed (ECEF) can be performed:
Similar to the receiver antenna, the phase centre variation is applied as a correction
to the geometric range to give the observed range
Here the azimuth counts clockwise from the j-axis towards the i-axis when looking in
the direction of the negative k-axis or towards deep space (Rothacher and Schmid,
2010).
In the set of antenna offsets used with the the IGS realisation of the ITRF2005
reference frame (IGS05), there are only published receiver antenna offsets and vari-
ations at the two GPS carrier frequencies, L1 and L2. There are no offsets and
variations given for the GLONASS L1 and L2 frequencies, therefore the calibrated
values for the nearest GPS frequency are used, both in orbit and clock generation
and subsequent PPP coordinate estimation, for GLONASS observations. The dif-
ference between the GPS and GLONASS specific corrections are in the range of a
few millimetres but may reach 10 mm when combined in the ionosphere-free com-
bination as shown in Figure 3.2. However when estimating float ambiguities the
average difference between the GPS and GLONASS corrections will be absorbed by
the ambiguity, reducing the 24 h static positioning error to 1 mm or less (Dach et al.,
2011b).
The increase in robot calibrations has increased the number of antennas for which
GLONASS specific phase centre offsets exist. With the change to the IGS realisation
of the ITRF2008 reference frame (IGS08) by IGS analysis centres on 17 April 2011,
an updated set of antenna calibrations has been introduced that include GLONASS
specific offsets.
38
Chapter 3. Error sources and mitigation in PPP
90
60 7.5
50 5.0
40 2.5
30
0.0
20
−2.5
10
−5.0
0
0 60 120 180 240 300 360
Azimuth (degrees)
Figure 3.2: Difference between the ionosphere free PCV for GLONASS
and GPS from GNSS specific calibrations as a function of azimuth
and elevation for the antenna JAV_RINGANT_G3T NONE using
the robot calibration values from [Link]
general/[Link]
For a fixed receiver the antenna orientation does not change, however, the satellites
rotate about their z axis in order to keep their solar panels directed at the Sun. As
the geometry changes due to the satellite’s orbit, there is also an apparent rotation
of the satellite with respect to the receiver.
39
Chapter 3. Error sources and mitigation in PPP
satellite body coordinate unit vectors (x̂0 , ŷ 0 , ẑ 0 ) and the local receiver unit vectors
(i.e. North, East, Up) denoted by (x̂, ŷ, ẑ):
As the satellite completes a full rotation, this must be recognised and the accumu-
lated rotations corrected for by adding full cycle terms of ±2π to the correction term
(Equation 3.25).
For mobile receivers the effect of the receiver antenna rotation is to change all the
measured carrier phase observables by the same angular amount. For GPS only
positioning using only one linear combination of carrier phase observables, this will
add the same range bias to the used linear combination for all satellites. Such a
bias is indistiguishable from a receiver clock bias, and will therefore add a bias to
the estimated receiver clock offset but will not affect positioning accuracy (Kouba,
2009a).
For multi constellation GNSS positioning where the used linear combinations of ob-
servables do not have the same wavelength, a different magnitude bias will be added
to each observable type, depending on the wavelength. In the case of GLONASS
carrier phase observables, where each satellite is on a different frequency, then the
bias added to each satellite’s carrier phase observable will differ.
A complete receiver antenna rotation would add a one cycle bias on both carrier
phase observables. The range effect of this on the ionosphere free combination used
in this thesis is
f 2 λ1 − f22 λ2
∆Φ = 1 2 (3.28)
f1 − f22
Comparing the effect at GPS frequencies to the effect at the furthest GLONASS
frequency (k = 7 in Equations 2.1 and 2.2) gives a difference of 1.9 mm and less than
0.5 mm over the range of GLONASS frequencies. The difference between the effect
for GPS satellites and the mean effect for an observed set of GLONASS satellites
can be considered as an addition to the GLONASS receiver clock bias, leaving only
the effect between GLONASS satellites biasing the observed range. For a spinning
receiver this difference would accumulate with every rotation but this is not a typical
user scenario. The use of a gyroscope or compass to track receiver rotations could
allow the receiver orientation to be determined and used in Equation 3.25 to correct
this effect. In this thesis, however, whilst satellite phase windup is modelled, the
40
Chapter 3. Error sources and mitigation in PPP
There are two situations when satellite orientation does not agree with the model
described in Section 3.6. If a satellite crosses between the Earth and the Sun, the
definition of its body frame requires it to rotate very rapidly; this is called a “noon
turn”. The satellites have a limited rotational speed, so during this rotation period
the true orientation of the satellite lags the modelled orientation. For Block IIR
GPS satellites a similar effect is seen when the Earth is between the Sun and the
satellite, known as a “midnight turn”.
The other case is when the satellite enters the Earth’s shadow. GPS satellites use
a solar sensor to orientate themselves towards the Sun. As the satellite enters the
Earth’s shadow the light sensor can no longer track the Sun. In the shadow, Block
II/IIA GPS satellites start yawing with a maximal hardware yaw rate of approxi-
mately 0.10 ◦ s−1 to 0.13 ◦ s−1 . The behaviour during the post-shadow recovery pe-
riod, which lasts about 30 min, is undefined and cannot be properly modelled. More
recent GPS satellites from Block IIR maintain the nominal yaw attitude during a
shadow crossing. The complete details of this behaviour are described by Bar-Sever
(1996) and a simplified model suitable for network solutions and PPP is given by
Kouba (2009b).
Dilssner et al. (2011) used the j component of the antenna offset of the GLONASS-M
satellites to determine their behaviour during Earth shadow crossing. They found
that using the nominal yaw attitude model can introduce range errors of up to
±19 cm and ±27 cm during noon and midnight turns respectively.
This special behaviour, during turns and shadow crossing, must either be correctly
modelled or observations from satellites during these periods must be removed. In
the case of early GPS and all GLONASS satellites, for which the behaviour during
Earth shadow crossings is undefined or unknown, there is no option but to remove
the observations.
In this thesis the satellites for which the behaviour during shadow crossings is un-
known or poorly defined are removed from the solution. Therefore this effect does
not contribute to the error budget but does reduce the number of useable satellites
when a satellite is performing a noon or midnight turn.
41
Chapter 3. Error sources and mitigation in PPP
The solid Earth is deformed by changing gravitational forces due to the movement
of the Sun, Moon and other planets. This is a periodic effect, the magnitude of
which depends on latitude. The Earth body tide deformation tidal variation can
be as large as 30 cm in the vertical and 5 cm in the horizontal (Kouba, 2009a).
Site displacements caused by tides of spherical harmonic degree and order (nm) are
characterised by the Love number hnm and the Shida number lnm . The values of
these numbers depend on latitude and tidal frequency (Wahr, 1981). The periodic
tidal effect is dominated by diurnal and semi-diurnal tides so can be largely averaged
out when performing static positioning over an entire day (Kouba and Héroux, 2001).
For 5 mm precision, only the second-degree tides and a height correction term are
necessary (Kouba, 2009a; McCarthy et al., 1989). This gives a site displacement
vector ∆r, in Cartesian coordinates of:
3
( " ! # )
X GMj r4 h i h2 2 h2
∆r = 3
3l2 (R̂j · r̂) R̂j + 3 − l2 (R̂j · r̂) − r̂
j=2 GM Rj 2 2 (3.29)
+ [−0.025m sin φ cos φ sin θg + λ] r̂
where GM , GMj are the gravitational parameters of the Earth, the Moon (j = 2)
and the Sun (j = 3); r, Rj are the geocentric state vectors of the station, the Moon
and the Sun with corresponding unit vectors r̂ and R̂j , respectively; l2 and h2 are
the nominal second degree Love and Shida dimensionless numbers (nominal values
0.6078 and 0.0847); φ, λ are the site latitude and longitude (positive east) and θg is
Greenwich Mean Sidereal Time.
A more complete model which takes into account the variability of the Love and
Shida numbers due to effects such as the Earth’s ellipticity, the Coriolis force due
to Earth rotation, and mantle anelasticity is given by McCarthy and Petit (2003).
This model is recommended when a positional precision of 1 mm is desired (Kouba,
2009a). A standard implementation of this model is given in the Fortran subroutine
42
Chapter 3. Error sources and mitigation in PPP
dehanttideinel.f (McCarthy and Petit, 2003, chap. 7), which is used in this thesis.
Ocean tide loading (OTL) has a similar deformation effect to the solid Earth tides
but is caused by the load of the ocean tides. The deformation caused by ocean
tide loading is almost an order of magnitude less than that caused by Earth body
tides (Kouba, 2009a). These tidal effects are much more localised, largely confined
to coastal regions. By convention there is no permanent part of tide loading. The
vertical site displacements due to OTL may reach values of several centimetres for
coastal sites and reduce with distance from the coast (Urschl et al., 2005). For single
epoch positioning to 5 cm accuracy, or sites more than 1000 km from the coast, OTL
may be largely discounted (Kouba, 2009a). The OTL effect is not modelled in this
thesis, in common with several PPP implementations (Abdel-Salam, 2005; Leandro
et al., 2011)
The weight of the atmosphere causes a load on the Earth’s surface. This load
varies with changes in atmospheric pressure causing both vertical and horizontal
displacements of the Earth’s surface. These displacements on average have an RMS
of 2.6 mm for the vertical component and 0.6 mm for the horizontal component, with
peak to peak variations as large as 20 mm in the vertical component and 3 mm in
the horizontal component (Petrov, 2004).
Petrov (2004) provides a model based on using input pressure field data from a
numerical weather model. In an analysis of the model performance based on Very-
Long-Baseline Interferometry (VLBI), the vertical and horizontal displacements
could be computed with errors less than 15%. Dach et al. (2011a) investigated
the impact of this model on the generation of orbit and clock products from a global
network. They found the effect of Atmospheric Pressure Loading (APL) to be clearly
visible and improvement in station repeatabilities of up to 20% when including the
Petrov (2004) correction at an observation level. Dach et al. (2011a) also found a
change in estimated satellite positions when including the APL correction. This
implies that if APL is ignored in the orbit generation then some of the effect is
absorbed into the estimated satellite orbits.
43
Chapter 3. Error sources and mitigation in PPP
Currently the IGS analysis centres do not correct for APL when generating orbits.
Therefore when using these products in PPP some of the effect of APL is already
included in the orbit product.
Urquhart (2009) showed that including APL corrections did not have a significant
effect on positioning given the current accuracy of PPP. Due to the small effect
and the lack of appropriate pressure data, this effect is ignored in this thesis. As
APL affects the receiver coordinates rather than biasing the satellite pseudoranges,
then ignoring APL will not add to the observation error budget, but will reduce the
resulting position accuracy by a few millimetres as given above.
The different code signals transmitted from a satellite all take different paths through
the satellite on their way to the antenna. This causes a slight bias between different
code signals such as C1 and P1. Such a bias is known as a Differential Code Bias
(DCB).
IGS products are aligned to the ionosphere free combination of P1 and P2, therefore
if using a different set of observations they must be corrected. For example, the Leica
1200 receiver records C1 and P2 observables, so for use with orbit and clock products
aligned to the ionosphere free combination of P1 and P2 (such as IGS products),
the C1 observations must be corrected to a P1 equivalent observation using P1-C1
bias values. The magnitude of the GPS P1-C1 bias can reach up to 0.6 m (Kouba,
2009a).
Gao et al. (2001) investigated the nature of the P1-C1 biases using a single an-
tenna split to multiple receivers. This allowed the estimation of constant satellite
and receiver dependent DCBs and a satellite independent time varying DCB. They
showed the constant satellite dependent P1-C1 biases can be estimated with accu-
racy at the level of a few centimetres. CODE produces a monthly set of GPS DCBs
for P1-C1 that can be used to correct for the bias in receivers observing C1 (Schaer
and Steigenberger, 2006).
The case for GLONASS is more complicated as each signal is on a different frequency
44
Chapter 3. Error sources and mitigation in PPP
and so takes a different route through both the satellite and the receiver. This
makes GLONASS DCBs specific to each satellite-receiver pair so corrections cannot
be distributed in the same way as for GPS. Using GLONASS code observations
requires estimating the DCB for each channel.
3.14 Multipath
The signal received at an antenna may have taken more than one path from the
satellite. Objects such as buildings, trees and even the ground can block or reflect
a signal so that the path travelled is increased; this is particularly true in urban
environments.
For code observations the magnitude of this effect is essentially unlimited. The
effect of multipath on carrier phases can be estimated by considering the interfer-
ence between a direct and indirect signal at the receiver. The two signals may be
represented by
AD = a cos φ, AR = βA cos φ + δφ (3.31)
with a the amplitude and φ the phase of the direct signal. The amplitude of the
reflected signal is affected by the damping factor β and the phase is delayed by the
base shift δφ due to the increased geometrical path. Then AD and AR are then the
instantaneous amplitude of the direct signal and the reflected signal respectively.
δφ
βM = 2 cos (3.33)
2
1
δφM = δφ (3.34)
2
45
Chapter 3. Error sources and mitigation in PPP
There is no simple model for multipath as it is highly dependent on the local en-
vironment. Choke ring or multiple antennas can be used to reduce the effect of
signals from the ground (Ray, 1999). For a static site the observation errors can be
tracked over time and due to the repeating nature of current GNSS satellite orbits,
a map of the multipath seen at the site can be generated and then used to correct
for permanent errors due to fixed objects (Ragheb et al., 2007; Wanninger and May,
2001).
Lau and Cross (2007) created a multipath model based on ray tracing signals in a
model of the local environment. Even for simple environments with one large reflec-
tor they found the accuracy of the model highly dependent on the exact specification
of the reflection source. For a moving receiver in an ever changing environment this
kind of correction is not possible.
The carrier to noise ratio is commonly used to describe the quality of a received
GNSS signal (Rost and Wanninger, 2009). This measure is a ratio of the received
signal power to the noise power. Collins and Langley (1999) compared the effect of
even weighting, elevation dependent weighting, and carrier to noise weighting. They
found both elevation dependent weighting and carrier to noise weighting improved
positioning accuracy compared to even weighting. In general there was little (less
than 1 mm) difference between the elevation weighting and carrier to noise weighting
coordinate estimates. This is largely due to the similarity of the two weighting
functions in a low multipath environment.
A cycle slip is the failure of the GNSS receiver to maintain lock on the carrier phase,
causing the number of elapsed cycles of the carrier wave to be miscounted (Blewitt,
1990). This can arise due to high receiver dynamics, ionospheric scintillation, or
46
Chapter 3. Error sources and mitigation in PPP
weak or obstructed signals. The loss of lock causes an apparent jump in the carrier
phase ambiguity of a whole number of cycles. Cycle slips may occur for just one
satellite or all visible satellites, known as total loss of lock.
If undetected, a cycle slip adds a constant bias to all subsequent carrier phase
observations for that satellite. Therefore cycle slip detection is essential if accurate
positioning is to be achieved. Ideally the integer value of the cycle slip is determined
and subsequent observations are corrected to remove this bias, known as cycle slip
fixing. If cycle slip fixing is not possible then the carrier phase ambiguity must be
re-estimated to include the new bias.
In assessments of the models used to account for the above error sources, the most
common practice in the literature is to report results in terms of quantifying the
effect on the GNSS estimated coordinates. This is due to the difficulty in identifying
the true value of the effect that is being modelled. Analysing the estimated coordi-
nates for biases, RMS scatter and seasonal variation is both simpler to achieve and
directly quantifies the development in terms of improved coordinates for potential
users of the model.
The magnitude of some of the error sources, such as atmospheric delays and multi-
path, are also time and location dependent. The variation can be large, for example
the difference in ionospheric activity during the day and at night. Multipath is highly
location dependent, and also affected by variations in receiver tracking bandwidth
at static or kinematic receivers.
Table 3.3 summarises the average contributions to the range error of the residual
error sources following modelling, elimination and reduction methods as described
above. The total error budget is calculated as the root sum of squares of the in-
dividual contributing factors. The effect of multipath on code phase and carrier
phase observations has been omitted as it is so dependent on the local environment,
however we can note that carrier phase multipath is at the centimetre level, whereas
code phase multipath is at the metre level. Atmospheric effects and multipath both
increase at lower elevation, this is often accounted for by elevation dependent ob-
servation weighting and an elevation cutoff threshold.
47
Chapter 3. Error sources and mitigation in PPP
GPS GLONASS
code phase code phase
Satellite orbit 2.5
Satellite clocks 2.5 2.5 > 2.5 > 2.5
Higher order ionosphere 4
Residual troposphere 5
Tracking noise 60 3 60 3
Multipath - - - -
Ocean tide loading 1.5
GLONASS receiver PCV 0.5
Other 2
Total 60 8.3 60 8.3
3.17 Summary
This chapter has considered the error sources that affect positioning accuracy in
PPP. Suitable corrections, models and mitigation strategies to address these error
sources have been identified.
The several metre range bias introduced by the ionosphere, and the lack of a suitably
accurate model, leads to the selection of the ionosphere free linear combinations of
code and phase observations to remove first order ionosphere effects. Models are
identified to account for solid Earth tides, phase centre variation and offset, phase
windup and relativistic effects. Given it is not possible to model multipath for
a moving receiver, observation weighting is used to account for the likely effect of
multipath at different elevations, and signal reception quality. This is also consistent
with the expected increase in unmodelled tropospheric and ionospheric delays at low
elevations.
Estimated parameters are required for the receiver clock offset (including hardware
bias), one for GPS and one for GLONASS. An estimated residual wet troposphere
zenith delay is required.
Finally, the remaining residual errors not accounted for by these methods, and
their contribution to the error budget have been considered. The models, estimated
parameters and residual range error budget form the basis for the PPP algorithm
laid out in the next Chapter.
48
Chapter 4 PPP method
This chapter presents the justification for the use of Kalman filtering as the esti-
mation technique used in this thesis, as well as describing the theoretical model of
Kalman filtering and how it is implemented in PPP.
Weighted least squares can be used to extend the least squares technique to mea-
surements with differing observational noise. In GNSS navigation, measurements
are at intervals and it is desirable to estimate a new position with each new set
of measurements. In classical least squares each set of observations creates a new
estimated position independent of all previous position estimates.
Multiple sets of measurements can be combined to reduce the effect of the obser-
vational noise. In batch least squares this is achieved by combining all the sets of
observations in the estimation process; each new set of observations adds to the size
of the set of equations to be solved.
The recursive nature of the Kalman filter means that the contribution of past obser-
vations is contained in the covariance estimates of the estimated parameters. This
makes the addition of a new set of observations numerically efficient compared to
the batch least squares method. In batch least squares the addition of a new set of
measurements increases the number of observations, and hence the size of the set of
linear equations that need to be solved.
The PPP observation model discussed is non-linear, therefore the Kalman filter must
49
Chapter 4. PPP method
be adapted. There are two approaches to linearising the Kalman filter, the linearised
Kalman filter and the Extended Kalman Filter (EKF). The linearised Kalman filter
starts from a nominal trajectory about which the Kalman filter equations are lin-
earised using a truncated Taylor series expansion evaluated at the nominal trajectory
(Grewal and Andrews, 2001). The EKF, sometimes known as the Schmidt-Kalman
filter (Schmidt, 1976), performs the linearisation about the estimated trajectory.
The linearised Kalman filter can reduce the real time computational burden as the
linearisation can be performed offline, whereas for the EKF it can only be computed
once the estimated position is known. The disadvantage of the linearised Kalman
filter is that the deviation of the actual trajectory from the nominal trajectory tends
to grow with time leading to increasing linearisation errors (Grewal and Andrews,
2001). The EKF is therefore more robust against nonlinear approximation errors,
as the filter only assumes linearity over the range of state estimation errors. Addi-
tionally, for kinematic positioning the nominal trajectory is often not known.
An alternative to the EKF is the Unscented Kalman Filter (UKF), which uses a de-
terministic sampling approach to capture the effects of the nonlinear system on the
mean and covariance of the state (Wan and Van Der Merwe, 2000). The carefully
chosen sample points completely capture the true mean and covariance of the Gaus-
sian Random Variable, and when propagated through the true nonlinear system,
capture the posterior mean and covariance accurately to the 3rd order (Taylor se-
ries expansion) for any nonlinearity. The EKF, in contrast, only achieves first-order
accuracy.
Other nonlinear filtering algorithms include particle filters and the numerical solu-
tion of the Fokker-Planck Equation (Daum, 2005). Whilst these can lead to im-
proved performance of estimation in nonlinear systems, this is often at the cost of
computational complexity.
Due to the ease of implementation, efficient run time and wide use within the liter-
ature the EKF was used in this thesis. The principles of the EKF together with the
implementation equations are described in detail below.
In simple terms, at a given observation epoch, the expected values of the measure-
ments are estimated based on the predicted state (position, receiver clock offset,
troposphere correction and ambiguities) using the precise orbit and clock products
50
Chapter 4. PPP method
and applying the models outlined in Chapter 3. The actual measurements and
the modelled measurements, together with the state and measurement covariances,
are combined in the Kalman filter to produce an updated state based on the new
information contained in the current set of observations (Figure 4.1).
The state of a system is the value of a set of parameters that are of interest at
a given time. Often these parameters may not be directly measurable but are
linked to variables that are observable through a measurement model. The dynamic
model of a system describes how the state at one time is related to the state at a
subsequent time. The EKF estimates the state of a system based on the predicted
state generated by the dynamic model and a set of measurements.
1) the Prediction Step where the time dependent state vector and its associated
covariances are updated based on the system dynamic model,
2) the Update Step where the predicted state is updated based on the measurements
at that point.
Full derivation is not given here as there are many examples in texts such as Grewal
and Andrews (2001) and Anderson and Moore (1979). Details of the models mak-
ing up the Kalman filter are given; the PPP specific versions of these models are
discussed later.
51
Chapter 4. PPP method
where xk is the state vector at epoch k, fk−1 the state transition function describing
the dynamic model, wk−1 the zero mean process noise which describes unmodelled
changes in x not contained in f that can be described as a normally distributed
random vector with associated covariance matrix Qk−1 .
where hk is the measurement model and vk the zero mean, normally distributed
measurement noise with covariance matrix Rk .
Given the dynamic model (Equation 4.1) and measurement model (Equation 4.2),
the predicted state and modelled set of observations (given that state) are defined
by the following nonlinear implementation equations:
52
Chapter 4. PPP method
The estimated state based on the predicted state and the measured observations is
calculated using the following steps:
the covariance matrix can be propagated forwards in time to give the a priori
covariance matrix:
[1] [1]T
Pk (−) = Fk−1 Pk−1 (+)Fk−1 + Qk−1 (4.6)
• The system dynamic model is linearised about the predicted state to give the
design matrix:
[1] ∂hk (x)
Hk ≈ (4.7)
∂x x=x̂−
k
which weights the influence of the observations on the state based on the
combined statistics of the predicted state and the observations.
• Finally the a posteriori covariance matrix of the estimated state can be com-
puted:
h i
[1]
Pk (+) = I − Kk Hk Pk (−) (4.10)
The derivation of the Kalman filter requires that the random noise vectors v and w
53
Chapter 4. PPP method
where E is the expectation operator. The implication is that any unmodelled biases
in the system noise or observation noise will propagate into the solution. Hence the
requirement for the precise models and corrections detailed in Chapter 3.
[1] T
h i h i
[1]
Pk (+) = I − Kk Hk Pk (−) I − Kk Hk + Kk Rk KkT (4.14)
This is the sum of a symmetric positive definite matrix and a symmetric non-negative
definite matrix, ensuring Pk (+) is symmetric positive definite (Bucy and Joseph,
1968).
This section describes the specifics and mathematical equations of the measurement
model used in this GPS/GLONASS positioning software.
The model for code phase Pi and carrier phase Φi GPS pseudorange observations is:
I
Pi = ρ + dtrop + cdt + + Pi (4.15)
fi2
I
Φi = λi φi = ρ + dtrop + cdt − 2 + λi bi + Φi (4.16)
fi
54
Chapter 4. PPP method
• ρ the satellite to receiver range corrected for earth deformation and antenna
effects,
• Pi and Φi respectively, the code and carrier phase noise including multipath
and receiver measurement error,
Here the satellite and receiver FCB discussed in Section 2.5.1 have been included in
the carrier phase bias bi .
The ionosphere free combination is formed for both code and phase observations to
eliminate the ionospheric delay:
f12 P1 − f22 P2
Pc =
f12 − f22
(4.17)
f 2 P − f22 P2
= ρ + dtrop + cdt + 1 21
f1 − f22
f12 Φ1 − f22 Φ2
Φc =
f12 − f22
f 2 λ1 b1 − f22 λ2 b2 f12 Φ1 − f22 Φ2 (4.18)
= ρ + dtrop + cdt + 1 +
f12 − f22 f12 − f22
= ρ + dtrop + cdt + bc + Φc
The model for GLONASS is similar, with the addition of an extra two terms to ac-
count for the difference between GPS and GLONASS time (dtR ), and the frequency
55
Chapter 4. PPP method
k
dependent differential code biases (Bc,R ) (Cai and Gao, 2007). The time difference
includes any offset between the time the receiver samples the GPS observations,
and the time the receiver samples the GLONASS observations. Only applicable to
the code observations is the differential code bias for each observation frequency.
Equations 4.19 and 4.20 are respectively the ionosphere free code and carrier phase
observation equations for GLONASS. The index k, used as both a subscript and a
superscript denotes the frequency dependence of the denoted terms. The subscript
R used to denote GLONASS specific observations and parameters is chosen to match
the use of the character R in RINEX observation format files to identify GLONASS
satellites.
2 2
k f1,k P1 − f2,k P2
Pc,R = 2 2
f1,k − f2,k
(4.19)
f 2 P1 − f2,k
2
P2
= ρ + dtrop + cdt + cdtR + 1,k 2 2
+ Bkc,R
f1,k − f2,k
2 2
f1,k Φ1 − f2,k Φ2
Φkc,R = 2 2
f1,k − f2,k
2 2 2 2
(4.20)
f1,k λ1,k b1 − f2,k λ2,k b2 f1,k P1 − f2,k P2
= ρ + dtrop + cdt + cdtR + 2 2
+ 2 2
f1,k − f2,k f1,k − f2,k
If using precise orbits referenced to a single time frame, the difference between GPS
time and GLONASS time is contained within the individual satellite clock correc-
tions so requires no further consideration. For example the ESA final GLONASS
clock product is aligned to GPS time (Springer, 2010). This does not however imply
that dtR ≡ 0 as there is no guarantee that the receiver clock offset between GPS and
GLONASS sampling time is zero, either due to sampling timing or different signal
propagation paths within the receiver. What it does mean is that the satellite clock
offsets correct both GPS and GLONASS satellite clocks to the same timescale and
that the ECEF satellite coordinates for both systems are referenced to the same
epoch.
From the measurement model the unknown parameters that are either of direct
interest or cannot be accurately modelled are:
56
Chapter 4. PPP method
where (xr )i are the three Cartesian coordinates of the receiver position within the
reference frame defined by the precise orbit products used, dZwet the wet tropospheric
zenith delay correction, cdt the receiver clock offset with respect to the time frame
of the precise clock products, cdtR the receiver GLONASS sampling time offset with
respect to its GPS sampling time, bi and bi,R respectively the GPS and GLONASS
carrier phase biases for the visible set of satellites. Here bi refers to the ionosphere
free bias bC (introduced in Equation 4.18) for GPS satellite i; similarly for the
GLONASS biases bi,R .
The exact order of the elements in the state vector is not important but for the
following description of design and noise matrices, and in the software produced,
the above order is used.
For each epoch the modelled value of all GNSS observations made at that epoch is
calculated based on the predicted state (Equation 4.3). First the predicted position
is adjusted for all Earth deformation effects, antenna offset and antenna phase centre
offset.
The range ρsr from receiver r to satellite s depends on the location of the satellite at
the time of transmission. The time of transmission is itself dependent on the time
of signal travel i.e. the range. Both precise and broadcast orbits are given in an
ECEF frame so the rotation of the frame between time of transmission ttrans and
time of reception trec must be taken into account. Combined, these effects give the
simultaneous equations:
where R(trec −ttrans ) is the rotation matrix that maps points in ECEF at time ttrans to
points in ECEF at time trec . There are two possible approaches to solving Equations
4.22 and 4.23. Either the pseudorange may be used as ρ in Equation 4.23, or the two
equations may be iterated starting from a nominal value of ρ (e.g. 20 000 km) which
converges rapidly. In this thesis the second method has been used as it enables
modelling of satellite observations without the corresponding observation data.
The satellite antenna offset, phase centre offset and variations are then applied
57
Chapter 4. PPP method
to determine the effective transmission point given the calculated centre of mass.
The additional corrections due to satellite and receiver clock offsets, phase windup,
tropospheric delay, and the relativistic correction are added to the range having
been calculated as described in Chapter 3. The ionosphere free combination is then
formed for both the code and carrier phase observations.
The design matrix is the Jacobian of the measurement model (Equation 4.7). The
full matrix form is shown:
∂h ∂h ∂h ∂h ∂h ∂h ∂h ∂h
... ...
∂(xr )1 ∂(xr )2 ∂(xr )3 ∂dZwet ∂dt ∂dtR ∂b1 ∂b1,R
(4.24)
.. .. ..
. . .
with each row being a first order differentiation of the observation equation for each
measurement, h, with respect to each element in the state vector. For the three
position elements:
1
2
∂h ∂ X s (xr − xs )i
= ((x )j − (xr )j )2 = (4.25)
∂(xr )i ∂(xr )i j ρ
where xs is the satellite transmission point, xr the receiver reception point and ( )i
denotes the ith component of a vector.
For the estimated correction to the ZTD, differentiating the wet tropospheric delay
from Equation 3.12 with respect to the estimated correction gives:
∂h
= mwet (4.26)
∂dZtrop
To avoid large differences in the magnitude of the values of H [1] which can lead to
problems due to numerical round off (Grewal and Andrews, 2001), the receiver clock
corrections and carrier phase biases are chosen to be estimated in metres, i.e. the
58
Chapter 4. PPP method
state vector contains cdt, cdtR , λc bi and λc,k(i) bi,R . This gives partial derivatives of:
∂h
=1 (4.27)
∂dct
∂h 0 GPS satellites
= (4.28)
∂dctR 1 GLONASS satellites
∂h 1
GPS satellite i carrier phase observations
= (4.29)
∂λc bi 0 code and other carrier phase observations
∂h 1
GLONASS satellite i carrier phase observations
= (4.30)
∂λc,k(i) bi,R 0
code and other carrier phase observations
For each set of observations, the observation covariance matrix R in Equation 4.9
must be constructed. The ionosphere free observations are modelled as independent
so off diagonal elements of R are zero. The code phase and carrier phase observations
are assumed to be uncorrelated.
The described set of GNSS observations are heterogeneous in nature, mixed code
and carrier phase, two GNSS systems, and have quite different measurement noise.
To optimally combine a mixed set of observations the stochastic properties of each
must be known so each observation is given the correct weighting in the filter.
Due to the wavelength of the pseudo-random noise signal that forms the code ob-
servation, there is a limit to how accurately a receiver is able to measure the code
pseudorange. A static receiver can use a narrow correlator when tracking the code
whereas, due to changing receiver dynamics, a kinematic receiver may have to use a
wide correlator reducing the accuracy of the measurements made. Due to the consid-
erably shorter wavelength of carrier phase observables and the frequency dependency
of multipath (Lachapelle, 1991), the carrier phase noise is considerably lower than
the code noise. Typical noise values for GPS are summarised in Table 4.1.
Elevation dependent weighting of the observations can be used to account for the
increased noise due to the higher level of multipath, tropospheric and ionospheric
delays experienced at lower elevations. Collins and Langley (1999) showed the most
appropriate model was based on the signal to noise ratio. Modelling the observa-
tions as independent gives the following equation for the diagonal elements of the
59
Chapter 4. PPP method
Measurement Noise
C code 10 cm to 300 cm
P code 10 cm to 30 cm
Carrier Phase 0.2 mm to 5 mm
SNR
Rii = Ci × 10− 10 (4.31)
where SNR is the signal to noise ratio and Ci a constant factor specific to the
receiver carrier loop tracking bandwidth. However, much archived RINEX data
does not include detailed signal to noise measurements; in this case the following
elevation dependant weighting is used
The off diagonal elements of R in the standard PPP model are zero.
Lau and Cross (2006) noted the orthogonal nature of the SNR and the carrier
phase multipath. Using linear combinations of the observations to determine the
pseudorange multipath, each set of observations from a satellite can be classified as
containing a multipath error or not. By modifying the SNR based stochastic model
when a multipath error is detected, up to 10% improvement in 3D RMS over the
unmodified SNR model was found in high multipath environments.
Hesselbarth and Wanninger (2008) phase interpolated ESA GNSS clock products
tabulated at 5 min intervals to a 1 Hz frequency and measured the GLONASS carrier
phase residuals during PPP. They found the noise of GLONASS carrier phase obser-
vations to be comparable to GPS observations. In light of the high-rate GLONASS
clock products now available from ESA, the relative noise of GPS and GLONASS
carrier phase is compared in Chapter 7.
f12 f22
σΦ2 c = ( )2 2
σ + ( )2 σ 2 (4.33)
f12 − f22 Φ1 f12 − f22 Φ2
60
Chapter 4. PPP method
where σΦ2 1 and σΦ2 2 are the carrier phase variances on L1 and L2 respectively.
In the software produced the elevation dependent model in Equation 4.32 is used.
The stochastic model of the estimated parameters is needed for Kalman filtering.
There are two things to consider: first the initial uncertainty which describes the
accuracy with which the initial parameters are known, and second the epoch to
epoch stochastics.
The initial covariance P0 that describes the initial state vector x̂0 is set depending
on the source of the initial state. Initially one may assume no correlation between
elements of the state vector, requiring only diagonal elements of P0 that represent
the variance of each element of x̂0 . If certain components of the initial state are
known, for example starting the processing from known coordinates then the initial
covariance may be provided. Another alternative is to use the method of Bancroft
(1985) to solve directly for an initial position with an accuracy of at least 20 m. An
initial value for the carrier phase biases may be found by differencing the carrier
phase observations with the pseudorange observations:
λc bi = Φc − Pc (4.34)
with initial variance given by the pseudorange observation noise stochastic model.
In the absence of a known dynamic model for the receiver coordinates, the unknown
parameters of the state vector may be modelled as random walk processes with
sufficient process noise to capture the epoch to epoch dynamics (Kouba and Héroux,
2001; Zumberge et al., 1997). The carrier phase bias parameters are, in the absence
of cycle slips, assumed to be constant over time.
where g(t) is the system rate dynamic model and w(t) the random white noise
defined by the spectral density matrix q(t). The values of the process noise matrix
Qk can be found, to a first order approximation, by the propagation of the system
61
Chapter 4. PPP method
A random walk process is well suited to receiver coordinate, troposphere delay and
clock offset estimation as the rate of change is independent of the current value. For
a random walk process, Equation 4.35 becomes:
[1]
Fk = I (4.38)
With the use of a random walk stochastic model Equation 4.36 simplifies to:
Z ∆tk
Qk = q(t)dt (4.39)
0
The values of the vector q(t) are highly dependent on the characteristics and dy-
namics of a particular receiver. Whilst the correct values may not be known the
spectral density may be chosen to reflect the expected rate of change of a parameter.
For example, at a static site the position is constant so qxi is chosen to be zero. In
a survey vessel operating at speeds of up to 10 m s−1 the spectral density can be
chosen as qxi = 102 m2 s−2 .
Tuning the spectral density to a particular application can improve the accuracy of
the resulting estimated positions. Over constraining by specifying too low a spectral
density results in a reduction in sensitivity to high rate changes in state whereas
under constraining leads to excessive noise in the estimated state.
Table 4.2 shows the stochastic model used for the estimated parameters in the EKF.
Obtaining a consistent stochastic model is critical to the estimation process but is
not straight forward (El-Mowafy, 2011). If there is too much process noise, the
estimated parameters will be less accurate due to the influence of observation noise.
If there is too little process noise for the estimated parameters, then the filter will
be very slow to converge and may not respond to the true dynamics of the system
62
Component Stochastic model Process noise standard deviation
Position random-walk process static 0, kinematic speed specific
Receiver clock error GPS random-walk process 100 m s−0.5
GPS-GLONASS system time offset random-walk 100 m s−0.5
63
Troposphere wet delay correction random-walk 10−3 m s−0.5
Carrier phase ambiguities constant 0
The correct process noise depends on the receiver dynamics, the variability of the
receiver clock and the variation in the troposphere delay at a given time. Therefore
standard values will not be correct all of the time, but the aim is that they are
relevant values for the majority of the time considered.
Due to the recursive nature of the Kalman filter, any error in the measurements used
not only corrupts the solution for that epoch but also future epochs. It is therefore
very important to identify any such errors and correct for or remove them. The
methods for dealing with uncorrectable outliers and potentially correctable cycle-
slips are described in the next section.
Due to multipath (see Section 3.14), undetected cycle slips and modelling errors,
some observations will contain erroneous data. To check if a set of observations
is consistent with the stochastic model associated with them, the Local Overall
Model test may be used (El-Mowafy, 2010; Teunissen, 1990). The prefit residuals vk
are defined as the difference between the modelled observations and the measured
observations:
vk = zk − Hk xk (−) (4.40)
where zk are the measured observations, Hk the design matrix and xk (−) the pre-
dicted state.
As the Kalman filter is based on the prefit residuals being normally distributed, a
Chi Squared test can be used to test if the prefit residuals are normally distributed.
describes the expected variance and correlation between residuals based on the com-
64
Chapter 4. PPP method
bination of each observation uncertainty and the uncertainty in the predicted state.
If the global test statistic detects an outlier then a local test is performed to identify
which observation is most likely to be the outlier. Each prefit residual (vk )i is
normalised:
(vk )i
wi = q (4.44)
(Cvk )ii
The index of the most likely outlier is given by wmax the largest absolute normalised
residual:
wmax = max |wi | (4.45)
The whole process is then repeated on the reduced set of observations until the global
model test passes. If no observation fails the local test then the set of observations
is inconsistent with the current state and the stochastic model. This could be due to
the measurement variances used in the filter being too small or the state being over
constrained through not adding enough process noise to describe the state dynamics.
Adding additional process noise to the state covariance, Pk = Pk−1 + cI for some
large c, effectively resets the filter, allowing it to re-converge.
If using a positioning algorithm such as the Kalman filter, for which each carrier
phase bias term is estimated over more than one epoch, cycle slips need to be
identified and ideally corrected for (see Section 3.15). There are several methods to
detect cycle slips:
65
Chapter 4. PPP method
• differencing the expected cycle count found by integrating the Doppler mea-
surement from the observation,
For a dual frequency receiver, a cycle slip is an integer change in one or both of the
carrier phase biases:
(∆n1 , ∆n2 ) = (b01 − b1 , b02 − b2 ) (4.47)
where b01 and b02 are the new values of the phase biases after a cycle slip. A cycle slip
can happen on one or both of the carrier frequencies at once, so any detection and
repair scheme must be sensitive to both cases, and able to separate out the effect of
the slip on each of the frequencies.
The wide-lane phase combination of the two carrier frequencies is given by:
f1 Φ1 − f2 Φ2
Φδ ≡ (4.48)
f1 − f2
f1 f2
=ρ+I 2 + λ δ bδ + M (4.49)
f1 − f22
To isolate the wide-lane bias, the following pseudorange combination may be sub-
tracted from Φδ :
f 1 P 1 + f 2 P2
Pδ ≡ (4.50)
f1 + f2
f1 f2
=ρ+I 2 +M (4.51)
f1 − f22
Then subtracting Equation 4.48 from Equation 4.50, the wide-lane bias, also known
as the Melbourne-Wübbena combination, is given by
1
bδ = (Φδ − Pδ ) (4.52)
λδ
66
Chapter 4. PPP method
1
hbδ ii = hbδ ii−1 + (bδi − hbδ ii−1 ) (4.53)
i
1 n o
σi2 = σi−1
2
+ (bδi − hbδ ii−1 )2 − σi−1
2
(4.54)
i
with the calculation of the variance being an approximation with error of O(1/i2 ),
where i is the number of epochs since the last cycle slip.
ΦI ≡ Φ1 − Φ2
= I + λ 1 b1 − λ 2 b2
(4.55)
= I + λ1 (b1 − b2 ) + (λ1 − λ2 )b2
= I + λ 1 bδ − λ I b 2
where the ionospheric wavelength λI ≡ (λ2 − λ1 ) ≈ 5.4 cm. Due to the short
wavelength of the ionospheric phase combination, the pseudorange cannot be used
to remove the ionospheric delay as the noise level of pseudorange observations is too
high. Instead the ionospheric combination is time differenced, removing many of
the highly time correlated errors. As the ionospheric combination is free from the
receiver–satellite geometry, what remains in the time difference is the rate of change
67
Chapter 4. PPP method
ΦI (tk ) − ΦI (tk−1 )
∆ΦI = (4.56)
tk − tk−1
The theory is that the rate of change of the ionosphere is limited to some bounding
value. If |∆IΦ | is greater than some threshold value then a cycle slip is detected.
Having detected a cycle slip, the next step is to find the integer values of the slip
allowing the estimated cycle slip to be corrected to match the current state of the
receiver. If the data is being post processed then the sections of cycle slip free carrier
phase observations, commonly know as “phase connected arcs”, may be identified.
For each pair of arcs, interrupted by a cycle slip, the two values of hbδ i are known.
The wide-lane cycle slip is then determined by subtracting the two hbδ i. How close
this value is to an integer, combined with the uncertainty in the two values of
hbδ i, allow the integer discontinuities to be identified with a high level of statistical
confidence. Similarly for the ionospheric combination, as the arc before and after
the slip are known, polynomial fitting may be performed. The polynomial may then
be extrapolated forwards and backwards across the disconnect allowing the cycle
slip to be reliably determined.
For processing in real-time mode these options are not available, so an instanta-
neous method to fix the slips is required. Banville and Langley (2010) proposed a
method based on the techniques used in relative positioning ambiguity fixing but
using time differenced measurements. Most of the errors affecting GNSS signals are
highly time correlated so over a short time period they may be mostly removed by
time differencing. The time difference is composed of the relative receiver/satellite
velocity, clock drift and any cycle slips:
The cycle slips detected, as well as the receiver velocity and clock drift, are solved for
in a least squares adjustment. In comparison to the real valued bias term in the PPP
observation model (Equation 4.16), the cycle slip bias term, ∆ni , is integer valued.
The techniques of integer least squares, such as the LAMBDA method (Teunissen,
1995) may now be applied to determining the cycle slips. The success rate of this
technique is dependent on the number of simultaneous cycle slips. In the case of
total loss of lock where there is a cycle slip on all channels, the pseudorange is used
to constrain the least squares solution. Thus, in this extreme case the success of the
68
Chapter 4. PPP method
If it is not possible to fix a cycle slip then the row and column in the covariance
matrix Pk (−) corresponding to the bias can be zeroed and then reinitialised from
the carrier phase - pseudorange difference (Equation 4.34).
Observations contain the range from satellite to receiver, path delays and advances
due to atmospheric effects, the error due to inaccurate values for satellite position
and clock offset, as well as relativistic effects. Also there are deformation effects
where a given location on the Earth moves with respect to the chosen reference
frame due to Earth body tides, ocean tide loading and atmospheric pressure loading.
Finally there is the accuracy with which the receiver is able to measure the received
signals and the effect of multipath. The combined residual effect of these phenomena
on the range after modelling and mitigation comprises the error budget (Chapter 3).
where H [1] is the linearised observation equation and P is the weight matrix.
The subscript capital X denotes that the coordinates are in an ECEF reference
frame.
69
Chapter 4. PPP method
the diagonal elements are used for the following DOP definitions for Geometric
Dilution Of Precision (GDOP), Positional Dilution Of Precision (PDOP) and Time
Dilution Of Precision (TDOP):
√
GDOP = qXX + qY Y + qZZ + qtt (4.60)
√
PDOP = qXX + qY Y + qZZ (4.61)
√
TDOP = qtt (4.62)
The positional elements of the global cofactor matrix QX can be transformed into
the local cofactor matrix Qx in the topocentric local coordinate system with axes
along the local north, east and up directions as:
q
nn
qne qnu
T
Qx = RQx R = qne qee qeu (4.63)
qnu qeu quu
where the rotation matrix RT = [neu] contains the axes of the local coordinate
system. Due to the invariance of the trace of a matrix with respect to rotation, the
PDOP value in the local system is identical to the value in the global system.
Two further DOP definitions that capture the dilution of precision in the local hori-
zontal and vertical, Horizontal Dilution Of Precision (HDOP) and Vertical Dilution
Of Precision (VDOP) are defined as:
√
HDOP = qnn + qee (4.64)
√
VDOP = quu (4.65)
It is then possible to map the expected range error due to the error budget discussed
in Section 3.16, to a theoretically achievable positional accuracy, as the product of
DOP and measurement accuracy.
70
Chapter 4. PPP method
Therefore given an error budget with standard deviation σrange , the positional accu-
racy would be PDOPσrange .
Many studies have shown the feasibility of PPP in achieving centimetre accurate
static PPP and decimetre level in kinematic mode (Bisnath, 2004; Bisnath and
Gao, 2008; Dixon, 2006; Gao et al., 2005; Kouba, 2009a; Muellerschoen et al., 2001).
Using the DOP to project the expected error from the measurement domain to the
position domain, it is possible to obtain a simplistic view of the expected positional
accuracy. The simplified calculation using DOP assumes that the UERE is the same
for every satellite and there is no correlation between the range errors. Multipath,
ionosphere, and troposphere delays all increase at low elevations, and are accounted
for using elevation dependent weighting. Errors in the satellite clock and orbit
products are correlated; Zumberge et al. (1997) note that combining orbit and clock
products from different sources results in degraded position repeatabilities, showing
the two not to be independent, so that the absolute measures of accuracy in both
orbits and clocks are reduced when combined to produce a range error bias.
Figure 4.2 shows the DOP values for a receiver in Newcastle upon Tyne, UK on
14/04/2011 using a 10◦ elevation cutoff. As is expected from the increased number
of satellites the DOP values when including GLONASS are reduced. Additionally
the DOP spikes present when only considering GPS are greatly reduced with a more
consistent DOP.
Using the average values from the above example and the residual error budget
for carrier phase observations from Chapter 3, gives the one sigma horizontal and
vertical positioning accuracy shown in Table 4.3.
GPS
HDOP 1.0 x 8.3 = 8.3 cm
VDOP 1.6 x 8.3 = 13.3 cm
GPS/GLONASS
HDOP 0.75 x 8.3 = 6.2 cm
VDOP 1.1 x 8.3 = 9.1 cm
71
Chapter 4. PPP method
Figure 4.2: GPS (top) and GPS/GLONASS (bottom) DOP values for
14/04/2011 (Data from Trimble GNSS Planning Online [Link]
[Link]/GNSSPlanningOnline)
The values in Table 4.3 are instantaneous accuracies assuming the ionosphere free
carrier phase ambiguity has been estimated with sufficient accuracy to not affect
the error budget (centimetre accuracy) i.e. after convergence. For static positioning
the accumulated average of many observations will remove the effect of any zero
mean biases on the estimated coordinates resulting in the higher accuracy seen in
the literature (Geng et al., 2010c).
Initially, the ambiguities are unknown and the accuracy is determined by the pseu-
dorange error budget and the dilution of precision. Most of the pseudorange noise is
due to tracking noise and multipath. As these are not constant bias parameters, time
averaging reduces the impact of the errors on the range error. This is the mechanism
behind convergence whereby sufficient observations are built up in the estimation
filter to reduce the uncertainty in the value of the carrier phase ambiguities, and
therefore increase the weight of the carrier phase observations in the position esti-
mation. Therefore the convergence time is largely defined by the pseudorange noise
and geometry.
72
Chapter 4. PPP method
4.10 Summary
This chapter has explained the principles of the Extended Kalman filter, as well
as the basis for choosing this technique by comparison with alternative estimation
methods. The implementation of the PPP method is described, including definition
of the observation model. The mathematical basis of the estimation process, as well
as the detailed methods for detecting outliers and cycle slips, are taken forwards
into the next chapter which describes the software developed to implement these
models.
73
Chapter 5 Software development
In order to achieve the research objectives, and investigate the impact of combining
GLONASS with GPS data on the performance of PPP as a navigation method, a
software package pppncl was created. This software was capable of performing PPP
in real time mode, with the ability to process data from GLONASS as well as GPS.
This software was created based on the theory described in the preceding two chap-
ters. The source code for Track (Chen, 1998) was used as a starting point, providing
some of the signal propagation delay models, planetary locations, Earth tide routines
and GPS-only RINEX file reading.
The following sections provide a brief description of the software design and ar-
chitecture used, as well as reflecting on the approach to developing this software
package.
Following definition of the research problem and consideration of the proposed ex-
periments to be undertaken, it was determined that pppncl needed the following
functionality in its design:
For the software to operate in real-time mode, and function as if it were operating
in real-time, there cannot be multiple passes through the input data. Therefore the
Kalman filter must operate in filtering mode, where the estimated state at a given
time is based only on observations available at or before that time. Additionally,
the cycle slip detection algorithms are designed such that they do not require future
data.
To allow processing of kinematic datasets the software was designed with no con-
straint, or assumption, that the predicted position is close to the previous position,
therefore all modelled parameters that depend on the receiver position are recalcu-
lated with every new observation epoch.
74
Chapter 5. Software development
The real-time mode operation forms the design for the high level architecture, which
is discussed in Section 5.3.
In Chapter 3 the models required to minimise the PPP error budget are described,
as not all of these existed in Track the following needed to be implemented:
• ANTEX format file parsing to read phase centre offsets and variations
Having made a software program capable of processing GPS data using the PPP
method, extending this software to include additional GLONASS observation data
did not require significant changes to the high level software design. GLONASS
is very similar to GPS; both systems are dual frequency with the two base carrier
frequencies at a similar frequency. As the same theory applies to both systems for
many of the physical effects that are modelled, the same processing strategy can be
applied.
75
Chapter 5. Software development
Another modification required was to update the RINEX reading routines to parse
GLONASS observation and navigation data, as well as GPS data.
The increased satellite numbers due to the inclusion of GLONASS require many of
the data storage arrays and matrices to be increased in size. Specifying the array
size parameters as a constant allows the same value to be reused in the many places
it is required whilst providing the ability to easily change all such instances by only
modifying the value in one place. The code was written such that the dimensions of
the storage arrays are calculated based upon the number of satellite systems and the
number of satellites within each system. This approach offers a simple solution to
updating this area of the software to include additional satellite networks in future.
The software program pppncl developed as part of this thesis is written largely in
Fortran 90, with certain reused components in Fortran 77. Fortran was chosen
due to the ease of reusing existing code from Track that contained many relevant
subroutines. Fortran 90 is a superset of the language specified in Fortran 77; the
additions include modules, array syntax and derived types.
Modules allow for grouping of functions, subroutines, data types and variables.
Grouping sections of code that are conceptually related helps to organise the code
during creation and further development. Additionally module interfaces allow the
compiler to check that functions and subroutines are being called with the correct
number and type of input and output variables passed; this helps catch a common
cause of programming error.
76
Chapter 5. Software development
Array syntax provides native language support for simple vector mathematics. This
allows for elimination of many loops that would be required in Fortran 77.
Derived types allow for the representation of abstract quantities in a single variable.
pppncl makes use of this to enable the use of object-oriented concepts in Fortran 90
as described by Decyk et al. (1997). For example the Melbourne-Wübenna cycle slip
detector described in Section 4.8.2 is implemented in this way. The running averages
and standard deviations, along with the required book-keeping, are encapsulated
within one variable.
Whilst Fortran provides good syntax for mathematical operations and produces fast
numerical code it was found during development of pppncl to be less suited to exper-
imentation than higher level more dynamic languages such as Python. Indeed the
absence of more complex data types and only a minimal standard library in Fortran
often lead to slower implementations due to the laborious nature of implementing
basic functionality from scratch.
The majority of variables in Track are shared through Fortran common blocks. The
use of such global variables is often considered to make programs more difficult to
understand (Wulf and Shaw, 1973). An attempt was therefore made to create pppncl
such that subroutines were made smaller and operated only on those variables that
were passed in.
An attempt was made to create an interface to the modelling and linear combination
of measurements, as outlined below, as a single function. Initially this seemed to
work well but as the software grew as additional models and GLONASS observations
were added, issues began to arise. As the modelling was separate from the code that
dealt with the estimated state vector, this resulted in many values being copied into
and out of the Kalman filter matrices and state vector. Had the modelling been
carried out as a distinct step, the prefit residuals and observation process noise
could have been calculated directly in the Kalman filter.
Due to the experimental nature of the developed software, it was sometimes found
that a calculated value that had been encapsulated within a function was needed at
a higher level in the call stack. The advantage of widespread use of shared variables
such as used in Track is that there is no need to predict this requirement in advance.
The structure of the program is described below with details of the implementation
given in Appendix A. During the design phase the potential future conversion to a
true real-time PPP positioning engine was considered, hence each loop starts with
77
Chapter 5. Software development
∗ Model carrier phase and code phase pseudoranges for each satellite
78
Chapter 5. Software development
There are several file formats designed and used by the wider GNSS community
to allow exchange of information. The use of open standards allows the use of ob-
servation data from the many receivers, and navigation data produced by analysis
centres which is freely available on the Internet. By implementing these input for-
mats in pppncl, this research was able to benefit from the work of the wider GNSS
community.
This section describes the types of input file that can be read by pppncl when
performing PPP. The different formats that can be used for each type of input data
are explained.
pppncl can read the RINEX 2 format for observation data and navigation data (or-
bits and clocks). RINEX 2 was designed with only GPS in mind, and later adapted
to allow mixed data from GPS and GLONASS satellites in the same file. It has be-
come apparent that it has shortcomings as a receiver observation exchange format
for files containing observation data of more than one satellite system, each one with
different observation types. To address this Gurtner and Estey (2007a) proposed
the RINEX 3 standard. Whilst it seems logical that RINEX 3 will replace RINEX
2 as the standard data exchange format, the not inconsiderable effort required to
update all the tools used by the wider GNSS community has meant that RINEX
3 adoption has been very slow. For that reason only a RINEX 2 reader has been
implemented in pppncl currently.
The common format for precise navigation data is the SP3 revision c format (Hilla,
2002). The precise navigation data produced by IGS analysis centres is distributed
in SP3 format, in daily (24 hour) files. For observation windows of longer than a day,
multiple SP3 files are required to span the observation interval. There are known
discontinuities at the day boundary in the generated precise orbits (Griffiths and
Ray, 2009), in pppncl interpolation is performed across this discontinuity. As there
is not a commonly used tool to concatenate multiple SP3 files together, the facility
to read directly concatenated SP3 files was added to pppncl. Multiple header blocks
are therefore allowed in pppncl SP3 input files.
High rate clock data, essential to high rate PPP can be read in RINEX 2 or 3 clock
format (Ray and Gurtner, 2010). As with the SP3 files the ability to read multiple
concatenated files was implemented to enable processing periods of more than one
day.
79
Chapter 5. Software development
Within pppncl the different navigation formats are abstracted behind a common
navigation data interface. The underlying navigation data may be a broadcast
ephemeris in RINEX navigation format or a precise ephemeris in SP3 format, pos-
sibly augmented with higher rate satellite clock corrections in RINEX clock format.
Within the modelling functions the underlying navigation data source is hidden
providing a simple “get satellite position at time t” subroutine.
For the mapping between GLONASS satellite number and GLONASS frequency
number, the broadcast GLONASS navigation data is always required even when
using precise orbits in SP3 format.
The standard format for distributing satellite and receiver antenna phase centre
offsets and variations is the ANTEX format (Rothacher and Mader, 2003; Rothacher
and Schmid, 2010). An ANTEX file parser was written to read ANTEX v1.4 format
files which contain the required satellite and receiver PCO and PCV as discussed in
Section 3.7.
DCBs are distributed in a tabular format for use by CC2NONCC (Romero, 2010). In
order to avoid having to pre-process RINEX observation files from receivers that
record C1 rather than P1, the ESOC updated version of the CC2NONCC software
routines were integrated into pppncl.
Within pppncl there are some parameters that are needed in order to run and others
for which there are a range of choices based on enabling specific functionality or a
numerical value describing a particular property. For example, the location of the
input and navigation files is required. The ability to enable optional functionality
such as using GLONASS satellites and change state vector process noise, for example
to correspond to the expected receiver dynamics, was also required.
The use of a configuration file allows for the input parameters to be adjusted easily.
Compared to command line options this provides a documented record of the input
options and allows for easily repeatable test runs. The syntax was kept simple to
ease creation and parsing with a keyword on each line identifying the parameter
that is being configured.
80
Chapter 5. Software development
The output files and formats used were created to satisfy the needs of the analysis
produced in the following chapters. They are all space separated tabular text files.
Plain text output files do not require special tools for manipulation and provide the
maximum compatibility for use with plotting tools.
To aid setup and help with running repetitive tasks a number of Python scripts have
been created. The Python programming language provides strong text manipulation
functionality along with a comprehensive standard library, including the ability to
easily download files from the Internet.
For every run of pppncl the following tasks must be performed: Acquiring and
preparing all the required input files then creating a configuration file that specifies
the prepared input files. To avoid having to do this manually when the desire is to
run many sets of data that share a largely common set of configuration parameters,
a set of Python scripts was created to facilitate the set up.
Scripts were written to address two common use cases of pppncl: estimate the
position of a receiver using observations from a given RINEX file, and estimate
the position of one or more named sites from the IGS network. The configuration
parameters not relating to the input files can be specified as a template to the
scripts. This facilitates estimating many sites using the same set of configuration
parameters.
The data processing involved in generating the results presented in Chapters 6 and
7 required many runs of pppncl. The run time of pppncl is approximately 4 minutes
for 24 hours of positioning at a 1Hz data rate on a 2.7 GHz Intel i5-2557M processor;
performing this operation hundreds or thousands of times on a single computer takes
a prohibitively long time.
81
Chapter 5. Software development
During use these scripts became the user facing interface to the program. Given that
the software has become a composite of two languages each with different strengths,
a more flexible way of constructing the program would be to create the split between
Python and Fortran at a lower level. In this way the models and processing steps
such as outlier detection, cycle slip detection and repair, and Kalman filtering would
become composable units. This would simplify the “glue” code that sets the program
up and calls each unit in turn. New ideas or approaches could be rapidly prototyped
in the high level language, only being re-implemented in Fortran if successful and
there was a need for reduced runtime.
5.8 Summary
This chapter has discussed the key areas that were considered in the creation of
pppncl from a software design and implementation perspective. This includes the
approach used to implement both a PPP software package and the extension to
include GLONASS. Related sections have reflected on the high level software archi-
tecture, specifically the importance of code structuring.
The needs for input and output data formats have been examined; the approach
taken enabled full use of available datasets. Further, the practical implications
of using pppncl to analyse large and varied datasets were addressed, with specific
attention to the tools used to configure the program, process multi-day data, as well
as the computer power required. The combination of all these factors has enabled
the successful creation of pppncl, the use of which has produced the results required
to answer the research questions investigated in this thesis within a practical time
frame.
For reference the specific implementation details of pppncl are given in Appendix A.
82
Chapter 6 PPP validation
Before using the created program, pppncl, to quantify the impact of the addition of
GLONASS observations on GPS PPP, it was necessary to verify the performance of
the PPP implementation using only GPS observations. This verification included
testing the performance of the program for both static and kinematic receivers.
In validation of both kinematic and static PPP implementation, the following char-
acteristics are desirable. The coordinates of the receiver must be known, to provide
the reference solution against which the estimated position is compared. In the ideal
case, there is no error in the reference coordinates. The apparent accuracy of the po-
sition estimation is affected by the accuracy of the truth coordinates. This reduces
the ability to determine the true accuracy of the estimated coordinates. For example
any correlation in the errors in the estimated and truth coordinates, will appear to
improve the accuracy of the estimated coordinates. A globally distributed set of
sites is desirable to capture the range of the parameters which change with location;
these include different satellite geometries, Earth tide and atmospheric variations.
Additionally, analysis of multiple sites also ensures that any conclusions drawn are
representative of the performance of the PPP implementation being validated.
It is preferable to cover the longest possible time period to account for temporal
variation. Ideally, this would be one long continuous dataset (>1 year), or multiple
shorter datasets dispersed throughout the year. There is a minimum useful time
length for a specific dataset, as it needs to be long enough to capture the conver-
gence period (order of hours). In the literature, static positioning accuracy is often
considered over a 24 hour period, therefore this minimum length is useful for cross
comparison of results. By contrast, in kinematic validation, the minimum useful
time period is the period for the solution to converge, approximately 4 hours based
on typical convergence times seen by Bisnath and Gao (2008). A longer period is
preferable, providing more data to analyse the performance metrics in a converged
state. In the ideal case, the validation scheme would cover a period of several years,
to include a variety of stages in the solar cycle, a longer period temporal parameter.
Whilst the handling of multipath and cycle slips is important in terms of under-
standing the achievable accuracy and integrity of PPP, in terms of validating the
83
Chapter 6. PPP validation
software, they complicate the assessment of best case accuracy and convergence
time. It is known that multipath degrades the positioning accuracy but as no new
techniques have been implemented in pppncl to mitigate it, this is not a focus area
of the validation. Similarly, cycle slips, satellite obstructions or loss of lock will
interrupt the filter convergence, and therefore reduce the utility of a dataset in as-
sessing convergence time. Therefore, a “clean” environment is desirable, with low
multipath and an uninterrupted sky view.
The other aspect of a validation scheme that should be considered is the statistical
methodology used to compare the results of the PPP program to the reference
solution, in particular accounting for uncertainties in the estimated parameters. For
static positioning, the statistical measures that should be compared to validate the
test solution versus reference solution are the mean, standard deviation, and the
RMS. For kinematic datasets, in addition to these measures, the convergence time
was also evaluated to show equivalent results to those published in literature. In
this thesis, the static and kinematic PPP performance results were compared to
the specific reference solution in each case, as well as in more general terms, to the
performance metrics given in literature.
Whilst the above characteristics describe the ideal validation scheme, there are prac-
tical limitations to executing any validation scheme. In terms of the reference so-
lution, exact coordinates do not exist. Any set of coordinates is measured with
respect to some reference frame; both the definition of the reference frame, and the
measurements introduce errors into the truth coordinates given. This complicates
the analysis of the PPP estimation accuracy; with no error, any difference between
the PPP and reference solution can be attributed to the PPP solution.
84
Chapter 6. PPP validation
Throughout this thesis, orbit and clock products are used which have been obtained
from IGS analysis centres and Veripos. Different products are used depending on
the requirements of the test, and to provide comparison of available accuracy be-
tween different products. The advantages of using IGS analysis centre clocks are
summarised by Kouba (2009a):
"The users of the IGS products in fact take full advantage of the IGS AC
global analyses, properly combined and quality checked, all in accordance
with the current international conventions and standards."
Use of the publicly available IGS analysis center products makes the tests repeatable.
Whilst the software has been formulated in a real-time mode, the majority of tests
are performed using final orbit and clock products which are not available until one
to two weeks after the date they relate to. In assessing static positioning, much of
the published work, for example Leandro et al. (2011) and Bisnath and Gao (2008)
make use of the combined IGS final orbit and clock products. In order to provide a
comparison with other published results in assessing the implementation of pppncl,
IGS final orbits were used in the static validation.
The individual orbits produced by the IGS analysis centers agree to within 1 cm
to 2.5 cm 1D weighted RMS (Griffiths and Ray, 2009). The error budget derived
in Section 3.16 can therefore be applied to position estimation using orbits and
clocks from individual IGS analysis centres. For the GPS only kinematic validation,
CODE orbits and clock products were used due to the availability of high-rate clocks
tabulated at 5 seconds that capture the short term variability of the satellite clocks,
reducing this error when performing 1 Hz positioning (Bock et al., 2009).
GPS/GLONASS PPP requires GPS and GLONASS orbit and clock products. Cur-
rently only ESA (Springer, 2010) and IAC (Oleynik et al., 2006) produce a GLONASS
precise clocks product. The ESA products are used in this thesis as their formula-
tion as a combined GNSS product removes the need for coordinate and time scale
transformations. Published accuracies for the ESA GLONASS clock products are
not available as with only one other analysis centre producing a clock product a
robust combination against which to compare the ESA product is not available.
The results obtained using the Veripos orbit and clock products are indicative of the
85
Chapter 6. PPP validation
potential performance of real time positioning as these are real time products; the
other results, using the IGS analysis centre post processed final products, demon-
strate the best available performance.
The initial stage was to assess the precision, accuracy and convergence time of the
software for static receivers. Different locations around the Earth experience dif-
ferent satellite geometries, and variations in tropospheric and ionospheric effects.
Therefore, a range of globally distributed sites was selected in order to give a com-
plete picture of performance. The IGS is a voluntary collaboration of over 200
organisations that maintain a network of more than 300 globally distributed GNSS
tracking stations (Kouba, 2009a). These sites are installed according to a set of best
practices in order to minimise multipath and provide high quality data. Additionally
relevant metadata for the sites is available including antenna offsets. All the sites
for which 24 h of data was available were used to provide the fullest global coverage.
Coordinates for the receivers in the IGS network are estimated weekly and released
by the IGS in SINEX format. The horizontal coordinate accuracy is 2 mm to 3 mm
and the vertical accuracy is 7 mm (Ferland and Piraszewski, 2009). This provides
a set of reference coordinates with which a PPP solution can be compared. The
location of the 281 sites used in this study is represented graphically in Figure 6.1.
Figure 6.1: Map of the sites used in static PPP GPS verification,
showing the global distribution
86
Chapter 6. PPP validation
The concentration of sites in Europe and the absence of sites in North-Central Africa
reflects the distribution of the underlying IGS network. The full list of station names
and coordinates may be found in the Appendix.
The Crustal Dynamics Data Information System (CDDIS) archive stores daily ob-
servation files for sites in the IGS network (Noll and Dube, 2001). Dual frequency
code and carrier phase observations are stored at 30 s observation intervals in com-
pressed RINEX format files. The observation data described was downloaded for
each of the stations for 10 October 2010 (GPS week 1605, day 0). Temporal changes
are therefore not considered in this test but are covered in Section 6.1.3.
24 h of observation data was processed for each site in pppncl using the following
options:
• The final satellite ephemerides from the IGS network solution in SP3 format
were downloaded for the day of study and the preceding and following day
from [Link] and concatenated into one file
in order that orbit interpolation could be performed at the beginning and end
of the studied day without being subject to boundary effects.
• Additional final IGS satellite clock offsets tabulated at the higher rate of 30 s
were also downloaded in RINEX clock format from the same location. The
clock products contain data from 00:00:00 to 23:59:30 so, due to the linear
interpolation used for the clock offsets, only data for the processed day is
needed. There is a small possibility of extrapolation for the final epoch of the
day at 23:59:30 but as most receivers maintain their clock to within 1 ms of
GPS time it would be a very small extrapolation.
• For receivers not recording the P1 code observation and only recording the
C1 code observation the CODE P1C1 weekly bias solutions, available at ftp:
//[Link]/aiub/CODE/p1c1bias.2000p, were used to convert the C1
observation into a P1 equivalent observation.
• To correct for antenna phase centre offsets and variations, the values from
the IGS antex file available at [Link]
general/[Link] were used.
• Initial receiver coordinates were taken from the RINEX file header for each
site.
87
Chapter 6. PPP validation
• An initial position uncertainty of 200 m was used as the RINEX header co-
ordinates are only approximate and to reflect the more general positioning
scenario where precise initial coordinates are not known. This also ensures
the PPP results reflect the accuracy of the method and are not simply the
results of tightly constrained initial values.
For each site the full 24 h of observation data was processed and the post fit state
vector saved to a file at each 30 s observation epoch. The observation data was pro-
cessed in a real-time mode reading and processing one epoch at a time, sequentially
starting at midnight. At each epoch the following were estimated:
• an ionosphere free carrier phase bias for each satellite, modelled as a constant.
Reference coordinates for each site were extracted from the SINEX file ftp://
[Link]/gps/products/1605/[Link] containing the IGS com-
bination solution for GPS week 1605. The final position estimate from the pppncl
solution after processing the complete 24 h was compared to the reference coordi-
nates from the IGS weekly solution. As no process noise is added to the coordinate
covariance at each epoch, the values of the estimated coordinates at each epoch are
the current estimate of the coordinates based on all the observation data incorpo-
rated up to that point. The final position estimate at 23:59:30 is therefore the single
estimated set of coordinates based on the complete day’s observations.
Figures 6.2, 6.3 and 6.4 show the distribution of the differences in the PPP coordi-
nates with respect to the reference solution in the North, East and Height directions
respectively, with the statistics summarised in Table 6.1. As would be expected the
precision is highest in the North component and lowest in the Height. The Height
component is always less accurate in GNSS due to all the satellites being above
the receiver and the high correlation between the Height, receiver clock offset and
zenith tropospheric delay. The reduced accuracy of the East component is a feature
commonly seen in carrier phase GNSS positioning with floating ambiguities. This
88
Chapter 6. PPP validation
is caused by the higher correlation of the phase biases with the East component, at
mid and low latitudes, due to the North-South ground track of GPS satellites at the
equator in the ECEF reference frame (Blewitt, 1989).
40
N (0, 0.003)
35
30
Number of sites
25
20
15
10
0
−0.015 −0.010 −0.005 0.000 0.005 0.010 0.015
North component error (m)
Figure 6.2: Static PPP error in the North component compared to the
IGS weekly solution. A Normal distribution with standard deviation of
the reference coordinates is overlaid.
70
N (0, 0.003)
60
50
Number of sites
40
30
20
10
0
−0.03 −0.02 −0.01 0.00 0.01 0.02 0.03
East component error (m)
Figure 6.3: Static PPP error in the East component compared to the
IGS weekly solution. A Normal distribution with standard deviation of
the reference coordinates is overlaid.
89
Chapter 6. PPP validation
60
N (0, 0.007)
50
Number of sites
40
30
20
10
0
−0.06 −0.04 −0.02 0.00 0.02 0.04 0.06
Height component error (m)
Figure 6.4: Static PPP error in the Height component compared to the
IGS weekly solution. A Normal distribution with standard deviation of
the reference coordinates is overlaid.
Table 6.1: Static PPP error compared to IGS weekly solution for 282
stations on 10 October 2010
The accuracy achieved in this test will represent the best static performance of
pppncl. This is due to the use of IGS sites which are of very high quality. The
use of high grade geodetic receivers minimises measurement noise; the common use
of choke ring antennas and the careful location of sites reduces receiver multipath.
Therefore two of the major contributions to the error budget are minimised in this
data set.
Kouba (2009a) shows the coordinate RMS achieved at 36 IGS Reference frame
stations during 24 h static positioning using GPS Pace PPP. For GPS week 1516
(January 25–31, 2009) using IGS Final orbit and clock products, RMS values of
3 mm, 5 mm and 14 mm in North, East and Height, respectively, were obtained.
Thus it can be concluded there is not a considerable difference between the pppncl
results and those from Pace PPP.
90
Chapter 6. PPP validation
A further comparison may be made for the estimated tropospheric ZTD. The IGS
generate a tropospheric product containing the total (hydrostatic + wet) ZTD at
5 min intervals (Byun and Bar-Sever, 2009). Using a post-processed PPP method,
daily observation files are processed for each site in the IGS network. The key
features of the processing approach are (Bar-Sever et al., 1998):
• Software: GIPSY
• A priori hydrostatic delay based on altitude (2.3 m at sea level), and 0.1 m for
the wet delay
Formal errors for the final tropospheric product are typically (one sigma) 1.5 mm
to 5 mm (Byun and Bar-Sever, 2009). Actual errors are known to be higher, due
to systematic errors in the GPS orbits and clock products. The difference between
the IGS final tropospheric product and the ZTD product, from each IGS analysis
centre, for all sites, has a daily bias of up to 5 mm. The standard deviation of this
measure for each day is 3 mm to 10 mm (Byun and Bar-Sever, 2009).
Using the IGS final ZTD product as a reference, the pppncl estimated ZTD can be
evaluated. The IGS final ZTD product was downloaded for all sites for 30 July 2010.
91
Chapter 6. PPP validation
The two sites showing the largest and smallest daily variation in ZTD were chosen
to highlight the performance at the two extremes. HNPT, located in North-East
America, showed the largest variation in ZTD with a range of 194.3 mm. MAUI,
located in Hawaii, had the smallest daily ZTD variation with a range of 6.9 mm.
The daily RINEX observation files for HNPT and MAUI were processed using ppp-
ncl. The same processing options as described in the generation of the IGS ZTD
product were used, with the exception of the GMF in place of the Niell mapping
function, and no estimation of atmospheric gradients. Observation data was pro-
cessed at 30 s and the total ZTD extracted at 5 min intervals matching the IGS
product.
2.70
IGS
2.65 ppptrack
Total zenith delay (m)
2.60
2.55
2.50
2.45
2.40
0 4 8 12 16 20 24
Time (h GPST)
Figure 6.5: Total ZTD for IGS station HNPT, 30 July 2010
Figures 6.5 and 6.6 show the estimated total ZTD from pppncl and the IGS final
product. The real-time mode nature of pppncl is seen during the first 30 min as
the filter converges. For HNPT the pppncl estimate follows the rapid drop in the
delay, although a slight lag can be seen. This is likely to be due to the forward only
nature of the pppncl Kalman filter combined with the low weighting of satellites at
low elevations.
The mean bias and standard deviation between the pppncl estimates and the IGS
product are not far from those reported in Byun and Bar-Sever (2009). For MAUI
with the low change in ZTD throughout the day, the difference between the PPP
92
Chapter 6. PPP validation
IGS
1.70
Total zenith delay (m) ppptrack
1.68
1.66
1.64
1.62
0 4 8 12 16 20 24
Time (h GPST)
Figure 6.6: Total ZTD for IGS station MAUI, 30 July 2010
Table 6.2: The daily bias and standard deviation for pppncl estimated
ZTD relative to the IGS ZTD product. The first hour of processing was
not included to remove errors due to convergence. All units are in mm.
estimated ZTD is comparable to the variability seen amongst IGS analysis centres.
HNPT, representing an extreme case, shows almost twice the bias seen at MAUI,
however much of this can be attributed to the lag in response of the filter estimate,
an artifact of real-time mode processing.
93
Chapter 6. PPP validation
MDVJ
NANO BRST
KOUR
REUN TOW2
CONZ
Figure 6.7: Map of the seven IGS stations used to characterise GPS
PPP convergence behaviour
Figures 6.8, 6.9 and 6.10 show the evolution of the difference between the PPP
solution and the IGS05 reference coordinates for each site on a log timescale. As
in the previous test, the North component shows the smallest error and also the
fastest convergence time. It should be noted that a larger scale has been used for
the Height plot due to the larger errors present.
0.2
BRST NANO
CONZ REUN
KOUR TOW2
0.1 MDVJ
North error (m)
0.0
−0.1
−0.2
Table 6.3 shows the mean, minimum and maximum of the magnitude of the position
errors after processing 1 h, 2 h, 3 h, 6 h, 12 h and 24 h of observation data. The
improvement in positioning accuracy with time is seen clearly in both the graphs and
94
Chapter 6. PPP validation
0.2
BRST NANO
CONZ REUN
KOUR TOW2
0.1 MDVJ
East error (m)
0.0
−0.1
−0.2
0.4
BRST NANO
CONZ REUN
KOUR TOW2
0.2 MDVJ
Height error (m)
0.0
−0.2
−0.4
95
Chapter 6. PPP validation
the table. 10 mm level accuracy is reached after 1 h for the North component, 3 h for
the East component and 6 h for the Height component. There is little improvement
in the accuracy of North and Height after 12 h, whereas the East component accuracy
does improve during the final 12 h.
Using the same procedure as above, 154 days of observation data were processed
using pppncl for the IGS station BRST. The days processed were in the period
30 July to 31 December 2010. The resulting coordinates estimated for every 30 s
observation interval were compared to the IGS weekly coordinate estimates for the
GPS week of that observation.
The resulting position errors in both horizontal and 3D were calculated. The hor-
izontal error is given by the distance of the PPP coordinates from the reference
coordinates in the horizontal plane defined by the North and East components.
The 3D error is given by the total distance of the estimated coordinates from the
reference coordinates. Given the North, East and Height errors dN , dE, and dH:
√
2D horizontal error = dN 2 + dE 2 (6.1)
and
√
3D error = dN 2 + dE 2 + dH 2 (6.2)
For each 30 s interval of the day this gave a set of 154 position errors, one from
each of the days processed. As the processing started at midnight each day the
coordinate error at any time during the day can be seen as the achievable accu-
racy after that length of PPP processing. By considering the 154 days as a whole,
96
Chapter 6. PPP validation
the statistical properties of pppncl static GPS PPP can be examined. The errors
for each epoch were sorted and the value corresponding to the 50th , 68th and 95th
percentile extracted. Figures 6.11 and 6.12 show how the value which bounds the
given percentile in the 2D horizontal and 3D coordinate errors evolves with time.
Horizontal accuracy of 10 cm is reached 95% of the time after 1 h, and 5 cm accuracy
is achieved within 2 h. For 3D positioning, in 95% of cases 10 cm accuracy is reached
within 2 h, and 5 cm within 5 h.
0.20
0.18 50.0 %
68.0 %
0.16
95.0 %
Horizontal error (m)
0.14
0.12
0.10
0.08
0.06
0.04
0.02
0.00
(3.75) (7.5) (15) (30) 1 2 4 8 24
Observation length, parenthesised values in minutes, others hours
Héroux et al. (2004) looked at PPP convergence times using IGS precise orbits
over a 7 day period. The convergence time showed a high degree of day-to-day
variability, both within a single site and between sites (Bisnath and Gao, 2008). In
their study, convergence to within 10 cm of the IGS solution took from 30 min to
4 h. The convergence times for pppncl are comparable to these values.
97
Chapter 6. PPP validation
0.20
0.18 50.0 %
68.0 %
0.16
95.0 %
0.14
3D error (m)
0.12
0.10
0.08
0.06
0.04
0.02
0.00
(3.75) (7.5) (15) (30) 1 2 4 8 24
Observation length, parenthesised values in minutes, others hours
Positioning a site for a whole year can reveal any errors that manifest as seasonal
trends and allows repeatability to be assessed. Leandro et al. (2011) processed the
observation data from the IGS station UNBJ for every day in 2008 using their GAPS
PPP software. 24 h static positioning was performed each day to create a time series
of coordinate errors. This test was repeated using pppncl.
For the site UNBJ the full year of 2008 was processed using the same configuration
as described in Section 6.1. Each daily RINEX observation file was used to estimate
a single static position for that day.
The continental plates making up the Earth’s crust are drifting at different veloci-
ties with respect to the ITRF reference frame. This means that a “stationary site”
is in fact moving at the plate’s velocity in ITRF. The velocities for station UNB1
(which was replaced by UNBJ in 2006) in ITRF2005 are Vx : −16.8 mm year−1 , Vy :
−1.3 mm year−1 and Vz : 5.1 mm year−1 (from [Link]
solutions/2005/doc/ITRF2005_GPS.[Link]). The velocities given in ITRF2005
are linear so are unable to describe any seasonal variability in the coordinates due
to, for example, atmospheric pressure loading effects.
98
Chapter 6. PPP validation
To allow for this drift, PPP estimated coordinates were grouped by GPS week
and compared to the IGS estimated coordinates for that week. Comparing to a
weekly solution effectively ignores the tectonic change in position coordinates during
a week. The maximum resulting error in considering reference coordinates at the
weekly level can be determined by the station’s ITRF2005 velocity. For UNBJ
this is 0.3 mm, an order of magnitude lower than the uncertainty in the reference
coordinates. Therefore ignoring the intraweek velocity will have negligible impact
on the resulting PPP error estimate.
The difference between the pppncl daily position estimates and the IGS weekly
reference coordinates is shown in Figures 6.13, 6.14 and 6.15. The coordinate errors
show constant variability throughout the year. The error in the North component
is consistently below 1 cm, with East and Height components at the 2 cm level. A
linear fit to the North and East components gives low annual velocity errors of
−0.3 mm year−1 and −1.3 mm year−1 .
0.015
0.010
North Error (m)
0.005
0.000
−0.005
−0.010
−0.015
0 50 100 150 200 250 300 350 400
Day of year (2008)
Figure 6.13: North 24 h pppncl position error for IGS station UNBJ
√ √
The horizontal ( dN 2 + dE 2 ) and 3D ( dN 2 + dE 2 + dH 2 ) error in the PPP so-
lution is tabulated in Table 6.4 and graphed in Figures 6.16 and 6.17. There is
no visible seasonal trend in the coordinate RMS. The horizontal coordinates agree
with the IGS published coordinates to within 1 cm and the 3D coordinates to within
2 cm.
99
Chapter 6. PPP validation
0.06
0.04
0.02
East Error (m)
0.00
−0.02
−0.04
−0.06
0 50 100 150 200 250 300 350 400
Day of year (2008)
Figure 6.14: East 24 h pppncl position error for IGS station UNBJ
0.06
0.04
Height Error (m)
0.02
0.00
−0.02
−0.04
−0.06
Figure 6.15: Height 24 h pppncl position error for IGS station UNBJ
2D 3D
Mean Error (mm) 5.5 12.6
100
Chapter 6. PPP validation
0.06
0.05
0.04
Error (m)
0.03
0.02
0.01
0.00
0 50 100 150 200 250 300 350 400
Day of year (2008)
0.08
0.07
0.06
0.05
Error (m)
0.04
0.03
0.02
0.01
0.00
0 50 100 150 200 250 300 350 400
Day of year (2008)
101
Chapter 6. PPP validation
Leandro et al. (2011) compared the coordinates estimated by GAPS PPP software
with IGS solutions for the year 2008. The reported average 2D and 3D errors of
7.7 mm and 16.3 mm are slightly worse than those obtained using pppncl (Table 6.5).
GAPS pppncl
Mean horizontal error (mm) 7.7 5.5
Mean 3D error (mm) 16.3 12.6
Table 6.5: Mean horizontal and 3D position error of GAPS and pppncl
compared to IGS weekly solutions for 2008
Due to the presence of trees around the perimeter of the airfield, the testing was
completed in the runway area to ensure the vehicle remained well clear of any sky
view obstructions. The runway used is 1.8 km long and approximately 0.1 km wide
orientated in a SW to NE direction. The site, shown in Figure 6.18, is located at
55.02◦ N, 1.87◦ W.
• Inertial navigation
102
Chapter 6. PPP validation
Post processed relative carrier phase GPS positioning is often used to provide ref-
erence coordinates (Armatys et al., 2003). This is in part due to the simplicity,
as it only requires a reference GNSS base station. Also, there are no issues with
offsets between the GNSS antenna and the location of the receiver reference point.
In other methods the relative position between the GNSS antenna reference point,
and the reference point of that technique, must be determined and accounted for.
In this thesis, relative carrier phase GPS positioning is used to provide truth coor-
dinates due to equipment availability and Track, a carrier phase double differencing
kinematic positioning software with ambiguity resolution (Chen, 1998).
The test took place from 08:00 to 14:30 UTC on 22 August 2008. A Leica GX1230
receiver with external battery was placed in a Landrover vehicle. A Leica AX1202
antenna was mounted on the vehicle roof by means of a roof rack. The receiver
was set to record code and carrier phase on both GPS frequencies at 1 Hz. The
vehicle was driven around the airfield for six and a half hours, along various routes
(Figure 6.19), and at a range of speeds from 2 m s−1 to 34 m s−1 (4.5 mph to 77 mph),
(Figure 6.20). The vehicle was kept stationary for one hour at the start, to allow the
103
Chapter 6. PPP validation
initial set of ambiguities to be found for the reference solution via static processing
if required. Then the vehicle was driven at a slow speed of 2 m s−1 to 3 m s−1 (4.5
to 7 mph) for 2 hours back and forth along the length of the runway. Next followed
a period of slow speed driving, with continual changes in direction. The final pe-
riod involved driving up and down the runway at increasing speeds up to 34 m s−1
(77 mph).
1200
1000
800
600
North (m)
400
200
−200
−400
−200 0 200 400 600 800 1000
East (m)
Figure 6.19: Kinematic route showing the distance of the vehicle from
the base station, marked by a red dot
35
30
25
Speed (m/s)
20
15
10
0
0 1 2 3 4 5 6 7
Time since start of test (hours)
Figure 6.20: Speed of the vehicle during the kinematic airfield test
104
Chapter 6. PPP validation
In this kinematic test the true coordinates of the rover were not known. A double
difference solution was used to provide a reference solution. The double differencing
solution requires a base station nearby, ideally as close as possible, and definitely
less than 10 km, as described in Section 2.4.1. A second Leica GX1230 receiver with
a Leica AX1202 antenna was therefore placed on a tripod at the side of the airstrip
to act as a base station.
The coordinates of the base station were determined using a double difference carrier
phase static baseline to the IGS station at MORP. The coordinates of MORP were
held fixed to the values determined in the IGS weekly average coordinate solution
of 3645667.799, -107277.178, 5215053.561 in metres X, Y, Z in the IGS05 reference
frame. The coordinates of the airfield base station were then estimated relative to the
MORP coordinates giving the base station coordinates in IGS05. This ensured the
location of the base station placed on the airfield, and hence the reference solution
coordinates, were in the same IGS05 frame as the orbits used in the kinematic PPP
solution.
To ensure the reference coordinates and the kinematic PPP estimated coordinates
are in the same reference frame and both locating the same point, the following
common elements were used for both sets of processing.
Final orbits from CODE (Dach et al., 2009) were used as they provide high rate
clock products tabulated at 5 s intervals suitable for 1 Hz positioning. This is the
highest rate available from any of the IGS analysis centres so therefore the best
suited to 1 Hz PPP positioning. The orbits are in SP3 format. The orbit and
clock products were downloaded from [Link]
1493/ (orbit filename cod14935.sp3).
Absolute phase centre corrections consistent with the orbits were used in the ANTEX
format. The offset file used is available from [Link]
station/general/[Link]
The Leica GX1230 records C1 and P2 observables. The C1 observable was cor-
rected to P1 using the CC2NONCC program available from [Link]
aiub/bcwg/cc2noncc/. The P1–C1 bias values from the file p1c1bias.2000 were
used which contains the monthly bias estimates generated by CODE as part of their
global network analysis.
105
Chapter 6. PPP validation
To provide a reference truth solution with which to compare the PPP solution, post
processed coordinates were estimated using MIT Track (Chen, 1998). Track is a
post processing double difference GPS processing software package. It allows for
static and kinematic positioning of an unknown site relative to a base station.
Track runs through the data several times as each step in the processing is per-
formed. First an approximate position is determined at each epoch using the code
observations. Then using the code solution as an a priori position estimate, a
Kalman filter is used to estimate the position and floating ambiguities using carrier
phase observations. The ambiguities are then resolved in one of two ways: either a
search is performed in ambiguity space or those ambiguities that are sufficiently well
defined by the Kalman filter are fixed to the closest integer. The Kalman filtering is
then repeated with the removal of the fixed ambiguities from the state vector. This
reduces the number of degrees of freedom and constrains the solution, potentially
allowing further ambiguities to be well enough defined to be fixed.
The Kalman filter can be operated in forward, backward or smoothing mode. The
forward and backward modes only use past or future data respectively to determine
the current estimate. The smoothing mode uses both past and future observations,
improving the quality of the solution at the expense of increased computation time.
The configuration used for Track is a combination of the short baseline static mode
(SHORT) and the high rate aircraft mode (AIR). Due to the short baseline lengths
involved, L1 and L2 observations were used separately. This reduces the observation
noise versus using the ionosphere free combination Lc , whilst relying on the iono-
spheric delay being sufficiently similar for both receivers to cancel out in the double
difference.
obs_file
base base.08o F
rove rover.08o K
site_pos
base 3662715.9801 -120118.3281 5202960.6893
106
Chapter 6. PPP validation
mode short
site_stats
base 4.7 4.7 4.7 4.7 4.7 4.7
ante_off
base 0.0000 0.0000 0.0000 LEIAX1202
rove 0.0000 0.0000 0.0000 LEIAX1202
cut_off 7
OUT_TYPE NEU
ANTMOD_FILE [Link]
USE_GPTGMF
BACK_TYPE SMOOTH
The smoothing mode of the filter was selected to avoid an initial period of reduced
accuracy due to the initialisation of the Kalman filter in Track. A low elevation cutoff
of 7◦ was used due to the open nature of the site allowing additional low elevation
satellite observations to be included without being subject to high multipath and
continual obstruction. The LEIAX1202 antenna has been calibrated by robot at
Geo++ GmbH so there are phase centre variation corrections available down to an
elevation of 0◦ in the IGS05 ANTEX file [Link].
As the baseline length was short, with a maximum length of just under 1.5 km (Fig-
ure 6.21) and all the ambiguities were resolved, the quality of the reference solution
was deemed acceptable. Formal errors of the reference solution were in the range
2 mm to 12 mm as shown in Figure 6.22 and Table 6.6. The East component shows
the lowest standard deviations. Both the horizontal and the vertical components
follow a similar trend with two periods of increased standard deviation occurring
together.
The PDOP is a measure of the strength of the satellite geometry when estimating
the position. Figure 6.23 shows that there are no PDOP spikes and consistently six
107
Chapter 6. PPP validation
1600
1200
1000
800
600
400
200
0
0 1 2 3 4 5 6 7
Time since start of test (hours)
Figure 6.21: Baseline length during kinematic airfield test
or more satellites visible above the elevation cutoff. The periods of slightly higher
standard deviations in the position solution correlate with the periods of reduced
satellite availability and increased PDOP.
108
Chapter 6. PPP validation
20
North σ (mm)
10
0
20
Height σ (mm) East σ (mm)
10
0
20
10
0
0 1 2 3 4 5 6 7
Time since start of test (hours)
Figure 6.22: Formal errors of the Track reference solution at Albemarle
12
Number of visible satellites
10
0
0 1 2 3 4 5 6 7
Time since start of test (hours)
Figure 6.23: Satellites above the elevation cut-off used in the reference
solution
109
Chapter 6. PPP validation
• The same 7◦ elevation cut off as used for the reference solution.
• Additional CODE final high rate clocks from cod14935.clk_5 were used.
These are tabulated at 5 s intervals in the RINEX clock format.
• The data was processed in a real-time mode, meaning that only data recorded
up to and including the time of each estimated position epoch is used to
produce that epoch.
North Error (m)
0.2
0.0
−0.2
0.2
East Error (m)
0.0
−0.2
Height Error (m)
0.2
0.0
−0.2
0 1 2 3 4 5 6 7
Time since start of test (hours)
Figure 6.24: Error of pppncl coordinates with respect to the double
difference Track solution at Albemarle Airfield
The graphs in Figures 6.24 and 6.25 show the North, East and Height error of the
pppncl solution with respect to the reference coordinates. The initial convergence
period can be seen in each component (approximately 1 h). In the North and East
110
Chapter 6. PPP validation
0.00
−0.05
0.05
East Error (m)
0.00
−0.05
Height Error (m)
0.05
0.00
−0.05
1 2 3 4 5 6 7
Time since start of test (hours)
Figure 6.25: Zoomed in view of the error of the pppncl coordinates with
respect to the double difference Track solution after 1 h at Albemarle
Airfield
Time (min)
Absolute Error (m) 0.2 0.1 0.05
North 0 3 34
East 0 40 62
Height 5 167 382
Table 6.8: Time taken until kinematic pppncl coordinates are within a
given difference from the reference coordinates, Albemarle Airfield
111
Chapter 6. PPP validation
Table 6.8 lists the time taken until all subsequent coordinate estimates are within
0.2 m, 0.1 m and 0.05 m of the reference coordinates for North, East and Height
respectively. Within 5 min, all coordinates converge to 20 cm accuracy, and 5 cm
accuracy in the horizontal is obtained after 1 h. This is a very stringent definition
of convergence with a 100% requirement on all subsequent coordinates to be within
the limit, compared to one often used of two sigma RMS.
6.3 Conclusion
The estimated zenith path delay due to tropospheric effects was compared to the
IGS final product. Centimetre level agreement was seen for sites with both low and
high variability in the zenith delay.
Convergence times for a subset of static sites were analysed. An horizontal accuracy
of 10 cm was reached 95% of the time after one hour, and 5 cm was achieved within
2 h. For 3D convergence, 10 cm accuracy was reached 95% of the time in 1.5 h, and
5 cm accuracy was reached within 5 h. These figures are comparable to those in
other published papers (Bisnath and Gao, 2008; Héroux et al., 2004).
The coordinate errors estimated using 24 h static PPP for IGS station UNBJ for
each day in the year 2008 were compared to those obtained by Leandro et al. (2011)
using GAPS PPP. The average coordinate errors obtained by pppncl of 5.5 mm in
112
Chapter 6. PPP validation
the horizontal and 12.6 mm in 3D were slightly better than those obtained using
GAPS PPP.
Based on the results obtained in comparison to other published work, pppncl has
been seen to perform at a similar level to other PPP implementations. It can
therefore be concluded that the PPP method has been successfully implemented
and provides a useful platform from which to perform additional research.
113
Chapter 7 Combined GPS/GLONASS PPP
This Chapter examines the effect of including GLONASS in addition to GPS satel-
lite observations into PPP. Tests were carried out to evaluate the impact on the
chosen performance metrics (Section 2.9), specifically convergence time, positioning
accuracy and continuity. Firstly, the relative magnitude of the GPS and GLONASS
range error budgets are assessed, in order to complete the observation stochastic
model used in the testing to assess the impact of GLONASS on positioning accu-
racy. The following performance metrics are evaluated: accuracy and convergence
time for a range of static sites, and then accuracy, convergence time and continuity
for a variety of kinematic positioning use cases; these were chosen to reflect some of
the most likely applications and environments of PPP use.
For the assessment of the range error budget, it is desirable to remove as many of
the estimated parameters as possible and replace them with accurate truth values;
this is an attempt to measure the actual error budget, rather than the residual error
after positioning estimation. Without this approach, the range error can bias the
estimated parameters, reducing the observation range residuals, which could cause
the error budget to be underestimated. On this basis, IGS static sites are used for
the estimation of the error budget, as they have known coordinates with a high level
of accuracy; these coordinates can be held fixed when calculating the residuals. The
other benefit of the IGS static sites is that they operate in a clear sky environment,
so can measure down to low elevations and at all azimuth angles, and therefore
have a low level of cycle slips. This is important as any undetected cycle slip will
bias the calculated residuals. This static site case would be expected to provide
a lower bound on the error budget; for a kinematic receiver, one would expect to
see increased multipath and higher receiver measurement noise, caused by increased
receiver dynamics.
The ideal test design would cover a geographically diverse range of areas, in order
to see all satellites with a range of elevation and azimuth angles, the corresponding
range of satellite nadir angles. This also allows the range errors, due to spatially
varying unmodelled error sources, to be sampled from a range of locations. In this
Chapter, 87 static sites spread across the land areas of the Earth were investigated.
Ideally, the tests would be done over long periods, as some of the error sources are
temporally variable, this also allows understanding of seasonal effects. At a mini-
mum, datasets are required to cover a long enough period to show the convergence
of the solution. The static datasets were evaluated for one week based on availabil-
ity. To aid the statistical significance of any conclusions drawn, it is desirable to
114
Chapter 7. Combined GPS/GLONASS PPP
have multiple datasets. There are constraints in collecting kinematic PPP datasets,
and limited published data available; in this thesis three kinematic datasets were
collected from a range of vehicles.
The first kinematic dataset was collected by a pickup truck driving on public roads
in the UK. This dataset was designed to represent the challenges posed by buildings
and other sky obstructions, which result in a lot of multipath, providing a stressed
case on which to assess any improvement in continuity of the solution with the
addition of GLONASS.
The next dataset considered was collected from an offshore survey vessel operating
off the North coast of Norway. This dataset was chosen as it represents a real use
case of PPP, in a commercial environment where there is a real desire to address
the issues of convergence time and continuity (Veripos, personal communication
July 2010). Additionally, the test was performed at a high latitude, where better
geometry would be expected from the GLONASS satellites due to their higher orbital
inclination angle. This dataset was the only dual frequency GPS/GLONASS data
available from Veripos.
The investigation of these first two datasets suffered from the lack of an accurate
reference solution. The final dataset chosen had a GNSS reference station within
2.5 km throughout. It was collected on board a small survey vessel operating in
inshore waters on the East coast of England.
In order to quantify the relative accuracy of the GLONASS and GPS carrier phase
observations, the observation residuals can be compared. Observation data for a
range of static sites was downloaded for GPS week 1605 and run through pppncl with
the position not estimated but held fixed at the value provided by the ESA weekly
solutions. There were 87 sites in the IGS network recording GPS and GLONASS
observations for week 1605.
115
Chapter 7. Combined GPS/GLONASS PPP
Figure 7.1 shows the resulting ionosphere free LC carrier phase residuals for GPS
and GLONASS respectively. The colouring indicates the logarithm of the proportion
of residuals occurring at the given elevation angle for a residual of the corresponding
size. The complete set of residuals from all 87 sites and all satellites observed during
the week are included in the plots. For both GPS and GLONASS, the majority of
all residuals fall in a narrow band centred about zero. Both satellite systems exhibit
an increase in observation noise at lower elevation angles, with GLONASS showing
a slightly more marked increase.
The null hypothesis was stated as the mean of the difference between the RMS of
the GPS and GLONASS ionosphere free carrier phase residuals being zero. Based
on the apparent different behaviour of GLONASS at elevations below 40◦ , the RMS
of the residuals is considered in two groups, those above 40◦ , and those below. A
two tailed paired t-test was used to assess the null hypothesis. The p-value for
elevation angles below 40◦ was 0.11, and the p-value for elevation angles above
40◦ was 0.42. Therefore the null hypothesis is accepted, indicating there is not
a significant difference in the mean difference between the GPS and GLONASS
ionosphere free carrier phase residuals, in either of the elevation groups considered.
An elevation dependent observation noise function was fitted to these RMS values.
Two distinct functions were fitted to match the features seen at different elevation
angles:
0.0063
E < 25◦
sin E
s(E) =
(7.1)
0.017 − 7.5 × 10−5
×E E > 25◦
This observation noise function was then used to generate an elevation dependent
weighting function:
1
w(E) = (7.2)
s(E)2
In the following sets of PPP coordinate estimates generated using pppncl, the car-
rier phase observations were weighted using the elevation dependent function in
Equation 7.2.
116
Chapter 7. Combined GPS/GLONASS PPP
−6.0
0.2
−7.5
LC residuals (m)
0.1
Log proportion
−9.0
0.0 −10.5
−0.1 −12.0
−13.5
−0.2
−15.0
0 10 20 30 40 50 60 70 80 90
Elevation angle (◦ )
(a) GPS
−6.0
0.2
−7.5
LC residuals (m)
0.1
Log proportion
−9.0
0.0 −10.5
−0.1 −12.0
−13.5
−0.2
−15.0
0 10 20 30 40 50 60 70 80 90
Elevation angle (◦ )
(b) GLONASS
Figure 7.1: Ionosphere free carrier phase residuals for 87 IGS stations
over 5 days using ESA orbit and clock products. The colouring represents
the normalised log frequency of residuals occurring at a given value.
117
Chapter 7. Combined GPS/GLONASS PPP
0.040
GPS
0.035
GLONASS
0.030
LC residual RMS (m)
0.025
0.020
0.015
0.010
0.005
0.000
10 20 30 40 50 60 70 80 90
Elevation angle (◦ )
dinates are fixed whilst determining the observation noise using PPP, the receiver
clock offset and the additional zenith tropospheric delay must still be estimated.
How the observations are weighted in the PPP that generates the residuals will
affect the values of those residuals. For example under-weighting the GLONASS
observations will cause them to under-influence the estimated parameters. This
then causes GPS observation noise, through incorrect estimation of the unknown
parameters, to appear in the GLONASS residuals.
Despite its limitations this simple method provides a starting point from which to
investigate the impact of GLONASS on PPP.
In comparison with the a priori error budget given in Section 3.16, the residuals
determined in this test are considerably smaller. At higher elevation angles (greater
than 50◦ ) the a priori error budget of 8.3 cm overbounds all the observed residuals.
This is likely to be due to the estimated parameters absorbing some of the range er-
ror. In particular, the high correlation between the troposphere delay, receiver clock
and height component may contribute to this (Witchayangkoon, 2000). Constant
biases are also absorbed in the estimated parameters as they are not observable and
so cannot be separated from the carrier phase ambiguity.
118
Chapter 7. Combined GPS/GLONASS PPP
The Veripos orbit and clock products were provided in SP3 format tabulated at 5 s
intervals. This is at a higher rate than the ESA and IGS orbit and clock products
which are tabulated at 15 min intervals for the orbits, and 30 s intervals for the clock
products.
The convergence time has been identified as a limitation of current PPP, in both
the float and ambiguity fixed methods (Section 2.5). The addition of GLONASS
satellites to GPS PPP stands to improve the DOP and has therefore been proposed
as a method to reduce the convergence time (Section 2.9).
Perform static PPP using pppncl with as many GPS/GLONASS sites as possible
for the 6 days over the period 10 to 15 October 2010, for which the Veripos Apex
G2 orbits were available. The RMS of the PPP estimated coordinates, with respect
to the weekly coordinate solutions, describes the repeatability of static PPP at that
site. To compare the impact of the different orbital products, static PPP coordinate
estimation was performed with IGS final, ESA final and Veripos Apex G2 real-time
orbit and clock products. For the ESA and Veripos products two processing runs
were performed, once using only the GPS satellites and once using both GPS and
GLONASS. The IGS clock products do not contain GLONASS satellite clock offset
corrections, so could only be used to perform GPS PPP, but were included to act
as the “best currently available” GPS only PPP benchmark.
To assess the influence of GLONASS on the convergence time for each of the dif-
119
Chapter 7. Combined GPS/GLONASS PPP
ferent orbit and clock products, the RMS of the coordinate errors was calculated
at each of the 30 s observation intervals for each of the orbit and GPS/GLONASS
combinations.
The observation data from all the IGS sites recording GPS and GLONASS observa-
tion data were downloaded from CDDIS (Noll, 2010). This gave a total of 87 sites
with a global distribution for the 6 day period (Figure 7.3). Observation data was
processed for 1 h, 2 h, 4 h, 8 h, 12 h and 24 h periods to quantify the improvement in
repeatability with increasing observation length. GPS only solutions were produced
using IGS, ESA and Veripos Apex G2 orbit and clock products.
Using the same processing strategy as described in Section 6.1, coordinates were
produced for all of the sites, for each orbit provider, and for all seven observation
lengths. For the ESA and Veripos orbit and clock products a second set of processing
was performed with GLONASS observations enabled. The GLONASS carrier phase
observation noise was set equal to the GPS carrier phase observation noise using
the values found in Section 7.1. All other parameters were kept the same as for the
GPS only processing.
A total of 133,110 runs of pppncl were required to process all the combinations of
orbit/clock provider, site and satellite system used. The resulting coordinates for
each site were then compared to the IGS weekly solution and then a per site coor-
dinate RMS was calculated for each 30 s position estimation interval. Figures 7.4,
7.5 and 7.6 show the mean RMS of the PPP estimated coordinates against time for
120
Chapter 7. Combined GPS/GLONASS PPP
1.0
ESA GPS
ESA GG
IGS GPS
VERIPOS GPS
North RMS (m)
VERIPOS GG
0.1
0.01
0 20 40 60 80 100 120
Processing length (minutes)
As would be expected from the reduced accuracy possible in real-time orbit and
clock product generation, the accuracy obtained with the Veripos Apex G2 orbit
and clock product is lower than that obtained with the corresponding ESA orbit
and clock products.
Table 7.1 gives the mean and standard deviation for the time until the RMS of the
site coordinate estimation error is less than 10 cm for the set of sites used in this test.
This table shows the reduction in convergence time due to the addition of GLONASS
was between 17% and 47%. A one-tailed paired t-test was performed on the null
hypothesis that the mean RMS convergence time of the GPS and GPS/GLONASS
121
Chapter 7. Combined GPS/GLONASS PPP
1.0
ESA GPS
ESA GG
IGS GPS
VERIPOS GPS
East RMS (m)
VERIPOS GG
0.1
0.01
0 20 40 60 80 100 120
Processing length (minutes)
Figure 7.5: Mean RMS of static PPP estimated East coordinates with
respect to IGS reference coordinates.
1.0
ESA GPS
ESA GG
IGS GPS
VERIPOS GPS
Height RMS (m)
VERIPOS GG
0.1
0.01
0 20 40 60 80 100 120
Processing length (minutes)
122
Chapter 7. Combined GPS/GLONASS PPP
Provider p value
ESA final N 9.5E-03
E 3.2E-11
H 4.4E-10
Veripos N 2.2E-07
E 2.5E-24
H 8.9E-15
Table 7.2: GPS and GPS/GLONASS one tailed paired t-test p values
for the comparison of PPP convergence time
PPP methods, for a given orbit and clock provider, were the same. Here the one
tailed t-test was selected, as the alternative hypothesis is that the convergence time
with GLONASS is lower than without. A paired t-test was used as both sets of
data sample the same sites, with the measurement being the time until the RMS
coordinate error was less than 10 cm. For both the ESA final orbits and the Veripos
real-time orbits, in each of the coordinate components (North, East and Height),
there was a statistically significant difference in the convergence time between GPS
and GPS/GLONASS PPP (Table 7.2). These results suggest that the addition of
GLONASS to PPP can reduce the convergence time to reach decimetre accuracy.
Martín et al. (2011) investigated the change in convergence time of eight sites when
including GLONASS observations using MagicGNSS PPP (Píriz et al., 2008). For
the day of 1 Jan 2010 they found no significant reduction in convergence time.
Between 1 Jan 2010 and the start of the test performed here (10 Oct 2010), nine
new GLONASS satellites were launched and two satellites were withdrawn from
service. The results presented in this thesis show that using pppncl with the current
GLONASS constellation, a clear reduction in convergence time is achieved with
GPS/GLONASS PPP over a much larger set of sites, in contrast to the findings of
123
Chapter 7. Combined GPS/GLONASS PPP
Table 7.3: Test statistics for the comparison of GPS and GPS/-
GLONASS PPP 3D coordinate RMS after 24 h
The accuracy of PPP estimated coordinates static coordinates after longer estima-
tion periods of up to 24 h is show in Figures 7.7, 7.8 and 7.8. At longer processing
lengths the difference between the RMS of the GPS and GPS/GLONASS (denoted
as GG in the figures) PPP solutions reduces. Comparing the 3D coordinate RMS
error between GPS only and GPS/GLONASS PPP, the null hypothesis that the
coordinate RMS between for the two methods after 24 h for a given orbit and clock
provider come from the same distribution. A Wilcoxon signed-rank test (Lowry,
2011) was performed for the comparison for both ESA and Veripos orbit and clock
products. From the values given in Table 7.3, for p<0.05 we accept the null hypoth-
esis for the ESA orbit and clock products and reject it for the Veripos products. We
therefore conclude that for the ESA Final orbit and clock product GPS/GLONASS
PPP does not produces statistically different coordinate RMS errors than GPS PPP.
For the Veripos orbits we conclude that there is a statistically significant difference,
with the median coordinate RMS error for GPS/GLONASS of 13.7 mm being larger
than that measured with GPS PPP.
This would imply that once the Kalman filter has converged and the carrier phase
ambiguities are accurately determined, the limiting factor in improved RMS is not
the number and geometry of the visible satellites, but either the accuracy of the
orbit and clock products or the models used in calculating the signal delay.
124
Chapter 7. Combined GPS/GLONASS PPP
0.07
ESA GPS
0.06 ESA GPS/GLONASS
IGS GPS
0.05
Veripos GPS
North RMS (m)
0.04
Veripos GPS/GLONASS
0.03
0.02
0.01
0.00
1 2 4 8 16 32
Processing length (hours)
Figure 7.7: RMS of the static PPP estimated North coordinates with
respect to IGS reference coordinates.
0.10
ESA GPS
ESA GPS/GLONASS
0.08
IGS GPS
Veripos GPS
East RMS (m)
0.04
0.02
0.00
1 2 4 8 16 32
Processing length (hours)
Figure 7.8: RMS of the static PPP estimated East coordinates with
respect to IGS reference coordinates.
125
Chapter 7. Combined GPS/GLONASS PPP
0.18
ESA GPS
0.16
ESA GPS/GLONASS
0.14 IGS GPS
Veripos GPS
Height RMS (m)
0.12
Veripos GPS/GLONASS
0.10
0.08
0.06
0.04
0.02
0.00
1 2 4 8 16 32
Processing length (hours)
Figure 7.9: RMS of the static PPP estimated Height coordinates with
respect to IGS reference coordinates.
The method used for cycle slip detection (Section 4.8.2) for each satellite is inde-
pendent of the other visible satellites as it only depends on linear combinations
of observations from that satellite. Therefore the addition of GLONASS does not
impact operation of the cycle slip detection step.
The fixing of cycle slips makes use of the carrier phase observations free from cycle
slips to estimate the cycle slips prior to integer fixing using the LAMBDA method.
The inclusion of GLONASS in addition to GPS provides more observations; if there
are cycle slips on the carrier phase measurements for several satellites, then with
more satellites, there are more likely to be sufficient unaffected satellites. This
enables the cycle slips to be resolved and therefore corrected, providing continuity
without requiring reinitialisation. An example of an extreme case of this where the
receiver loses lock on all GPS satellites is described in the inshore survey vessel
dataset in Section 7.5.3.
The method used for outlier detection is described in Section 4.8.1. The addition of
GLONASS satellites affects the outliers in two different ways, the ability to detect
outliers is impacted, as well the positioning accuracy in the case of undetected
126
Chapter 7. Combined GPS/GLONASS PPP
outliers. The statistic used in the initial overall model test is derived from all of the
observations made at a particular epoch.
This test was completed by considering the IGS station, HERS, a dual frequency
GPS/GLONASS receiver logging data at 1 Hz. The same processing strategy was
used as described in Section 6.1, using the ESA orbit and clock products.
To test the addition of the step type to the carrier phase, a single cycle was added to
the L1 and L2 carrier phase measurements of GPS satellite ID 3 for the period from
30 min after the start of the dataset until 60 min after the start. This combination of
slips in L1 and L2 was chosen as the change is not measurable in the wide lane linear
combination, making it harder to detect. Both the GPS and the GPS/GLONASS
processes detected the cycle slips immediately, both at addition and subtraction of
the step bias, and both were able to successfully identify the correct cycle slip. This
step resulted in no impact to the positioning solution.
To test the addition of the step type in code phase, for the same time period as the
above test, a 25 m bias was added to the C/A code phase measurement. In the GPS
only solution this bias caused the overall model’s test statistic to be larger than the
threshold value. The ionosphere free code phase linear combination was identified as
an outlier and removed from the estimation process for the entire period of the bias
addition. The graph of the test statistic and the thresholds is shown in Figure 7.10,
the effect on positioning performance is shown in Figure 7.11. For the GPS only
solution, the 25 m bias was detected and removed for the entire period, it therefore
did not impact the position apart from the effect of having one less code phase
measurement in use. For the GPS/GLONASS PPP solution, this bias did not cause
the overall model test statistic to exceed the threshold until 2887 s (approximately
16 min after the addition of the bias), the test statistic was not consistently above
the threshold until 3189 s. After this time it was identified and removed from the
estimation process. The late detection was due to the increased threshold value
calculated in the outlier detection process.
Despite the failure of the GPS/GLONASS PPP to remove the step bias for the
127
Chapter 7. Combined GPS/GLONASS PPP
majority of the addition period the 3D coordinate error remained lower than that
of the GPS/GLONASS solution at all times. Also visible in Figure 7.10 is the large
increase in the test statistic for both methods after the removal of the step bias at
3600 s. This is erroneously detected as a cycle slip using the Melbourne-Wübenna
method, the slip resolution method then “successfully” fixes the apparent cycle slip
to the wrong value. This is then detected by the outlier detection and the ambiguity
is reset.
100
GPS test statistic
GPS threshold
80
GG test statistic
GG threshold
60
40
20
0
0 1000 2000 3000 4000 5000 6000 7000 8000
Time (s)
Figure 7.10: Overall model test statistics for GPS and GPS/-
GLONASS PPP with code phase step bias, vertical black lines show
start and end of bias addition at 1800 s and 3600 s respectively
To test the effect of a slowly increasing bias in the carrier phase, a ramp bias was
added increasing at 1 cm every second. In the GPS only solution this bias was
initially detected after 12 s, and in the GPS/GLONASS PPP solution, after 13 s
(Figure 7.12). The carrier phase observation is then excluded from the solution, but
not reset for up to 5 epochs. This approach to carrier phase outlier detection is an
attempt to avoid unnecessarily resetting ambiguities due to any short term outliers.
In both methods, the carrier phase ambiguity is reinitialised for this satellite, but
as the bias continues to grow this triggers repeated outlier detections, causing the
carrier phase bias to be repeatedly reset. The repeated outlier detections are shown
in Figure 7.12 by the test statistic repeatedly crossing the threshold. The reset
interval for GPS is lower, at an average of 21 s, than for GPS/GLONASS with
an average of 28 s. The GPS only solution is able to detect the bias at a smaller
128
Chapter 7. Combined GPS/GLONASS PPP
0.25
GPS
GPS+GLONASS
0.20
3D position error (m)
0.15
0.10
0.05
0.00
0 1000 2000 3000 4000 5000 6000 7000 8000
Time (s)
Figure 7.11: 3D RMS coordinate error for GPS and GPS/GLONASS
PPP with code phase step bias, vertical black lines show start and end
of bias addition
magnitude due to the lower threshold as for the code phase step bias detection.
Despite the earlier detection in GPS only PPP, there is a clear reduction in the po-
sition accuracy in the GPS only solution at the start of the bias introduction (Fig-
ure 7.13). With fewer observations the GPS only PPP solution is negatively influ-
enced by the initial period when the bias is undetected, whereas the GPS/GLONASS
solution does not suffer from an obvious discontinuity in coordinate accuracy.
To test the addition of the ramp bias to the code phase, a ramp bias was added at
time 1800 s, increasing at 1 m every second. This causes the test statistic to slowly
grow as the bias grows (Figure 7.14). In the GPS only solution this bias was detected
after 23 s, and in the GPS/GLONASS PPP solution, after 28 s. Unlike in the case
for the carrier phase ramp bias, once the bias grows above the detection threshold,
it is then detected for the rest of the period it is added. The low weight of the
code phase measurement, due to the expected multipath and receiver measurement
noise contribution to the error budget, cause the effect of the undetected period
of bias to be fairly small. Figure 7.15 shows the impact on 3D positioning accu-
racy with reference to the IGS weekly coordinates, there is no noticeable decrease
in the GPS/GLONASS coordinate accuracy during the 28 s for which the bias is
undetected.
129
Chapter 7. Combined GPS/GLONASS PPP
100
GPS test statistic
GPS threshold
80
GG test statistic
GG threshold
60
40
20
0
1750 1800 1850 1900 1950 2000 2050 2100
Time (s)
Figure 7.12: Overall model test statistics for GPS and GPS/-
GLONASS PPP with a carrier ramp bias, vertical black lines show start
and end of bias addition
0.25
GPS
GPS+GLONASS
0.20
3D position error (m)
0.15
0.10
0.05
0.00
0 1000 2000 3000 4000 5000 6000 7000 8000
Time (s)
Figure 7.13: 3D RMS coordinate error for GPS and GPS/GLONASS
PPP with a carrier ramp bias, vertical black lines show start and end of
bias addition
130
Chapter 7. Combined GPS/GLONASS PPP
100
GPS test statistic
GPS threshold
80
GG test statistic
GG threshold
60
40
20
0
1700 1750 1800 1850 1900
Time (s)
Figure 7.14: Overall model test statistics for GPS and GPS/-
GLONASS PPP with a code phase ramp bias, vertical black lines show
start and end of bias addition
0.25
GPS
GPS+GLONASS
0.20
3D position error (m)
0.15
0.10
0.05
0.00
0 1000 2000 3000 4000 5000 6000 7000 8000
Time (s)
Figure 7.15: 3D RMS coordinate error for GPS and GPS/GLONASS
PPP with a code phase ramp bias, vertical black lines show start and
end of bias addition
131
Chapter 7. Combined GPS/GLONASS PPP
To summarise this outlier detection test, the GPS only solution appears to be able
to detect both smaller code phase biases and smaller carrier phase biases. In the
GPS/GLONASS solution, the additional degrees of freedom, due to the additional
observations on which outliers could occur, increases the threshold at a given signif-
icance, therefore increasing the minimal detectable bias.
To assess the impact of using GLONASS combined with GPS on kinematic PPP,
a variety of kinematic datasets were collected using a range of transport modes as
described below.
This test involved driving a vehicle on public roads. A Leica GS10 receiver was
mounted with an AS10 antenna on a pickup truck above the main cab to give the
best clear sky view possible. The vehicle was driven for approximately 6 h on 20 Jan
2011 from Newcastle-upon-Tyne to Carlisle around the Lake District area and to
Southwest Scotland, in the UK. The route passed through villages and hilly areas,
representative of the UK countryside. As well as driving at a range of speeds, there
were three periods during the journey when the vehicle was stationary. Data was
collected at 1 Hz, recording all GPS/GLONASS satellites in view. Observation data
from static sites along the vehicle’s route that form part of the Leica Smartnet
Network RTK service were downloaded via the Leica Smartnet website (Burbidge,
2006). For subscribed users of the Smartnet service it is possible to access the 1 Hz
observation data from the reference network, with the past 30 days available for
download. This provided a set of reference stations to use in estimating a set of
“truth” coordinates using the double-differenced post-process carrier phase method.
132
Chapter 7. Combined GPS/GLONASS PPP
MIT Track was used to produce a set of reference coordinates with which to compare
the GPS and GPS/GLONASS PPP solutions. The route passed near three base
stations: CARL, KIRK, DRUM, as shown in Figure 7.16. As these base stations are
run by the Ordnance Survey (Cruddace, 2006), the known coordinates are given in
ETRS89 and are therefore not compatible with the ESA GPS/GLONASS orbit and
clock products which are given in ITRF2005. For each of the three base stations,
Dr. Nigel Penna generated GIPSY PPP-based daily static coordinates using the
ESA Final orbits. The weighted mean coordinates from each day of the surrounding
week (composed of seven 24 h periods) were used as the reference locations. Table
7.4 shows the mean coordinates produced using GIPSY and the standard error of
the mean. The sub-millimetre standard errors in all coordinate components are well
below the achievable centimetre-level kinematic PPP accuracy, so are taken as truth
in the following analysis.
All three reference stations were included in the Track coordinate estimation mean-
ing multiple baselines can be formed. This means that the shortest baseline available
is always included in generating the reference coordinates, to avoid having to process
in multiple short discontinuous batches. The reference stations were held fixed at
the calculated coordinates. All three of the reference stations used are equipped with
the LEIAR25 LEIT Leica antenna for which phase centre corrections are available in
the IGS ANTEX file ([Link]).
All the receivers observe C1 and P2 code observations so the monthly estimated
P1-C1 biases from CODE were used to correct the C1 observable to a P1 equivalent
observable.
Using the following configuration file, double difference coordinates for the receiver
mounted on the pickup truck were produced using Track:
OBS_FILE
carl data/carl020k00_p1.11o F
kirk data/kirk020k00_p1.11o K
133
Chapter 7. Combined GPS/GLONASS PPP
drum data/drum020k00_p1.11o K
rove data/leicapickup_p1.11o K
SITE_POS
carl 3671344.2931 -188441.0622 5194774.1971
drum 3601680.6339 -235186.8396 5241367.8543
kirk 3672068.2812 -259826.9475 5191172.8759
rove 3693831 -176914 5179434
MODE AIR
SITE_STATS
all 0.005 0.006 0.007 0.0 0.0 0.0
rove 200 200 200 4.7 4.7 4.7
ATM_STATS
all 0.1 0.0001
ante_off
all 0.0 0.0 0.0 LEIAR25 LEIT
rove 0.0 0.0 0.0 LEIAS10 NONE
ANTMOD_FILE [Link]
CUT_OFF 10
OUT_TYPE GEOD+NEU+XYZ
INTERVAL 1
USE_GPTGMF
BACK_TYPE SMOOTH
The kinematic dataset was very challenging to process with regular complete loss
of lock due to the vehicle passing under bridges and behind other obstructions such
as forests. Due to periods with a low number of visible satellites, high dilution of
precision and short continuous observation arcs, the reference solution has many
holes and periods where confidence in the result is low due to Track being unable
to fix the carrier phase ambiguities to integer values.
134
Chapter 7. Combined GPS/GLONASS PPP
Figure 7.16 shows the sections of the route that Track was able to process are
shown in relation to the base station locations. The gaps show periods for which
the coordinate uncertainties produced by Track were greater than a threshold of 1 m
for coordinate output. From Figure 7.17 it can be seen that the problematic periods
occur whilst the vehicle is moving. If the formal errors for the estimated coordinates
are taken into account, it can be seen that for all periods during which the vehicle
is in motion, the estimated coordinates are only accurate to 20 cm (one sigma) at
best.
−5˚ −4˚ −3˚ −2˚
DRUM
55˚ 55˚
KIRK CARL
135
Chapter 7. Combined GPS/GLONASS PPP
45
Distance from nearest base station (km) 40
35
30
25
20
15
10
0
11:00:00 12:00:00 13:00:00 14:00:00 15:00:00 16:00:00
Time (GPST)
Figure 7.17: Distance of the pickup truck from the nearest base station
1.0
North σ (m)
0.8
0.6
0.4
0.2
0.0
1.0
East σ (m)
0.8
0.6
0.4
0.2
0.0
1.0
Height σ (m)
0.8
0.6
0.4
0.2
0.0
11:00:00 12:00:00 13:00:00 14:00:00 15:00:00 16:00:00
Time
Figure 7.18: Pickup truck test Track-estimated reference coordinate
uncertainty
136
Chapter 7. Combined GPS/GLONASS PPP
fields.
The pickup truck data was processed with pppncl twice, once using only GPS ob-
servations, and once using both GPS and GLONASS observations. The following
processing strategy was used:
The estimated position solutions were then differenced from the reference solution
to produce a time series.
As for the reference solution, the PPP solution also struggled to converge (Fig-
ures 7.19 and 7.20) due to the repeated complete loss of lock which forced reinitiali-
sation of the carrier phase ambiguity estimates. The only periods where continuous
satellite tracking was possible coincided with periods when the vehicle was station-
ary.
As shown in Figure 7.19, during the three periods from 10:00–11:00, 11:45–12:40 and
14:30–15:45, when the vehicle was stationary, the PPP solution converges towards
the reference solution. In the remaining periods, the receiver did not maintain lock
on the satellites for long enough for the carrier phase ambiguities to be sufficiently
well determined in the Kalman filter for accurate positioning. Figure 7.18 shows
that the accuracy of the reference solution obtained in this period was worse than
40 cm, so is anyway not useful as truth coordinates.
After the receiver resumes tracking satellites after loss of lock, pppncl attempts to
determine the introduced cycle slip in the carrier phase observations in order to avoid
137
Chapter 7. Combined GPS/GLONASS PPP
0.2
0.0
−0.2
−0.4
11:00 12:00 13:00 14:00 15:00 16:00
Height (m)
0.4
0.2
0.0
−0.2
−0.4
11:00 12:00 13:00 14:00 15:00 16:00
0.4
0.2
0.0
−0.2
−0.4
11:00 12:00 13:00 14:00 15:00 16:00
0.4
East (m)
0.2
0.0
−0.2
−0.4
11:00 12:00 13:00 14:00 15:00 16:00
Height (m)
0.4
0.2
0.0
−0.2
−0.4
11:00 12:00 13:00 14:00 15:00 16:00
138
Chapter 7. Combined GPS/GLONASS PPP
reinitialising the estimates of the carrier phase ambiguities. If there has been total
loss of lock, the cycle slip fixing method relies on the code observations to determine
the discontinuity in the carrier phases. In many cases, cycle slip fixing is not possible
after the complete loss of lock due to the low accuracy of the pseudorange during the
initial period after the receiver resumes tracking of the satellites. This can be seen
in the Melbourne-Wübbena combination of code and carrier phase (Equation 4.52)
which, being invariant of the receiver satellite geometry, only contains the wide-lane
carrier phase ambiguity and any measurement noise in the signals. Excluding the
measurement noise, the Melbourne-Wübbena combination should be constant. The
high variability seen during epochs 10–30 in Figure 7.21, with the nominally constant
wide-lane bias changing by approximately 1.6 m, shows the poor accuracy of the
code and carrier tracking. This is visible in the Melbourne-Wübbena combination
for GPS satellite 2 after the receiver resumes tracking the satellite following total
loss of lock.
24
Melbourne-Wubbena
Widelane wavelengths (∼86.2 cm)
Running average
22 Detection thresold
20
18
16
0 20 40 60 80 100 120
Observation epoch +3.488129×108
139
Chapter 7. Combined GPS/GLONASS PPP
of pppncl during the period 13:00-14:30. This dataset highlights the importance of
a short convergence time in challenging environments where loss of lock is common.
An attempt was made to process this dataset via alternative PPP programs using
the PPP Software Centre (Banville et al., 2009). This was unsuccessful with none
of the online PPP services returning a solution, reinforcing that the issue seen with
this case study was the poor quality of data rather than the implementation of the
PPP method in pppncl.
The provided dataset contains observations recorded at 1 Hz. The only nearby
reference station which archives data at a rate of 1 Hz or higher is NYA2 in Svalbard.
There are no other stations in the region providing data at the required 1 Hz rate
(NYA1 in Svalbard, TRO1 in Tromso, Norway and VAR2 in Northeast Norway all
archive data at 30 s observation intervals).
NYA2 is a proposed IGS site, and therefore has high quality monumentation; it
can thus be assumed that the receiver and installation are suitable for use as a
reference station. However, as can be seen from Figure 7.23 the ship’s route ranged
from between 70 kilometres to over 1000 kilometres from the base station. This is
a considerably larger distance than the ∼10 kilometres for which the single baseline
140
Chapter 7. Combined GPS/GLONASS PPP
141
Chapter 7. Combined GPS/GLONASS PPP
1400
Distance from base station (km)
1200
1000
800
600
400
200
0
13 14 15 16 17
Jun
2010
Date
relative carrier phase positioning would be expected to provide high quality results.
However, given the remote location, and the limited sea based datasets available,
this non ideal base station was used to generate the reference solution to allow some
comparison to be made.
Precise coordinates for NYA2 are not available, however it is only 167 m North West
of the IGS station NYA1. A precise weekly coordinate solution is available for NYA1
from ESA in the same, ITRF2005, reference frame as the ESA orbit and clock solu-
tions to be used in the PPP positioning. Using Track in short baseline static mode,
NYA2 was positioned relative to NYA1, providing coordinates for NYA2 given in
Table 7.5. Due to the short base line, the L1 and L2 carrier phase observations could
be used directly resulting in formal errors of less than a millimetre for the estimated
NYA1-NYA2 baseline. The coordinate uncertainty for NYA2 was calculated using
variance propagation of the combined uncertainty in the coordinates of NYA1 and
the baseline uncertainty. The resulting estimated coordinates for NYA2 have 1 mm
formal errors in the North and East directions.
A reference solution was created using Track with NYA2 as the reference station
using the previously estimated coordinates. The achieved accuracy of the reference
solution is shown in Figure 7.24. The average formal error of the reference solution
142
Chapter 7. Combined GPS/GLONASS PPP
NYA1 NYA2
Coordinate Position (m) σ (mm) Position (m) σ (mm)
x 1202433.7652 1.2 1202379.4810 1.3
y 252632.3309 1.1 252474.5700 1.2
z 6237772.6440 5.2 6237786.3375 5.4
0.2
North (m)
0.1
0.0
0.2 13 14 15 16 17
Jun
East (m)
2010
0.1
0.0
0.2 13 14 15 16 17
Height (m)
Jun
2010
0.1
0.0
13 14 15 16 17
Jun
2010
Date
The entire five day period was processed using pppncl in one continuous session.
ESA final orbit and clock products were used due to the availability of precise
GLONASS satellite clock corrections. The orbit and clock products for the obser-
vation period were concatenated to provide one file of continuous corrections. The
143
Chapter 7. Combined GPS/GLONASS PPP
• Processed once using only GPS and again with GPS and GLONASS
The resulting position solutions were then differenced from the reference solution
coordinates to produce a time series for both GPS and GPS/GLONASS PPP.
Figures 7.25 and 7.26 show the coordinate errors in the North, East and Height com-
ponents with respect to the reference solution for GPS PPP and GPS/GLONASS
PPP solutions respectively. In both solutions, for GPS and GPS/GLONASS com-
bined, the period from 16 June 2010 shows considerably worse accuracy with respect
to the reference coordinates than the initial period from 13 to 15 June 2010.
Looking at the distance of the ship from the base station during the period con-
sidered (Figure 7.23), it is clear that the apparent period of poor accuracy from 16
June correlates with the period when the ship was a long distance from the base
station. Despite the consistency of the formal errors for the reference coordinates
reported by Track (Figure 7.24), it would seem that the true accuracy of the ref-
erence coordinates visibly deteriorates, as would be expected over such long base
line lengths (greater than 500 kilometres). Therefore to ensure sufficient accuracy of
the reference solution, the following analysis is restricted to the period during which
the ship was less than 110 kilometres from the reference station. This resulted in
approximately two days from the overall dataset, for which usable reference coor-
dinates can be used to make a comparison of GPS versus GPS/GLONASS PPP
positioning.
Figures 7.27 and 7.28 show the GPS PPP and GPS/GLONASS PPP solutions re-
spectively by comparison with the reference solution, for the period 05:00 13 June to
24:00 14 June with reliable reference coordinates. Immediately obvious is the large
discontinuity at 10:00 on 13 June 2010 in the GPS solution which is absent in the
144
Chapter 7. Combined GPS/GLONASS PPP
0.4
0.2 2010
0.0
−0.2
−0.4
13 14 15 16 17
Height (m)
0.4 Jun
0.2 2010
0.0
−0.2
−0.4
13 14 15 16 17
Jun
2010
Figure 7.25: Difference between the GPS PPP coordinates and the
reference coordinates, for the offshore survey vessel. Red lines show the
one sigma uncertainty in the reference solution.
North (m)
0.4
0.2
0.0
−0.2
−0.4
13 14 15 16 17
Jun
East (m)
0.4
0.2 2010
0.0
−0.2
−0.4
13 14 15 16 17
Height (m)
0.4 Jun
0.2 2010
0.0
−0.2
−0.4
13 14 15 16 17
Jun
2010
145
Chapter 7. Combined GPS/GLONASS PPP
0.4
0.2 2010
0.0
−0.2
−0.4
06:00 12:00 18:00 00:00 06:00 12:00 18:00
Height (m)
0.4 14-Jun
0.2 2010
0.0
−0.2
−0.4
06:00 12:00 18:00 00:00 06:00 12:00 18:00
14-Jun
2010
Figure 7.27: Difference between the GPS PPP coordinates and the ref-
erence coordinates, for the offshore survey vessel when less than 110 km
from the base station.
North (m)
0.4
0.2
0.0
−0.2
−0.4
06:00 12:00 18:00 00:00 06:00 12:00 18:00
14-Jun
East (m)
0.4
0.2 2010
0.0
−0.2
−0.4
06:00 12:00 18:00 00:00 06:00 12:00 18:00
Height (m)
0.4 14-Jun
0.2 2010
0.0
−0.2
−0.4
06:00 12:00 18:00 00:00 06:00 12:00 18:00
14-Jun
2010
146
Chapter 7. Combined GPS/GLONASS PPP
GPS GPS/GLONASS
North 6.7 6.3
East 6.8 4.8
Height 18.1 12.9
GPS/GLONASS solution for the same period. This indicates that the additional
GLONASS satellites help to make the positioning method more robust. Addition-
ally, there are a number of jumps seen in the North and East GPS solutions which
are not present in the more robust GPS/GLONASS solution. These events are very
short (one to two epochs), and indicative of bad data below the MDB having a
larger impact on the estimated position in the GPS solution.
The reported formal errors in the Height component of the reference coordinates are
5 cm to 15 cm (Figure 7.24). This leads to apparently noisy PPP Height solutions,
and as such the difference between GPS and GPS/GLONASS PPP is less apparent.
Table 7.6 shows the standard deviations of the coordinate differences with respect
to the reference coordinates of the two solutions. There is a clear improvement
with the addition of GLONASS in the East and Height components. One of the
major factors in this reduction is the difference between the two solutions at 10:00
13 June. In the GPS only pppncl solution there is a clear discontinuity at this time,
followed by a period of re-convergence that, in the case of the Height component,
lasts approximately 1 h. The GPS/GLONASS solution maintains enough continuous
observations over this period to avoid reinitialisation, highlighting the increased
reliability of the dual system solution.
Comparing the Height component of the GPS and GPS/GLONASS PPP solutions
(Figure 7.29) shows the addition of the extra GLONASS observations has a mea-
surable effect on the coordinates. The RMS of the difference over the whole period
is 26 cm. Referring to the difference with respect to the base station (Figures 7.25
147
Chapter 7. Combined GPS/GLONASS PPP
and 7.26) it appears that metre level difference between the two PPP solutions is
down to errors in the GPS only solution. However as the accuracy of the reference
coordinates is so low during this period this is not a conclusive finding.
Figure 7.30 shows the height of the GPS/GLONASS pppncl PPP solution compared
with the Veripos Ultra solution. For much of the central part of the dataset, from
approximately midday 13 June to midday 15 June, the Veripos solution is using
broadcast orbit and clock corrections. This was due to the high latitude of the vessel
which prevented the Veripos Ultra orbit and clock corrections being received from
the communications satellite. Unfortunately this period, for which the Ultra orbit
corrections were unavailable, coincides with the only usable part of the reference
solution (05:00 13 June to 24:00 14 June as mentioned earlier). Therefore with
148
Chapter 7. Combined GPS/GLONASS PPP
this dataset it is not possible to make a quantitative comparison of the pppncl and
Veripos Ultra solutions, as the true coordinates are not known.
Looking at Figure 7.30 for the period 15 to 16 June, what appears as 2 m level noise
reflects the real height change as the ship moves over the waves. Towards the end
of 16 June there is a spike in the Veripos Ultra Height solution followed by a clear
discontinuity. This feature is highly unlikely to be correct given the ship is floating
on the sea. The pppncl GPS/GLONASS solution does not show this discontinuity,
but the extent to which this is an improvement over the Veripos Ultra solution
unfortunately cannot be determined. With the absence of a high quality reference
solution there is little more that can be drawn from comparing these two datasets.
This test involved comparison of GPS and GPS/GLONASS PPP methods to po-
sition a moving vessel in an inshore sea environment. A Leica GS10 receiver with
an AS10 antenna was mounted on the foredeck of the Port of Tyne vessel, Lynceus,
during underwater surveying operations in the mouth of the River Tyne, North East
England. The River Tyne has high ground surrounding it, and is a working port
with large ships and oil rigs moored in the river, or transiting through. The land
149
Chapter 7. Combined GPS/GLONASS PPP
shape combined with other shipping can potentially lead to a noisy environment
from a data collection point of view.
Data was collected during 24 March 2011, from 09:12 until 14:50 (GPS time), includ-
ing a 10 min gap in data recording at 11:30 (GPS time) due to unplanned receiver
shutdown. During this period the boat made several trips up and down the river
and performed a detailed survey of the river entrance; the ground track followed is
shown in Figure 7.31. The boat was always less than 2.5 km from the base station
used (Figure 7.32); this is within the range over which a double differencing car-
rier phase solution using one base station would be expected to provide centimetre
accurate coordinates.
Figure 7.31: Ground track of the survey vessel Lynceus. The red
triangle marks the location of the GNSS site NSLG used as a base station
The ITRF2005 coordinates of the reference station NSLG were estimated relative to
the IGS station MORP. Using the IGS weekly coordinates for MORP to provide ac-
cess to the ITRF2005 reference frame, a static baseline GPS solution performed with
150
Chapter 7. Combined GPS/GLONASS PPP
2.5
1.5
1.0
0.5
0.0
10:00:00 11:00:00 12:00:00 13:00:00 14:00:00
Time (GPST)
Figure 7.32: Distance of the survey vessel Lynceus from the base sta-
tion NSLG
Track was used to estimate the coordinates of NSLG. Table 7.7 gives the coordinates
used for MORP and the estimated coordinates of NSLG. The one sigma uncertainty
for the NSLG coordinates was obtained by the method of variance propagation.
MORP NSLG
Coordinate Position (m) σ (mm) Position (m) σ (mm)
x 3645667.762 4.1 3664667.331 4.2
y -107277.138 1.5 -91646.772 1.6
z 5215053.596 4.8 5201998.175 4.9
Apart from two short satellite arcs (one at the start and the other at the end, of
151
Chapter 7. Combined GPS/GLONASS PPP
less than 40 s) all ambiguities were fixed to integer values. The formal errors of the
coordinates produced by Track are shown in Figure 7.33. Both the North and East
components have been estimated with reported formal errors of 1 cm to 2 cm, and
the Height with 3 cm to 4 cm formal errors.
0.10
North (m)
0.08
0.06
0.04
0.02
0.00
10:00 11:00 12:00 13:00 14:00
0.10
East (m)
0.08
0.06
0.04
0.02
0.00
10:00 11:00 12:00 13:00 14:00
0.10
Height (m)
0.08
0.06
0.04
0.02
0.00
10:00 11:00 12:00 13:00 14:00
Time (GPST)
Figure 7.33: Tyne survey vessel Track reference coordinate formal
errors
Again the data was processed using pppncl, once with only GPS observations and
once with GPS/GLONASS observations combined. All other aspects of the PPP
strategy for the two solutions were identical as follows:
152
Chapter 7. Combined GPS/GLONASS PPP
The position solutions were then differenced from the reference solution to produce
a time series of coordinate errors with respect to the reference solution.
Figure 7.34 shows the difference between the GPS only PPP estimated coordinates
and the reference coordinates. It can be seen that although there is slow initial
convergence, the GPS only PPP is capable of achieving coordinates accurate to
within the uncertainty in the reference solution (1 cm to 2 cm in the North and
East components, and 5 cm in the Height component). However, looking at both
the initial period, and after the receiver reset at 11:30 (GPS time), the convergence
time to the highest achieved accuracies is up to 90 min (Figure 7.36).
At 10:38 (GPS time) cycle slips are detected on 7 of the 10 visible GPS satellites
due to the receiver losing lock on the GPS carrier phase for 10 s. After the receiver
reestablished carrier phase tracking, the values of the cycle slips were not successfully
identified and corrected for using the method described in Section 4.8.2. The affected
satellites were GPS 17, 32, 31, 25, 13, 24 and 2.
After unsuccessful cycle slip fixing, the carrier phase ambiguities for these satellites
are reset and then re-estimated in pppncl. This causes the reduction in accuracy
visible in Figure 7.34 during the following 25 min period whilst the solution re-
converges.
Figure 7.35 shows the difference between the GPS/GLONASS PPP estimated coor-
dinates and the reference coordinates. Comparing Figures 7.34 and 7.35, it can be
seen that there is a clear reduction in convergence time when GPS and GLONASS
are combined. From Figures 7.36 and 7.37, it can be seen that the convergence
time with GLONASS is reduced to approximately half an hour, one hour less than
for the GPS only PPP. Here convergence is defined as the time until 2 cm RMS is
achieved in the North and East components and 4 cm RMS is achieved in the Height
component.
With the GPS/GLONASS combined observations, the issue seen in the GPS only
observations at 10:38 due to cycle slips is not present, as there are no cycle slips
on the GLONASS satellites at this time. This demonstrates the benefit of using
GLONASS in addition to GPS to provide a more robust positioning solution.
153
Chapter 7. Combined GPS/GLONASS PPP
North (m)
0.4
0.2
0.0
−0.2
−0.4
10:00 11:00 12:00 13:00 14:00
0.4
East (m)
0.2
0.0
−0.2
−0.4
10:00 11:00 12:00 13:00 14:00
Height (m)
0.4
0.2
0.0
−0.2
−0.4
10:00 11:00 12:00 13:00 14:00
0.4
0.2
0.0
−0.2
−0.4
10:00 11:00 12:00 13:00 14:00
0.4
East (m)
0.2
0.0
−0.2
−0.4
10:00 11:00 12:00 13:00 14:00
Height (m)
0.4
0.2
0.0
−0.2
−0.4
10:00 11:00 12:00 13:00 14:00
154
Chapter 7. Combined GPS/GLONASS PPP
0.05
0.00
−0.05
−0.10
10:00 11:00 12:00 13:00 14:00
0.10
Height (m)
0.05
0.00
−0.05
−0.10
10:00 11:00 12:00 13:00 14:00
Figure 7.36: Large scale difference between GPS PPP and reference
coordinates for the Tyne survey vessel. Red lines mark the one sigma
uncertainty in the reference solution.
0.10
North (m)
0.05
0.00
−0.05
−0.10
10:00 11:00 12:00 13:00 14:00
0.10
East (m)
0.05
0.00
−0.05
−0.10
10:00 11:00 12:00 13:00 14:00
0.10
Height (m)
0.05
0.00
−0.05
−0.10
10:00 11:00 12:00 13:00 14:00
155
Chapter 7. Combined GPS/GLONASS PPP
not bring higher ultimate accuracy than GPS only coordinate estimates, instead the
benefit is faster convergence time and a more robust positioning solution.
Table 7.8 shows the overall RMS values for the full time period considered. The
improvement in RMS seen is largely due to the improved convergence time of the
GPS/GLONASS solution.
Table 7.8: Tyne survey vessel PPP coordinate RMS with respect to
the reference coordinates for the complete survey period. Units for all
values are cm.
Considering the final 2 h period from 13:00 (GPS time) after both solutions have
converged, there is little difference in the accuracy of the two solutions (Table 7.9).
The error in both the PPP solutions is comparable to the uncertainty in the double
differencing reference coordinates so it is not possible to determine the underlying
accuracy of the PPP solutions.
Table 7.9: Tyne survey vessel converged PPP coordinate RMS with
respect to the reference coordinates for the final 2 h. Units for all values
are cm.
Figure 7.38 reinforces the above observations: the improvement brought about by
the use of GLONASS in addition to GPS is seen primarily during the initial conver-
gence period and in the case when there are problems with the cycle slips of GPS
satellites. Secondly, after convergence, the accuracy achieved by GPS/GLONASS
combined is comparable to that achieved by GPS only coordinate estimation.
156
Chapter 7. Combined GPS/GLONASS PPP
1.4
1.2
Reduction in error (m)
1.0
0.8
0.6
0.4
0.2
0.0
−0.2
10:00 11:00 12:00 13:00 14:00
Time (GPST)
Figure 7.38: Reduction in horizontal error of the GPS/GLONASS
PPP solution compared with the GPS only solution with respect to the
reference coordinates
ambiguities with a high degree of confidence. Figure 7.39 shows the estimated value
of carrier phase ambiguities varies by up to 2 m during this re-estimation period for
the cycle slip event at 10:38 (GPS time). There is some curvature in the value of
the estimated ambiguities; while this is happening the value of the carrier phase
ambiguity is incorrect, leading to the reduced accuracy seen in Figure 7.36 and
Figure 7.38.
With the eight additional GLONASS satellites visible in the GPS/GLONASS pro-
cessing, the GLONASS carrier phase observations, with their well determined am-
biguity estimates, provide precise observations allowing for accurate determination
of the receiver coordinates. With the coordinates well defined, the initial estimates
of the post cycle slip GPS carrier phase ambiguities in the Kalman filter are more
accurate. This can be seen in the reduced variability in the GPS carrier phases
after the cycle slip event (Figure 7.40). This is coupled with reduced a posteriori
variance of the GPS carrier phase ambiguities. This causes the GPS carrier phase
observations to carry more weight in the Kalman Gain matrix, thus contributing
more to the following position estimates than in the GPS only case. These effects
combined not only provide continuous high quality positioning solutions through
such an event, but also create a system that rapidly returns to full accuracy.
157
Chapter 7. Combined GPS/GLONASS PPP
Figure 7.39: GPS carrier phase ambiguity estimates after a loss of lock
event during GPS only PPP. The legend identifies the GPS satellites
involved.
Carrier phase ambiguity estimate (m)
G02
5
G04
G10
G13
0 G17
G20
G23
−5 G24
G25
G31
G32
−10
10:37:00 10:39:00 10:41:00 10:43:00 10:45:00
Time of day (GPST on 24 Mar 2011)
Figure 7.40: GPS carrier phase estimates after loss of lock event during
combined GPS/GLONASS PPP. The legend identifies the GPS satellites
involved.
158
Chapter 7. Combined GPS/GLONASS PPP
A frequency analysis was made of the height component for the final 2 h converged
period of this dataset as the GPS and GPS/GLONASS PPP solutions show quite dif-
ferent behaviour (Figures 7.36 and 7.37). The Lomb-Scargle periodogram (Scargle,
1982) in Figure 7.41 shows how the amplitude of the signal varies at different wave-
lengths. The peak at 600 s in the GPS solution is reduced in the GPS/GLONASS
solution, however there is more power in the GPS/GLONASS solution at short
wavelengths corresponding to high frequency noise. This high frequency noise is
most likely due to either receiver measurement errors or high frequency noise in the
GLONASS clocks.
10-1
10-2
10-3
Amplitude (m)
10-4
10-5
10-6 0
10 101 102 103 104
Wavelength (s)
Figure 7.41: Frequency analysis of the height component timeseries
for the final 2 h of the Tyne survey vessel dataset
7.6 Conclusion
The GLONASS measured carrier phase range residuals for GPS week 1605 were
compared to the GPS carrier phase range residuals for the same period using static
PPP at 87 IGS sites. Despite an apparent larger GLONASS RMS range residual
at elevation angles lower than 40◦ , the mean difference in RMS range residuals
159
Chapter 7. Combined GPS/GLONASS PPP
between GPS and GLONASS was not found to be statistically significantly (p <
0.05) different to zero. This is consistent with the results found by Hesselbarth
and Wanninger (2008). A single elevation weighting function was fitted to the
carrier phase range residuals and used in the subsequent comparison of GPS and
GPS/GLONASS PPP. This result would benefit from a longer measurement time
period, which would allow the seasonal effects of the residual unmodelled error
sources to be evaluated.
Previous studies have been inconclusive as to the benefit of the inclusion of GLONASS.
Cai and Gao (2007) suggested that the lack of significant impact on positioning re-
sults found in their 2007 investigation was due to the incomplete GLONASS network
of 12 satellites. Since then the GLONASS network has increased to 23 satellites;
all the testing in this thesis has been done with this enlarged network. It would
be possible to assess the impact of the number of satellites in combined networks
on positioning accuracy. This could be done either by use of historical datasets
gathered at stages during the development of the GLONASS network, or simulated
by removal of specific satellites from the dataset. However, this simulated approach
would not capture the effect of the improvement in the accuracy of the GLONASS
orbit and clock products due to the enlarged network. In doing this investigation, an
understanding may be gained as to what portion of a network needs to be included to
see a significant improvement in positioning performance. This could predict when
improved positioning would be obtained from the planned additional networks, for
example, Galileo or BEIDOU.
In kinematic mode, the GPS/GLONASS PPP was shown to perform with the same
160
Chapter 7. Combined GPS/GLONASS PPP
ultimate accuracy as GPS only PPP with a total position error RMS of 3 cm. A
clear improvement is seen in the robustness of the GPS/GLONASS PPP coordinates
resulting in reduced coordinate RMS errors over whole sessions. The additional
GLONASS satellites can maintain centimetre accurate positioning in the presence
of GPS only cycle slips. Whilst the inclusion of GLONASS satellites reduces the
power of the outlier detection tests, the number and size of the jumps introduced by
undetected outliers is reduced when GLONASS. Convergence time is seen to reduce
with GLONASS by up to a factor of three, as in the case of the boat on the Tyne.
Analysis of the pickup truck showed the difficulty in obtaining a high quality dataset
for a land based route in the United Kingdom. Although this may be representative
of real world usage, with poor quality datasets, it was not possible to produce high
quality reference coordinates and therefore it was hard to draw conclusions from the
comparison made between any test solutions. An improvement in the convergence
time was seen during the three periods whilst the pickup truck was stationary.
The testing in this chapter has reinforced how important it is to have a quality
reference solution when comparing PPP performance. With hindsight, using carrier
phase GPS as the reference solution limits the ability to perform the tests in remote
locations, or in challenging environments such as obstructed sky-view or in the
presence of high multi-path. In these situations, the GPS reference solution suffers
from degraded accuracy due to reduced satellite numbers and multipath effects
caused by the local environment that are not removed by the double differencing.
In assessing PPP, the advantage of using GPS is the shared antenna, avoiding any
offset in sensors that would have to be accounted for in the other options proposed. A
disadvantage of using relative carrier phase GPS as the reference solution is that it is
161
Chapter 7. Combined GPS/GLONASS PPP
still influenced by multipath. This can produce both a biased GPS reference solution
and biased PPP solution, as they can be biased by the same error; this potentially
reduces the reported error. This potential correlation has not been accounted for in
this thesis.
162
Chapter 8 Conclusion
8.1 Introduction
The research work contained in this thesis focused on investigating the impact of
using GLONASS satellite observations in addition to GPS observations in PPP.
In the offshore industry and work in other remote areas, there is a demand for
shortened convergence times and improved robustness of positioning techniques.
The investigation has been completed in several main parts:
2. Creation of the program pppncl, using the starting point of Track to build a
program to perform PPP, then validation in both static and kinematic modes.
Section 8.2 summarises the conclusions made from these investigations, and recom-
mendations for further investigations are given in Section 8.3.
8.2 Conclusions
Based on the investigations described in the previous chapters, the following sections
summarise the conclusions that can be drawn.
• A PPP software program, pppncl, has been created which in GPS only mode is
capable of providing positioning comparable to other published software using
this method.
• pppncl can read RINEX 2.11 observation data of unlimited length. Naviga-
tion data (satellite orbit and clock products) can be provided in either SP3 or
RINEX 2 navigation format. Additional high rate clock products in RINEX
clock format can be read and incorporated. Receivers providing C1 code ob-
servations instead of the P1 observations required for compatibility with IGS
163
Chapter 8. Conclusion
clock products can be used with the addition of code P1C1 bias corrections
provided in the CC2NONCC format.
• With the addition of the mappings from GLONASS satellite number to trans-
mission frequency provided in RINEX 2 GLONASS navigation format, it is
possible to include GLONASS carrier phase observations alongside GPS in the
coordinate estimation.
• pppncl produced zenith wet tropospheric delays show good agreement at the
centimetre level with the IGS published troposphere product, even during
periods of high variability.
• pppncl was tested in kinematic mode using a vehicle on an open airfield, driving
at a range of speeds, over a period of several hours, and found to be capable
of producing coordinates with an accuracy of 1.5 cm RMS in the North and
East components, and 3.4 cm in the Height component.
• The accuracy seen with pppncl was consistent throughout the test and not
adversely affected by the increased velocity or varying dynamics of the vehicle.
164
Chapter 8. Conclusion
• The relative noise in the GPS and GLONASS carrier phase observations was
quantified by analysis of the carrier phase residuals during static PPP. No
significant difference between the GPS and GLONASS residuals was found.
• Using two different sets of GPS/GLONASS orbit products, the ESA Final
product and the experimental Veripos Ultra real-time product, the effect
of the addition of GLONASS on static PPP repeatability was assessed at
a range of processing lengths. The improvement in repeatability seen for
GPS/GLONASS PPP compared with GPS only after 1 h of static positioning
with 30 s observation intervals was 20%, 2% and 24% for the North, East and
Height components respectively using ESA Final orbit and clock products.
When Veripos Apex G2 orbit and clock products were used, the improvement
after the same time was 24%, 10% and 14% for the North, East and Height
components respectively.
Positioning quality was compared for several routes covering a range of transport
modes to understand the impact of using GPS/GLONASS observations versus GPS
only observations in kinematic PPP. In each case a reference solution was generated,
using the dual frequency double differencing method as implemented in Track, and
used as truth coordinates.
• The challenges facing PPP were clearly demonstrated in a test with a road
vehicle on public roads, where a continual loss of lock on satellites was seen.
The resulting short continuous phase connected arcs kept the PPP solution
continuously returning to the initialisation phase. This also significantly re-
duced the accuracy of the double difference reference solution due to the short
165
Chapter 8. Conclusion
time lengths without loss of lock, meaning the carrier phase ambiguities could
not be well enough defined to be fixed at their integer values. With such a
poor reference solution, no meaningful quantitative comparison between the
GPS only and GPS/GLONASS PPP solutions could be made. Qualitatively,
however, during the three 1 h periods for which a usable reference solution was
available, the GPS/GLONASS solution showed faster convergence and hence
more accurate positioning than the GPS only solution.
• With the data collected during the boat’s route on the River Tyne, it can
be seen that there is a clear reduction in convergence time when GPS and
GLONASS are combined. The convergence time with GLONASS was reduced
to approximately half an hour, three times less than for the GPS only PPP. Ad-
ditionally, with the GPS/GLONASS combined observations, continuous high
accuracy positioning was maintained during a period where the receiver lost
lock on all GPS carrier phase observations. This demonstrated the benefit
of using GLONASS in addition to GPS to provide a more robust positioning
solution.
• The accuracy of the fully converged GPS and GPS/GLONASS PPP solutions
was similar with 3 cm 3D RMS. However the session horizontal RMS was
almost halved due to the increased robustness and faster convergence time of
the GPS/GLONASS solution.
• The total length of convergence time varies between the datasets, depending
on the specific environment, although an overall improvement is seen with the
addition of GLONASS across all the datasets studied. This applies not only
to the time to the first usable measurements after powering up the receiver,
but also to the time to recover after catastrophic events such as high levels of
interference, accidental receiver reset, or total loss of lock due to obstructions
such as bridges or dense tree cover.
• With ideal observations, when the receiver maintains lock for a long enough
period (greater than 1 h) the RMS of GPS only PPP and GPS/GLONASS
PPP are comparable. However, in the less than ideal case, often found in the
166
Chapter 8. Conclusion
• The use of the Condor distributed computing platform allowed for rapid pro-
cessing of the many discrete units of work such as when processing the 133,110
167
Chapter 8. Conclusion
The following are suggestions for further research into this area:
• With the recent change to the IGS08 reference frame in April 2011, the IGS
introduced GLONASS specific receiver antenna phase centre offsets and vari-
ations. The impact of this on PPP positioning could be investigated, with a
possibility of increased accuracy.
• One could investigate the use of auto correlation and adaptive Kalman filtering
techniques for automatic quantification of the received observation noise.
• In October 2011, the first launch of Galileo in-orbit validation satellites takes
place. This presents an opportunity for further PPP positioning improvement
with the addition of Galileo to GPS/GLONASS satellite observations. The
inclusion of Galileo would bring improvements in satellite geometry due to
increased number of satellites. The Galileo system uses CDMA, like GPS, so
avoids the complication of each satellite transmitting on a separate frequency,
and therefore there are no inter-frequency biases on the same signal. Galileo
should also be less susceptible to multipath due to the code signal used.
168
Chapter 8. Conclusion
Galileo time scales present in the precise clock products, and any biases intro-
duced between Galileo and the other utilised GNSS in the receiver. Galileo
will transmit on E1 and E5a frequency bands, in common with GPS to im-
prove compatibility with GPS (Hofmann-Wellenhof et al., 2008). Therefore
the same ionosphere free combination can be used for Galileo as for GPS.
• The error budget for Galileo is expected to be lower than that of GPS due
to the more stable atomic clocks in the satellites and the more modern rang-
ing codes that are less susceptible to multipath. Therefore, considering GPS
combined with Galileo, it would be expected to provide improved performance
in terms of the metrics described in Section 2.9, by comparison with GPS/-
GLONASS.
• GLONASS and GPS systems are being modernised, with GLONASS adding
CDMA encoding, the addition of the L5 frequency to GPS, and additional
civilian signals such as L1C and L2C. Further work could be done to make
best use of these new features in positioning accuracy and robustness.
169
Chapter 8. Conclusion
• One of the key benefits of combined GPS/GLONASS PPP is the reduced ini-
tialisation time. The recent developments in ambiguity fixed PPP still require
an initial convergence period of up to 90 min until the integer valued biases
are sufficiently well determined (Laurichesse et al., 2009). The addition of
GLONASS observations to ambiguity fixed PPP could allow for a reduced
initialisation time until the integer biases are identified.
170
References
Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz,
J., Greenbaum, A., Hammarling, S., McKenney, A., and Sorensen, D. (1999).
LAPACK Users’ Guide. Philadelphia, PA: Society for Industrial and Applied
Mathematics, 3rd ed.
Armatys, M., Muellerschoen, R., Bar-Sever, Y., and Meyer, R. (2003). Demonstra-
tion of decimeter-level real-time positioning of an airborne platform. In Proceed-
ings of ION NTM . Anaheim, CA, USA.
Arroyo-Suarez, E., Riley, J., Glang, G., and Mabey, D. (2005). Evaluating a global
differential GPS system for hydrographic surveying. In Proceedings of OCEANS
2005 MTS/IEEE, vol. 3, (pp. 2557–2563). Institute of Electrical and Electronics
Engineers. doi:10.1109/OCEANS.2005.1640155.
Banville, S., and Langley, R. (2010). Instantaneous cycle-slip correction for real-
time PPP applications. Navigation, Journal of the Institute of Navigation, 57 (4),
325–334.
Banville, S., Langley, R. B., and Santos, M. C. (2009). The Precise Point Positioning
Software Centre: An Insight Into Online PPP Services. In Poster presented at the
IAG 2009 Meeting. Buenos Aires, Argentina.
Bar-Sever, Y. E. (1996). A new model for GPS yaw attitude. Journal of Geodesy,
70 , 714–723. doi:10.1007/BF00867149.
171
REFERENCES
Barker, R., Lapucha, D. D., and Wood, T. (2002). The impact of high performance
GPS on the offshore marine survey, navigation and positioning industry. In Pro-
ceedings of the Offshore Technology Conference, 6-9 May 2002, Houston, Texas.
doi:10.4043/14195-MS.
Belehaki, A., Cander, L., Zolesi, B., Bremer, J., Juren, C., Stanislawska, I., Dialetis,
D., and Hatzopoulos, M. (2006). Monitoring and forecasting the ionosphere over
Europe: The DIAS project. Space Weather, 4 . doi:10.1029/2006SW000270.
Beutler, G., Davidson, D., Langley, R., Santerre, R., Vanicek, P., and Wells, D.
(1984). Some theoretical and practical aspects of geodetic positioning using carrier
phase difference observations of GPS satellites. Technical Report 109, University
of New Brunswick, Canada.
Bisnath, S. (2004). Precise orbit determination of low earth orbiters with a sin-
gle GPS receiver based, geometric strategy. Tech. Rep. 220, University of New
Brunswick.
Bisnath, S., and Gao, Y. (2008). Current state of precise point positioning and
future prospects and limitations. In F. Sansò, and M. G. Sideris (Eds.) Observing
our Changing Earth, vol. 133 of International Association of Geodesy Symposia,
(pp. 615–623). Springer. doi:10.1007/978-3-540-85426-5_71.
Blewitt, G. (1989). Carrier phase ambiguity resolution for the global positioning sys-
tem applied to geodetic baselines up to 2000 km. Journal of Geophysical Research,
94 (B8), 10187–10203. doi:10.1029/JB094iB08p10187.
Bock, H., Dach, R., Jäggi, A., and Beutler, G. (2009). High-rate GPS clock correc-
tions from CODE: support of 1Hz applications. Journal of Geodesy, 83 , 1083–
1094. doi:10.1007/s00190-009-0326-1.
172
REFERENCES
Böhm, J., Heinkelmann, R., and Schuh, H. (2007). Short note: A global model
of pressure and temperature for geodetic applications. Journal of Geodesy, 81 ,
679–683. doi:10.1007/s00190-007-0135-3.
Böhm, J., Niell, A., Tregoning, P., and Schuh, H. (2006a). Global Mapping Function
(GMF): A new empirical mapping function based on numerical weather model
data. Geophysical Research Letters, 33 (7), L07304. doi:10.1029/2005GL025546.
Böhm, J., Werl, B., and Schuh, H. (2006b). Troposphere mapping functions for GPS
and very long baseline interferometry from European Centre for Medium-Range
Weather Forecasts operational analysis data. J. Geophys. Res., 111:B02406 . doi:
10.1029/2005JB003629.
Boucher, C., and Altamimi, Z. (2001). ITRS, PZ-90 and WGS 84: current real-
izations and the related transformation parameters. Journal of Geodesy, 75 (11),
613–619. doi:10.1007/s001900100208.
Braasch, M., and Van Dierendonck, A. J. (1999). GPS receiver architectures and
measurements. In Proceedings of the IEEE, vol. 87, (pp. 48–64). doi:10.1109/
5.736341.
Bucy, R. S., and Joseph, P. D. (1968). Filtering for Stochastic Processes, with
Applications to Guidance. New York: Wiley.
Byun, S., and Bar-Sever, Y. (2009). A new type of troposphere zenith path delay
product of the International GNSS Service. Journal of Geodesy, 83 (3), 1–7. doi:
10.1007/s00190-008-0288-8.
Cai, C., and Gao, Y. (2007). Precise point positioning using combined GPS and
GLONASS observations. Journal of Global Positioning Systems, 6 (1), 13–22.
doi:10.5081/jgps.6.1.13.
Chen, G. (1998). GPS Kinematic Positioning for the Airborne Laser Altimetry at
Long Valley, California. Ph.D. thesis, Massachusetts Institute of Technology.
Department of Earth, Atmospheric, and Planetary Sciences.
Collins, J., and Langley, R. (1999). Possible weighting schemes for GPS carrier
phase observations in the presence of multipath. Final contract report for the US
173
REFERENCES
Colombo, O., Sutter, A. W., and Evans, A. G. (2004). Evaluation of precise, kine-
matic GPS point positioning. In Proceedings of the ION GPS GNSS 2004 Meeting,
(pp. 21–24).
Dach, R., Böhm, J., Lutz, S., Steigenberger, P., and Beutler, G. (2011a). Evaluation
of the impact of atmospheric pressure loading modeling on GNSS data analysis.
Journal of Geodesy, 85 (2), 75–91. doi:10.1007/s00190-010-0417-z.
Dach, R., Brockmann, E., Schaer, S., Beutler, G., Meindl, M., Prange, L., Bock,
H., Jäggi, A., and Ostini, L. (2009). GNSS processing at CODE: status report.
Journal of Geodesy, 83 , 353–365. doi:10.1007/s00190-008-0281-2.
Dach, R., Schmid, R., Schmitz, M., Thaller, D., Schaer, S., Lutz, S., Steigenberger,
P., Wübbena, G., and Beutler, G. (2011b). Improved antenna phase center models
for GLONASS. GPS Solutions, 15 , 49–65. doi:10.1007/s10291-010-0169-5.
Dahlquist, G., and Björck, A. (2003). Numerical methods. Dover Publications Inc.
Daum, F. (2005). Nonlinear filters: beyond the Kalman filter. IEEE Aerospace and
Electronic Systems Magazine, 20 (8), 57–69. doi:10.1109/MAES.2005.1499276.
Davis, J., Herring, T., Shapiro, I., Rogers, A., and Elgered, G. (1985). Geodesy by
radio interferometry: Effects of atmospheric modeling errors on estimates of base-
line length. Radio Science, 20 (6), 1593–1607. doi:10.1029/RS020i006p01593.
Dilssner, F., Springer, T., Gienger, G., and Dow, J. (2011). The GLONASS-M
satellite yaw-attitude model. Advances in Space Research, 47 (1), 160–171. doi:
10.1016/[Link].2010.09.007.
174
REFERENCES
Dixon, K. (2006). StarFire TM: A global SBAS for sub-decimetre precise point
positioning. In Proceedings of ION GNSS , (pp. 2286–2296). Fort Worth, Texas.
Dodd, D. (2007). Utility of Ionosphere and Troposphere Models for Extending the
Range of High-accuracy GPS.. Ph.D. thesis, The University of Southern Missis-
sippi.
Doherty, P. H., Klobuchar, J. A., and Kunches, J. M. (2000). Eye on the ionosphere:
The correlation between solar 10.7 cm radio flux and ionospheric range delay. GPS
Solutions, 3 , 75–79. doi:10.1007/PL00012820.
Dong, D.-N., and Bock, Y. (1989). Global Positioning System Network analysis with
phase ambiguity resolution applied to crustal deformation studies in California.
Journal of Geophysical Research, 94 (B4), 3949. doi:10.1029/JB094iB04p03949.
Dow, J., Neilan, R. E., and Rizos, C. (2009). The International GNSS Service in a
changing landscape of Global Navigation Satellite Systems. Journal of Geodesy,
83 , 191–198. doi:10.1007/s00190-008-0300-3.
Edwards, S., Clarke, P., Penna, N., and Goebell, S. (2010). An examination of
Network RTK GPS services in Great Britain. Survey Review, 42 (316), 107–121.
doi:10.1179/003962610X12572516251529.
Estey, L. H., and Meertens, C. M. (1999). TEQC: The multi-purpose toolkit for
GPS/GLONASS data. GPS Solutions, 3 , 42–49. doi:10.1007/PL00012778.
Euler, H., Zebhauser, B., Townsend, B., and Wübbena, G. (2002). Comparison of
different proposals for reference station network information distribution formats.
In Proceedings of ION GPS , (pp. 2349–2358). Portland, OR.
Feng, S., Ochieng, W., Moore, T., Hill, C., and Hide, C. (2009). Carrier phase-based
integrity monitoring for high-accuracy positioning. GPS Solutions, 13 (1), 13–22.
doi:10.1007/s10291-008-0093-0.
175
REFERENCES
Gao, Y., Lahaye, F., Héroux, P., Liao, X., Beck, N., and Olynik, M. (2001). Modeling
and estimation of C1–P1 bias in GPS receivers. Journal of Geodesy, 74 (9), 621–
626. doi:10.1007/s001900000117.
Gao, Y., and Shen, X. (2002). A new method for carrier-phase-based precise point
positioning. Navigation, 49 (2), 109–116.
Gao, Y., and Wang, M. (2008). Real-time kinematic OTF positioning using a single
GPS receiver. In M. G. Sideris (Ed.) Observing our Changing Earth, vol. 133 of
International Association of Geodesy Symposia, (pp. 655–668). Berlin, Heidelberg:
Springer. doi:10.1007/978-3-540-85426-5_76.
Gao, Y., Wojciechowski, A., and Chen, K. (2005). Airborne kinematic positioning
using precise point positioning methodology. Geomatica, 59 (1), 275–282.
Ge, M., Gendt, G., Rothacher, M., Shi, C., and Liu, J. (2008). Resolution of GPS
carrier-phase ambiguities in precise point positioning (PPP) with daily observa-
tions. Journal of Geodesy, 82 (7), 401–401. doi:10.1007/s00190-007-0208-3.
Geng, J., Meng, X., Dodson, A., Ge, M., and Teferle, F. (2010a). Rapid re-
convergences to ambiguity-fixed solutions in precise point positioning. Journal
of Geodesy, 84 (12), 705–714. doi:10.1007/s00190-010-0404-4.
Geng, J., Meng, X., Dodson, A., and Teferle, F. (2010b). Integer ambiguity resolu-
tion in precise point positioning: method comparison. Journal of Geodesy, (pp.
1–13). doi:10.1007/s00190-010-0399-x.
Geng, J., Meng, X., Teferle, F., and Dodson, A. (2010c). Performance of precise
point positioning with ambiguity resolution for 1 to 4 hour observation periods.
Survey Review, 42 (316), 155–165. doi:10.1179/003962610X12572516251682.
176
REFERENCES
Görres, B., Campbell, J., Becker, M., and Siemes, M. (2006). Absolute calibration of
GPS antennas: laboratory results and comparison with field and robot techniques.
GPS Solutions, 10 (2), 136–145. doi:10.1007/s10291-005-0015-3.
Grewal, M. S., and Andrews, A. P. (2001). Kalman Filtering: Theory and Practice
Using MATLAB. John Wiley & Sons, Inc., 2nd ed.
Griffiths, J., and Ray, J. R. (2009). On the precision and accuracy of IGS orbits.
Journal of Geodesy, 83 (3-4), 277–287. doi:10.1007/s00190-008-0237-6.
Grinter, T., and Roberts, C. (2011). Precise point positioning: where are we now.
In Proceedings of IGNSS Symp.. Sydney, Australia.
Gurtner, W., and Estey, L. (2007a). RINEX – the receiver independent ex-
change format-version 3.00. [Link]
[Link].
Gurtner, W., and Estey, L. (2007b). RINEX: The receiver independent exchange for-
mat version 2.10. [Link]
txt.
Hatanaka, Y. (2008). A compression format and tools for GNSS observation data.
Bulletin of the Geographical Survey Institute, 55 , 21–30.
Hegarty, C., and Chatre, E. (2008). Evolution of the global navigation satellite
system (GNSS). Proceedings of the IEEE, 96 (12), 1902–1917. doi:10.1109/
JPROC.2008.2006090.
Héroux, P., Y. Gao, J. K., Lahaye, F., Mireault, Y., Collins, P., Macleod, K.,
Tétreault, P., and Chen, K. (2004). Products and applications for precise point
positioning - moving towards real-time. In Proceedings of ION GNSS 2004 , (pp.
1832–1843). Long Beach, CA.
Herring, T., King, R., and McClusky, S. (2006). GPS analysis at MIT, GAMIT
reference manual, release 10.3. Massachusetts Institute of Technology.
177
REFERENCES
Hilla, S. (2002). Extending the standard product 3 (SP3) orbit format. In Pro-
ceedings of the International GPS Service Network, Data, and Analysis Center
Workshop. Ottawa, Canada.
Keshin, M., Le, A., and van der Marel, H. (2006). Single and dual-frequency precise
point positioning: approaches and performances. In NAVITEC 2006 . Noordwijk,
The Netherlands.
Kouba, J., and Héroux, P. (2001). Precise point positioning using IGS orbit and
clock products. GPS solutions, 5 (2), 12–28. doi:10.1007/PL00012883.
Lachapelle, G. (1991). GPS observables and error sources for kinematic positioning.
In Proceedings of IAG International Symposium, vol. 107, (pp. 17–26). New York:
Springer Verlag. doi:10.1007/978-1-4612-3102-8_2.
178
REFERENCES
Lau, L., and Cross, P. (2006). A new signal-to-noise-ratio based stochastic model for
GNSS high-precision carrier phase data processing algorithms in the presence of
multipath errors. In Institute of Navigation - 19th International Technical Meeting
of the Satellite Division, ION GNSS 2006 , vol. 1, (pp. 276–285). Fort Worth, TX.
Lau, L., and Cross, P. (2007). Development and testing of a new ray-tracing ap-
proach to GNSS carrier-phase multipath modelling. Journal of Geodesy, 81 (11),
713–732. doi:10.1007/s00190-007-0139-z.
Laurichesse, D. (2011). The CNES real-time PPP with undifferenced integer ambi-
guity resolution demonstrator. In Proceedings of the ION GNSS 2011 . Portland,
Oregon.
Laurichesse, D., Mercier, F., Berthias, J. P., Broca, P., and Cerri, L. (2009). Integer
ambiguity resolution on undifferenced GPS phase measurements and its applica-
tion to PPP and satellite precise orbit determination. Navigation, 56 (2), 135–149.
Leandro, R., Santos, M., and Langley, R. (2011). Analyzing GNSS data in
precise point positioning software. GPS Solutions, 15 , 1–13. doi:10.1007/
s10291-010-0173-9.
Leick, A. (2004). GPS satellite surveying. John Wiley & Sons, 3 ed.
Li, X., Zhang, X., and Ge, M. (2011). Regional reference network augmented precise
point positioning for instantaneous ambiguity resolution. Journal of Geodesy,
85 (3), 151–158. doi:10.1007/s00190-010-0424-0.
Martín, A., Anquela, A. B., Capilla, R., and Berné, J. L. (2011). PPP technique
analysis based on time convergence, repeatability, IGS products, different software
processing and GPS+Glonass constellation. Journal of Surveying Engineering,
137 (3), 99–108. doi:10.1061/(ASCE)SU.1943-5428.0000047.
McCarthy, D. D., Boucher, C., Eanes, R., Fukushima, T., Herring, T., Lieske, J.,
Ma, C., Montag, H., Pâquet, P., Reigber, C., et al. (1989). IERS Standards
179
REFERENCES
McCarthy, D. D., and Petit, G. (2003). IERS Conventions (2003). IERS Technical
Note 32, Verlag des Bundesamts für Kartographie und Geodäsie, Frankfurt am
Main.
Mendes, V., and Langley, R. (1998). Tropospheric zenith delay prediction accuracy
for airborne GPS high-precision positioning. In Proceedings of The Institute of
Navigation 54th Annual Meeting, (pp. 337–347).
Metcalf, M., and Reid, J. K. (1999). Fortran 90/95 Explained. Oxford University
Press.
Monteiro, L. S., Moore, T., and Hill, C. (1999). What is the accuracy of DGPS?
Journal of Navigation, 58 (2), 207–225. doi:10.1017/S037346330500322X.
Muellerschoen, R. J., Reichert, A., Kuang, D., Heflin, M., Bertiger, W. I., and Bar-
Sever, Y. E. (2001). Orbit determination with NASA’s high accuracy real-time
global differential GPS system. In Proceedings of ION GPS-2001 . Salt Lake City,
UT.
Newcomb, S. (1895). Tables of the motion of the Earth on its axis and around the
Sun. Astronomical papers prepared for the use of the American ephemeris and
nautical almanac , [Washington, Bureau of Equipment, Navy Dept., 1895], 6 (6),
169.
Niell, A. E. (1996). Global mapping functions for the atmosphere delay at radio
wavelengths. Journal of Geophysical Research, 101 (B2), 3227–3246. doi:10.
1029/95JB03048.
Noll, C., and Dube, M. (2001). The IGS global data center at the CDDIS — an
update. Physics and Chemistry of the Earth, Part A: Solid Earth and Geodesy,
26 (6-8), 603–604. Proceedings of the First COST Action 716 Workshop Towards
180
REFERENCES
Operational GPS Meteorology and the Second Network Workshop of the Interna-
tional GPS Service (IGS). doi:10.1016/S1464-1895(01)00108-9.
Ochieng, W. Y., Sauer, K., Walsh, D., Brodin, G., Griffin, S., and Denney, M. (2003).
GPS integrity and potential impact on aviation safety. Journal of Navigation,
56 (1), 51–65. doi:10.1017/S0373463302002096.
Oleynik, E., Mitrikas, V., Revnivykh, S., Serdukov, A., Dutov, E., and Shiriaev, V.
(2006). High-accurate GLONASS Orbit and Clock Determination for the Assess-
ment of System Performance. In ION GNSS , (pp. 2065–2079).
Petrie, E., Hernández-Pajares, M., Spalla, P., Moore, P., and King, M. (2011).
A review of higher order ionospheric refraction effects on dual frequency GPS.
Surveys in Geophysics, 32 , 197–253. doi:10.1007/s10712-010-9105-z.
Petrov, L. (2004). Study of the atmospheric pressure loading signal in very long base-
line interferometry observations. Journal of Geophysical Research, 109 , B03405.
doi:10.1029/2003JB002500.
Píriz, R., Calle, D., Mozo, A., Navarro, P., Rodríguez, D., and Tobías, G. (2009).
Orbits and clocks for GLONASS precise-point-positioning. In Proceedings of the
22nd International Technical Meeting of The Satellite Division of the Institute of
Navigation (ION GNSS 2009), (pp. 2415–2424).
Píriz, R., Mozo, A., Navarro, P., and Rodriguez, D. (2008). magicGNSS: Precise
GNSS products out of the box. In Proceedings of the 21st International Technical
Meeting of the Satellite Division of The Institute of Navigation (ION GNSS 2008),
(pp. 1242–1251).
Press, W., Teukolsky, S., Vetterling, W., Flannery, B., and Metcalf, M. (1996).
Numerical recipes in Fortran 90: the art of parallel scientific computing. Volume
2 of Fortran numerical recipes. Cambridge University Press.
Ragheb, A., Clarke, P., and Edwards, S. (2007). GPS sidereal filtering: coordinate-
and carrier-phase-level strategies. Journal of Geodesy, 81 (5), 325–335. doi:
10.1007/s00190-006-0113-1.
Ray, J., and Gurtner, W. (2010). RINEX extensions to handle clock information.
[Link]
181
REFERENCES
Rocken, C. (2005). Atmospheric water vapor and geoid measurements in the open
ocean with GPS. Geophysical Research Letters, 32 , L12813. doi:10.1029/
2005GL022573.
Rocken, C., Mervart, L., Johnson, J., Lukes, Z., Springer, T., Iwabuchi, T., and
Cummins, S. (2011). New real-time global GPS and GLONASS precise position-
ing correction service: Apex. In Proceedings of the 24th International Technical
Meeting of The Satellite Division of the Institute of Navigation (ION GNSS 2011),
(pp. 1825–1838). Portland, OR.
Rost, C., and Wanninger, L. (2009). Carrier phase multipath mitigation based on
GNSS signal quality measurements. Journal of Applied Geodesy, 3 (2), 81–87.
doi:10.1515/JAG.2009.009.
Rothacher, M., and Mader, G. (2003). Receiver and satellite antenna phase center
offsets and variations. In P. Tétreault, R. Neilan, and K. Gowey (Eds.) Proceedings
of the Network, Data and Analysis Centre 2002 Workshop, (pp. 141–152). Ottawa.
[Link]
Rothacher, M., and Schmid, R. (2010). ANTEX: The antenna exchange format,
version 1.4. [Link]
txt.
182
REFERENCES
Schaer, S., and Steigenberger, P. (2006). Determination and use of GPS differential
code bias values. In IGS workshop in Darmstadt.
Thain, D., Tannenbaum, T., and Livny, M. (2005). Distributed computing in prac-
tice: the condor experience. Concurrency - Practice and Experience, 17 (2-4),
323–356. doi:10.1002/cpe.938.
Thomas, I. D., King, M. A., Bentley, M. J., Whitehouse, P. L., Penna, N. T.,
Williams, S. D. P., Riva, R. E. M., Lavallee, D. A., Clarke, P. J., and King,
E. C. (2011). Widespread low rates of Antarctic glacial isostatic adjustment
revealed by GPS observations. Geophysical Research Letters, 38 , L22302. doi:
10.1029/2011GL049277.
Urquhart, L. (2009). Atmospheric pressure loading and its effects on precise point
positioning. In Proceedings of the 22nd International Technical Meeting of The
183
REFERENCES
Satellite Division of the Institute of Navigation (ION GNSS 2009), (pp. 658–667).
Savannah, GA.
Urschl, C., Dach, R., Hugentobler, U., Schaer, S., and Beutler, G. (2005). Validating
ocean tide loading models using GPS. Journal of Geodesy, 78 (10), 616–625.
doi:10.1007/s00190-004-0427-9.
Vermeer, M. (1997). The precision of geodetic GPS and one way of improving it.
Journal of Geodesy, 71 (4), 240–245. doi:10.1007/s001900050091.
Wahr, J. (1981). The forced nutation of an elliptical, rotating, elastic, and ocean less
Earth. Geophys. J. Roy. Astron. Soc., 64 , 705–727. doi:10.1111/j.1365-246X.
1981.tb02691.x.
Wan, E., and Van Der Merwe, R. (2000). The unscented Kalman filter for non-
linear estimation. In Proceedings of the IEEE 2000 Adaptive Systems for Signal
Processing, Communications, and Control Symposium, (pp. 153–158). Institute
of Electrical and Electronics Engineers. doi:10.1109/ASSPCC.2000.882463.
Wang, C., Feng, Y., Higgins, M., and Cowie, B. (2010). Assessment of commer-
cial network RTK user positioning performance over long inter-station distances.
Journal of Global Positioning Systems, 9 (1), 78–89. doi:10.5081/jgps.9.1.78.
Weber, G. (2006). Streaming real-time IGS data and products using NTRIP. In
Proceedings of the IGS2006 Workshop. ESOC Darmstadt, Germany.
Weber, G., Mervart, L., Lukes, Z., Rocken, C., and Dousa, J. (2007). Real-time clock
and orbit corrections for improved point positioning via NTRIP. In Proceedings
of ION GNSS 2007 , (pp. 1992–1998). Fort Worth, Texas: the 20th International
Technical Meeting of the Satellite Division of The Institute of Navigation.
Welch, G., and Bishop, G. (1995). An introduction to the Kalman filter. TR 95-041,
University of North Carolina at Chapel Hill, Chapel Hill, NC 27599-3175. URL:
[Link]
184
REFERENCES
Wu, J., Wu, S., Hajj, G., Bertiger, W., and Lichten, S. (1993). Effects of antenna
orientation on GPS carrier phase. Manuscripta Geodaetica, 18 , 91–98.
Wulf, W., and Shaw, M. (1973). Global variable considered harmful. ACM SIG-
PLAN Notices, 8 , 28–34. doi:10.1145/953353.953355.
Zumberge, J., Heflin, M., Jefferson, D., Watkins, M., and Webb, F. (1997). Precise
point positioning for the efficient and robust analysis of GPS data from large
networks. Journal of Geophysical Research, 102 (B3), 5005–5017. doi:10.1029/
96JB03860.
185
Appendix A Software reference
This appendix provides reference information on the software created as part of this
thesis. The pppncl configuration file format is described in A.2. Output formats
are described. Utility scripts functionality and command line options are given.
Overview of the key functions and data types used in the pppncl source code.
This section describes the types of input file that are needed by pppncl when per-
forming PPP. The different formats that can be used for each type of input data are
explained.
Observation data can be read in RINEX 2.1 format (Gurtner and Estey, 2007b).
As the file is read sequentially epoch by epoch, there is no limit to the length of
the observation file that can be used. Files containing up to 14 observation types
at each epoch can be read. The software places no restriction on the observation
type identifiers in the RINEX file. The limit of 14 observations allows reading of
files containing the currently complete following set of observations:
• P code: P1, P2
Legacy GPS only files with no satellite system identifier or mixed GNSS system files
can be read.
RINEX 2 was designed with only GPS in mind, and later adapted to allow mixed
data from GPS and GLONASS satellites in the same file. It has become apparent
that it has shortcomings as a reciever observation exchange format for files con-
taining observation data of more than one satellite system, each one with different
observation types. To address this Gurtner and Estey (2007a) proposed the RINEX
3 standard. Whilst it seems logical that RINEX 3 will replace RINEX 2 as the
186
Appendix A. Software reference
standard data exchange format the not inconsiderable effort required to update all
the tools used by the wider GNSS community has meant that RINEX 3 adoption
has been very slow. For that reason only a RINEX 2 reader has been implemented
in pppncl at this time.
pppncl can use broadcast orbits in RINEX format or precise orbits in SP3 format,
with optional high rate clocks in RINEX clock format. For the RINEX GPS navi-
gation format and the SP3 format, the routines from Track were used as a starting
point. The SP3 routines were enhanced to allow reading of files containing multiple
GNSS systems so that GLONASS data could be used.
The RINEX format contains a specification for the file format containing the broad-
cast ephemeris data; this can be read by pppncl.
Precise orbits and clocks can be read in SP3 format (Hilla, 2002). Repeated header
blocks are ignored so for observation periods of longer than a day, or to avoid
interpolation errors at the beginning and the end of the day, multiple SP3 files may
be concatenated into a single file. Subsequent header blocks will be ignored.
High rate clock data, essential to high rate PPP can be read in RINEX 2.00 or 3.00
clock format (Ray and Gurtner, 2010). Only the satellite clock offset records (lines
beginning AS) are read in by pppncl, other lines are skipped.
For the mapping between GLONASS satellite number and GLONASS frequency
number, GLONASS navigation data is required. This can be read in RINEX 2.10
GLONASS navigation format. The CDDIS produce a merged navigation file con-
taining the union of all broadcast navigation information received across all sites
stored in the CDDIS archive.
The IGS maintain a list of satellite and receiver antenna phase centre offsets and
variations in ANTEX format. An ANTEX file parser was written to read ANTEX
v1.4 format files (Rothacher and Mader, 2003; Rothacher and Schmid, 2010). The
list of receiver antennas that are contained in the IGS list can be found at http:
//[Link]/igscb/station/general/rcvr_ant.tab.
187
Appendix A. Software reference
For positioning using older orbit and clock products that were generated using rela-
tive phase centre offsets, the GAMIT (Herring et al., 2006) GPS satellite information
format [Link] may be used. Note this only provides block mean satellite phase
centre offsets and no phase centre variations so accuracy will be limited. This file
can either be obtained from a GAMIT installation in ~/gg/tables/[Link] or
online at [Link]
pppncl integrates the ESOC updated version of the CC2NONCC software routines
(Romero, 2010) for correcting C1 observables to P1 like observables as described in
Section 3.13. This allows reading CC2NONCC correction format P1C1 biases. His-
torical monthly average values are available from CODE at [Link]
aiub/CODE/p1c1bias.2000 in one file. This is only updated periodically, so for
the latest values or if values are only needed for a particular month they may be
obtained from [Link] where YYYY
denotes the four digit year, YY the two digit year and MM the two digit month. E.g.
[Link] for P1C1 biases from July
2011.
The pppncl processing options are specified in a configuration file, largely based on
the configuration file format for Track. All lines that do not start with a space are
ignored as comment lines. The commands are case insensitive. Here follows a list
of keywords used in the file, with an explanation of their effects and the format of
any data they expect.
In the following descriptions of the configuration options, values that are required by
a given command are enclosed in angled brackets (< >), optional values are enclosed
in square brackets ([ ]). These should be replaced including the brackets with the
required value or in the case of optional items ommited. For clarity, required spaces
are shown using the character ␣ and newlines with the ¶ character.
Where a filename is required as an option the value should be either the relative
path and filename or an absolute path and filename of the input file that should
be used. Relative paths are interpreted relative to the directory pppncl is invoked
188
Appendix A. Software reference
List of commands:
␣obs_file¶
␣␣<site name>␣<RINEX observation file>¶
This required command defines the RINEX observation file that is used as input
data. <site name> is a four letter identifier describing the name of the site and
may be chosen freely. <RINEX observation file> is the path and filename of the
observation file to be used, there is a limit of 256 characters.
␣nav_file␣<navigation file>␣[SP3]¶
This required command defines the navigation file used for orbits and clock data.
<navigation file> is the filename of the navigation file, and may be a maximum
length of 256 characters. By default, this is assumed to be a RINEX GPS broadcast
navigation file. If SP3 format navigation data is to be used, the optional SP3 marker
should follow the name of the navigation file. This is done rather than explicitly
auto detecting the format as some providers use an alternative file extension to .sp3.
This command is required for use of GLONASS observations. It specifies the file-
name of the GLONASS broadcast navigation file used to map the GLONASS satel-
lite number to transmitting frequency. The maximum length of the filename is 256
characters.
␣clk_file␣<clock file>¶
This optional command is used to specify that high rate clock corrections in RINEX
clock format should be used. <clock file> is the filename of the clock data to be
used, with a limit of 256 characters. If this command is supplied, the navigation file
must be in SP3 format.
␣site_pos¶
␣␣<site name>␣<X>␣<Y>␣<Z>¶
This command is used to specify the initial coordinates of the receiver in Earth
centred Earth fixed coordinates, in the reference frame of the provided navigation
data. The site name must match the site name given in the obs_file command.
The units for <X>, <Y> and <Z> are metres. These values are also used as the
189
Appendix A. Software reference
coordinates to which the North, East and Height relative coordinate output file,
[Link] is referenced.
␣site_stats¶
␣␣<site name>␣<APR X>␣<APR Y>␣<APR Z>␣<noise X>␣<noise Y>␣<noise Z>¶
This required command specifies the initial coordinate uncertainty and the epoch
to epoch process noise. <site name> must match the name given in the obs_file
command. <APR X>, <APR Y> and <APR Z> are the X, Y and Z initial coordinate
standard deviations, units are in metres. <noise X>, <noise Y> and <noise Z> are
the X, Y and Z process noise added at each epoch, units are in metres per second.
As the coordinates are modelled as random walk, the variance added at each epoch
is dt × noise2 , where dt is the observation step interval.
For static positioning, the <noise> values should be set to zero. For kinematic po-
sitioning the <noise> values should match the expected dynamics of the receiver.
Specifying a large process noise will allow the solution to respond faster to high
dynamic situations but will also increase the effect of observation noise on the esti-
mated coordinates.
␣atm_stats¶
␣␣all␣<APR ATM>␣<noise ATM>¶
This required command specifies the initial uncertainty and noise in the estimated
correction to the wet zenith tropospheric delay. <APR ATM> is the initial zenith
tropospheric delay correction standard deviation, units are in metres. <noise ATM>
is the process noise added at each epoch, units are in metres per second. As the
wet zenith tropospheric delay correction is modelled as random walk, the variance
added at each epoch is dt × noise2 . Setting both values to zero will disable the wet
tropospheric correction estimation.
␣clk_stats¶
␣␣all␣<initial std dev>␣<per epoch process noise>¶
This required command specifies the initial uncertainty and noise in the estimated
receiver clock offset. <initial std dev> is the initial estimated receiver clock offset
standard deviation, units are in metres. <per epoch process noise> is the process
noise added at each epoch, units are in metres per second. As the receiver clock
offset is modelled as random walk, the variance added at each epoch is dt × noise2 .
190
Appendix A. Software reference
␣ante_off¶
␣␣<station>␣<ARP offset>␣[L1 PCO]␣[L2 PCO]¶
If absolute antenna offsets are being used with a provided ANTEX file, then this
command is not required unless the ARP offset provided in the RINEX file needs
to be overridden, in which case the North, East and Up offsets should be given in
units of metres.
If relative phase centre offsets are being used then additionally, the L1 and L2 PCOs
may be given as North, East and Height offsets in units of metres.
␣data_noise␣<L1>␣<L2>␣<P1>␣<P2>¶
This optional command is used to specify the standard deviation of the zenith
observation process noise for each of the GPS L1, L2, P1 and P2 observables, with
units in metres. If the C1 observable is being used in place of P1, then the process
noise specified for P1 will be used for the C1 observable.
␣data_type␣<data type>¶
• P1
• P2
• UC alternative ionosphere free code and phase Gao and Shen (2002)
This optional command specifies an elevation cutoff angle in degrees below which
observations are not used. If this command is not given, the default value of 10◦ is
used.
␣gnss_type␣<GNSS types>¶
This optional command specifies the GNSS types to be used, with the same letters
191
Appendix A. Software reference
as in the RINEX file. The usable options are “G” for GPS only positioning, and
“GR” for combined GPS and GLONASS positioning. If this option is not specified,
the default is GPS only. Note, for the GR option to work, a GLONASS navigation
file must be specified with the g_nav_file command.
␣block_file␣<block file>¶
This optional command specifies the GAMIT style satellite information file. It is
only required when using relative phase centre offsets. This command is used to
determine which block of GPS satellites (for example GPS 2a or 2b) each GPS
satellite is from. This allows the block specific satellite phase centre offset to be
applied.
␣absolute_antenna␣<true | false>¶
This command is used to enable absolute antenna offsets if it is set to true. If abso-
lute antenna offsets are to be used, an ANTEX format file containing the required
antenna offsets must be specified using the antex_file command. The receiver
antenna is read from the header of the RINEX observation file.
This command specifies the filename of the ANTEX file containing antenna phase
centre offset variation data. For consistency, the same set of antenna offsets should
be used as those used in the generation of the orbit and clock data.
This optional command specifies the filename of the CC2NONCC format, P1C1 bias
file. This file is only required if the RINEX observation file does not include P1
observations, or if the RINEX file contains both P1 and C1 observations and it is
desired to use the C1 observation as a fall back for missing P1 observations.
␣code_type␣<C1 | P1>¶
192
Appendix A. Software reference
This optional command specifies what should be done when a cycle slip is detected
but is not able to be fixed. The options are described here:
• If the option noise is chosen, then a large amount of process noise will be
added to the diagonal element of the a priori covariance matrix corresponding
to the ambiguity for that satellite,
• If the reset option is chosen ambiguity is reinitialised from scratch using the
same method as if the satellite had just appeared above the horizon. All
correlations in the state vector covariance matrix are zeroed.
If this command is not given, then the reset option is used by default.
␣start_time␣<YYYY>␣<MM>␣<DD>␣<HH>␣<mm>␣<SS.S>¶
This command specifies the start time for processing observations. Observations
in the RINEX file before this point will be skipped. <YYYY> is the four digit year,
<MM> is the two digit month, <DD> is the two digit day of month, <HH> is the two
digit hours, <mm> is the two digit minutes, and <SS.S> is a floating point value for
seconds.
␣stop_time␣<YYYY>␣<MM>␣<DD>␣<HH>␣<mm>␣<SS.S>¶
This command specifies the end time for processing observations. Observations in
the RINEX file after this point will not be used. <YYYY> is the four digit year,
<MM> is the two digit month, <DD> is the two digit day of month, <HH> is the two
digit hours, <mm> is the two digit minutes, and <SS.S> is a floating point value for
seconds.
193
Appendix A. Software reference
A.3.1 [Link]
This file contains the modelled and estimated zenith tropospheric delay. The format
is as follows:
where:
Est ZWD is the estimated correction to the zenith wet delay in metres.
Sigma Est ZWD is the formal error of the estimated correction to the zenith wet
delay, one sigma.
53405.00␣0.0519␣2.2908␣0.0761␣0.0073¶
A.3.2 [Link]
This set of files contain the state of the Melbourne-Wübenna cycle slip detector
(Section 4.8.2). One file is generated for each satellite with the format mw_SXX.out
with SXX replaced by the satellite number (eg G23) . The file format is as follows:
where:
hbδ i is the running average estimate of the wide-lane bias (Equation 4.53),
σhbδ i is the approximate standard deviation of the running average (Equation 4.54).
354229978.0␣-24.47839984␣-24.45920839␣0.6682412131¶
194
Appendix A. Software reference
A.3.3 [Link]
This file records outlying observations detected using the method described in Sec-
tion 4.8.1. It specifies the observation count at which an outlier was detected, the
satellite number and GNSS constellation, and the number of consecutive times the
observation from that satellite has been detected as an outlier.
A.3.4 [Link]
This file contains a record of all the cycle slips detected by the two different meth-
ods as described in Section 4.8.2. An example is given below, with a Melbourne-
Wübenna cycle slip detected for GPS satellite 24 at time 33,170 seconds since start
of the day.
A.3.5 lc_residuals.out
This file contains the LC carrier phase residuals at each epoch. Column 1 contains
the time in seconds since the start of the day, columns 2 to 33 give the GPS LC
residuals in metres for satellites 1 to 32 in order. If a satellite is not visible at that
epoch, NaN is printed instead of a number. Columns 34 to 65 contain the GLONASS
LC residuals in a similar manner.
A.3.6 pc_residuals.out
In a similar manner to the lc_residuals file, this file contains the pseudorange iono-
sphere free residuals. The format is the same as that described for the lc_residuals
file.
195
Appendix A. Software reference
A.3.7 [Link]
This file contains the North, East and Up difference (in metres) between the esti-
mated PPP coordinates at each epoch and the initial coordinates for the site given
in the command file. This is useful in static positioning when it contains the error
in the PPP solution with respect to the known coordinates. The format of the file
is as follows:
33179.0␣-1433.0455␣-1039.5154␣1.2711¶
A.3.8 [Link]
This file contains the dilution of precision given the used observations at each epoch.
Satellites below the elevation cutoff or observations detected as outliers are not
included in the calculation of dilution of precision. The entries are as follows:
<PDOP>␣<HDOP>␣<VDOP>¶
where PDOP is the positional dilution of precision, HDOP is the horizontal dilution of
precision and VDOP the vertical dilution of precision.
A.3.9 residuals_by_pos.out
This file contains the carrier phase residuals for each satellite visible at every epoch,
including the elevation and azimuth angles of the satellite at that epoch.
<YY>␣<DD>␣<MM>␣<HH>␣<mm>␣<SS.S>␣<Sat>␣<Az>␣<Elev>␣<LC residual>¶
The <YY> is the two digit year, <MM> is the two digit month, <DD> is the two digit day
of month, <HH> is the two digit hour, <mm> is the two digit minute, and <SS.S> is a
floating point value for seconds. <Sat> is of the format: one character specifying the
network (G for GPS, R for GLONASS), followed by two digit satellite ID number.
<Az> and <Elev> are given in degrees, and <Az> is counted clockwise from North
with a range 0 to 360 degrees. <LC residual> is given in metres.
196
Appendix A. Software reference
11␣3␣24␣9␣12␣51.0␣R02␣242.15␣77.20␣2.37224480947872962E-003¶
A.3.10 slip_count.out
This file contains information about the cycle slips detected and how many could
be repaired.
<epoch>␣<#WL>␣<#NL>␣<#Sats>␣<#WL_fix>␣<#NL_fix>␣<#Sats_fix>¶
where <epoch> is simply the integer tally of observation epochs processed so far,
<#WL> is the number of wide-lane cycle slips detected, <#NL> is the number of
narrow-lane cycle slips detected, <#Sats> is the number of satellites with cycle slips,
<#WL_fix> is the number of wide-lane cycle slips that could be fixed, <#NL_fix> is
the number of narrow lane cycle slips that could be fixed and <#Sats_fix> is the
number of satellites for which both the wide-lane and narrow-lane cycle slips could
be fixed.
353␣1␣1␣1␣0␣0␣0¶
A.3.11 [Link]
This file contains the entire estimated state vector at every epoch. Column 1 is the
seconds since the start of the day. The order of the subsequent columns is the same
as that specified in the state vector in Section 4.3 . The final column contains the
number of observations used to produce the estimate at that epoch. The units of
all the estimated states are in metres.
A.3.12 [Link]
This file contains the ECEF estimated coordinates of the receiver at each epoch in
the reference frame defined by the orbit and clock products provided. The format
of the file is as follows:
197
Appendix A. Software reference
where <Est ZWD> is the estimated correction to the zenith wet delay, <dt> is the
receiver clock offset and <GLONASS dt> is the receiver GPS-GLONASS hardware
bias. The units of all values apart from the <seconds of day> are in metres.
179.0␣3665815.4628␣-92718.0128␣5201177.2662␣-1.783E-003␣-13.0870␣0.0¶
where <command file> is the path and name of a command file. The command
file format is described in Section A.2. Any relative paths in the command file
are evaluated as relative to the current working directory not the location of the
command file. All output is created in the current working directory.
A.5 Scripts
To aid setup and help with running repetitive tasks a number of Python scripts
have been created. Some of the scripts require teqc (Estey and Meertens, 1999) for
RINEX file concatenation and crz2rnx for compressed RINEX file decompression
(Hatanaka, 2008).
A.5.1 [Link]
This script can be used to download files from the CDDIS archive. The script
downloads the given file type for a chosen year and day of year, decompresses it
if necessary and then prints the location of the downloaded file to stdout. For all
options the placeholder <year> should be replaced with the 4 digit year and <doy>
replaced with the day of year for which the data is required. The types of file that
can be downloaded and the command line format for that file type are:
198
Appendix A. Software reference
A.5.2 [Link]
A.5.3 run_directories.py
This script is designed to run every pppncl command file in a directory tree using a
Condor cluster to parallelise the execution of multiple command files. A command
file is taken to be any file with extension “.cmd”. It is called as:
199
Appendix A. Software reference
command file found the given executable will be run with the command file passed
as a command line argument to the executable. The executable will be invoked in
the directory containing the command file. For running on Condor the command
file and all input files mentioned in the command file must be visible to the cluster.
-l, --local
Run the executable on the local machine rather than using Condor.
-b, --batch
Submit to Condor in batches of 3000. This is useful if the number of jobs
submitted overwhelms the Condor dispatcher.
A.5.4 run_files.py
This script runs an executable multiple times, once for each command file specified.
The syntax is:
where EXE is the path to a copy of pppncl and FILES is a space separated list of one
or more command files. The output for each command file is placed in the directory
containing the command file so each command file should be in a separate directory.
A.5.5 setup_rinex.py
This program creates a basic pppncl command file given a RINEX observation file.
It also downloads all the necessary navigation data. This includes an SP3 orbit file,
high rate clock corrections, antenna phase centre offsets in ANTEX format, P1C1
bias data and GLONASS broadcast navigation data if required. It is called using:
where rinexfile is the GNSS observation file for which the pppncl configuration
file is to be generated.
200
Appendix A. Software reference
--data=<archive>
The location into which to download the navigation data. If the required files
are already present in the archive they will not be re-downloaded.
--center=<provider>
The navigation products provider to use, options are:
esa ESA products which include GLONASS orbit and clock corrections
cod CODE products which offer clock corrections tabulated every 5 seconds
The initial coordinates of the site are set to the values in the header of the RINEX
observation file.
A.5.6 setup_site.py
This program creates a basic pppncl command file for static positioning for any sta-
tion in the IGS network. It includes downloading and decompressing the observation
data and all navigation data required. The script looks in three places for a set of
initial coordinates for the site. The first two use the given station name to look up
the coordinates within the given file. In priority order the locations are: the .ssc
SINEX summary file for the chosen navigation data provider, the RINEX clock file
header from the chosen navigation data provider and finally the coordinates in the
RINEX header file.
-y <year>, --year=<year>
The year to use.
-d <doy>, --doy=<doy>
The day of year to use.
201
Appendix A. Software reference
--highrate
Download and use 1 Hz observation data. The individual hourly high rate files
are downloaded and concatenated using teqc (Estey and Meertens, 1999).
--data=<archive>
The location into which to download the navigation data. If the required files
are already present in the archive they will not be re-downloaded.
--center=<provider>
The navigation products provider to use, options are:
esa ESA products which include GLONASS orbit and clock corrections
cod CODE products which offer clock corrections tabulated every 5 seconds
A.5.7 multi_position.py
This program creates a set of pppncl command files for a list of IGS stations over a
list of days. With the default settings it:
• downloads and decompresses all the required observation and navigation data,
30 s observation data is used;
• creates a three day SP3 file spanning the requested day to avoid orbit inter-
polation boundary effects at the start and end of the 24 h period;
• creates a directory named after the current time stamp in which all the output
is placed;
• creates a pppncl command file for each site and day with settings for a 24 h
static PPP solution;
202
Appendix A. Software reference
• runs all of the command files using the pppncl executable version found in the
path using the Condor compute cluster if it is not very busy, otherwise run
pppncl locally;
• generates a set of plots and statistics on the coordinate accuracy and carrier
phase residuals.
--site-file=<sitefile>
A file containing a list of four character site names, one per line. Lines starting
with ’#’ are skipped as comment lines.
-y <year>, --year=<year>
The year to use
-d <doy>, --doy=<doy>
The day of year to use. This may consist of a comma separated list of values
to specify multiple days. Each day listed is setup as a separate PPP run. If
a range of days is wanted, this may be specified using the format begin-end.
The following example shows how to process days 1, 3 through 10 inclusive
and 27 as individual 24 hour sessions:
--doy=1,3-10,27
--data=<archive>
The location into which to download the navigation data. If the required files
are already present in the archive they will not be re-downloaded.
--center=<provider>
The navigation products provider to use, options are:
esa ESA products which include GLONASS orbit and clock corrections
cod CODE products which offer clock corrections tabulated every 5 seconds
-o <dir>, --o=<dir>
Directory in which to place the output. Use in conjunction with --no-datedir
to place the output directly in <dir>. Note this directory should be visible
from the Condor cluster.
203
Appendix A. Software reference
--ppp <exe>
The location of the pppncl executable to use. Note this should be a path that
is accessible from the Condor cluster.
--continue-on-error
Continue running even if it was not possible to setup one of the requested
sites.
--single-sp3
Do not create a three day spanning SP3 file, use only the SP3 file for the
relevant day.
--setup-only
Stop after creating the command files and downloading all the required data.
--no-datedir
Do not create an extra sub directory with the current date containing the out-
put, place the output directly in the current directory or the location specified
in --output.
--highrate
Download and use 1 Hz observation data. Creates a single file for each day
using teqc.
The command file template is a command file with certain values set as placeholders.
The placeholders are then replaced with the relevant value by the multi_position.py
program. The format is based on the Python [Link] functionality.
The following placeholders will be replaced by the values specific to each site and
day. The values that can be substituted are:
$station
The name of the GNSS station.
$rinex_file
The full filename of the RINEX observation file.
204
Appendix A. Software reference
$combined_sp3
The full filename of the SP3 navigation file.
$clk
The full filename of the high rate RINEX clock file.
$x, $y, $z
The receiver coordinates from the relevant IGS analysis centre weekly SINEX
file.
$absolute
The correct value for the ABSOLUTE_ANTENNA command based on the selected
analysis centre and date.
$antex_file
The full filename of the latest available ANTEX antenna phase centre file.
$block_file
The full filename of the latest available GAMIT satellite information file
[Link].
$p1c1bias_file
The full filename of the latest available P1C1 bias file.
$gnav_file
The full filename of the GLONASS navigation file.
As many or as few as desired of the above placeholders may be placed in the template
file. If a placeholder appears multiple time all copies of it will be substituted.
205
Appendix A. Software reference
␣obs_file¶
␣␣$station␣$rinex_file¶
¶
␣nav_file␣$combined_sp3␣SP3¶
¶
␣clk_file␣$clk¶
¶
␣site_pos¶
␣␣$station␣$x␣$y␣$z¶
¶
␣site_stats¶
␣␣$station␣20.0␣20.0␣20.0␣20.0␣20.0␣20.0¶
¶
␣clk_stats¶
␣␣all␣60000␣1000¶
¶
␣data_type␣KH¶
¶
␣cut_off␣5¶
¶
␣gnss_type␣G¶
¶
␣absolute_antenna␣true¶
¶
␣antex_file␣$antex_file¶
¶
␣code_type␣C1¶
The software program pppncl developed as part of this thesis is written largely in
Fortran 90, with certain reused components in Fortran 77. Each set of functionality
is packed in a Fortran module. Other modules or functions that make use of a
particular module must explicitly state this with the use statement. This allows
inter-module dependencies to be easily discovered when reading the source code.
The structure of the program and the major modules and key functions are described
below. Except where stated these modules and functions consist of new code written
206
Appendix A. Software reference
The command strings for the command file are specified in batch_file_cmd.txt.
This file is converted into a Fortran 90 source file batch_file_cmd.f90 by the Perl
script create_batch_cmds.pl.
batch_file.f90 actually performs the reading of the command file, matching com-
mands in the command file to those given in batch_file_cmd.txt. It then reads
the data associated with that command and sets the relevant variables in the module
batch_file_data.
The module driver contains the main loop which forms the core of the software. The
subroutine loops over all the epochs in an observation file performing the following
operation for each observation:
• Using calc_model_data:
207
Appendix A. Software reference
– Model observations
Rinex reading and writing functionality is contained within the Fortran module
rinexlib in libraries/rinex.
Observation files
The following types and functions are used to open a RINEX observation file and
read the information in the header:
Data types
rinex_file_t:
Represents a RINEX observation file. Encapsulates the path and Fortran unit
to a RINEX observation file as well as the RINEX header information for that
file.
Functions
open_obs_file():
Opens a RINEX observation file for reading.
208
Appendix A. Software reference
rinex_file_data_type_present():
Returns .true. if a given observation type is present (e.g. “L1”) in the RINEX
file.
rinex_cc2noncc_applied():
Returns .true. if P1C1 biases have already been applied to this RINEX ob-
servation file using CC2NONCC.
antenna_delta():
Returns the antenna offset from the RINEX observation header. Values are
Height, East and North displacements of the base of the antenna relative to
the monument in metres.
antenna_type():
Returns the antenna type from the RINEX observation header.
approx_pos():
Returns the approximate receiver coordinates from the RINEX observation
header.
Once the file has been opened, the observations in the RINEX file are read in one
at a time and queried using the following:
Data types
observation:
Represents a single epoch in a RINEX observation file.
Functions
zero():
Clears an observation.
obs_sats():
Returns an array of the satellites present in an observation.
get_next_obs():
Reads the next observation in a RINEX file. Returns .true. on success,
.false. if there are no more observations to read. The satellite networks for
which observations will be read in are selected by passing a string containing
the RINEX format letter for each of the required networks. For example the
string “GR” would read in GPS and GLONASS observations, other satellite
networks present in the file would be skipped.
209
Appendix A. Software reference
read_data_1 and read_data_2 perform the actual reading and parsing of the
RINEX observation.
get_data_type():
Fills a two dimensional array with the requested data type e.g. “L1”, “P2” for
all satellites in an observation. The first dimension of the array is the satellite
corresponding to the satellites in the observation returned by obs_sats(). The
second dimension is the requested data types in the order they were passed.
data_type_present():
Similar to get_data_type() but returns an array of logical values indicating
if each of the observation types was present for each of the satellites in the
given observation epoch.
Data types
nav_record:
A data structure containing all the values from a RINEX GPS navigation data
entry.
Functions
open_nav_file():
Opens a given RINEX GPS navigation file, reads the header and the data
entries.
read_header():
Reads header information from a RINEX GPS navigation file.
read_data():
Reads the data records of a RINEX GPS navigation file.
get_nav_for_prn():
Used by get_sat_clk_offset and get_sat_pos to find the broadcast ephemeris
data for a given GPS satellite at a given time.
get_sat_clk_offset():
Calculates the clock offset for a given GPS satellite at a given time using the
broadcast ephemeris. Clock offset is relative to GPS time.
210
Appendix A. Software reference
get_sat_pos():
Calculates the coordinates for a given GPS satellite at a given time using the
broadcast ephemeris. The coordinates are x, y, z in metres in the WGS84
Earth centred reference frame and may be Earth fixed or inertial depending
on the chosen option.
Data types
nav_g_t:
Contains a set of GLONASS broadcast navigation records.
g_nav_record_t:
A GLONASS broadcast navigation record.
Functions
open_nav_g_file():
Opens a given RINEX GLONASS navigation file, reads the header and data
entries.
get_sat_slot():
Looks up the slot number for a particular satellite number at a given time.
The slot number can then be used to determine the transmitting frequencies
of the satellite. The returned value corresponds to k in Equations 2.1 and 2.2.
Internal functions
read_header():
Reads header information from a RINEX GLONASS navigation file.
read_record():
Reads an individual data record of a RINEX GLONASS navigation file into a
g_nav_record_t.
store_record():
Stores a record in a nav_g_t.
211
Appendix A. Software reference
The module sp3 contains functions for parsing SP3 format (Hilla, 2002) files and in-
terpolating the data to calculate satellite coordinates and clock offsets. The compile
time constants max_sat, max_systems and max_sp3 defined in the module config
affect the size of the arrays used to store SP3 data.
When interpolating the tabulated SP3 values there must be at least two entries
before and two entries after the required time, to avoid the interpolation errors
that appear at the edge of the interpolation interval known as Runge’s phenomenon
(Dahlquist and Björck, 2003).
Data types
sp3_header_t:
Stores the information contained in an SP3 file header.
sp3_data_t:
Stores the tabulated satellite coordinates and clock offsets from an SP3 file.
Functions
read_header():
Parses an SP3 file header.
read_sp3():
Opens and parses an SP3 file, multiple satellite systems are supported.
eph_to_xyz():
Calculates the coordinates for a given satellite in either an ECEF or Earth
Centred Inertial (ECI) coordinate system at a given time. Lagrange interpo-
lation is used to interpolate the tabulated values to the required time (Press
et al., 1996).
eph_to_xyz_v():
As eph_to_xyz but also returns the satellite velocity by evaluating the deriva-
tive of the interpolating Lagrange polynomial.
comp_svs_clk():
Calculates receiver clock offset for a particular satellite at a given time.
212
Appendix A. Software reference
Internal functions
find_index():
Finds the index in the SP3 data of start and end points for an interpolation
of a given order for a particular satellite at a particular time. If possible the
returned range will evenly span the required time. Missing data or proximity
to the start and end points of the parsed SP3 data will cause the returned
range to not be centred on the required time. If the requested time is outside
the data available or the found span would result in interpolation near the
edge of the span (resulting in problems as described above), a start index of
−1 is returned.
The Fortran module clk contains functions for reading and interpolating RINEX 3.0
clock files (Ray and Gurtner, 2010). The compile time constants max_clk, max_sat
and max_systems from the module config affect the size of the arrays used for
storing data read in from a RINEX clock file.
Data types
clk_data_t:
Stores the parsed clock offset data.
Functions
read_clk():
Opens and reads a RINEX 3.0 clock format file.
comp_svs_clk():
Linearly interpolates the clock offset for a given satellite at a particular time.
The offset is given in seconds.
The module navigation_data provides an abstraction over the used satellite naviga-
tion type. The underlying navigation data may be a broadcast ephemeris in RINEX
navigation format or a precise ephemeris in SP3 format, possibly augmented with
higher rate satellite clock corrections in RINEX clock format.
213
Appendix A. Software reference
Functions
setup_navigation_data():
Reads in and parses the selected navigation files.
load_glonass_nav():
Loads an additional RINEX GLONASS navigation file to obtain the mapping
from satellite number to GLONASS frequency number.
get_sat_pos():
Calculates the position (in metres X, Y, Z) of a given satellite at a particular
time using the corresponding method for the underlying navigation type. The
returned coordinates are in the reference frame of the provided navigation
data.
get_sat_pos_v():
Calculates the position (in metres X, Y, Z) and velocity (in m s−1 ) of a given
satellite at a particular time using the corresponding method for the underlying
navigation type.
get_sat_clk_offset():
Calculates the clock offset in seconds of a given satellite at a particular time
using the corresponding method from the underlying navigation type.
get_sat_frequency():
Returns the two carrier frequencies of the given satellite.
The module cc2noncc provides functions to read and apply differential code biases.
The file parsing uses code from CC2NONCC v6.3 provided by ESOC (Romero, 2010).
Data types
cc2noncc_t:
Stores a set of DCB.
Functions
cc2noncc_new():
Reads in a set of DCB biases in the CC2NONCC format (Section A.1.4) for a
given month.
214
Appendix A. Software reference
cc2noncc_applyc1p1():
Applies a set of C1-P1 biases to a set of C1 observations given the observations
and the corresponding satellites (Equation 3.30).
The module glonass_f implements Equations 2.1 and 2.2 to determine the GLONASS
frequency for a given channel number.
Functions
glonass_frequency_for():
Calculates the L1 and L2 carrier frequencies for a given GLONASS channel
number.
Data types
sat_phase_center_t:
Used to store the satellite phase centre data.
rcv_antenna_t:
Used to store the receiver phase centre offsets and variations.
Functions
phase_center_from_block_list():
Loads satellite phase centre offset block mean values in the GAMIT [Link]
format into a sat_phase_center_t.
phase_center_from_antex():
Loads satellite phase centre offset block mean values from an ANTEX file into
a sat_phase_center_t.
215
Appendix A. Software reference
get_sat_phase_centre_offset():
Gets the phase centre offset for a particular satellite at a given time given the
satellite coordinates. The satellite body frame offset is converted to ECEF
coordinates at the given time (Section 3.7). The units of the offset are metres.
phase_center_variation():
Returns the satellite phase centre variation given a satellite and the nadir
angle it is viewed from. The returned value is in metres.
rcv_antenna_new_manual():
Initialises a rcv_antenna_t using a given phase centre offset for L1 and L2
carrier frequencies.
rcv_antenna_new_antex():
Initialises a rcv_antenna_t using antenna phase centre offset and variation
values from an antenna in an ANTEX file. The antenna can be searched for
in an antex file using find_matching_rcv_antenna from the module antex.
rcv_antenna_offset():
Gets the receiver antenna offset in metres for the L1 and L2 carrier frequencies.
rcv_antenna_variation():
Gets the receiver antenna variation at a given elevation angle in metres for the
L1 and L2 carrier frequencies.
The modules antex and svs contain analogous functions to those given above that
implement the loading, look-up and linear interpolation of the values for a particular
format.
The module theory contains functions to model the expected GNSS observations
(Section 4.4), calculate the partial derivatives with respect to the state vector (Sec-
tion 4.5), determine the observation stochastics (Section 4.6) and compute the prefit
residual vector (Equation 4.40).
Data types
theory_t:
Contains modelling options such as the antenna phase centre model and the
atmospheric model to use.
216
Appendix A. Software reference
model_out_t:
Contains the output of the modelling including the prefit residuals, the obser-
vation partial derivatives and the observation stochastics for a given observa-
tion combination for a particular satellite.
Functions
calc_model_data():
Given a predicted state, a set of observations and the time of observation,
models the requested observation combination filling a model_out_t for each
satellite. These steps correspond to Equations 4.4 and 4.7 of the Kalman filter.
Internal functions
predict_ranges():
Computes the predicted code and carrier phase observations (Equations 4.15
and 4.16) on the L1 and L2 frequencies for a set of GPS and GLONASS satel-
lites at a given time. A state vector (Section 4.3) of receiver coordinates,
receiver clock offset, wet tropospheric delay correction and estimated ambigu-
ities are needed as inputs.
get_transmit_time():
The transmit time from each satellite is found by iteration as described in
Section 4.4 so the computation may be performed without requiring any real
observations. The satellite coordinates and clock offset are found using the
selected navigation data (Section A.6.6).
o_minus_c():
Calculates the observed minus computed measurements for the chosen set of
observation combinations e.g. LC and PC. Also computes the partial deriva-
tives of the computed combinations with respect to the state vector (Sec-
tion 4.3) and the modelled variance of the observation combinations.
form_partials():
Calculates the partial derivatives of a range measurement to a satellite with
respect to the receiver coordinates and the tropospheric correction (Equations
4.25 and 4.26).
fill_partials():
Fills in the partial derivatives array of a model_out_t data type.
217
Appendix A. Software reference
check_millisecond_jump():
Given a set of prefit residuals, determines if the receiver has inserted a mil-
lisecond jump in its clock. Certain receivers, for example the Leica GS10,
compensate for internal clock drift by inserting millisecond jumps whenever
the receiver clock is offset from GPS time by more than half a millisecond.
Failure to detect this can lead to errors in calculating the satellite position
and clock offset.
Tropospheric modelling
Data types
atm_zenith_model_t:
An atmospheric zenith model.
atm_mapping_funtion_t:
An atmospheric mapping function.
Functions
model_zenith():
Models the tropospheric zenith path delay using the Saastamoinen zenith delay
model at a given time, latitude and height (see Section 3.3, Equations 3.4 and
3.3).
atm_zenith_model_hydrostatic_delay():
Gets the value of the modelled hydrostatic delay from an atm_zenith_model_t
in metres.
atm_zenith_model_wet_delay():
Gets the value of the modelled wet delay from an atm_zenith_model_t in
metres.
atm_mapping_init():
Initialises an atm_mapping_funtion_t to use either:
The code for the Niell mapping functions was used from Track, the GMF
mapping function uses the routine provided at [Link]
218
Appendix A. Software reference
at/~ecmwf1/gmf.f_hu.
atm_mapping_function_update():
Updates the mapping function state based on the time, receiver coordinates
and satellite elevation.
atm_mapping_function_hydrostatic():
Returns the value of the hydrostatic mapping function.
atm_mapping_function_wet():
Returns the value of the wet mapping function.
tropospheric_delay():
Calculates the total tropospheric path delay due to the modelled wet and dry
zenith delays and the wet and dry mapping functions. A correction to the
zenith wet delay may be passed in (see Equation 3.12).
The module svs contains models for the nominal yaw attitude that allow phase
windup determination and mapping satellite body frame antenna offsets to an ECEF
frame.
Functions
SUN20():
Calculates a unit vector and range to the Sun at a given time. The code from
Track was used which is based on the tables of the Sun by Newcomb (1895).
windup():
Calculates the phase windup in radians for a given satellite and receiver loca-
tion and time (Section 3.8). Keeps track of complete rotations to avoid adding
single cycle discontinuities.
Relativistic correction
The module relativity contains functions that calculate the relativistic corrections
due to orbital eccentricity and gravitational bending.
219
Appendix A. Software reference
Functions
relativity_eccentricity_correction():
Calculates the relativistic correction due to satellite orbit eccentricity (Equa-
tion 3.18).
relativity_gravitational_correction():
Calculates the propagation correction for a satellite signal due to space-time
curvature induced by the Earth’s gravitational field (Equation 3.19).
The module earth_tide contains routines from Track for calculating the distortion
at a given point and time due to solid Earth tides (Section 3.10).
Functions
earth_tide_dXYZ():
Calculates the X, Y, Z coordinate shift in metres due to solid Earth tides at a
given position and time.
Data types
mwdata_t:
Stores the values of the running mean and variance from Equations 4.53 and
4.54 for a single satellite.
mw_detector_t:
Contains a set of mwdata_t, one for each satellite and the chosen number of
standard deviations outside of which a cycle slip is detected.
csdetector_ion_history_t:
Contains ΦI (tk−1 ) and tk−1 for use in Equation 4.56 at the next epoch.
220
Appendix A. Software reference
csdetector_ion_t:
Contains a set of csdetector_ion_history_t, one for each satellite and the
absolute ionosphere combination rate of change outside of which a cycle slip
is detected.
Functions
csdetector_mw_new():
Initialises a mw_detector_t with a chosen cycle slip detection threshold in
units of number of standard deviations.
csdetector_mw_detect():
Given a set of carrier phase observations, returns an array of logical values
indicating if a cycle slip was detected for each satellite using the Melbourne-
Wübbena combination method described in Section 4.8.2.
csdetector_mw_force_slip():
If a cycle slip has been detected using some other method which would in-
fluence the wide-lane bias, this function should be used to reset the running
average used in the cycle slip detection.
csdetector_ion_new():
Initialises a csdetector_ion_t with a chosen detection threshold given in
metres.
csdetector_ion_detect():
Given a set of carrier phase observations, returns an array of logical values
indicating if a cycle slip was detected for each satellite using the ionosphere
combination method described in Section 4.8.2.
The module slip_fix implements cycle slip fixing using the integer least squares
time differenced method of Banville and Langley (2010) as described in Section 4.8.2.
Functions
fix_slips():
Attempts to fix the cycle slips previously identified.
ls():
A wrapper around the LAPACK least squares function DGELS that performs
least squares estimation using QR factorisation (Anderson et al., 1999). Used
by fix_slips to perform the initial cycle slip estimation.
221
Appendix A. Software reference
The module lambda is a Fortran 90 interface wrapper to the LAMBDA integer least
squares implementation available from Delft University of Technology. The function
lambda4 returns a set of integer valued estimates given a set of float estimates and
their associated covariance matrix. All calculations are performed using the Delft
code; lambda4 makes the use of the method simpler by allocating the space required
for work vectors.
The module blunder implements the outlier detection method described in Sec-
tion 4.8.1.
Functions
LOM_test():
Given a set of prefit residuals and the prefit residual covariances found using
residual_covariance, first identifies if the set of prefit observation residuals
are consistent with their given covariance and if so identifies the most likely
outlier. If the most likely outlier does not fail the local test (Equation 4.46),
−1 is returned as the index of the most likely outlier.
log_failure():
Used to add an entry to [Link] when an observation is detected as an
outlier.
If they are available the LAPACK routines dgetrf and dgetri are used to perform
matrix inversion using LU decomposition (Anderson et al., 1999).
Functions
ekf():
Implements the extended Kalman filter update step given the predicted state
222
Appendix A. Software reference
[1]
x̂k (−), prefit residuals (zk − ẑk ), linearised system dynamic model Hk , obser-
vation covariance matrix Rk and a priori state covariance matrix Pk (−).
Elements in the state vector for which the corresponding row in the a priori
state covariance matrix Pk (−) is zero will not be changed by the update step
(Equations 4.8 and 4.9). To avoid performing multiple matrix multiplications
and inversions on unnecessarily large matrices, the rows of Pk (−) are checked
at the start of this subroutine. Using Fortran 90 array indexing only those
rows and columns of the input vectors and matrices are selected and passed
to ekf_ that actually performs the filtering.
residual_covariance():
Calculates the expected uncertainties in the pre-fit observation residuals given
the observation stochastic model, the design matrix and the predicted state
vector covariance. Returns Cvk from Equation 4.42.
Internal functions
ekf_():
Called by ekf to perform the filtering once the active elements of the input
matrices have been extracted.
kalman_gain():
Calculates the Kalman gain matrix (Equation 4.9).
condition_on_measurement():
Adjusts the predicted state based on the measured observations (Equation 4.8).
a_posteriori_covariance():
Calculates the a posteriori covariance matrix (Equation 4.10).
a_posteriori_covariance_joseph():
Calculates the a posteriori covariance matrix using the alternative Joseph form
(Equation 4.14).
223
Appendix B IGS sites and reference coordinates
The below table lists the IGS sites and their coordinates, for GPS week 1605, used
in static PPP accuracy assessment. Coordinates are given in ITRF2005 using the
GRS80 ellipsoid.
224
Appendix B. IGS sites and reference coordinates
225
Appendix B. IGS sites and reference coordinates
226
Appendix B. IGS sites and reference coordinates
227
Appendix B. IGS sites and reference coordinates
228
Appendix B. IGS sites and reference coordinates
229
Appendix B. IGS sites and reference coordinates
230
Appendix B. IGS sites and reference coordinates
231
Appendix B. IGS sites and reference coordinates
Table B.1: List of sites and their coordinates, for GPS week 1605, used
in static PPP validation. Coordinates are given in ITRF2005 using the
GRS80 ellipsoid
232









