GPS/INS Operation in Shadowed Environments
Except where reference is made to the work of others, the work described in this thesis is
my own or was done in collaboration with my advisory committee. This thesis does not
include proprietary or classi ed information.
Benjamin J. Clark
Certi cate of Approval:
George T. Flowers
Professor
Mechanical Engineering
David M. Bevly, Chair
Assistant Professor
Mechanical Engineering
Andrew J. Sinclair
Assistant Professor
Aerospace Engineering
Joe F. Pittman
Interim Dean
Graduate School
GPS/INS Operation in Shadowed Environments
Benjamin J. Clark
A Thesis
Submitted to
the Graduate Faculty of
Auburn University
in Partial Ful llment of the
Requirements for the
Degree of
Master of Science
Auburn, Alabama
August 9, 2008
GPS/INS Operation in Shadowed Environments
Benjamin J. Clark
Permission is granted to Auburn University to make copies of this thesis at its
discretion, upon the request of individuals or institutions and at
their expense. The author reserves all publication rights.
Signature of Author
Date of Graduation
iii
Vita
Benjamin J. Clark was born in Athens, AL in February 1982. After living with his
parents in Neil, AL, the Clark family returned to Athens where his parents still live. After
graduating from East Limestone High School in 2000, he began attending Freed-Hardeman
University where he met Follin Childs. Benjamin was awarded a B.S. in Mathematics in 2004
from Freed-Hardeman University. Just after their graduation, the couple were married and
moved to Opelika, AL where Benjamin began work toward a B.S. in Mechanical Engineering
at Auburn University. During these two years, he worked with Dr. David Bevly in the
GPS and Vehicle Dynamics Laboratory. Upon completion of his undergraduate education,
Benjamin entered graduate school in the summer of 2006.
iv
Thesis Abstract
GPS/INS Operation in Shadowed Environments
Benjamin J. Clark
Master of Science, August 9, 2008
(B.S., Auburn University, 2006)
(B.S., Freed-Hardeman University, 2004)
111 Typed Pages
Directed by David Bevly
This thesis presents the analysis techniques developed for monitoring GPS signals in
harsh shadowed environments such as heavy foliage. It also details a method selected
for improved performance by combining raw GPS information with an Inertial Navigation
System (INS). Normal GPS operation in shadowed areas su er from position jumps of tens
to hundreds of meters. The developed analysis reveals that these errors are due to the
quickly changing local errors that cause a GPS receiver to report erroneous position spikes.
Monitoring variables for the signal strength and change in multipath are employed to keep
track of the environmental e ects on the GPS measurements. A new visualization technique
is also developed to qualitatively monitor the environmental e ects.
From the visualization technique, the e ects of the shadowing environment are shown
to simultaneously a ect the signal strength and multipath. It is shown that foliage cover
causes these e ects to occur spontaneously as a signal travels through and around obsta-
cles. To mitigate these errors, a GPS/INS closely coupled system is implemented which
uses inertial sensors to smooth the erroneous GPS jumps. Introduction of alternative sen-
sors allows for integrity monitoring in the form of GPS outlier measurement rejection so
v
that local environmental e ects can be detected and removed from the navigation solution.
The resulting implementation reveals that it is possible to operate in these GPS harsh en-
vironments without su ering from the position jumps of tens or hundreds of meters. This
implementation allows for navigation in foliage cover comparable to the under ten meter
accuracy of standard GPS in clear environments.
vi
Acknowledgments
First and foremost, I am thankful to God for so many opportunities. I am thankful
for my wife, Follin, who has put me through school here at Auburn University. I want
to acknowledge my family in helping me to make learning a habit and my professors in
equiping me with tools to follow that pursuit.
vii
Style manual or journal used Journal of Approximation Theory (together with the style
known as \aums"). Bibliograpy follows van Leunen?s A Handbook for Scholars.
Computer software used The document preparation package TEX (speci cally LATEX)
together with the departmental style- le aums.sty.
viii
Table of Contents
List of Figures xii
List of Tables xiv
1 Introduction 1
1.1 Objectives and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Prior Art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Outline of Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2 Overview of GPS Position Computation 6
2.1 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Positioning Fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Satellites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3.1 Signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3.2 GPS Control Segment . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3.3 Ephemeris . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4 Positioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3 Overview of GPS Errors 27
3.1 Satellite Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.1.1 Errors in Ephemeris Data . . . . . . . . . . . . . . . . . . . . . . . . 28
3.1.2 Satellite Clock Errors . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.2 Atmospheric Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2.1 Ionospheric Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2.2 Tropospheric Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3 Local Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3.1 Multipath Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3.2 Receiver Measurement Error . . . . . . . . . . . . . . . . . . . . . . 36
3.4 GPS Observable Error Models . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4 Shadowed Environments 38
4.1 Environmental Error Monitoring Parameters . . . . . . . . . . . . . . . . . 38
4.1.1 Signal Attenuation Monitoring Parameter . . . . . . . . . . . . . . . 39
4.1.2 Multipath Monitoring Parameter . . . . . . . . . . . . . . . . . . . . 41
ix
4.1.3 Error Visualization for a Single Satellite . . . . . . . . . . . . . . . . 43
4.1.4 Error Visualization for All Satellites in View . . . . . . . . . . . . . 44
4.1.5 Analysis Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5 Error Mitigation 48
5.1 Di erential Corrections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.1.1 Position Di erencing . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.1.2 Code Di erencing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.1.3 Carrier Smoothed Code Di erencing . . . . . . . . . . . . . . . . . . 54
5.1.4 Static Comparison of Di erential Methods . . . . . . . . . . . . . . . 55
5.2 Inertial Aiding Sensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.2.1 Coupling Architectures . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.2 INS Mechanization Equations . . . . . . . . . . . . . . . . . . . . . . 58
5.2.3 State Vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2.4 State Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.2.5 State Measurement Relations . . . . . . . . . . . . . . . . . . . . . . 63
5.3 Fault Detection and Exclusion . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.3.1 Normalized Innovation Parameter . . . . . . . . . . . . . . . . . . . 68
5.3.2 FDE Threshold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
6 Results 71
6.1 Pseudorange Calculated Position . . . . . . . . . . . . . . . . . . . . . . . . 71
6.2 Closely Coupled Implementation Results . . . . . . . . . . . . . . . . . . . . 75
6.3 Closely Coupled with FDE Implementation Results . . . . . . . . . . . . . . 76
6.4 Comparison of Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7 Conclusion 82
7.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
Bibliography 85
Nomenclature 88
A Coordinate Frames and Transformations 90
A.1 Coordinate Frames . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
A.1.1 Earth-Centered Inertial (ECI) Frame . . . . . . . . . . . . . . . . . . 90
A.1.2 Earth-Centered Earth-Fixed (ECEF) Frame . . . . . . . . . . . . . . 91
A.1.3 Navigation Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
A.1.4 Body Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
A.1.5 Orbit Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
x
A.1.6 Satellite Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
A.2 Coordinate Transformations . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
A.2.1 Coordinate Transformation Matrix . . . . . . . . . . . . . . . . . . . 95
A.3 ECEF and Geodetic Coordinates . . . . . . . . . . . . . . . . . . . . . . . . 96
xi
List of Figures
1.1 NovAtel Reported Position in Shadowed Environment . . . . . . . . . . . . 2
2.1 ECEF User Position . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 GPS Orbital Planes [20] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Orbit Angles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4 Position Using Least Squares Estimate . . . . . . . . . . . . . . . . . . . . . 23
2.5 Position Using Weighted Least Squares Estimate . . . . . . . . . . . . . . . 25
3.1 GPS Ephemeris Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.2 Atmospheric Layers [3] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.3 Typical Klobuchar Ionosphere Model [25] . . . . . . . . . . . . . . . . . . . 33
3.4 Ionospheric Obliquity Factor Model [25] . . . . . . . . . . . . . . . . . . . . 33
3.5 Three Multipath Types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.1 Illustration of GPS Signal Shadowing Due to Foliage . . . . . . . . . . . . . 39
4.2 Static Data Collection Location . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.3 Single Satellite Error Visualization . . . . . . . . . . . . . . . . . . . . . . . 45
4.4 Aerial View of Static Data Location . . . . . . . . . . . . . . . . . . . . . . 46
4.5 All-In-View Satellite Error Visualization . . . . . . . . . . . . . . . . . . . . 46
5.1 Common Satellite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.2 Single Di erence Diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.3 Two Common Satellites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
xii
5.4 Method Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.5 Coupling Architectures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.6 Fault Detection and Exclusion Threshold Selection . . . . . . . . . . . . . . 69
6.1 Dynamic Path . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
6.2 Clear Sky and Heavy Foliage Environments Along Dynamic Path . . . . . . 73
6.3 Pseudorange Calculated Positions . . . . . . . . . . . . . . . . . . . . . . . . 74
6.4 Zoom on Pseudorange Calculated Positions . . . . . . . . . . . . . . . . . . 74
6.5 Zoom on Closely Coupled Positions . . . . . . . . . . . . . . . . . . . . . . . 76
6.6 Comparison of Closely Coupled and Pseudorange Solutions . . . . . . . . . 77
6.7 Comparison of Closely Coupled with FDE and Pseudorange Solutions . . . 77
6.8 Normalized Innovations Over Entire Path . . . . . . . . . . . . . . . . . . . 78
6.9 Normalized Innovations Over Zoom Segment . . . . . . . . . . . . . . . . . 78
6.10 Overlay of Positioning Methods In Heavy Foliage . . . . . . . . . . . . . . . 79
6.11 Overlay of Positioning Methods In Clear Sky . . . . . . . . . . . . . . . . . 80
6.12 Overlay of Positioning Methods In Light Foliage . . . . . . . . . . . . . . . 80
A.1 Earth-Centered Inertial (ECI) Frame . . . . . . . . . . . . . . . . . . . . . . 91
A.2 Earth-Centered Earth-Fixed (ECEF) Frame . . . . . . . . . . . . . . . . . . 92
A.3 Navigation Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
A.4 Body Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
A.5 Orbital Plane With Respect To Equatorial Plane . . . . . . . . . . . . . . . 94
A.6 Orbit Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
xiii
List of Tables
2.1 Navigation Message Subframes . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2 Clock and Ephemeris Parameters . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Receiver Parameters from [24] . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1 GPS Standard Error Budget [27] . . . . . . . . . . . . . . . . . . . . . . . . 27
4.1 Frequency-Lock-Loop Parameters from [36] . . . . . . . . . . . . . . . . . . 40
5.1 System Noise Covariance Matrix Values . . . . . . . . . . . . . . . . . . . . 63
xiv
Chapter 1
Introduction
1.1 Objectives and Motivation
The Global Positioning System (GPS) is gaining increased attention since reaching its
full operational capability (FOC) in 1995 [21]. Although the system was designed for a user
to navigate with a budget under $10,000, user equipment prices have currently dropped
below the $100 mark and are continuing to drop as more capable computing resources are
developed. Currently GPS navigation capabilities are found on aircraft, automobiles, cell
phones, pets, and most any mobile user. These increasing navigation capabilities are leading
to further advances in other  elds. One such  eld is autonomous vehicle control, as can be
seen by the use of GPS in the Darpa Grand Challenges [33]. With the global nature and
long-term stability of GPS, such applications bene t from its use in being able to locate
and orient a user accurately, consistently, and cheaply.
However, even with its prevalence and wide acceptance, there are certain applications
for which stand-alone GPS does not meet the necessary performance requirements. One of
these is navigation in signal-harsh environments such as areas around many tall buildings
(urban canyon) or heavy tree cover. Much work has been done in investigating GPS oper-
ation in urban canyon environments, but signi cantly less has been done in evaluation of
GPS operation in foliage cover. These foliage environments are described in the literature
as shadowed environments and pose some unique di culties to GPS operation [18]. An
example of this operation is shown in Figure 1.1. This plot shows the position solutions
reported by a NovAtel DL-3V GPS receiver in a rural subdivision with moderate foliage
1
cover. As can be observed, inaccurate position jumps of tens of meters can be seen as the
receiver tries to operate in this signal-harsh environment.
Figure 1.1: NovAtel Reported Position in Shadowed Environment
It is noted that such operation could lead to serious failures in a control system that
uses such positions as the total navigation solution. Thus, the objective of this research is
to analyze the cause of this performance degradation as well as ways to lessen its impact.
The goal is a highly reliable system that performs well even in the aforementioned di cult
environments.
1.2 Prior Art
Analysis of GPS operation has been investigated by several sources. Each of these
studies had its own individual focus and target application. The use of GPS in measuring
lengths of movement data was investigated in [10]. Here it was observed that the forest
2
canopy e ects on GPS caused more position jumps and therefore a longer traveled distance
was reported. Positioning errors were attributed to the canopy?s interference with the satel-
lite pseudoranges. Various forest cover types were studied in [38] with respect to their e ect
on GPS accuracy. In [16] a signal interruption probability (SIP) parameter was investigated
in di erent forest conditions which was found to be related to the accuracy of carrier-phase
di erential operation in these environments. The link between integer ambiguity resolution
and the SIP parameter was the main focus while the main e ect monitored was the total
loss of satellite tracking. In [30] the position dilution of precision (PDOP) value was investi-
gated as a monitoring variable for operation in foliage environments. However, it was found
that this geometry-based variable was not a good indicator of performance since a satellite
that provided a good PDOP could be erroneous enough to make a solution worse than if
it were ignored. In work done by Lachapelle, the multipath e ect was studied as the main
error source in foliage environments [22]. This work compares performance in foliage areas
due to narrow or wide correlator spacing. This spacing corresponds to a receiver?s ability
to reject multipath, which is a large error contributor in these environments. A method for
visualizing the multipath as a function of a satellite?s azimuth and elevation was given in
[17]. This technique yields impressive images describing the correlation of multipath over
various regions around the receiver.
In order to operate in the foliage conditions, sensor integration has been investigated.
A method proposed in [5] used a loosely coupled sensor fusion approach with GPS, gyro-
scope, encoders, and compass. This operation was implemented on a system with a human
operator. To provide an increase in accuracy, a type of terrain-referencing was used by the
operator of the robot. Steps were also taken to adjust the weighting on the measurements
3
through the measurement covariance matrix when jumps in the GPS solution occurred using
a two meter threshold for detection.
1.3 Contributions
This thesis develops a detailed analysis of errors experienced by GPS users in shad-
owed environments. These errors include not only loss of tracking on the satellites but
also their degradation even while still being available. This analysis includes the use of
error monitoring variables capable of accurately describing the harsh e ects of the foliage
environment. A new visualization mechanism is also presented that allows for correlations
between these monitoring variables to be qualitatively observed. With this analysis, system
design conclusions are drawn for simplifying the implementation of a system to navigate in
shadowed environments.
Also given are the details of an implemented system which combines GPS/INS cou-
pling with an outlier rejection scheme. The coupling allows for continuous operation and
more consistent results than stand-alone GPS. The rejection algorithm is a method of fault
detection and exclusion (FDE) that increases the robustness of the navigation solution.
1.4 Outline of Thesis
Chapter 2 describes some background in GPS operation used throughout this thesis.
Details are presented of how measurements are generated and how these measurements are
used to calculate user position. Chapter 3 describes common errors that a ect this solution
experienced by all GPS users. These common errors arise from various sources and con-
tribute to the overall error in the GPS measurements. Chapter 4 describes the shadowed
4
environment and the errors speci c to such environments. The error monitoring variables
are described and derived for analyzing the e ectiveness of GPS in these environments.
Then the visualization scheme is described and utilized for the improved navigation system
design. Chapter 5 describes several error mitigation algorithms including di erential opera-
tion and coupling with an inertial navigation system (INS). This chapter also describes the
FDE algorithm used for systematic outlier rejection. Chapter 6 presents the results of the
GPS/INS and FDE system. Finally, conclusions are presented in Chapter 7.
5
Chapter 2
Overview of GPS Position Computation
Before delving into the details of the errors that a ect GPS, some foundational infor-
mation is presented. This chapter describes some of the basics a receiver uses to calculate
a position using GPS observables. This includes the calculation of the satellite?s position
from broadcast ephemeris parameters as well as using these positions with measured ranges
to the satellites. Also described are the various parts of the overall system that make posi-
tioning possible. This includes the satellites and GPS control segment, which oversees the
satellite signals and position in space.
2.1 Nomenclature
In this thesis, a superscript represents the coordinate frame in which the components
are reported. If a variable refers to a measure between two frame origins, a subscript is used
with the frame the measure is from listed  rst and the frame the measure is to listed second.
The frames mentioned in this paper are the inertial frame, i, the Earth-Centered-Earth-
Fixed (ECEF) frame, e, the navigation frame, n, the body frame, b, the orbit frame, o, and
the satellite frame, s. Detailed descriptions of these frames are provided in Appendix A and
[14]. For example, the position of a user is taken to be the center of both the navigation and
the body frames (which coincide but have di erent orientations). Therefore, expression of
the user position in the ECEF frame would be expressed as reeb. This is read as: position of
the body frame origin with respect to the ECEF origin expressed in ECEF coordinates. The
magnitude of a vector quantity is given without a coordinate frame in which the quantity
6
is reported in (no superscript). For example, the range from a user to a satellite is given
as rbs. Vector quantities that do not correspond with a single set of reference frames are
denoted with an overbar. That is, the state vector would be labeled as  x.
2.2 Positioning Fundamentals
GPS is built upon several concepts, the  rst of which is time-ranging. If a radio signal
is traveling at the known speed of light, c, and leaves its satellite source at a known time, ts,
a user can determine the range, rsb, from the source by timing when the signal is measured,
tm, by
rsb = c(tm ts) (2.1)
GPS position solutions can be calculated in an Earth-Centered Earth-Fixed (ECEF)
frame, described in Appendix A. In this coordinate system, a user?s position can be given
by a longitude angle  , a latitude angle  , and a range reb. This situation is illustrated in
Figure 2.1 and details of this transformation are given in Appendix A.
The second concept on which GPS is built is Euclidean distance. To calculate the
geometric range between a user at reeb = (xeb;yeb;zeb) and a source (satellite) at rees =
(xes;yes;zes), a user would apply the Euclidean distance formula
rbs =
q
(xes xeb)2 + (yes yeb)2 + (zes zeb)2 (2.2)
7
Figure 2.1: ECEF User Position
In the case of a GPS user, the range and the satellite position are known. Thus, mul-
tiple ranges are needed from multiple satellites in order to determine the user?s position.
Assuming the time is perfectly known by all satellites and the user, exactly three ranges
are needed to solve for the user?s ECEF position. However, the assumption that the user?s
clock is perfectly synchronized with GPS time used by the satellites is a poor one. Such
systems would require expensive atomic clocks continuously running to keep track of this
time. In order to reduce the required cost of user equipment, a clock bias state is also
estimated along with the user?s position. Thus a receiver needs to solve for four unknowns
at a time to determine the correct position. Due to this fact, the range measurement is
corrupted by the receiver?s clock bias and various errors such as atmospheric delay and is
therefore called a pseudorange measurement,  . For a satellite s, the overall pseudorange
measurement model therefore becomes
 s =
q
(xes xeb)2 + (yes yeb)2 + (zes zeb)2 +c tu + s (2.3)
8
where  tu is the user?s clock o set error, c tu is the user?s clock bias in the units of meters,
and  s is the pseudorange error for satellite s.
2.3 Satellites
There are currently 31 GPS satellites in orbit in 6 orbital planes [34]. The GPS constel-
lation was originally designed with four satellites in each of these orbital planes with slots
for spares [21]. The orbital planes described by the original design are shown in Figure 2.2.
This design decision was to provide global availability and integrity, the details of which can
be found in [20]. The constellation was designed to ensure users have the required minimum
of 4 satellites in view for all users, neglecting local blockages.
Figure 2.2: GPS Orbital Planes [20]
2.3.1 Signals
Each GPS satellite operates as a signal transmitter orbiting the earth every half of
a sidereal day (11 hours 58 minutes) [21]. This implies that a user will see a satellite in
approximately the same location at approximately the same time (4 minutes earlier) each
9
day. At a single full orbit, the observer will be on the opposite side of the Earth. This fact is
an important consideration when taking into account repeatability in GPS experimentation.
Each satellite uses a frequency synthesizer to drive subsystems generating three components
of the signal: carrier, code, and navigation data message.
Carrier
The carrier subsystem generates the sinusoidal signals that center the frequency of
the overall GPS signal. Currently all satellites (Block II/IIA/IIR/IIR-M) broadcast in the
L-band at the L1 (1575.42 MHz) and the L2 (1227.6 MHz) frequencies. The satellites use
other frequencies in the L-band and future implementations will include access on the L5
band (1176.45 MHz) [21]. The main design purpose of the carrier signals was to provide a
means to e ectively transmit the messages in the other subsystems. However, carrier phase
measurements are also the basis of highly accurate positioning techniques when some of the
inherent di culties of the method are overcome. Since a receiver is capable of accurately
measuring the carrier phase to a small fraction (about 1%) of the carrier wavelength, on
the L1 signal this implies an accuracy of approximately 1.9 mm. This measurement is
di cult to use since the whole number of wavelengths between the user and the satellite is
not known. This problem is known as integer ambiguity. Techniques which solve a portion
of this ambiguity between two receivers has had considerable research to utilize the carrier
measurement for relative positioning between the receivers [23]. This technique is commonly
referred to as Real-Time Kinematic (RTK) operation.
10
Code
The code subsystem includes two distinct set of code generators. These codes are used
by the receiver to identify what satellites are available and also to determine the time of
the signal arrival. Both codes are based on the same theory of maximal length codes and
are chosen for their excellent autocorrelation and crosscorrelation properties. They are
binary, deterministic, pseudo-random noise (prn) sequences that are distinct between any
two satellites. The coarse-acquisition (C/A) code is modulated using binary phase shift
keying (BPSK) onto the L1 carrier. It is a code that is 1023 bits (called chips for code
sequences) long and repeats at a rate of 1.023 Mchips/s, or once every millisecond. The
P(Y) code is a much longer code that is truncated after every 7-day period for a length of
6.1871x1012 chips at a rate of 10.23 Mchips/s. The increased frequency of the P(Y) code
allows for more accurate measure of the signal time of arrival. This code is modulated onto
both the L1 and L2 carriers but phase shifted 90 from the C/A code to avoid interference.
The P(Y) code is currently encrypted for use only by US-government approved receivers
and is not considered further in this thesis.
The special properties of these codes allow for the tracking of multiple satellites on the
same carrier frequency, a property called code-division multiple access (CDMA). A receiver
can generate a local replica of the code and multiply the incoming signal with this local
replica to track and determine time of arrival.
Navigation Message
The last signal subsystem generates the navigation message. This is a 50 bit/s message
that includes all the data parameters needed by the receiver to calculate a satellite?s position
11
(called ephemeris), health, clock correction, and other items as well as the GPS week, time,
etc. The navigation message was designed to transmit a large amount of data e ciently
while maintaining a relatively slow bit rate. Each frame of the navigation message is broken
into  ve subframes, the contents of which are shown in Table 2.1 [19].
Table 2.1: Navigation Message Subframes
Subframe Contents
1 Clock Correction & Satellite Quality
2 Ephemeris
3 Ephemeris
4 Almanac & Ionosphere & UTC Correction
5 Almanac
2.3.2 GPS Control Segment
The GPS control segment tracks the progress of the satellites and uses their measure-
ments to estimate health, clock errors, path, etc. Parameters are then formed and an uplink
loads the desired parameters back onto the satellites to be broadcasted to the users through
the navigation message.
The GPS control segment consists of 4 control stations capable of relaying updated
information to the satellites, 6 monitor stations to track the satellite?s progress, and 10 Na-
tional Geospatial-Intelligence Agency (NGA) tracking stations [21]. The increasing number
of sites in use has several bene ts to the overall system. First, the control segment is better
capable of detecting errors in the individual satellites which leads to quicker response to
gross failure detection and user warning. This warning can be issued by the control segment
by switching the satellite?s status from healthy to not recommended for use. Secondly, in-
creased tracking capabilities have led to improved path estimation for the satellites. With
12
better knowledge of the path, more accurate ephemeris parameters can be transmitted to
the users. This improvement has led to the general user?s improved positioning accuracy as
well [8].
2.3.3 Ephemeris
For the user to calculate his position, he must have measured ranges from at least four
satellites for reasons explained in Section 2.2. The receiver must also be able to determine
the position of the range source (i.e. the satellite) at the time of position calculation.
This is accomplished by reading the ephemeris parameters from the navigation message
transmitted by the satellites.
The ephemeris parameters are actually estimated Keplerian orbit parameters along
with other parameters that de ne corrections to Kepler?s motion model. These parameters
are listed and described in Table 2.2, which is taken from [11].
Satellite Clock Correction
Part of the ephemeris parameters are used by the receiver to calculate the satellite?s
clock error. Although the atomic clocks aboard the space vehicles are highly accurate and
stable, using them without correction would lead to a large degradation in accuracy. By
using Equation (2.1), it is seen that the transit time, tt, (time the signal took to go from
the satellite to the receiver) and the transmission time, ts, are
tt = tm ts = ruc   uc (2.4)
ts = tm  uc (2.5)
13
Table 2.2: Clock and Ephemeris Parameters
Parameter Description Units
TGD Group Delay s
IODC Issue of data, clock -
toc Clock data reference time s
af2 Second-order correction to the satellite clock s/s2
af1 First-order correction to the satellite clock s/s
af0 Constant correction to the satellite clock s
M0 Mean anomaly at reference time semicircles
 n Mean motion di erence from computed value semicircles
e Eccentricity -p
A Square root of the semimajor axis pm
  0 Right ascension at reference time semicircles
i0 Inclination angle at reference time semicircles
! Argument of perigee semicircles
 _ Rate of right ascension semicircles/s
 _i Rate of inclination angle semicircles/s
Cuc Cosine harmonic correction term to the argument of latitude rad
Cus Sine harmonic correction term to the argument of latitude rad
Crc Cosine harmonic correction term to the orbit radius m
Crs Sine harmonic correction term to the orbit radius m
Cic Cosine harmonic correction term to the angle of inclination rad
Cis Sine harmonic correction term to the angle of inclination rad
toe Ephemeris reference time s
IODE Issue of data, ephemeris -
All of the above variables are either known or measured. The satellite clock error is modeled
as a quadratic function with parameters af0, af1, and af2. They are available from the
navigation message in subframe 1 shown in Table 2.1. It is used in function form with time
referenced from the clock correction reference time, toc. The clock correction is calculated
as
dtsv = af0 +af1 (ts toc) +af2 (ts toc)2 +dtr TGD (2.6)
14
where TGD is the group delay clock correction for single frequency users (L1 or L2 only).
This is necessary since the clock delay cannot be corrected without dual frequency use [7].
Also, dtr is a relativistic correction term calculated as
dtr = FepAsin (Ek) (2.7)
where the orbit parameters, e, A, and Ek, are described in the satellite position section and
F is a constant whose value is  4:442807633x10 10 s/pm [7]. The user can calculate the
corrected transmission time with the satellite clock correction and therefore the time from
the ephemeris reference epoch (toe) as
t = ts dtsv (2.8)
tk = t toe (2.9)
The time from the ephemeris reference epoch is de ned to account for beginning or ending
of week crossovers. The logic to account for the crossover is that if the value is greater than
302,400 s, then 604,800 s are subtracted from it. Likewise, if the value is less than -302,400
s, then 604,800 s are added to it. This approach tracks the time di erence since the value
of t is reset at the start of a new GPS week [7].
Satellite Position
Determining where a satellite is in its orbit amounts to  nding the angle it is from
the orbit?s perigee (point of closest approach to the Earth). This angle is the eccentric
anomaly and is shown in Figure 2.3 as E. The mean anomaly, M, is used due to the highly
15
nonlinear change of speed of a body in orbit and the corresponding di culty in calculating
the angle traveled over time. The mean anomaly has a range between 0 and 2 . This value
is also shown in Figure 2.3 and allows for a linear relation between angles along an auxiliary
circle versus time and is calculated by time-scaling the mean motion parameter n, shown
in Equation (2.11).
Figure 2.3: Orbit Angles
Using the ephemeris parameters described in Table 2.2, the satellite position can be
calculated. The mean motion parameter is calculated as
n =  n +  !is (2.10)
where  !is = p =a3 and  is the Earth?s gravitational constant, which has a value of
3:986004418 1014m3=s2 according to the WGS84 standard [14]. The mean anomaly is an
16
angular measure that is linear with time and is calculated as
M = M0 +tkn (2.11)
With this value, the eccentric anomaly, E, is solved iteratively until convergence using
Kepler?s equation
E = M +esin (E) (2.12)
Similarly, E will have a value within the range of 0 and 2 . With the eccentric anomaly,
the true anomaly,  , can be found as
 = tan 1
0
B@
p1 e2 sin (E)
1 ecos (E)
cos (E) e
1 ecos (E)
1
CA (2.13)
where tan 1 refers to the four-quadrant inverse tangent. Next the argument of latitude,  ,
can be calculated to lie within the range of 0 and 2 as
 =  +! (2.14)
Using the C terms from the ephemeris parameters in Table 2.2, corrections for the argument
of latitude, radius, and inclination can be calculated as
 uoos = Cus sin (2 ) +Cuc cos (2 ) (2.15)
 roos = Crs sin (2 ) +Crc cos (2 ) (2.16)
 i = Cis sin (2 ) +Cic cos (2 ) (2.17)
17
uoos =  + uoos (2.18)
roos = A(1 ecos (E)) + roos (2.19)
i = i0 + i+  _itk (2.20)
Now the 2D position in the orbital plane can be found as
xoos = roos cos (uoos) (2.21)
yoos = roos sin (uoos) (2.22)
The correction to the ascending node  , which is the angle determining where the orbit
crosses the Earth, can then be calculated to be within 0 to 2 as
 =  0 +
 _
   e
 
tk  e (toe +tt +dtsv) (2.23)
With these values, the ECEF position of the satellite at time ts can be calculated as
xees = xoos cos ( ) yoos cos (i) sin ( ) (2.24)
yees = xoos sin ( ) +yoos cos (i) cos ( ) (2.25)
zees = yoos sin (i) (2.26)
Satellite Velocity
For some processing procedures, the satellite velocity is also required. A numerical
approximation to this value can be generated by calculating the satellite positions at two
successive times and dividing by the time di erence in the two epoch. Assuming the time
18
di erence is a su ciently small value, this value can be used as the satellite velocity. Al-
ternatively, the orbit Equations 2.10 thru 2.26 can be di erentiated to get an analytical
solution for satellite velocity [14].
_E =  !is +  n
1 ecos (E) (2.27)
_ = sin ( )
sin (E
_E (2.28)
_roos = Aesin (E) _E + 2 (Crs cos (2 ) Crc sin (2 )) _ (2.29)
_uoos = (1 + 2Cus cos (2 ) 2Cuc sin (2 )) _ (2.30)
Therefore the time derivative of the satellite?s orbit coordinates are
_xoos = _roos cos (uoos) roos _uoos sin (uoos) (2.31)
_yoos = _roos sin (uoos) roos _uoos cos (uoos) (2.32)
The changes in the ascending node and inclination are calculated as
_ = __   e (2.33)
_i = _ _i+ 2 (Cis cos (2 ) Cic sin (2 )) _ (2.34)
19
Thus the satellite ECEF velocity is found as
_xees = _xoos cos ( ) _yoos cos (i) sin ( ) + _iyoos sin (i) sin ( ) (2.35)
+
 
 e __ 
 
(xoos sin ( ) +yoos cos (i) cos ( ))
_yees = _xoos sin ( ) + _yoos cos (i) cos ( ) _iyoos sin (i) cos ( ) (2.36)
+
 
 e __ 
 
( xoos cos ( ) +yoos cos (i) sin ( ))
_zees = _yoos sin (i) _iyoos cos (i) (2.37)
2.4 Positioning
Four satellites are required to determine a user?s position with reasonable accuracy
with the four unknowns of user ECEF position and clock bias, (xu;yu;zu;c tu). Additional
satellites allow for an overdetermined solution and other techniques, to be discussed in
subsequent chapters, make use of this situation. Several hand-held GPS receivers allow an
additional constraint of being on the Earth?s surface to remove altitude as an unknown.
This approach only requires three satellites to determine two position coordinates along
with the receiver?s clock bias [13]. However, this is often a poor assumption as local height
changes relative to the Earth geoid are di cult to account for, and the user may violate
this constraint.
20
Considering a single time epoch, the user?s state vector  X is made up of the receiver?s
ECEF position reeb and clock bias in range form c tu.
 X =
2
66
66
66
66
4
xeb
yeb
zeb
c tu
3
77
77
77
77
5
(2.38)
The measurement vector at this epoch is generated by placing all available pseudorange
measurements in a column vector for all m satellites being tracked,
 Y =
2
66
66
4
 1
...
 m
3
77
77
5
(2.39)
The observation-state relationships are found by placing the right sides of Equation (2.3)
into vector form so that
 G =
2
66
66
4
q
(x1 xu)2 + (y1 yu)2 + (z1 zu)2 +c tu
...
q
(xm xu)2 + (ym yu)2 + (zm zu)2 +c tu
3
77
77
5
(2.40)
To e ciently solve this set of equations for the state, linearization is used to change the
problem into a recursively solved system. The linearization is taken about point  X and a
state correction vector is de ned as
 x =  X  X (2.41)
21
Using Taylor series expansion on Equation (2.40) and neglecting higher order non-linear
terms yields
 Y   G   X + @  G
@  X
  
  
  X  x
=  Y +H x
=
2
66
66
4
  1
...
  m
3
77
77
5
+
2
66
66
4
 ( e u1)T 1
... ...
 ( e um)T 1
3
77
77
5
 x (2.42)
where   represents the estimated pseudorange using the linearization point as the user?s
ECEF position and clock bias and  e us represents the unit vector from the user to the
satellite. It is determined by
 e us =
2
66
66
4
x u xs
r s
y u ys
r s
z u zs
r s
3
77
77
5
(2.43)
Therefore a measurement correction vector is calculated as
 y =  Y   Y = H x (2.44)
A least-squares solution to Equation (2.44) minimizes the error between the measurement
correction and the transformed state correction according to a cost function, where the error
is de ned as
 =  y H x (2.45)
22
In the case of least-squares, this cost function is
JLS =  T (2.46)
Thus, it minimizes the sum of the squared error. The solution of this problem is well
documented in several sources [31, 32] and is calculated as the solution to
 
HTH
 
^x = HT  y (2.47)
Methods to solve this system include the pseudoinverse, decomposition, etc. The result of
this least-squares solution method on a typical set of static data is shown below in Figure
2.4. Notice in this  gure that this method yields jumps approximately tens of meters from
the mean position. The source of these errors will be discussed in Chapter 3.
Figure 2.4: Position Using Least Squares Estimate
An improved solution can be found by using a weighted least squares method if variance
information is available for each measurement [32]. In this method, the variances for each
pseudorange measurement are placed along the diagonal of a covariance matrix R. Using
23
these weights changes the cost function to be
JWLS =  TR 1 (2.48)
The resulting minimization gives the relation
 
HTR 1H
  1
HT ^x = R 1 y (2.49)
Similar methods used to solve Equation 2.47 can be used to solve Equation 2.49.
As will be described later in Chapter 3, many of the errors that a ect a GPS receiver
have been well studied and their variances described in the literature [27]. Major error
sources a ecting a well placed static receiver are atmosphere errors and receiver errors.
The overall pseudorange variance can be approximated as
 2 =  2tDLL + 2atm (2.50)
where  2tDLL is the thermal noise variance of the receiver?s delay locked loop, which generates
the pseudorange measurements, and  2atm is the variance of the atmospheric e ects on the
pseudorange measurement. Taking the total e ect of the ionospheric and tropospheric
error gives a variance approximated as  2atm = 5:22 [24]. The thermal noise variance of
the receiver?s delay locked loop is shown as a function of the carrier to noise ratio C=N0 as
shown in Equation 2.51 [24].
 tDLL =  c
s
4F1d2Bn
C=N0
 
2(1 d) + 4F2dTC=N
0
 
(2.51)
24
The description of the receiver speci c terms are given in Table 2.3 along with values taken
from [24]. Using this variance, an example set of position solutions is shown in Figure 2.5.
Table 2.3: Receiver Parameters from [24]
Parameter Description Value Units
F1 discriminator correlator factor 1 -
F2 discriminator type factor 1 -
d correlator spacing 1/2 chips
Bn code loop noise bandwidth 2 Hz
T prediction integration time 2 ms
 c code chipping rate 293.05 chips/s
Compared to solutions shown in Figure 2.4, this method produces a less noisy and more
accurate position solution by de-weighting pseudoranges that are comparatively weak (low
C=N0 values).
Figure 2.5: Position Using Weighted Least Squares Estimate
2.5 Summary
This chapter describes the basics necessary for a GPS receiver to calculate its position
using the ephemeris parameters from the GPS data message and the pseudoranges found
25
from determining the di erence in transmission time, ts, and the user?s measurement time
tu. Chapter 3 will build on this background and show the e ects of various error sources
on GPS operation.
26
Chapter 3
Overview of GPS Errors
Much work has been done in analyzing the errors in GPS since its use began. A central
source of information could be considered [27], which details the major error sources and
includes a commonly quoted error budget, shown in Table 3.1. Several error modes have
been heavily studied when it comes to the use of GPS receivers. They are typically grouped
into three di erent error regions: satellite, atmospheric, and local. The satellite errors are
the ephemeris data and satellite clock errors. The atmospheric errors are the ionosphere and
troposphere delays. Finally, the local errors are the multipath and receiver measurement
errors.
Table 3.1: GPS Standard Error Budget [27]
One-sigma error, m
Error source Bias Random Total
Ephemeris data 2.1 0.0 2.1
Satellite clock 2.0 0.7 2.1
Ionosphere 4.0 0.5 4.0
Troposphere 0.5 0.5 0.7
Multipath 1.0 1.0 1.4
Receiver measurement 0.5 0.2 0.5
UERE, rms 5.1 1.4 5.3
Filtered UERE, rms 5.1 0.4 5.1
3.1 Satellite Errors
The GPS constellation was engineered to provide high accuracy position using ranges
from the satellites. As with all implementations, there are still errors that must be accounted
for. The errors in this section deal with the design and implementation of the satellites.
27
Both of these error sources have broadcast correction models that signi cantly decrease their
e ect, but there are still unmodeled e ects that a ect the user in navigation calculations.
These errors are the same for all users of GPS and therefore are exactly correlated between
any two locations. This fact can be used in certain error mitigation techniques described
later.
3.1.1 Errors in Ephemeris Data
The satellite positions are not perfectly known to the user. To determine a satellite?s
position in ECEF coordinates, the user receives the satellite?s ephemeris parameters in the
navigation message. These parameters de ne the satellite?s path estimated by the GPS
control segment. This path is evaluated at the time the signal leaves the satellite?s antenna.
Since this path is an estimate, there is a di erence between the true satellite path and
the path described by the ephemeris parameters. These errors are typically 3-5 m rms [25].
However, most of this error is in the cross-track direction which means it a ects the direction
to the satellite rather than the range, as shown in Figure 3.1. Errors in the direction to
the satellite (which arise in the unit vector from the user to the satellite) have a much
smaller impact than along-track errors [27]. The resulting ranging error from the satellite?s
position is therefore on the order of 1.5 m. The orbit model is updated by the GPS control
segment periodically (anywhere from about eight hours to daily). At the time of update,
the ephemeris error is nearly zero since the time is near the reference time for the orbit. As
time progresses between ephemeris updates, the satellite position error grows.
28
Figure 3.1: GPS Ephemeris Errors
In order to compensate for these positioning errors, precise ephemerides are available
post-process (typically one to a few days later, depending on accuracy desired). The Na-
tional Geodetic Survey (NGS) posts sets of these orbit parameters at [2]. More complicated
models are used to generate these solutions. The model described in Section 2.3.3 was
chosen so its parameters would  t compactly in the navigation message. Increases in model
complexity would alter this message and therefore make signi cant changes to the GPS
signal structure.
3.1.2 Satellite Clock Errors
The satellites in the GPS use highly stable atomic clocks so the assumption of syn-
chronized signal transmission can be used. However, there are errors in the satellite clocks
due to drift. Most of these errors can be compensated for by using a clock error model. In
the navigation message,there is a clock model to use that is also estimated by the control
segment, described in Section 2.3.3. The error from satellite clock drift is also on the order
of 1.5 m [25]. This error is reduced from that found in Table 3.1 due to increases in control
29
segment capabilities and the use of more stable atomic clocks in more recent GPS satellite
blocks [21].
3.2 Atmospheric Errors
As the GPS signals travel through the Earth?s atmosphere, they are subject to various
changes in speed. The majority of these delays come from two layers of the atmosphere:
the ionosphere and troposphere. Due to the variation in composition in these layers, the
amount and means of the e ect are di erent. The various layers of the atmosphere are
shown in Figure 3.2, from [3]. Since these e ects are caused by the satellite signal passing
through a certain portion of the atmosphere, these errors are highly correlated with respect
to location. That is, if two users are in close proximity (usually taken to be in the kilometer
range), the atmospheric errors on a particular satellite signal are approximately equal. The
farther apart two receivers are, the less correlation between the atmospheric errors.
Figure 3.2: Atmospheric Layers [3]
30
3.2.1 Ionospheric Errors
As seen from Figure 3.2, the ionosphere is a portion of the Earth?s atmosphere extending
from around 85 km to over 600 km above the Earth?s surface. It is characterized by ionized
partiles due to excitation by the Sun?s ultraviolet radiation [27]. The magnitude of the
e ects the ionosphere has on the GPS signals is a function of the total electron content,
which can change quickly as time progresses. This e ect is dispersive, meaning its e ect
on the signal is a function of the signal frequency. This fact led, at least in part, to the
use of multiple frequencies on the GPS satellites. By using measurements on multiple
frequencies, the ionospheric error can be estimated and subsequently removed [6]. However,
this approach requires the use of dual frequency measurements. These measurements are
not often available on low-cost GPS receivers. This is due to the fact that the code on the L2
frequency is encrypted for use by authorized users only, such as the military. However, all
new satellites transmit a civilian code on the L2 frequency commonly called L2C [35]. Once
older satellites are phased out and replaced with updated ones, the ionosphere correction
will be available to all users.
Another e ect the ionosphere has on the satellite signal is that it a ects the code and
carrier signals di erently. The dispersive media in the ionosphere delays the code, called
group delay, while advancing the carrier phase, called phase advance. These two e ects are
equal in magnitude but opposite in sign [25]. This arises due to the fact that dispersive
media make signals of di erent frequencies travel at di erent rates through the media. Since
the code and carrier are signals of signi cantly di erent frequencies, their speeds are altered.
For a satellite s, once the L1 ionospheric delay (Is1) has been calculated, the L2 delay (Is2)
31
is also found by
Is2 =  Is1 (3.1)
 =
 f
L1
fL2
 2
(3.2)
Included in the satellite?s broadcast navigation message are parameters for an ionospheric
delay model called the Klobuchar model, shown in Equation (3.3). This model gives the L1
zenith ionospheric delay, ^Iz;L1, as a function of local time, t, as
^Iz;L1
c =
8
<
:
A1 +A2 cos
 2 (t A
3)
A4
 
jt A3j< A44
A1 otherwise
(3.3)
where A1 is the nighttime value of the zenith delay ( xed at 5 10 9 s), A2 is the amplitude
for the daytime values, A3 is the phase peak time ( xed at 14 hours, local time), and A4 is
the period of the cosine function. Using some typical values of these parameters (A2 = 7
m, A4 = 28 hours) yields corrections shown in Figure 3.3 [25]. Since the ionospheric e ect
is a function of the amount of ultraviolet rays absorbed by the atmosphere, the model has
to compensate more during the middle of the day when the Sun is directly overhead. An
obliquity factor, OF, is used to scale the zenith delay, ^Iz;L1, as a function of the elevation
angle to the satellite, el, as
OF (el) = 1:0 + 16:0 (0:53 el)3 (3.4)
This factor is needed since the amount of atmosphere the signal travels through is a function
of the elevation angle from the horizon. At low elevation (near the horizon), the signal travels
32
through more of the atmosphere. Thus, the ionospheric delay is more and the model must
be scaled. A plot of the obliquity factor is shown in Figure 3.4 [25].
Figure 3.3: Typical Klobuchar Ionosphere Model [25]
Figure 3.4: Ionospheric Obliquity Factor Model [25]
33
3.2.2 Tropospheric Errors
The troposphere also causes a delay in the signals passing through the atmosphere. Un-
like the ionosphere, however, the delay is not a function of the frequency, so dual frequency
measurements can not be used to compensate for the tropospheric error. This error mode
is often broken into wet and dry components to be modeled. The wet model is dependent
on the amount of moisture in the air and is therefore highly variable. However, it makes
up only approximately 10% of the tropospheric error when compared to the delaying na-
ture of the dry component. In order to compensate for this error (which is often neglected
compared to other errors), a model is required [25].
3.3 Local Errors
Local errors are those that a ect users on an individual basis. While the other error
groupings are correlated between users, these local e ects are either very weakly correlated
(in the case of multipath) or uncorrelated (with regards to receiver measurement errors).
Often these errors are best dealt with in receiver and antenna design as well as the user
choosing a good location for the antenna. In the shadowed environments, it is the local
errors that change drastically and are analyzed in this research.
3.3.1 Multipath Error
Multipath error occurs when the signal broadcast by the satellites bounces o surfaces
before reaching the receiver antenna. This has a twofold e ect on the signal. It  rst delays
the signal since it takes a longer path to the receiver antenna than the line-of-sight signal.
It also attenuates the signal since some of the energy is re ected in directions away from
34
the receiver. Research on this error source has increased in the last decade due to the
e ectiveness of mitigating other error sources such as the ionospheric delay, satellite errors,
and the removal of Selective Availability (SA). Due to these advances, one of the next major
sources of error comes from multipath [4].
Due to the fact that multipath is by nature highly dependent on the receiver antenna?s
immediate surroundings, it cannot be removed by di erential GPS processing methods
[9]. Multipath has been broken down into three di erent categories, each resulting from a
di erent type of impact from surrounding objects and each having di erent e ects on the
receiver [26]. The  rst is specular re ection which arises from re ection o of a smooth
surface. Due to the re ection, the amplitude of the signal is reduced, and the path from
satellite antenna to user antenna is longer. This results in a waveform that is both lower in
amplitude and out of phase with the line-of-sight signal [12]. This category greatly a ects
the carrier phase and is thus the major multipath e ect in precise applications. The second
multipath category is di raction. This e ect appears due to re ections around edges and
corners of objects in the path of the signal. Lastly, when a signal bounces o a rough surface,
many specular re ections occur. This is referred to as di use multipath. These three types
of multipath are illustrated in Figure 3.5. The e ects of multipath will be described more
in Section 4.1. It is also noted that the multipath e ect on the carrier phase measurement
is much less than that on the pseudorange measurement. This will be used as a simplifying
assumption in Section 4.1.
35
Figure 3.5: Three Multipath Types
3.3.2 Receiver Measurement Error
Several errors arise when a GPS receiver generates a measurement and brings the infor-
mation it is receiving from the satellites into the discrete domain for processing. Receivers
are subject to thermal noise in which the receiver tracking loops are attempting to drive an
error signal to zero. This error signal is the di erence between the incoming signal and the
local replica of the signal. As the noise on the measurement increases, as is the case with a
temperature increase, the uncertainty in the measurements is increased. Also, discretization
of the signal and its subsequent processing leads to greater error. However, as advances are
made in computing capabilities, these issues become less and less severe.
36
3.4 GPS Observable Error Models
Each of these error modes can be treated as a delay in the signal that corresponds
to an additive factor in the pseudorange and carrier phase observable models. Therefore,
for satellite s and L-band frequency l where l = L1, L2, the GPS observable models are
extended from Equation (2.3) to include the errors described in this chapter as
 sl = rbs +c( tu + ts) +Isl +Ts +M sl + sl (3.5)
 sl = rbs +c( tu + ts) +Nsl l Isl +Ts +M sl + sl (3.6)
where c ts is the satellite clock error, Isl is the ionospheric error, Ts is the tropospheric
error, M sl is the pseudorange multipath error, M sl is the carrier phase multipath error,
Nsl is the carrier phase integer ambiguity,  l is the wavelength for frequency l, and  sl
contains the other errors on the carrier phase measurement.
3.5 Summary
This chapter has discussed the basic common errors experienced by all users of GPS
receivers. As was mentioned, the three error groupings (satellite, atmospheric, and local)
also roughly correspond to the amount of correlation in these errors between users. It is the
local errors that are signi cantly a ected by the surrounding environment, and it is these
that are of concern in this thesis and will be analyzed in Chapter 4.
37
Chapter 4
Shadowed Environments
Shadowed environments such as heavy foliage are distinct in the errors a ecting GPS
operation. Although all users of GPS su er from error modes mentioned in Chapter 3,
the presence of a large number of trees also causes changes in these modes that must be
considered for a user to operate under these conditions. This and other environments like
urban canyons create signal blockages where a receiver is unable to track a satellite due to
an obstacle blocking satellite visibility. While this adversely a ects performance, there are
cases when a signal is degraded so much its presence in the solution increases the error more
than if it were not included at all. These environments are subject to some sort of signal
blockage that degrades receiver performance but does not cause the receiver to lose lock on
the signal. This degradation comes from quickly changing multipath due to a large number
of objects in close proximity as well as signal attenuation due to interfering objects, shown
in Figure 4.1. These are the main e ects considered in this research.
4.1 Environmental Error Monitoring Parameters
In order to study how these error sources occur in heavy foliage environments, the
e ects considered must be determined among the general errors occurring in the GPS mea-
surements. While monitoring the pseudorange and pseudorange rate errors is ideal, the true
calculation of these metrics requires the use of a highly accurate truth system. Therefore,
alternative monitoring methods were investigated. To determine signal attenuation e ects
38
Figure 4.1: Illustration of GPS Signal Shadowing Due to Foliage
the carrier to noise ratio (C=N0) is monitored. This gives information as to when a sig-
nal is degraded but loss-of-lock has not occurred. To monitor multipath, observation of a
parameter M1 has been shown to be related to a receiver?s pseudorange multipath [29].
4.1.1 Signal Attenuation Monitoring Parameter
Signal attenuation arises when the signal passes through foliage cover. These obstacles
cause both a delay and a reduced signal amplitude which increases the di culty of tracking
the signal. This reduced amplitude e ect can be seen in the C=N0 that is reported by the
GPS receiver. Decreases in the C=N0 increase the receiver?s delay-locked-loop uncertainty,
 tDLL, shown in Equation (2.51). The pseudorange variance,  2 , is the combination of the
tracking error and the atmospheric errors shown in Equation (2.50). Similarly the C=N0
also shows the increase in the pseudorange rate variance. The pseudorange rate is related
39
to the receiver?s measurement of the Doppler frequency as
_ s = cf
L1
fD;s (4.1)
The variance of the Doppler frequency is related to the noise on the Frequency-Lock-
Loop (FLL) that tracks the frequency of the incoming signal. The thermal variance of the
FLL is therefore taken to be [36]
 tFLL =  L2 T
s
4FBn
C=N0
 
1 + 1TC=N
0
 
(4.2)
Again the values of the parameters in Equation 4.2 are assumed from typical values found
in [36] and shown in Table 4.1.
Table 4.1: Frequency-Lock-Loop Parameters from [36]
Parameter Description Value Units
F loop factor 1 -
Bn frequency loop noise bandwidth 2 Hz
T predetection integration time 2 ms
 L carrier wavelength 0.19 m
However, the FLL (and therefore the pseudorange rate measurement) is also a function
of the dynamic stress error on the tracking loop, fe. These stresses are due to the dynamics
of the receiver antenna. Assuming mild dynamics with a dynamic stress of fe = 3Hz
and taking the FLL variance to be the pseudorange rate variance, i.e.  2_ =  2FLL, the
pseudorange rate variance is therefore [36]
 2_ =
 
 tFLL + fe3
 2
(4.3)
40
Since the variances on the pseudorange and pseudorange rate are both functions of the
C=N0, this parameter is chosen for monitoring the errors in shadowed environments.
4.1.2 Multipath Monitoring Parameter
The other dominant e ect to monitor is the changing multipath. As was discussed in
Section 3.3.1, objects near the GPS antenna can alter the propagated signal by re ection.
Depending on the local surroundings, these objects may re ect toward the antenna multiple
delayed copies of the same signal. A receiver attempts to track the combination of all copies
of the signal, which leads to decreased performance. Since the re ected path is always longer
than the direct path, the multipath signals are always delayed. The re ected signals are
also typically weaker than the line-of-sight signal depending on the nature of the re ection.
As seen in Figure 4.1, multipath often causes a large change in the line-of-sight to the
signal source. This error mode also a ects the receiver?s range rate measurement since the
antenna?s velocity projected onto the line-of-sight vector to the satellite is di erent from its
projection onto the re ected path.
In order to quantify these multipath e ects, a multipath monitoring variable, labeled
M1 in this thesis, is used. For a given satellite, three measurements are used to derive the
M1 parameter: the L1 pseudorange and L1 and L2 carrier phase. These measurements
are shown in their common form with the error sources being added to the true range in
Equations (3.5) and (3.6). It should be noted that the multipath on the carrier phase, M i,
is signi cantly less than the multipath on the pseudorange, M i [29], i.e. M i <<M i.
41
Using the carrier phase model in Equation (3.6) and the dual frequency ionosphere
relation in Equation (3.1) the L1 ionosphere delay can be found as
 1  2 = N1 1 N2 2 I1 +I2 +M 1 M 2
) I1 =  1  2  1  N1 1 N2 2  1  M 1 M 2  1 (4.4)
where the noise term is dropped. Since these parameters are using values from only one
satellite, s is also removed from the notation. Combining Equation (4.4) with the L1 carrier
phase in Equation (3.6) removes the ionosphere delay as
 1 +  1  2  1 = rs +c( tu + ts) +N1 1 + N1 1 N2 2  1
+M 1 + M 1 M 2  1 +T1 (4.5)
Then combining Equation (3.5) with Equations (4.4) and (4.5) gives
M1 =  1 
 
1 + 2  1
 
 1 +
 2
  1
 
 2 (4.6)
=  N1 1 2N1 1 N2 2  1 +M 1 M 1 2M 1 M 2  1 (4.7)
 B1 +M 1 (4.8)
where the carrier phase multipath variables have been neglected and the bias term
B1 = N1 1 2N1 1 N2 2  1 (4.9)
is a function of the integer ambiguities on the dual frequency carrier phase measurements.
42
This bias term, B1 can be removed by averaging or taking the  rst M1 measurement as
the bias term [29]. For this work, the change in multipath is being monitored so calculating
the di erence between successive M1 variables removes the bias term. This approximates
the change in L1 pseudorange multipath which is the e ect to be monitored. Since B1 is
a function of the integer ambiguities a cycle slip will change the bias. This would appear
as a spike in the M1 variable. This e ect is neglected in the analysis but is a source of
additional error that could be considered in future work.
Since L1 pseudorange and L1 and L2 carrier phase measurements are required for this
analysis, static dual-frequency data was collected using a NovAtel DL-V3 receiver with a
NovAtel dual-frequency antenna. The photo in Figure 4.2 shows the location of the static
analysis since it provided a good contrast between clear view and heavy foliage. This allowed
for clear distinctions between the two environments to be drawn. A visualization technique
was then used to evaluate these e ects in heavy foliage.
4.1.3 Error Visualization for a Single Satellite
In order to determine the total e ect of foliage cover on the satellite signals, a visual-
ization technique for a single satellite was developed. Both monitoring variables (M1 and
C=N0) can be plotted together to visualize the signal error at a point in time. For a single
satellite, this is shown in Figure 4.3. The change in MP1 is shown as the magnitude in
this error plot. The C=N0 is shown as the shade in the plot. These monitoring variables
together reveal the e ects of the environment on the signal. Figure 4.3 shows these mon-
itoring variables for a single satellite signal as it passes from heavy foliage to a clear sky
43
Figure 4.2: Static Data Collection Location
view. This diagram shows that the signal strength drops and multipath jumps correspond
with each other, although not directly.
4.1.4 Error Visualization for All Satellites in View
The single satellite visualization technique was extended to produce a visualization
technique for all satellites in view. For a better view of what happens in tree cover, the
static location shown in Figure 4.2 was again used. This location had heavy tree cover on the
west side and clear sky on the east side. The aerial view of this area is shown in Figure 4.4.
The combination of all satellite?s monitoring variables are shown in Figure 4.5. This plot is
similar to Figure 4.3 except that each of the M1 and C=N0 plots are bent to follow the path
that the satellite travels across the sky. This allows spatially correlated e ects to be shown.
Since there is a relatively distinct dividing line between the two regions, the distinction can
44
Figure 4.3: Single Satellite Error Visualization
also be seen in the sky plot of Figure 4.5. This plot shows the e ect that foliage has in
quickly changing the multipath as well as degrading the signal. Low-elevation e ects are also
shown as satellites get close to the outer circle. This plotting scheme is extended from [17]
where only the actual multipath is plotted with respect to satellite azimuth and elevation.
The visualization developed for this thesis shows the changing e ects of the multipath as
well as the signal strength for analysis in changing conditions.
4.1.5 Analysis Results
From the monitoring variable plots shown in Figures 4.3 and 4.5, the shadowing e ects
of heavy foliage are shown to degrade the satellite signals. Even when a signal is still being
tracked and is available to the system, low signal power (shown through low C=N0 values)
45
Figure 4.4: Aerial View of Static Data Location
Figure 4.5: All-In-View Satellite Error Visualization
46
and large jumps in the multipath (seen through M1 variable) indicate that accuracy is
degraded.
4.2 Summary
This chapter discusses the error modes speci c to shadowed environments and the
technique developed for their analysis. Monitoring variables were chosen to visualize the
e ects of foliage blocking a satellite signal. Although these e ects usually cause a signal
to drop out of view due to blockage, the analysis shows that even when satellites are
still available, their accuracy is compromised due to the shadowing e ects. This chapter
motivates the need to develop a solution for operation in these environments, which is
presented in Chapter 5.
47
Chapter 5
Error Mitigation
The analysis presented in Chapter 4 motivates the need for methods to mitigate errors
speci c to shadowed environments. These errors include signal attenuation and quickly
changing multipath. Both these error modes cause an increase in the uncertainty of the
GPS receiver measurements and thus degraded navigation performance. Error mitigation
is necessary to obtain accurate and reliable navigation solutions for operation in shadowed
environments. From the design requirements in Chapter 1, the system must be low cost,
operate continuously, and avoid solution spikes so it can be used in conjunction with vehicle
navigation and control systems. In Chapter 4 it was shown that a stand-alone GPS receiver
will not meet these requirements. Two methods of improved operation are presented here:
di erential corrections and sensor coupling. Although the use of di erential corrections
is very common and improves the accuracy of GPS solutions relative to a base receiver,
the shadowing speci c errors presented in Chapter 4 cause di erential methods to also fail
since these e ects are not compensated in the correcting algorithms. Therefore the solution
implemented in this thesis is integration with an inertial sensor in a closely coupled archi-
tecture with fault detection and exclusion (FDE) performed by an innovation monitoring
technique [14].
5.1 Di erential Corrections
GPS accuracy can be improved with the use of other local GPS receivers [28]. This is
accomplished by operating on GPS measurements to calculate the position of one receiver
48
relative to the other. With this as the goal, there are methods that can reduce or even
eliminate some of the aforementioned error types. The di erential research presented here
was conducted as an attempt to improve relative positioning between multiple receivers.
This relative position can be used to determine an accurate global position for a mobile
receiver if the base receiver is placed at a known, surveyed location.
5.1.1 Position Di erencing
The most basic of di erential positioning methods is to take the di erence of two inde-
pendently calculated positions. This method assumes the correlated range errors detailed
in the pseudorange and carrier phase Equations (3.5) and (3.6) result in a correlated error
in position. The di erencing removes position error that is common to the two receivers.
However, even though the error sources may be correlated between two receivers, their
combined e ects on the position solution changes the amount of correlation. This method
tends to remove some errors common to both calculated positions, but the combination is
not a linear di erencing of the errors. This is due to the fact that the di erencing is done
in the position domain rather than the range domain. In this case, the individual solu-
tions are calculated using biased and noisy measurements as described in Equations (3.5)
and (3.6) to generate a solution based on methods similar to the one described in Section
2.4. These independently calculated positions can also have corrections applied using com-
mercial or freely available corrections such as Omnistar or the Wide Area Augmentation
System (WAAS), respectively. This typically reduces the global position error of each of
the independent solutions to around three meters. After di erencing these two corrected
49
positions, the relative positioning error is on the order of one meter for a commercial re-
ceiver. Note that these di erences are subject to large jumps since errors can vary quickly
between receivers.
5.1.2 Code Di erencing
If the measurements of two receivers are temporally and spatially correlated, several of
the delay errors in the range measurement can be removed if the actual measurements are
di erenced before the position solution is calculated. This error correlation was described
previously in Chapter 3. If two receivers, A and B, can view the same satellite, j, as shown
in Figure 5.1, di erencing these two ranges removes the errors common between them. This
is called single di erencing, and it removes the common satellite clock error,  tj.
Figure 5.1: Common Satellite
50
The resulting single di erence models are shown in Equations (5.1) and (5.2) where
the L1 frequency is assumed.
 j;AB = rAB +c tAB +Ij;AB +Tj;AB +M j;AB + s;AB (5.1)
 j;AB = rAB +c tAB +Nj;AB  Ij;AB +Tj;AB +M j;AB + s;AB (5.2)
Note in both equations, the components A and B are combined into similar terms represent-
ing the di erence in each component between the receivers for the satellite j, noted as AB.
Also note the di erencing of the measurements includes the noise from each measurement,
thus increasing the noise of the new combined measurement. The single di erence measure-
ments can be useful when considered along with Figure 5.2. Here it is shown that the single
di erenced measurement is related to the baseline by the unit vector to the satellite, j.
Figure 5.2: Single Di erence Diagram
Assuming the distance to the satellite is much greater than the baseline length, the
unit vectors from receivers A and B are approximately the same. As an example, assuming
51
a baseline of 10 km with the average satellite altitude of 20,200 km above the receivers A
and B with the satellite directly above receiver A, the angle of separation,  , between the
two receivers as viewed from the satellite would be 20200 kmsin( 
2 rad)
= 10 kmsin( ) ) ? 0:000495 rad.
This angle also becomes the di erence in elevation angle from satellite B to the satellite
and therefore the error that would be included in the direction of the unit vector.
From Figure 5.2, Equation (5.3) can be formed.
SDj;AB = ujA bAB (5.3)
=jujAjjbABjcos ( )
Where  is the angle between the unit vector to the satellite, ujA and the baseline vector,
bAB. Since this angle is assumed to be  2 rad with an error of  at the point of consideration,
the equation expands to
SDj;AB =jejAjjbABjcos (   )
= 10 km sin ( )
?1:225mm
Therefore, even under these poor conditions, the error is insigni cant compared to the noise
on the single di erence measurement. Also note if the receivers are spatially correlated, their
atmospheric errors are approximately equal, thus Ij;AB?Tj;AB?0.
Using two single di erences for satellites j and k, these two measurements can be
combined to form double di erences, shown in Equations (5.4) and (5.5) and illustrated in
52
Figure 5.3
 jk;AB = rjk;AB + jk;AB (5.4)
 jk;AB = rjk;AB +Njk;AB + jk;AB (5.5)
This di erencing removes the user clock bias,  tAB, which was common to both single
Figure 5.3: Two Common Satellites
di erences. This removes the requirement that the clock bias must be estimated. However,
this comes at the expense of further increasing the noise on these new quantities. Also
notice that the integer ambiguities in the carrier phases (Nj;A, Nj;B, Nk;A, and Nk;B)
are lumped into a single ambiguity (Njk;AB) that must be estimated rather than the four
integer ambiguities in the original measurements shown in Equation(3.6). Similar to the
single di erence, the double di erence can be related to the baseline as
DDjk;AB = (ujA ukA) bAB (5.6)
53
Placing all available double di erences with each di erenced about satellite j?s single
di erence into matrix form yields a relation to estimate the baseline in Equation (5.7)
2
66
66
4
ujA ukA
...
ujA umA
3
77
77
5
 bAB =
2
66
66
4
DDjk;AB
...
DDjm;AB
3
77
77
5
(5.7)
5.1.3 Carrier Smoothed Code Di erencing
The methods described so far do not make use of the highly accurate carrier phase
measurement due to the di culty of estimating the double di erence integer ambiguities,
Njk;AB. However, these ambiguities are constant over time, assuming no cycle slips have
occurred. A smoothing method is used in which a double di erence range, DD, is estimated
using the carrier phase double di erence ( jk;AB) and the pseudorange double di erence
( jk;AB). By using the carrier phase double di erences at two time epochs ( njk;AB and
 n 1jk;AB) to determine the change in the double di erence range (DD), a low noise measure-
ment is obtained of the change in the double di erence over one time step. By combining
these measurements in a Kalman  lter, the double di erence range (DD) can be calculated
which uses the low noise carrier phase double di erences ( jk;AB) to remove the noise in-
herent in the code double di erence ( jk;AB). The resulting measurement is a very smooth
double di erence that lies along the mean of the code double di erence measurement. To
calculate the smoothed di erence, Equations (5.8) through (5.12) are employed.
DD = DD+ + 
 
 njk;AB  n 1jk;AB
 
(5.8)
P = P +Q (5.9)
54
Figure 5.4: Method Comparison
K = P (P +R) 1 (5.10)
DD+ = DD +K
 
DD + njk;AB
 
(5.11)
P = (I K)P (5.12)
Here P is the properly initialized covariance matrix, Q is the process noise covariance, R is
the measurement covariance, and K is the resulting Kalman gain.
5.1.4 Static Comparison of Di erential Methods
The three di erential methods are shown in Figure 5.4 for a static case. In this com-
parison the code double di erence performed the worst with a 3D standard deviation of 1.2
m. The carrier smoothed and position di erence estimates were fairly close to the same
performance with deviations of 0.54 m and 0.67 m, respectively.
55
The noise performance of the smoothed double di erence is slightly better than that
of the raw position di erencing due to the low noise of the carrier phase measurements.
However, both su er from drift about the true relative position magnitude (baseline). This
drift results in errors of approximately one meter from the actual baseline. Assuming a well
surveyed base position, the global position will be accurate to about one meter. This is an
improvement over the global accuracy of a single receiver. However, the double di erence
processing techniques do not provide a signi cant advantage over the raw position di er-
encing. Another di culty in applying these methods to roving vehicles is the requirement
for a surveyed base station and data link between the receivers. Also, as the baseline grows,
the di erential methods su er a drop in performance. This is due to the fact that the error
correlation decreases and each receiver may be tracking di erent satellites. Therefore an
onboard solution is desired and described in this thesis.
5.2 Inertial Aiding Sensors
An alternative to coupling multiple GPS receivers is to use alternative sensors, particu-
larly an inertial measurement unit (IMU). Since inertial sensors have high output rates and
low short-term noise, they are well suited for navigation and control over short periods of
time. The length of this time is dependent on the grade of IMU used [37] since these sensors
tend to drift over time. To keep the system low-cost, an automotive grade IMU, the Xbow
440, was used. The long-term stability and bounded nature of the GPS navigation solution
tends to make it a perfect sensor to couple with an IMU in a way that takes advantage of
the strengths of both systems while reducing the e ects of their weaknesses [14].
56
5.2.1 Coupling Architectures
For this system, three coupling architectures found in the literature were considered
[15]. These are shown in Figure 5.5. The IMU mechanization equations propagated the
desired states. For general operation, a six degree of freedom order was chosen to include
three body axes of both linear and rotational motion. The di erences in the considered
architectures stems from which measurements are used by the coupling  lter and whether
or not the GPS receiver is aided by the  lter. In a loosely coupled approach, the receiver
calculates a solution using at least four satellites, and its position and velocity are used as
measurements. This solution is independent of the coupling solution. One major disadvan-
tage of this method is that when the receiver cannot track enough satellites to compute the
independent position, velocity, and time (PVT) solution, no GPS aiding measurements are
available. Therefore the coupled solution accuracy is determined by the drift characteristics
of the low-cost IMU. In the closely coupled architecture, pseudorange and pseudorange rate
are used as measurements from any available satellite measurement. Therefore the updat-
ing is continued with less than 4 satellites. In this case, the GPS position calculation is
moved into the coupling  lter rather than being performed independently. Finally, a tightly
coupled system also uses pseudorange and pseudorange rate measurements, but the IMU
is used to aid the receiver in tracking GPS satellite signals. This approach increases the
capabilities of tracking in GPS harsh environments, particularly dynamic maneuvers. How-
ever, it requires the use of a GPS receiver that can accept the aiding measurements. The
nomenclature for these methods has changed over time but these labels e ectively delineate
the architecture being used. For this work, closely coupled integration is chosen since it
57
continuously updates errors but does not require GPS hardware to make use of the IMU
aiding information.
Figure 5.5: Coupling Architectures
5.2.2 INS Mechanization Equations
Since the coupling algorithm corrects solutions generated by the IMU, a navigation
solution is given by these corrected measurements from the IMU mechanization equations.
The IMU combined with its mechanization is called an inertial navigation system (INS). This
system includes the raw IMU measurements and the technique used to convert these mea-
sures into navigation states. The IMU measurements of angular rates, !bib = [gx gy gz ]T,
and speci c force, fbib = [fx fy fz ]T, are used to drive these states as shown in Equa-
tions (5.14 - 5.17). To simplify notation,  bib is the skew-symmetric form of the angular rate
58
vector !bib such that
 bib =
2
66
66
4
0  gz gy
gz 0  gx
 gy gx 0
3
77
77
5
(5.13)
Taking  t as the IMU measurement time step, the coordinate transformation matrix, Ceb,
from the body frame, b, to the ECEF frame, e, is propagated as [14]
Ceb = Ceb 
 
I3 +  bib t
 
 ( eie Ceb)  t (5.14)
The speci c force vector is translated to the ECEF frame by
feib = Ceb fbib (5.15)
With the acceleration due to gravity and earth rotation given as geb, the velocity state
is propagated as
veeb = veeb + (feib +geb 2 eie veeb)  t (5.16)
The position vector can then be updated as
reeb = reeb +veeb t (5.17)
59
Simple Euler numerical integration schemes were chosen for ease of computation. Al-
though the error due to this numerical approximation would tend to accumulate over long
time periods, with GPS corrections this e ect is reduced.
Equations (5.14 - 5.17) are processed whenever a new measurement is available from
the IMU. As correction measurements (GPS observables) are available, the position (reeb),
velocity (veeb), and coordinate transformation matrix (Ceb) are updated to include these
corrections. This correction process is detailed in subsequent sections. By periodically
applying the corrections, the INS solution represents the system?s best estimate of the
antenna?s navigation state at the current time. In the implementation described in Sections
5.2.4 and 5.2.5, this INS solution is used to linearized the state and measurement equations.
5.2.3 State Vector
The GPS/INS integration scheme includes 17 states used to keep track of the vehicle?s
motion, shown in Equation (5.18).
x =
2
66
66
66
66
66
66
66
66
66
66
4
  eeb
 veeb
 reeb
ba
bg
cdtu
_cdtu
3
77
77
77
77
77
77
77
77
77
77
5
(5.18)
These states are partitioned into six groupings of similar variables. The state vector in-
cludes three components of attitude error,  eeb, which are used to update the coordinate
60
transformation matrix from body to ECEF frames, Ceb; three components of velocity error
in the ECEF frame, veeb; and three components of position error in the ECEF frame, reeb.
The three accelerometer biases, ba, and three gyro biases, bg, are included and modeled as
constants with process noise. The receiver clock bias, cdtu, and drift, _cdtu, are included
with the clock drift modeled as a constant.
5.2.4 State Dynamics
Since the INS provides a set of nonlinear dynamic equations to propagate the state
vector in Equation (5.18) and the GPS measurements are nonlinear functions of these
states (to be described in Section 5.2.5), an Extended Kalman Filter (EKF) can be used in
the integration scheme [14]. In order to use the coupling architecture in an EKF form, the
equations given in Section 5.2.2 are linearized for propagation of the state errors ( x) and
its error covariance matrix (P). The state vector includes errors of the attitude, velocity,
and position states but keeps track of the actual IMU bias and receiver clock states rather
than their errors. The form of the state equation is
_ x = F x+ws (5.19)
61
where
F =
2
66
66
66
66
66
66
66
66
4
  eie O3 O3 O3 Ceb O2
F21  2 eie F23 Ceb O3 O2
O3 I3 O3 O3 O3 O2
O3 O3 O3 O3 O3 O2
O3 O3 O3 O3 O3 O2
OT2 OT2 OT2 OT2 OT2 F66
3
77
77
77
77
77
77
77
77
5
(5.20)
F21 =
2
66
66
4
0  fz fy
fz 0  fx
 fy fx 0
3
77
77
5
(5.21)
F23 = 2g0re2
en
reebreTeb 
 reeb  (5.22)
F66 =
2
4
0 1
0 0
3
5 (5.23)
where Im represents an identity matrix of size m m, O3 represents a 3 3 null matrix,
O2 is a 3 2 null matrix, and ws is the process noise included in the system derivation.
The state transition matrix at the current time step, k, can be approximated as
 k = I17 +F t (5.24)
This allows for the propagation of the state correction as
 x k =  k x+k 1 (5.25)
62
With the state transition matrix, the covariance matrix can be updated to
P k =  kP+k 1 Tk +Q (5.26)
Here Q is the system noise covariance matrix which is assumed to be a diagonal matrix
with entries given in Table 5.1. These values were tuned by trial and error to achieve the
desired performance.
Table 5.1: System Noise Covariance Matrix Values
State Value Units
Attitude Errors 0.01 (rad/s)2
Velocity Errors 0.01 (m/s)2
Position Errors 0.0005 (m)2
Accelerometer Bias Errors 0.000001 (m/s2)2
Gyro Bias Errors 0.000001 (rad/s)2
Clock Bias Error 0.1 (m)2
Clock Drift Error 0.01 (m/s)2
5.2.5 State Measurement Relations
The GPS/INS integration scheme uses measurements from all available satellites to
update the state corrections as
z = H x+wm (5.27)
where H is the state-measurement matrix and wm is the measurement error included for
derivation.
63
For satellite s, pseudorange,  s, and pseudorange rate, _ s, measurements are provided
by the receiver. The rate measurement comes from the receiver?s Doppler frequency mea-
surement, fDs. This measure is transformed to a rate measurement as
_ s = cf
L1
fDs (5.28)
where c is the speed of light and fL1 is the L1 carrier frequency. For each satellite, s, these
quantities ( s and _ s) make up the measurements available to the navigation system at each
measurement update.
Estimates of these measurements are calculated from the corrected INS solution, which
is the system?s best estimate of the state just before the measurement, and the GPS posi-
tion and velocity described in Sections 2.3.3 and 2.3.3. These estimated pseudorange and
pseudorange rate measurements are calculated as
^ s =
q 
reeb rees T  reeb rees +cdtu (5.29)
^_ s = (reeb rees)T (veeb vees)q 
reeb rees T  reeb rees 
+ _cdtu (5.30)
The di erences in the measurements ( s and _ s) and estimated measurements ( s and
_ s) give the measurement innovations. These di erences are used both to correct the state
estimates and to run the fault detection and exclusion algorithm, described in Section 5.3.
These innovations are calculated as
 z s =  s ^ s (5.31)
 z_ s = _ s ^_ s (5.32)
64
The innovations represent the new information provided to the navigation EKF. If the
estimated measurement is close to the actual measurement, the innovation corresponding to
this measurement will be low. In this case, the portion of the navigation  lter corresponding
to this innovation will have little e ect on altering the state correction ( x). Alternatively
if the estimated state has drifted signi cantly, the estimated and actual measurements will
diverge and the corresponding innovation will be relatively large. Thus the navigation  lter
will tend to alter the state correction to drive the innovation closer to zero.
Placing all of the innovations in a column vector forms the measurement vector as
z =
2
66
66
66
66
4
 z 1
...
 z_ 1
...
3
77
77
77
77
5
(5.33)
The integration equations use a standard Extended Kalman Filter algorithm with the
states described in Equation (5.18). The system is placed in linear form by di erentiating
the estimated measurements in Equations (5.29) and (5.30) with respect to each of the
states. Approximations for this measurement state relationship matrix are made in [14]
and are used here. For each row corresponding to a pseudorange, the resulting H matrix
entry is
H s = [O3 O3 usu O3 O3 1 0 ] (5.34)
65
where usu is the unit vector from the satellite s to the user. Similarly, for each row corre-
sponding to a pseudorange rate, the resulting H matrix entry is
H_ s = [O3 usu O3 O3 O3 0 1 ] (5.35)
Assuming uncorrelated measurements, the variances of these measurements are used
as entries along the diagonal of the measurement covariance matrix, R. The variance for a
pseudorange measurement is described in Equation (2.50). The variance for a pseudorange
rate measurement is assumed to be a function of the frequency-lock-loop (FLL) that gen-
erates the Doppler frequency measurement used in Equation (5.28). It is therefore taken
from values given in [21] where the FLL thermal variance is taken to be
 tFLL =  L2 T
s
4FBn
C=N0
 
1 + 1TC=N
0
 
(5.36)
with parameters de ned by the tracking loop architecture, shown in Table 4.1. The Doppler
frequency variance is then
 _ =
s
 2tFLL + f
2e
9 (5.37)
where fe is the dynamic stress error, taken to be 3 m/s. At a given time step, if a set of
measurements are available, the measurement innovation vector, z, is generated and the
Kalman gain calculated as
K = P k HT
 
HP k HT +R
  1
(5.38)
66
where the measurement covariance matrix (R) is a diagonal matrix with elements corre-
sponding to the pseudorange and pseudorange rate variances, described in Equations (2.50)
and (5.36), respectively. Both of these variances for each satellite signal are functions of the
signal?s C=N0, which represents the strength of the signal at the antenna. Lastly, the state
error and covariance matrix are updated using the measurements.
 x+k =  x k +Kz (5.39)
P+k = (I17 KH)P k (5.40)
5.3 Fault Detection and Exclusion
Fault detection refers to the ability of a system to determine if there is a signi cant
degradation in the GPS solution. The user is alerted if the system cannot guarantee solution
accuracy to a certain level [14]. Fault exclusion provides a means to remove the fault
and continue operation with the desired accuracy without having to alert the user. These
schemes are used for critical applications where GPS availability and accuracy are necessary,
such as aircraft landings. Traditionally these methods are used to detect and correct satellite
failures. In this thesis, the errors being studied cause sporadic measurement outliers due
to tree cover. These outliers are caused by large instantaneous jumps in multipath and
corresponding drops in the signal power at the antenna. They are considered outliers
since these measurements contain errors that cause the range or range rate to be outside
the expected accuracy of the GPS observables. Therefore a technique is used to exclude
outlying measurements and provide a more robust solution.
67
5.3.1 Normalized Innovation Parameter
Since the errors in multipath and signal attenuation a ect the receiver measurements
and these errors enter the solution through the innovation (described in Section 5.2.5),
the innovation was chosen as the detection parameter. However, since the innovation is a
function of measurement noise and not just these errors, jumps might not necessarily imply
a fault. Therefore a normalization routine was implemented to provide more consistent
operation. Each innovation is normalized by the square root of the C matrix diagonal,
where C = HPHT +R. The C matrix is calculated as part of Equation (5.38) and therefore
does not greatly increase the computational burden. Thus, the normalized innovation is
yi = zipC
ii
(5.41)
5.3.2 FDE Threshold
The normalized innovation parameters in Equation (5.41) are compared to a previously
determined threshold, yt. The measurement and its corresponding rows in the H and R
matrix are removed if it is faulty, i.e. yi > yt. The threshold value is determined using
the assumption that the normalized innovations are normally distributed with zero mean
and unit variance [14]. The FDE requires the designer to set a false alarm probability.
This is the probability that a non-erroneous measurement is rejected. Thus information
that could have been useful is rejected. From this probability, a threshold is calculated
by selecting a scaling on the normal distribution standard deviation to exclude the false
alarm probability as shown in Figure 5.6. By decreasing the false alarm probability, the
threshold increases and fewer outlier measurements will be rejected. The value used in this
68
work was 3, suggested from [14] to include 99.73 % of genuine measurements. This value
corresponds to a false alarm probability of 0.27 %, which is relatively high. However, since
this is a snapshot method, the e ect of rejecting a good measurement would only a ect a
single epoch.
Figure 5.6: Fault Detection and Exclusion Threshold Selection
5.4 Summary
This chapter details the methods investigated to compensate for the errors inherent
in general GPS operation and in operation in shadowed environments. For normal GPS
operation, the code di erential and carrier smoothed methods were presented. However, for
short baselines the performance of the code di erential algorithm is not greatly improved
over position di erencing due to the magni cation of noise and unmodeled errors in the
double di erence generation. The carrier smoothed code method was presented as an im-
proved positioning method for general operation. Also, this smoothing technique can be
69
adapted to other cases where pseudorange measurements are used. This decreases the noise
in the range measurement for improved performance. However, these techniques are not
viable in shadowed environments since satellite visibility is compromised and large multi-
path changes cannot be di erentially corrected. In heavy foliage, the likelihood of both
rover and base receivers accurately tracking at least  ve common satellites is low. As an
alternative to di erential operation, the low-level coupling architectures were described and
compared. The loosely coupled method has disadvantages similar to the di erential tech-
niques; namely, satellite visibility can cause the corrections to be unavailable. The closely
coupled method overcomes this di culty by using any raw GPS observables from visible
satellites. Therefore it is capable of supplying corrections even in low satellite visibility
situations. Alternatively, a tightly coupled method was described but its requirement for
special GPS equipment make it more di cult to implement, and therefore it is not explored
further in this thesis. In addition to coupling, a fault detection and exclusion method was
used to increase the integrity of the system?s solution. This is accomplished by rejecting
measurements considered to be outliers by comparing the normalized innovation parameter
to a pre-determined threshold.
70
Chapter 6
Results
To test the system, a path was chosen to include clear and heavy foliage areas, shown
in Figure 6.1. The route was located in Auburn, AL and began along a relatively clear
roadway, shown in the top of Figure 6.2. This section allows for good initialization of the
coupling  lter. It also gives a segment for comparing all methods without dealing with
the e ects of shadowed environments. Near the southern portion of the path, the test
included heavy foliage, shown at the bottom of Figure 6.2. This segment degraded signal
accuracy and satellite visibility, even dropping to only one satellite being reported by the
GPS receiver. This path was chosen to test the system in varying conditions and to compare
the system with stand-alone GPS solutions.
6.1 Pseudorange Calculated Position
As a baseline comparison, only the reported L1 pseudoranges and satellite positions
from broadcast ephemerides were used to calculate the receiver position. This method was
described in detail in Section 2.4. Due to the foliage conditions, there were points at which
no solution could be generated since too few satellites were reported. The overall position is
shown in Figure 6.3, and a zoomed segment is shown in Figure 6.4 detailing the operational
di culties of using GPS in foliage environments. As is clear from Figures 6.3 and 6.4,
this solution is unacceptable for position critical vehicle operations since position jumps
of approximately 100 meters occur in several places. Notice these jumps also occur in the
shadowed portion of the track located along the southern straight section.
71
Figure 6.1: Dynamic Path
72
Figure 6.2: Clear Sky and Heavy Foliage Environments Along Dynamic Path
73
Figure 6.3: Pseudorange Calculated Positions
Figure 6.4: Zoom on Pseudorange Calculated Positions
74
The zoomed region shown in Figure 6.4 shows one of the characteristic e ects of GPS
operation in shadowed environments. As a heavily shadowed satellite drops in and out of
visibility, it causes distinct spikes to be generated in a consistent direction. The magnitude
of these jumps is large indicating the large amount of shadowing a ecting the receiver in
this area.
6.2 Closely Coupled Implementation Results
Using the closely coupled algorithm (described in Section 5.2) by including the IMU
allows the navigation system to bridge GPS signal degradation. For this implementation,
an automotive-grade IMU (Crossbow 440) was used. This method also increases the system
solution rate to 50Hz, which is bene cial for vehicle control. The solution improves over
the pseudorange position since the e ect of jumps due to foliage is  ltered out, as seen in
Figure 6.5. The inclusion of the IMU in the  lter allows the navigation system to remain
on the path even though the shadowing e ects would cause the solution to jump away from
the path.
However, the solution still su ers from position and velocity jumps due to tree cover
errors on the GPS measurements. These appear due to the innovation outliers, and their
e ect can be seen along the straight southern portion of the path shown in Figure 6.6.
At the western jump, a satellite is reported at a single time epoch, and its innovation is
large compared to the other measurements (magnitude of  7 compared to < 2). In this
case the pseudorange measurement does not introduce the jump. The error occurs in the
pseudorange rate measurement. This fault causes erroneous corrections to be made to the
velocity in the navigation  lter. This velocity then propagates to the jump in position
75
Figure 6.5: Zoom on Closely Coupled Positions
solution which is corrected over time. This is shown as the stepping in Figure 6.6 where
the introduced velocity error causes the position solution to diverge over time before it is
corrected.
6.3 Closely Coupled with FDE Implementation Results
By including the normalized innovation monitoring described in Section 5.3, these
erroneous jumps are removed as shown in Figure 6.7. This implies that the correct outlier
was detected and removed from the solution. As can be seen, this system is able to bridge
parital signal loss while still estimating state corrections. In this case, the pseudorange
rate innovation causes a spike in the normalized innovation parameter corresponding to a
shadowed signal. The e ects of this spike are removed, and operation is continued without
that measurement. These results show the method is quite robust for operation in foliage
environments like those encountered in the test path.
76
Figure 6.6: Comparison of Closely Coupled and Pseudorange Solutions
Figure 6.7: Comparison of Closely Coupled with FDE and Pseudorange Solutions
77
The normalized innovations are shown for the entire path as well as this portion of the
path in Figures 6.8 and 6.9.
Figure 6.8: Normalized Innovations Over Entire Path
Figure 6.9: Normalized Innovations Over Zoom Segment
These plots show the threshold and the distribution of the normalized innovations. It
is noted that during the middle portion of Figure 6.8, the shadowing e ects are seen as a
78
spreading of the normalized innovations. The innovation thrown out at the beginning of
Figure 6.9 is the one that causes the jump at the left part of Figure 6.6.
6.4 Comparison of Methods
To compare the methods of operation in shadowed environments, the paths for each
method are plotted together above aerial photography via [1] in Figures 6.10 thru 6.12. As
can be seen, the pseudorange positions are subject to large jumps in solution in foliage areas
and even some jumps in clear sky areas. The GPS/INS closely coupled solution  lters these
jumps but also introduces jumps of lower magnitude since poor GPS measurements are still
used in the navigation EKF. However, these jumps are removed with the FDE technique as
seen in Figure 6.7.
Figure 6.10: Overlay of Positioning Methods In Heavy Foliage
79
Figure 6.11: Overlay of Positioning Methods In Clear Sky
Figure 6.12: Overlay of Positioning Methods In Light Foliage
80
6.5 Summary
In this chapter, implementation results are presented for the closely coupled system
with and without fault detection and exclusion described in Chapter 5. Comparison is made
to stand-alone GPS operation by calculating pseudorange positions at every measurement
update, wherever enough information was available. It is shown that the FDE method aids
in navigation due to its removal of error-causing outliers in the measurement domain.
81
Chapter 7
Conclusion
This thesis has shown that GPS operation in shadowed environments such as heavy
foliage poses diverse challenges. Although GPS receivers are subject to various error sources
from satellite, atmospheric, and local e ects, shadowed environments add complexity to the
signal degradation. The work presented in this thesis develops a method of analysis for the
errors speci c to shadowed environments, namely the rapidly changing multipath and the
signal attenuation caused by signal degradation. These sources of error a ect the ability
of a receiver?s tracking loops to e ectively track the incoming signal and thus negatively
impact the receiver?s performace in generating measurements. By using the error monitoring
parameters (M1 and C=N0) described in Chapter 4, the e ects of these error sources can be
tracked and visualized. This method makes error analysis more e cient than attempting
to accurately determine the exact range error e ects to which the GPS receiver is subject.
A visualization technique for both a single satellite signal and for all satellites in view was
developed to aid in the analysis of GPS operation in heavy foliage environments. As was
shown in Chapter 4, the presence of foliage between the GPS receiver and the broadcasting
satellite causes quickly varying and spontaneous drops in signal strength and also large
changes in the signal multipath. The occurrence of these errors was seen to correspond to
each other. This analysis led to the conclusion that foliage e ects on the GPS signal was a
quickly changing error mode and thus a snapshot method of fault detection could be used.
82
With insight into the errors a ecting navigation operation, methods were detailed to
handle such cases in Chapter 5. These methods included di erential operation and iner-
tial sensor integration. The performance of code-phase di erential methods were compared
and their shortcomings in shadowed environment operation described. Three integration
methods (loosely coupled, closely coupled, and tightly coupled) were described and their
advantages and shortcomings compared for operation in foliage. The closely coupled im-
plementation was chosen since this architecture makes use of all available GPS observables
along with the outputs of an IMU mechanization. This level of integration allows for oper-
ation in lower satellite visibility, which is common in shadowed environments. Comparison
of the closely coupled navigation system with raw GPS processing (independent of other
sensors) showed that integration reduced the large erroneous jumps in the navigation solu-
tion by combining GPS measurements with an INS calculated solution. It was also shown
in Chapter 5 that the closely coupled architecture alone could not operate e ectively in a
shadowed environment. Due to the nature of the foliage-speci c errors, a snapshot inno-
vation  ltering scheme was presented that allowed for normalized outlier rejection to add
integrity to the navigation solution. Inclusion of the fault detection and exclusion (FDE)
method further improved the operation of the navigation system during particularly dense
foliage cover.
7.1 Future Work
The next step in this research is to package the algorithms and hardware into a deploy-
able module for use on actual systems. This platform has been investigated and advances
will be forthcoming.
83
Further analysis of the multipath e ects could improve the processing by allowing
the navigation processor to estimate these errors. Therefore rejection could be avoided
by removal of these error sources. Improved multipath monitoring could also be extended
to other applications such as more advanced di erential processing methods like real-time
kinematic (RTK). This improvement would be bene cial since local e ects are errors that
di erential methods cannot intrinsically remove.
The navigation methods presented in this thesis can also be improved upon by adding
more sophisticated integrity monitoring schemes to the processing. With these statistical
techniques, more explicit bounds on the system integrity would be available. Therefore the
bounding thresholds could be set as a function of desired performance. Further performance
improvement would come by using a receiver capable of accepting aiding for signal tracking.
This could be developed as part of a software receiver platform as this functionality is not
available to receivers used in Auburn?s GPS and Vehicle Dynamics Laboratory.
84
Bibliography
[1] Gps visualizer, http://www.gpsvisualizer.com/.
[2] National geodetic survey precise gps orbits, http://www.ngs.noaa.gov/orbits/.
[3] Bhamer, Ionosphere, http://en.wikipedia.org/wiki/Ionospherevg, May 2007.
[4] Sunil B. Bisnath and Richard B. Langley, Pseudorange multipath mitigation by means
of multipath monitoring and de-weighting, Proceedings of KIS 2001 (Ban , Alberta),
June 2001.
[5] Michael H. Bruch, G.A. Gilbreath, J.W. Muelhauser, and J.Q. Lum, Accurate waypoint
navigation using non-di erential gps, Tech. report, Space and Naval Warfare Systems
Center, 2002.
[6] Rob Conley, Ronald Cosentino, Christopher J. Hegarty, Elliott D. Kaplan, Joseph L.
Leva, Maarten Uijt de Haag, and Karen Van Dyke, Performance of stand-alone gps,
Understanding GPS: Principles and Applications (Elliott D. Kaplan and Christopher J.
Hegarty, eds.), Artech House, 2006.
[7] Arinc Research Corporation, Navstar interface control document,
http://www.navcen.uscg.gov/gps/geninfo/IS-GPS-200D.pdf, April 2000.
[8] Tom Creel, Arthur J. Dorsey, Philip J. Mendicki, Jon Little, Richard G. Mach, and
Brent A. Renfro, New, improved gps: The legacy accuracy improvement initiative, GPS
World (2006).
[9] Donghai Dai, Todd Walter, C.J. Comp, Y.J. Tsai, P.Y. Ko, Per Enge, and J. David
Powell, High integrity multipath mitigation techniques for ground reference stations,
Proceedings of the 1997 ION GPS Conference, 1997.
[10] Nicholas J. DeCesare, John R. Squires, and Jay A. Kolbe, E ect of forest canopy on
gps-based movement data, Wildlife Society Bulletin 33 (2005), no. 3, 935{941.
[11] Jay A. Farrell and Matthew Barth, The global positioning system & inertial navigation,
McGraw-Hill, 1998.
[12] Julio C. Farret and Marcelo C. Santos, An alternative method for detection and mit-
igation of static multipath in l1 carrier phase measurements, Proceedings of the ION
NTM 2001, 2001.
[13] Garmin, Gps glossary, http://www8.garmin.com/aboutGPS/glossary.html, 2008.
85
[14] Paul D. Groves, Principles of gnss, inertial, and multisensor integrated navigation
systems, Artech House, 2008.
[15] Christopher R. Hamm, IV Warren S. Flenniken, David M. Bevly, and Daniel E. Law-
erence, Comparative performance analysis of aided carrier tracking loop algorithms in
high noise/high dynamic environments, Proceedings of the ION GNSS (Long Beach,
CA), ION, September 2004.
[16] Hisashi Hasegawa and Tetsuhiko Yoshimura, Estimation of gps positional accuracy un-
der di erent forest conditions using signal interruption probability, Journal of Forestry
Research 12 (2007), 1{7.
[17] Stephen Hilla, Plotting pseudorange multipath with respect to satellite azimuth and
elevation, GPS Solutions 8 (2004), 44{48.
[18] J. J. Spilker Jr., Foliage attenuation for land mobile users, Global Positioning System:
Theory and Application, Volume 1, American Institute of Aeronautics and Astronau-
tics, Inc., 1996.
[19] , Gps navigation data, Global Positioning System: Theory and Applications,
Volume 1 (Bradford W. Parkinson and James J. Spilker Jr., eds.), American Institute
of Aeronautics and Astronautics, Inc., 1996.
[20] , Satellite constallation and geometric dilution of precision, Global Positioning
System: Theory and Application, Volume 1, American Institute of Aeronautics and
Astronautics, Inc., 1996.
[21] Elliott D. Kaplan and Christopher J. Hegarty (eds.), Understanding gps: Principles
and applications, Artech House, 2006.
[22] G. Lachapelle, J. Hendriksen, and T. Melgara, Seasonal e ect of tree foliage on gps
signal availability and accuracy for vehicular navigation, Proceedings of the 7th In-
ternational Technical Meeting of The Satellite Division of the Institute of Navigation
(Salt Lake City, UT), September 1994, pp. 527{532.
[23] David Gary Lawrence, Aircraft landing using gps, Phd dissertation, Stanford Univer-
sity, September 1996.
[24] Helena Leppakoski, Heidi Kuusniemi, and Jarmo Takala, Raim and complementary
kalman  ltering for gnss reliability enhancement, Proceedings of the PLANS Confer-
ence, 2006.
[25] Pratap Misra and Per Enge, Global positioning system: Signals, measurements, and
performance, second ed., Ganga-Jamuna Press, 2006.
[26] Dignus-Jan Moelker, Multiple antennas for advanced gnss multipath mitigation and
multipath direction  nding, Proceedings of ION GPS-97, September 1997.
86
[27] Bradford W. Parkinson, Gps error analysis, Global Positioning System: Theory and
Application, Volume 1 (Bradford W. Parkinson and James J. Spilker Jr., eds.), Amer-
ican Institute of Aeronautics and Astronautics, Inc., 1996.
[28] Bradford W. Parkinson and Per K. Enge, Di erential gps, Global Positioning System:
Theory and Application, Volume 2 (Bradford W. Parkinson and James J. Spilker Jr.,
eds.), American Institute of Aeronautics and Astronautics, Inc., 1996.
[29] C. Rocken, C. Meertens, B. Stephens, J. Braun, T. VanHove, S. Perry, O. Rudd, and
M. McCallumand J. Richardson, Unavco academic research infrastructure (ari) receiver
and antenna test report, Tech. report, UNAVCO, 1996.
[30] P. Sigrist, P. Coppin, and M. Hermy, Impact of forest canopy on quality and accuracy
of gps measurements, Int. J. Remote Sensing 20 (1999), no. 18, 3595{3610.
[31] Robert F. Stengel, Optimal control and estimation, Dover Publications, 1994.
[32] Byron D. Tapley, Bob E. Schutz, and George H. Born, Statistical orbit determination,
Elsevier Academic Press, 2004.
[33] William Travis, Robert Daily, David M. Bevly, Kevin Knoedler, Reinhold Behringer,
Hannes Hemetsberger, Jrgen Kogler, Wilfried Kubinger, and Bram Alefs, Sciautonics-
auburn engineering?s low-cost high-speed atv for the 2005 darpa grand challenge, Jour-
nal of Field Robotics 23 (2006), no. 8, 579{597.
[34] United States Naval Observatory (USNO), Gps constellation status,
ftp://tycho.usno.navy.mil/pub/gps/gpstd.txt.
[35] Phillip W. Ward, John W. Betz, and Christopher J. Hegarty, Gps satellite signal
characteristics, Understanding GPS: Principles and Applications (Elliott D. Kaplan
and Christopher J. Hegarty, eds.), Artech House, 2006.
[36] , Satellite signal acquisition, tracking, and data demodulation, Understanding
GPS: Principles and Applications (Elliott D. Kaplan and Christopher J. Hegarty, eds.),
Artech House, 2006.
[37] IV Warren S. Flenniken, Modeling inertial measurement units and analyzing the e ect
of their errors in navigation applications, Master?s thesis, Auburn University, December
2005.
[38] Hayati Zengin and Ahmet Yesil, Comparing the performances of real-time kinematic
gps and a handheld gps receiver under forest cover, Tech. report, Istanbul University,
Faculty of Forestry, Department of Forest Management Planning, 2004.
87
Nomenclature
 x generic vector
 s pseudorange error for satellite s
 sl pseudorange error for satellite s and L-band freqency l
 sl carrier phase error for satellite s and L-band freqency l
 longitude angle
 l carrier wavelength for L-band freqency l
 latitude angle
 sl carrier phase measurement for satellite s and L-band frequency l
 s pseudorange from user to satellite s
c tu;s clock o set of user, satellite
c speed of light
e Earth-Centered Earth-Fixed (ECEF) coordinate frame
e body coordinate frame
i inertial coordinate frame
Is1;2 ionospheric delay on L1, L2 frequency for satellite s
M sl carrier phase multipath for satellite s and L-band freqency l
88
M sl pseudorange multipath delay for satellite s and L-band frequency l
n navigation coordinate frame
Nsl carrier phase integer ambiguity for satellite s and L-band frequency l
reeb;es ECEF position vector of user, satellite
reb;es user, satellite distance from ECEF origin
reb user range to ECEF origin
rsb range from satellite to user
s satellite coordinate frame
tm time a signal is measured
Ts tropospheric delay for satellite s
ts time a signal leaves the satellite
x   x vector from  origin to  origin expressed in  coordinates
x  magnitude of x vector from  orign to  origin
xeb;es ECEF x coordinate of user, satellite
yeb;es ECEF y coordinate of user, satellite
zeb;es ECEF z coordinate of user, satellite
89
Appendix A
Coordinate Frames and Transformations
Part of navigation solution presentation is a clear understanding of what coordinate
frames are being used. Also, transformations between these coordinate systems make the
combination of di erent system types (like GPS and INS) compatible. A description of the
frames used in this work is given followed by conversions between them used in this work.
A.1 Coordinate Frames
It is necessary to clearly delineate what coordinate frame system is being used and
what frame measures are related to. This is particularly important in combining systems
that give values of the same quantity in di erent frames.
A.1.1 Earth-Centered Inertial (ECI) Frame
The ECI frame is approximately an inertial frame, meaning it can be su ciently ap-
proximated as a non-rotating coordinate frame. This is not strictly true but for the purposes
of the calculations used, it can be approximated as such. The necessity of describing this
frame is that this is the reference system to which inertial sensors measure. The origin
of the ECI frame is located at the center of the Earth. Its z axis is coincident with the
Earth?s axis of rotation. The x axis points in a non-rotating direction, usually taken to be
the vernal equinox [14]. The y axis completes the right-handed coordinate system. This
frame is illustrated in Figure A.1.
90
Figure A.1: Earth-Centered Inertial (ECI) Frame
A.1.2 Earth-Centered Earth-Fixed (ECEF) Frame
The ECEF frame is a frame in which the origin is the center of the earth, the z axis is
coincident with the Earth?s axis of rotation, the x axis points along 0 longitude, and the
y axis completes the right-handed coordinate system by pointing along 90 longitude. This
frame is illustrated in Figure A.2.
A.1.3 Navigation Frame
The local navigation frame is de ned with an origin on the user body (usually center of
mass or sensor coordinate system center). It is the point for which the navigation solution
is sought [14]. One axis points toward the local North direction, co-planar with the ECEF
z axis. The Down direction is de ned to be perpendicular to the Earth ellipsoid in the
direction of the center of the Earth. The East direction completes the right-hand coordinate
91
Figure A.2: Earth-Centered Earth-Fixed (ECEF) Frame
system. This frame is illustrated in Figure A.3. Note that there are singularities for this
frame at each pole because of the nature of the North direction de nition.
A.1.4 Body Frame
The body frame is the frame rigidly attached to the user body (assumed to be coincident
with the sensor coordinate system). The x axis points in the forward direction of the body,
the z axis downward on the body, and the y axis completes the right-hand coordinate
system. This frame is illustrated in Figure A.4.
A.1.5 Orbit Frame
The orbit frame is a planar system with the origin at the Earth?s center and inclined
from the polar axis by an inclination angle, i, as shown in Figure A.5. This frame is also
rotated about the polar axis by  , the longitude of the ascending node. The ascending node
92
Figure A.3: Navigation Frame
Figure A.4: Body Frame
93
is the point at which a satellite in the orbit crosses the equatorial plane in the positive z
direction. This angle is measured from the IERS reference meridian as shown in Figure A.5
[14]. Within the orbit plane, the x axis is de ned from the Earth center to the ascending
Figure A.5: Orbital Plane With Respect To Equatorial Plane
node and the y axis completes the right-hand coordinate system as shown in Figure A.6.
The orientation of the plane is de ned by the argument of perigee, !. The perigee is the
point of closest approach in the orbit.
Figure A.6: Orbit Frame
94
A.1.6 Satellite Frame
The satellite frame is de ned for satellitesas the coordinate frame centered at the phase
center of the satellite antenna. Due to its use in this work, directions of the coordinate axes
are not important since the origin location is the only important set of coordinates to be
determined.
A.2 Coordinate Transformations
Transformation between sets of coordinate frames is necessary since various quantities
used in the implementation are measured or given in di erent coordinate frames. Also,
the user may need the resulting coordinates in a more usable form for their application.
For example, ECEF coordinates are not as helpful as navigation coordinates for vehicle
operation. Also, most mapping applications use geodetic coordinates since the origin of the
navigation coordinates is di cult to correlate with positions on the Earth?s surface.
A.2.1 Coordinate Transformation Matrix
In order to relate the various coordinate frames, a transformation matrix is used. To
take a vector quantity between two coordinate frames  and  that is expressed in coordinate
frame  , x   , to be expressed in coordinate frame  , x   can be accomplished by
x   = C  x   (A.1)
95
The inverse operation is accomplished by using the transpose of the transformation
matrix, that is
C  =
 
C  
 T
(A.2)
If an angular rate vector is expressed in skew-symmetric form then it can be transformed
to another coordinate frame as
    = C      C  (A.3)
The rate of change of the transformation matrix is found as
_C  = C      (A.4)
A.3 ECEF and Geodetic Coordinates
To translate between ECEF coordinates and geodetic coordinates, the Earth?s ellip-
soidal model must be used. For these transformations, the transverse radius of curvature,
RE, is needed. This is the radius of curvature of constant latitude lines calculated as
RE (Lb) R0p1 e2 sin2L (A.5)
where R0 is the equatorial radius and Lb is the geodetic latitude. Therefore ECEF coordi-
nates can be calculated from geodetic coordinates as
xeeb = (RE (Lb) +hb) cosLb cos b (A.6)
96
yeeb = (RE (Lb) +hb) cosLb sin b (A.7)
zeeb =
h 
1 e2
 
RE (Lb) +hb
i
sinLb (A.8)
and geodetic coordinates can be calculated from ECEF coordinates as
sinLb = z
eeb
(1 e2)RE (Lb) +hb
tan b = y
eeb
xeeb
hb =
q 
xeeb 2 + yeeb 2
cosLb  RE (Lb) (A.9)
This solution has to be iterated through since RE is a function of latitude, and latitude and
height must be solved for iteratively.
97