zotero/storage/266N6YDD/.zotero-ft-cache

fcT'L INST. OF STAND & TEf" f.:]f.
IIIDS blfiT17

NIST PUBLICATIONS

Nisr United States Department of Commerce Technology Administration National Institute of Standards and Technology

N/ST Technical Note 1385

NIST Technical Note 1385
Global Position System Receivers
and Relativity
Neil Ashby Marc Weiss
Time and Frequency Division Physics Laboratory National Institute of Standards and Technology 325 Broadway Boulder, Colorado 80303-3328
March 1999
'4rES
U.S. DEPARTMENT OF COMMERCE, William M. Daley, Secretary
TECHNOLOGY ADMINISTRATION, Gary R. Bachula, Acting Under Secretary for Technology NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, Raymond G. Kamnier, Director

National Institute of Standards and Technology Technical Note Natl. Inst. Stand. Technol., Tech. Note 1385, 52 pages (March 1999)
CODEN:NTNOEF
U.S. GOVERNMENT PRINTING OFFICE WASHINGTON: 1999
For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, DC 20402-9325

.

Contents

1 Introduction

1

2. Theory

2

2.1 Pseudorange

2

2.2 Phase Time

2

2.3 Relevant Relativity

3

2.4 Common Misunderstandings

4

2.5 User Corrections

6

2.6 Relativistic Doppler Effect and Use of the GPS Carrier

7

2.7 Rotation Matrix

8

3. Time Tagging at Receiver

9

— 3.1 Time Tagging at Receiver General Prescription

9

3.2 Examples

11

Example 3.21 Multi- Channel Receiver, Time Tagging at Receiver

12

Example 3.211 Erroneous Use of ECEF Coordinates

16

Example 3.212 Iteration With a Succession of Inertial Frames

16

Example 3.213 The Benefit of Better Initialization

17

4. Time Tagging at Transmitter

18

4.1 Code Boundary Measurements

18

4.2 Phase Time

19

4.3 Accounting for Receiver Motion

20

— 4.4 Time Tagging at Transmitter General Prescription

22

4.5 Examples

24

Example 4.51 Time Tagging at Transmitter

24

Example 4.52 Time Tagging at Transmitter, Neglecting Doppler Corrections 29

5. Time Transfer Receiver

29

Example 5.1 Earth-Fixed Single Channel Timing Receiver

32

6. References

32

Appendix A. The Eccentricity Correction

35

Appendix B. Linearization Algorithm for Time Tagging at Multichannel Receiver

37

Appendix C. Linearization Algorithm for Time Tagging at Transmitters

39

Appendix D. Equivalence of Integrated Doppler and One- Way Light Time

42

Appendix E. Numerical Check on First-Order Doppler Correction

45

111

GLOBAL POSITIONING SYSTEM RECEIVERS AND RELATIVITY
Marc Weiss National Institute of Standards and Technology
Time and Frequency Division
325 Broadway, Boulder, CO 80303-3324
Neil Ashby
Department of Physics, Campus Box 390
University of Colorado
Boulder, CO 80309-0390
ABSTRACT
We illustrate the general methods for applying relativistic corrections needed by a GPS receiver in providing time or position to a user. We focus on estimating the geometric range delay, the time GPS signals require to propagate from the transmitter to the receiver in vacuum. We discuss several common uses of GPS including time tagging at the receiver, time tagging at the transmitter, and earth-fixed timing receivers. We include examples
which illustrate these cases numerically.
KEYWORDS
GPS receivers; relativity; navigation; gravitational redshift; second-order
Doppler shift; pseudorange.
1. INTRODUCTION
A GPS receiver must make two corrections that are related to relativity in order to provide time or position to a user. We discuss these corrections and focus mostly on estimating the geometric range delay to, the time for GPS signals to propagate from the
transmitter to the receiver in vacuum. Proper estimation of t^ is essential for solving for position or time. This is an application of the relativistic principle of the constancy of the speed of light, which states that electromagnetic signals travel in Euclidean straight lines
with speed c relative to an inertial reference frame. We present a few cases which apply to many common uses of GPS and illustrate them numerically.
We present the theory behind corrections with references given for any derivations not
done here. Through the derivations, we show that tlu> Interfac-e Control Document (ICDGPS-200) specifications, as issued by the Joint Program Office of the Global Positioning

System [1], consistently cover the requirements of relativity down to within 1 ns or less.
We show that the ICD specifications include relativity corrections with enough accuracy

for most applications [2,3]. In particular, we discuss the relativistic Doppler effect, the
formula for its instantaneous magnitude, and its relationship with typical GPS receiver operation. We also address the use of carrier-phase measurements, which is not discussed

in the ICD.

Relativistic effects which are not usually modelled in GPS are: (a) a small effect on

satellite clock rates due to the earth's oblateness and (b) a small time delay due to the

slowing of electromagnetic signal propagation in the earth's gravitational field. The earth's

oblateness contributes a very small constant rate correction for satellite clocks, which is so

small (less than 50 ps per day) that it can be neglected. It also causes a periodic variation in
the GPS satellite clock time having twice the orbital period and a peak-to-peak amplitude

less than 200 ps. For a GPS satellite-to-user Hnk the gravitational time delay is less than

200 ps [2].

Our goal is to give simple recipes, so that GPS users may understand how to im-

We plement receiver designs consistently with the requirements of relativity.

then back

up these prescriptions with full explanations and derivations, so that those interested can

understand where the prescriptions come from and what approximations are involved. The

result is that this document provides a brief tutorial on how relativity is apphed for users
of the GPS broadcast ephemerides.

2. THEORY

2.1 Pseudorange
Generally, a GPS navigation receiver measures the arrival times, on a local clock, of
signals from four or more different satellites. The receiver then uses these measurements
along with information in the navigation messages to solve for four unknowns: user position
{x,y,z) and the receiver clock offset from GPS time. The signals from the sateUites can
be thought of effectively as continuous timing signals arriving at the receiver from the
satellite clocks. The receiver has its own local clock for comparison. The user measures
either (1) the times in the received timing signals at a specific local clock time, or (2) the
arrival times on the local clock of a specific time tag in the received timing signals. The difference between the reception time (according to the local clock) and the transmission
time at one satellite (according to GPS time), multiplied by the speed of fight, is called
the pseudoTU'fige. These pseudoranges are used to solve for user position and time.

2.2 Phase Time In addition to times measured on real clocks, we can think of a particular phase point
in the wave train between transmitter and receiver as nvnnerically equal to some time r,

which follows the wave as it propagates. This "phase time" r is transmitted from the
satellite at SV clock time tsv and propagates to the receiver where it arrives at receiver clock time t^. Thus the pseudorange is c{tpi — tsv)-

2.3 Relevant Relativity
Several relativistic effects have already been incorporated into the GPS system so,

for the ordinary user of broadcast ephemerides, only two relativistic corrections must be

considered. First, the receiver must apply a correction to the transmitted time to account

for relativistic effects arising from orbit eccentricity of the transmitting satellite. This is

the Atr term defined in the ICD. Second, the finite and universally constant speed c of

signals propagating in a vacuum from a transmitter to a receiver, relative to an inertial

frame (the geometric path delay), must be accounted for. The receiver must also account

for ionospheric and tropospheric delay corrections, which we do not consider here.

Three relativistic effects are germane to GPS. Rates of clocks in GPS are adjusted

(as they are for International Atomic Time) to match the rate that clocks would run on

the geoid of the earth. The geoid is a surface of constant gravitational potential in the

rotating frame in which the effects (2) and (3) described below add to a constant value. The three relativity effects are as follows.

(1) Constancy of the speed of light and relativity of synchronization. GPS time is

defined using the principle of the constancy of c to synchronize an imagined system of

clocks everywhere in space in the neighborhood of the earth (this is called Einstein syn-
chronization). GPS satellite clocks are in principle adjusted to agree with this imagined system of clocks. This network of synchronized GPS clocks realizes a coordinate time, a

system of self-consistent time markers with which to label events. This definition of GPS

time requires a locally inertial coordinate system. GPS time is thus defined relative to an

earth-centered inertial coordinate system (an ECIF), but the rate is set to match the rate

at which clocks would run on the geoid. An ECIF is also used to simplify the paths of

signals propagating from satellites, since, with sufficient accuracy for GPS, light travels in

Euclidean straight lines at the speed c in vacuum relative to such inertial frames [2].
(2) Second-Older Doppler shift. A clock moving with respect to an ECIF runs slower
relative to coordinate time in that ECIF than if it were at rest in the ECIF. This is the

time dilation effect due to the magnitude of the relative velocity, sometimes called the

second-order Doppler effect.

GPS For satellites in

orbits, the fractional frequency offset

needed to compensate for this is approximately -1-8.3 x 10~^^ relative to the rate of clocks

on the earth's geoid.

A (3) Gravitational frequency shift.

clock at rest in a lower gravitational potcMitial

runs slower relative to coordinate time than if it were at rest in a higher potential. This

is called the gravitational red shift. Thus, standard clocks closer to the earth nin slowt^-

than standard clocks farther away, since the gravitational potential heconu^s nu^re negative

closer to the earth. Clocks on GPS satellites run faster than clocks at rest on the earth's surface. Thus GPS satellite clock frequencies need to be adjusted by a fraction of about
—5.3 X 10~^° relative to the earth's geoid, to compensate for this effect.

To compensate for the relativistic effects described in paragraphs (2) and (3) above
for circular orbits, and as a consequence of the requirement that GPS satellite clocks run at the rate that a standard clock on the geoid would run, atomic clocks in GPS satellites
are given a fixed fractional frequency offset of —4.4645 x 10~^^. The user does not have

to be concerned about these rate corrections.

These three relativistic effects, however, explain why the user must still apply two

corrections. The first relativistic correction is the Atr term defined in the ICD. This term

corrects the satellite vehicle (SV) clock offset due to any eccentricity in the vehicle orbit.

Eccentricity produces periodic time excursions from the nominal fixed fractional rate offset

of —4.4645 X 10~^°. The second correction applies to users' estimates of the geometric range

delay, the delay from the transmitter to the receiver if the signal traveled in vacuum. Since

ECEF an

frame rotates relative to any inertial frame, light does not travel in a Euclidean

ECEF straight line in an

frame. It is incorrect to calculate the geometric range in an

ECEF frame using the magnitude of the vector difference between the receiver's position

at reception and the transmitter position at transmission. The distance is most easily calculated in an ECIF where signals do travel in Euclidean straight lines at the speed of

light, the constant c.
Many ECI coordinate systems, differing by constant spatial rotations from each other,

may serve these purposes. All frames that have the same origin at the earth's center and

do not rotate with respect to the "fixed" stars will define simultaneity in the same way.

All such frames are equivalent for determining the propagation delay. Yet users usually

need to refer their positions to the earth. For example they may need their coordinates in
the WGS-84 coordinate system, the ECEF frame upon which the broadcast ephemerides are based. Users fixed on the earth who want time from GPS often know their coordinates as constants in the ECEF frame. In fact, most users would hke to know their coordinates in the ECEF frame.

2.4 Common Misunderstandings There has been confusion in the GPS community over a number of issues related to
relativity [5,6]. For example, arguments have been presented that we need to define a coordinate system in the rest frame of the receiver or perhaps in that of the transmitter. Since measurements are made in the receiver, it might at first appear necessary to use a coordinate system comoving with the receiver. However, such arguments overlook the fact
— — that relativistic effects indeed, any physical effects can be described in any coordinate
system.
Another argument that has been advanced is that we need to compute corrections in

the rotating coordinate system of the ECEF frame, since that is the coordinate system in
which the user ultimately wishes to find his/her position. The relativity of simultaneity

has sometimes been presented as a principle necessitating such constructions.

Although measurements are made with a particular instrument and an instantaneous

rest frame is associated with that instrument, it is not necessary to use that rest frame
to analyze the data. This is a common misunderstanding underlying many of the errors

in the literature: the belief that we must always use reference frames which are comoving with the transmitter or receiver. The user does not need a coordinate system whose origin

coincides with the user's trajectory in spacetime to analyze effects local to the user.
It is a common practice in physics to choose a convenient coordinate system to solve a

problem. This is particularly germane to relativity, where the principle of covariance states

that, if a problem is understood in one coordinate system, it can be understood in any

other by making appropriate transformations. Thus, even though neither the user nor the

transmitter are at rest in an ECIF, the user may choose an ECIF to simplify computations
of position, then transform those coordinates to the ECEF frame, the WGS-84 system.

The pseudorange that the user measures is the same regardless of which coordinate system

the user chooses. It is simpler to use an ECIF for computations than to try to make all

corrections in the rotating ECEF coordinate system.

The fundamental measurement of a GPS receiver is the pseudorange. This is a mea-

surement of the difference between the time associated with the phase of a GPS signal

We and the time on the local clock at the arrival of that phase.

will discuss why this

measurement has the same value regardless of which coordinate system the user chooses

to do calculations.

The phase of an electromagnetic wave is a relativistic scalar. This follows simply

because the phase is the relativistic scalar product of two fourvectors: the wave propagation

fourvector and the position fourvector. This means the phase has the same numerical value

in all reference frames. It does not matter whether this phase is evaluated in the earth-

centered inertial frame, the rest frame of the transmitter, or the rest frame of the receiver.
For GPS signals, GPS time is encoded in the transmitted pseudorandom noise (PR.N)
signal. The bits in the PRN code are bi-phase modulated on the L-band carrier. Thus

when the time arrives for new bit of the code, the phase of the carrier either reverses 180°

or it does not, corresponding to whether the next bit is one or zero. The receiver measures
the time of arrival of a GPS signal by comparing it with a locally generated replica of that
signal. The receiver correlates the received signal against the; local replica and examines

the correlator output to find the correlation maximum, then locks the local code replica

to the transmitted code. Maintaining this lock in general requires adjusting the frequency

and time offset of the code replica in the correlator relative to the local clock. Tlu^ riMciviM-

system then determines the time difference between the received phase tunc and the lime
on the local clock. It is as if the wave front of tlu^ GPS signal luus an unambiguous time

mai'ker on it. Thus, the user can measure the phase of the received signal, hence the time
as received, against the time of the receiver's local oscillator. The difference between the phase time as received and the time of the receiver's local oscillator, after multiplication by c to give a quantity with units of length, is the pseudorange.
We can see why the pseudorange is independent of which coordinate system the re-
ceiver uses to solve for position and/or time. This follows since the phase is a relativistic
scalar and the measurement is made when the receiver clock is co-located with that GPS
signal phase. The measurement is essentially the time difference between the local clock
time at a particular point on that clock's space-time path, and the GPS signal phase time
at the same spacetime point.
The navigation problem could be completely solved using only the ECEF frame. It is complicated in the ECEF frame because hght and hence GPS signals no longer travel in
R Euclidean straight lines. The time for a signal to travel from transmitter at to receiver
— at r is not |r R|/c. Similarly we could set up a coordinate system in either the rest frame
of the receiver or of the transmitter. However, we emphasize that the problem is greatly simplified by analyzing the problem in an ECIF, then transforming results appropriately
to the ECEF frame. This is simplified by the fact that the time coordinate in ECIF's as we have defined them is identical to that in the ECEF frame. Thus, for example, the
relativity of simultaneity is not an issue because time, and simultaneity, are defined the same way in all these coordinate systems.
2.5 User Corrections
At any arbitrarily chosen instant the ECEF frame coincides with an ECI frame having identical x, y, z axes, but not rotating. Removing the rotation of the earth from the ECEF
frame defines a coordinate system that is close enough to an inertial frame to serve for
estimating the path delay. As time passes, the ECEF system rotates while this ECI system
remains behind, so to speak. Any such fixed coordinate system may serve as an ECIF for
estimating geometric range, using the Euclidean distance between the coordinates in the
ECIF of the satellite at transmission and the coordinates of the receiver at reception. The time for the signal to travel this path, the geometric path delay tf). may be estimated to high accuracy as the geometric range divided by the defined speed of light c. The slowing of
electromagnetic signal speed in earth's gravitational field (the Shapiro time delay) amounts
to less than 200 ps delay for GPS users near earth. [2]
The principle of the constancy of the speed of light in an inertial frame requires that an ECIF be used for geometric path delay if it is calculated by dividing the Euclidean
distance by c. Using such an ECIF greatly simplifies the problem of solving for a GPS user's position or time. Whereas GPS is intended to provide users with their position or time in the ECEF system, if we use ECEF coordinates for geometric range with the
simple Euclidean distance metric and divide by c to obtain geometric path delay, we would

make significant eiiois since we would be ignoring motion of the receiver during the signal
propagation time. In general, a navigation solution requires pseudorange measurements from at least
four satellites. These are used to obtain geometric range delay estimates, which in turn are used to solve for position. To use the simplicity of the Euclidean distance metric that comes with an ECIF, all satellite positions at the transmission epochs must be transformed
into the common ECIF. The user may then solve for the receiver position in this coordinate frame, and for the GPS time t corresponding to this solution. Finally, the user must find
the receiver coordinates in the ECEF frame at the GPS time t by accounting for the rotation of the ECEF frame between the chosen moment tc at which an ECI frame is defined, and the GPS time t.
With a receiver using four satellites we generally have five different times of interest
as natural candidates for tc'- either a single GPS transmission time and four different
reception times (time tagging at transmission time), or a single reception time and four
diflPerent transmission times (time tagging at reception time). In any case an ECI frame
can always be found which coincides with the ECEF frame at a chosen instant of time, and in this sense the ECEF frame can determine an inertial coordinate system. This answers
the need for a coordinate system in which light travels in straight lines with speed c during
the transmission of the signals, and in which all the GPS satellite clocks can in principle be synchronized. We will give a number of examples below for using this general prescription. Our focus here will remain on relativistic corrections for GPS users. We do not consider
corrections for nonrelativistic delays: ionospheric plasma delays, tropospheric delays due to water vapor, multipath interference, or receiver system delays. Nonrelativistic delays must certainly be accounted for; however, here our topic is relativity.
2.6 Relativistic Doppler Effect and Use of the GPS Carrier Many receivers that make use of the GPS carrier use only the accumulated phase to
smooth code measurements. In this case the instantaneous relativistic Doppler shift need not be considered. As we show in Appendix D, the integrated carrier as received over one cycle of the code is equivalent to the change in geometric range over that interval. Thus if a receiver using accumulated Doppler shifts keeps track of geometric range, then it automatically keeps track of the relativistic Doppler shift. Nevertheless some receivers
may use the instantaneous frequency of the received signals for some purpose. Such a
receiver must account for the relativistic Doppler effect.
A receiver system that uses the instantaneous Doppler shift of the received cai'rier
signal must also correct for the frequency shift of the received signal within the framework of relativity. As the true range between the satelUtc and receiver changes duo to relative motion, the carrier frequency changes due to the Doppler effect. The relationship between the received frequency /, and the transmitted frequency, or proper frequency F. is derived

D in Appendix with accuracy sufficient for GPS. This relationship is [2,3,4]

NV / (1 + A^{rR)/c^ - v2/2c2)
. Nv

~

F{1 + A^{rT)/c^ - y^/2c^)
1

'

c

{^-^^J

where

V = transmitter velocity in ECI coordinates,
V = receiver velocity in ECI coordinates,
N = ECI unit vector in the direction of signal propagation,
A$ = gravitational potential difference between

a clock and a reference clock on the geoid.

If gravitational potential differences can be ignored, the relationship between received and transmitted frequency can be written without approximation as

/

F

7(v)(l-^) 7(V)(1-^)'

(1)

where

7(^) =

(2)

, 1- ^

The change in frequency is closely related to the rate of change of range. In fact, in
Appendix D we show that the relativistic Doppler frequency shift equation can be derived by differentiating the geometric path delay [4,6]. In Appendix D we have included for
completeness the effect of the gravitational frequency shift to order c~^, which is sufficient for users of the GPS. This proves the conceptual equivalence, within the framework of Special and General Relativity, of the methods of pseudorange (to obtain an instantaneous
A position solution) and integrated Doppler frequency (to obtain changes in position).
user of integrated Doppler frequency accounts for changes in the geometric path delay simply by integrating the relativistic Doppler shift and relating this observable quantity to the change in geometric range. This is discussed in detail in Appendix D.

2.7 The Rotation Matrix Ahnost always the user of the broadcast ephemerides will want to use coordinate
transformations which correspond to rotations of the coordinate system about the WGS-
84 z axis. For convenience we write out here an example of such a rotation matrix.
Consider an ECIF with z axis which coincides with the defined WGS-84 axis. Let the position coordinates of some point of interest in these ECI coordinates be denoted by

xeci

^ECI = Veci

(3)

zeci

Suppose that the ECI axes coincide with the ECEF axes at the time tc- To calculate the ECEF coordinates of the point at time t, note that the angle through which the ECEF coordinates have rotated is G{t) = fte{t-tc), where l^e = 7.292 115 146 7 x 10~^ rad/s
is the WGS-84 value of the earth rotation rate. Then the coordinates of the point in the
ECEF frame are given by

cos0(t) sin0(t)

XeCI

XeCEF = R {t)^ECI — -sin0(t) cos0(i)

Veci

(4)

1

ZECI

— This is sometimes called a "passive" rotation that is, a rotation of the coordinate axes

keeping physical position vectors unchanged. In eq (4) the time interval t — tc must be

small since the earth rotation rate actually varies. With t — tc as large as 3 s, errors no

mm larger than 0.05

are introduced.

3. TIME TAGGING AT THE RECEIVER

Generally, a GPS navigation user measures the arrival times of signals from four
different satellites and uses them to solve for four unknowns: position x. y, z and time
t. A common way of doing so is to make the four measurements at one instant at the
receiver ("time tagging at the receiver"). Another method is to use signals which left
the satellites at a common GPS time ("time tagging at the transmitters"). The latter
method is more complex, both because the signals are not received simultaneously, and
because GPS times on the satellite clocks must be obtained by applying transmitted clock
correction coefficients. That is, the clocks on the satellites are synchronized only if the
user applies transmitted clock corrections. The GPS satellite clocks themselves may differ from GPS system time by up to 1 ms in the underlying ECI reference frame [1. Sections
20.3.3.3.1.9-20.3.3.3.3.1]. Thus in the case of time tagging at the transmitters we must account both for the motion of the user during the intervals between transmission and
reception, and the differences between SV clock time and GPS system time. On the other
hand, as we shall show, time tagging at the transmitter can be used for navigation while avoiding most rotations.
— 3.1 Time Tagging at the Receiver General Prescription
(1) At a chosen reception time, measure the transmission times using the received pseudo-random noise (PRN) codes from four satellites. This gives us the times
tsv according to the SV clocks at transmission for the signals received at the
chosen reception time.
We need the GPS time of transmission for each satellite mid the
common reception time according to the (possibly biiuscd) local clock.
We obtain the time of reception from our local clock. Tlu^ transmission

time is encoded in the received signal. Since we are locked to the code, we

PRN can determine the offset of the locally generated

sequence required

to maintain lock.

(2) Apply the prescription of Table 20- IV of the ICD-200, for corrections to obtain
GPS system time i, at transmission.
This includes the eccentricity correction At^, described in the ICD200, section 20.3.3.3.3.1, and is an effect due to relativity. At^ is a
correction which applies to the SV clock; it is the same correction no
matter where the receiver is or how the receiver is moving. The value of
Atr will generally be different for the clocks in different satellites. For example, for a clock in a satellite of orbital eccentricity e, transmitting
at an instant when the eccentric anomaly is E, the correction to be
subtracted from the transmitted time is

Atr = Fx 10-^°ev^sin E ,

(5)

where

F = -4.442 807 633 s/m^/^

(6)

A is the semi-major axis of the satellite in meters, e is the orbital eccenE tricity, and is the eccentric anomaly. Values of these quantities are

transmitted in the navigation message.

(3) Compute each satellite's position in the ECEF frame at its transmission time.
This is a straightforward application of the equations of Table 20-IV

of the ICD-200 [1] for determining the a:, y, z coordinates of the satellite

at the instant of transmission. Note that the broadcast message gives the

ECEF satellite positions at the instant of transmission in

coordinates

specifically in the WGS-84 reference frame.

(4) Choose an ECIF for computation of the path delays.
This choice is arbitrary, but some choices are more convenient than
others. Simplification may sometimes occur if the choice is appropriately
made. A natural choice for the case of time tagging at the receiver is to freeze the ECEF frame at a time close to the expected GPS time of
reception.
(5) Transform the ECEF coordinates of each SV obtained in step (3) into the
chosen ECIF. This will normally require at least three, and perhaps four, of such
ECEF position vectors to be rotated. If we freeze the ECEF frame at
tlie instant of trHnsmission of on(^ of the satellites, only three of these

10

rotations will be needed. Alternatively, we may freeze the ECEF frame

at any instant to define an earth-centered inertial system. If we choose

the GPS time equal to the local clock's reception time, the later correc-

tions may be small.

GPS If the local clock is sufficiently close to

time

at the instant of reception, then the last step below, step (7), would not

need to be done.

(6) Solve the path delay equations for the receiver's position and time.

This can be done by linearizing the propagation delay equations

and solving them iteratively. Here we solve simultaneously for geomet-

ric range delay t]j and user position. The user's motion through the

ECIF during propagation is included in these to estimates. Note that

other contributions to path delay which we do not consider here, such

as ionospheric and tropospheric delays, should be incorporated in this

We process.

use an initial estimate of position to solve for the linearized

corrections to receiver position and the time offset. If these corrections
are not small enough for the user, they may be used to obtain a next

estimate of position and time. This, in turn, may be used to obtain

the next linearized correction for position and time offset. This iterative

process converges rapidly.

(7) Rotate the user's position coordinates into the ECEF reference frame.

After finding the user's position in the chosen ECI coordinate sys-
tem and the correct GPS system time, the receiver's ECI position coordinates are rotated into the ECEF reference frame at the instant cor-

responding to the measured reception time. If the ECI frame chosen

ECEF in step (4) is sufficiently close to the

frame at the instant of

reception, this step may be omitted.

3.2 EXAMPLES
In these examples we use orbital data representing the GPS constellation as it ap-
proximately was in March of 1995. Values of the semi-major axes, eccentricities, and inchnations have been randomized within very small, but reasonable limits. Keplerian
orbit algorithms are used to propagate the SV orbits forward in time and to accurately compute SV eccentric anomalies. Values of Atsvj are randomized using a uniform distribution within a 1 ms range. The SV identification numbers used are fictitious. However, the actual transmitter and receiver data used in these examples are selected from among many possibihties resulting from a simulation of the entire GPS constellation. The data are
thus obtained from a simulation where ''truth" is accurately known and can \yc c t)mpai'ed
with navigation solutions.
11

Example 3.21 Mult i- Channel Receiver, Time Tagging at the Receiver
We choose four satellites from this constellation in view from our receiver location with
elevations above the horizon greater than 20°. Suppose the receiver is truly at geocentric
m R ECEF latitude and longitude 35° N, 0° E, and distance = Q 378 136.3 from earth's
center. Suppose the user wants to make a measurement of position and GPS time at his/her clock time t = 37 239.999 422 365 6 s of the week, and that this just happens to = be GPS time t 37240 seconds exactly, for purposes of illustration.
The actual receiver ECEF coordinates at this instant, which we give here for compar-
ison with the navigation solution, will be

i? cos 35°

m 5 224 663.389

WGS84

(7)

E sin 35°

m 3 658 348.690

Such accuracy is not justified, but we are specifying positions witiiin 1 nun so that the
convergence of the algorithms can be checked.
We follow the step numbering from the prescription given above.
(1) Suppose we already have the measured SV clock times tsvi for the signals
received at the chosen reception time.
= The reception time is specified on the local clock, which is t 37 239.999 422 365 6 s. This just happens to be at GPS coordinate time i = 37 240.000 000 000 s exactly and corresponds to a pseudorange c X 0.075 ms for SV 1. The GPS time should be recovered during the
computations.

We (2)

apply the prescription according to Sect. 20.3.3.3.3.1 from the ICD-

GPS-200 describing the user algorithm. We then obtain the system time t^,

GPS the

time of transmission, for each satellite i. The relativistic eccentricity

correction is part of this calculation. (The subscript i is not used in the ICD,

but is added here for clarity.) The data we actually use in the present example

were generated from a computer simulation which propagated orbits and clocks

for the entire GPS constellation forward in time. In this simulation the offsets

^isvi were randomly generated.

(3) Using the broadcast ephemeris, we obtain the ECEF coordinates of the
satellites from the prescriptions given in Table 20-IV of the ICD-GPS-200. From
the GPS time ti, the time interval tk = ti — toe from ephemeris reference epoch
is calculated. Then tk may then be used in the algorithm for computation of
ephemerides. This gives the x, y, z coordinates of the sateUite, in the ECEF frame,
at its transmission epoch. Table 1 gives the results after these steps. (These data
are actually generated from a computer simulation of the GPS constellation.)

12

Table 1. GPS satellite positions in ECEF coordinates.

SV#
1
2 3 4

Transmission epoch ti {s)
37 239.924 422 365 6 37 239.920 713 391 8 37 239.925 307 870 37 239.929 346 353 9

X^ (m)
13 005 878.255 20 451 225.952 20 983 704.633 13 798 849.321

y^ (m)
18 996 947.213 16 359 086.310 15 906 974.416
-8 706 113.822

Zj (m)
13 246 718.721
-4 436 309.875
3 486 595.546 20 959 777.407

In Table 1, GPS system times ti rather than the time intervals

tk

—

— ti

toe

from

ephemeris reference epoch,

are

used

to label

the

events. The subscript i varying from 1 to 4 labels data from the different

satellites. Table 1 assumes that the eccentricity correction discussed in
Appendix A has already been applied to obtain the transmission epoch

data given in the second column of the table. The position coordinates

ECEF of each satellite are given in the

frame, which is turning in inertial

space.

(4) To rotate from ECEF to ECIF coordinates, the inverse of the rotation matrix of eq (4) is required. Thus we need {R^{t^ — ic))~\ which must be
calculated and applied to the coordinates of each satellite individually. That is,

ti will be different for each satellite, but tc is the same for all satellites. For

this example, we choose an ECIF by fixing the rotation of the ECEF frame at

We — tc = 37 239.000 000 000 s.

stress that this choice is arbitrary better choices

are usually available. The purpose of this choice is to illustrate its arbitrariness

and also to illustrate how rapidly the iteration algorithm converges. With this

choice, for example, the rotation matrix for satellite 3 is

cosCte{t3 -tc) -sinpeitd -tc)

{R^{t3-tc)y

sintteitz -tc) cos Clf,{t3 -tc)

(8)

1

and fteits - tc) = 0.925 307 870 fte = 6.747 451 534 x 10"^ Corresponding

transformations are applied to each of the transmitted data sets.

(5) We transform the coordinates in Table 1 into this ECIF. The rotation

We C matrix is different for each transmission epoch.

use a subscript to indicate

that the inertial system has been arbitrarily chosen. The results from transform-

ing to the ECI frame are given in Table 2.

(6) The propagation delay equations to be solved arc discussed in Api)cndi,\

B. We give a quick summary here. The notation is cis follows. The four GPS

GPS satellites, at

times ti, t2, ^3, t^ send out signals from the ECIF locations

13

Table 2. Transmitted data transformed to the chosen ECI system

SV #
1
2 3 4

Transmission Epoch t^ (s)
37 239.924 422 365 6 37 239.920 713 391 8 37 239.925 307 870 37 239.929 346 353 9

Xa (m)

ya (m)

ZCi (ni)

13 004 597.642 20 450 127.566 20 982 631.270 13 799 439.294

18 997 823.895 16 360 459.358 15 908 390.245
-8 705 178.668

13 246 718.721
-4 436 309.875
3 486 595.546 20 959 777.407

ri5

1*2 5

1*3 5

r4 5

respectively.

(The subscripts now label data from the different

satellites, after transforming to the chosen ECI reference frame. The rotations do

not affect the times.) These four signals are received simultaneously at time i by a

receiver at position re- The problem is to determine t and r^ at the receiver. The velocity of the receiver does not have to be considered. The receiver's position

re at the time of reception t will be determined by the solution of four equations

which express the condition that the speed of propagation is c. The four equations

to be solved are

(re - Tjf - c\t - tjf =0; 3 = 1,2,3,4.

(9)

However, if the position yc and time t are known approximately, the equations can
be reduced to a system of four linear equations which can be solved by standard matrix inversion techniques. These equations are derived in Appendix B, and are

_ Ar - - = - - - (r(^)

r^.) .

c^{t^'^ tj)At

- (^c^{t'^'^

tjf

(r^''^

r,)^) .

(10)

In eq (10), the quantity c{t^^^ —tj) is the ith estimate of the pseudorange from the
(i)
receiver to the jth satellite, and r^ is the zth estimate of the receiver position in
ECI coordinates. Equation (10) is a system of Hnear inhomogeneous equations in the corrections Ar, At, which can be solved by matrix inversion. The matrix of coefficients of the unknowns Ar and At will usually be nonsingular, unless the
configuration of satellites is so unfavorable that eqs (9) do not have a solution
(such a condition is possible.) The solutions obtained will be approximate, but can be used to obtain new trial values; the iteration can be repeated as many
times as necessary to obtain the accuracy required.
To illustrate this process in the present case, as our initial guess at receiver position we take a worst case and assume the receiver is at the center of the earth. Also, we do not know the receiver clock bias. In this example we assume the clock bias is such that the reception event
14

occurs at t = ti +0.07500 s. Then Table 3 gives the results at each stage
of the iteration.
m The calculation converges to within 1 after three iterations. The
time at reception on the receiver clock, determined by the solution, is
tR = 37 240.000 000 000 Os

(7) These results must lastly be transformed into the WGS84 system by applying the rotation R^itu — tc)- That is, to obtain WGS84 coordinates at
the instant of reception, we must use a sidereal rotation corresponding to that instant. The necessary rotation matrix is

COSfteitR -tc) sniVt^{tR-tc)

R^tR-tc) = - sin ft,,{tR - tc] COS^e{tR-tc)

(11)

1

and Cte{tR — tc) — 7.292 115 146 7 x 10~^ because in this case, tR — tc happens to be exactly 1 s. Upon applying this rotation to the position coordinates given
in the last line of Table 3, the measured values of the receiver position are

m 5 224 663.388

WGSBA

m • 0.000

(12)

m 3 658 348.689

These agree with the true position coordinates to within a inillimeter.

Table 3. Results of iterative solution of propagation delay equations.

# Trial
(start) 1 2
3 4
5

User's X position
m m 5 057 363.392 m 5 226 931.552
m 5 224 663.780 m 5 224 663.374
m 5 224 663.374

User's y position
m m 2 355.126 m 354.224 m 380.983 380.988 m 380.988 m

User's z position
m m 3 541 092.792 3 659 938.391 m
3 658 348.973 ui 3 658 348.689 ni
3 658 348.689 m

User's clock time
37 239.999 422 365 6 s 37 239.997 532 305 8 s 37 240.000 033 455 9 s 37 240.000 000 006 s 37 240.000 000 000 s 37 240.000 000 000 (1 s

15

Table 4. Iterative solution ignoring earth rotation.

# Trial

User's x position (m) User's y position (ni) User's z position (m) User's clock time (s)

(staxt)
1
1 2
3 4 5

5 057 359.485 5 226 927.086 5 224 659.325 5 224 658.919 5 224 658.919

2 012.579 00.194 27.107 27.112 27.112

m 3 541 090.386 m 3 659 935.704 m 3 658 346.291 m 3 658 346.008 m 3 658 346.008

37 239.999 422 365 6 37 239.997 532 323 37 240.000 033 469 3 37 240.000 000 019 6 37 240.000 000 013 6 37 240.000 000 013 6

Example 3.211 Erroneous Use of ECEF Coordinates

One might think that effects due to earth rotation, during the very small signal prop-

agation times from transmitters to receiver, were sufficiently small that transforming to

ECI coordinates was not necessary. Let us therefore test this with the example data given

We in Table 1.

again initialize the iteration as in Example 3.21, above. The results of the

iteration, using the same algorithm as before are given in Table 4.

In this case, the solutions are in error by a significant amount. The receiver time is

off by 14 ns and the position is in error by almost 30 m. Thus, it is an error to ignore

earth rotation during the propagation of signals from transmitters to receiver.

However, these results suggest an approach which could result in a computationally ef-

ficient strategy for the receiver. Suppose that the receiver is continually monitoring signals

from four or more satellites and keeping track of its approximate position by neglecting

earth rotation, as in Example 3.211 above, but that an accurate position is not needed

until some particular moment. At that moment, the erroneous position and time could be

used as a starting point for a correct iterative solution of the one-way light time equations.

Since the starting point is not far off, this iteration should converge very rapidly. Thus,
transformations to the chosen inertial frame C would be made only when an accurate

position and time were needed, eliminating much computation at other times.

Example 3.212 Iteration with a Succession of Inertial Frames Another approach is to consider a succession of inertial frames in the iteration process,
each of which is aligned with the WGS84 frame at the assumed reception time. This would
have the nice advantage that at the end of the iteration, the ECIF would be exactly ahgned
with the WGS84 frame, so no final rotation would be necessary. The disadvantage is that
the ECIF coordinates of the transmission events must be recomputed at each stage of the iteration. Thus in the intermediate stages of iteration more computation is required. During actual use, we might have very good estimates of the receiver position available
to start the iteration. Thus, very few iterations could yield receiver position to sufficient
16

accuracy, making this approach favorable. To illustrate this approach, we shall use a succession of equations of the form

Xc = R^{tc){R^r'^WGSS4 - {R^{t - tc)r'XwGS84

(13)

We to obtain inertial coordinates of the transmission events.

will first make a reasoned

guess for the receiver time tn and then set tc — iR- One iteration will yield a more accurate

We estimate for Ir.

will then choose tc to be this new value and recompute satellite

transmission event data for this new tc-, and repeat the process as often as necessary. The

results are given in Table 5.

Example 3.213 The Benefit of Better Initialization It pays to be more careful with our initial guesses. To illustrate this, let us take as
our initial time the latest transmission epoch plus 75 ms and take the initial position to be
the subsatellite point of the latest transmission epoch. Again using a succession of inertial
frames as in the above example, the results are given in Table 6. What has been achieved
is more rapid convergence, compared with the previous example.

Table 5. Iterative solution with a succession of inertial frames.

# Trial

User's x position (m) User's y position (m) User's z position (m) User's clock time tc (s)

(start) 1 2
3 4

5 057 363.550 5 226 931.564 5 224 663.794 5 224 663.388

1 986.551 -25.990
-0.018
0.000

3 541 092.792 3 659 938.391 3 658 348.973 3 658 348.689

37 239.999 422 365 6 37 239.997 532 305 8 37 240.000 033 455 9 37 240.000 000 006 37 240.000 000 000

Table 6. Effect of better initialization.

# Trial

Uiser's X position (m) Uiser's y position (m) Uiser's z position (m) User's clock time tc (s)

(start) 1 2
3 4

3 313 469.280 5 190 759.859 5 224 663.389 5 224 663.389 5 224 663.389

-2 090 568.570
400.888 -0.892
-0.000
0.000

5 032 997.819 3 634 586.629 3 658 412.593 3 658 348.689 3 658 348.689

37 240.004 346 353 9 37 239.999 499 918 6 37 240.000 000 000 37 240.000 000 000 37 240.000 000 000

17

4. TIME TAGGING AT THE TRANSMITTER

Let us now consider time tagging at the transmitters with four sateUites. This ap-
proach requires a prior estimate of receiver velocity and perhaps acceleration over the interval of signal reception. If all signals leave satellites simultaneously with respect to the
ECEF frame, most users will receive these signals within 20 ms of each other. For most
users acceleration will contribute negligibly.
We address two separate issues for this type of receiver. The most important from
a conceptual point of view is that, during the 20 ms or so which may be required for the receiver to make four time-of-arrival measurements, the receiver will generally move
through the inertial frame. This motion must be accounted for in order to achieve accurate
position determination.
A second, more practical comphcation, is that the C/A code repeats every millisecond, and it might be easiest to choose the GPS system time tag tc for the measurements to be identical to one of the instants at which the C/A code begins to repeat. We call this a "code boundary." In the present example we shall illustrate use of C/A code in which the receiver only measures the arrival times of code boundaries on its local clock. We shall first discuss C/A code boundary measurements and the necessary time measurement
corrections. Then we shall discuss how to account for receiver motion.

4.1 Code Boundary Measurements
First we discuss the choice of the time tc at which the ECEF frame is to be "frozen" in place to create the ECIF. The time tc denotes a GPS system time. Unfortunately,

because of the clock and relativistic corrections represented by the quantity Atsv in the
ICD-GPS-200, the GPS system time t at a satellite differs from the satellite clock time,

tsv according to the equation

t^tsv-Atsv-

(14)

Thus, transmission of a a particular system time tc would occur at satellite clock time
which is numerically Atsv later (or earlier if Atsv is negative) by

tsv ^tc + Atsv.

(15)

The correction Atsvj foi' a particular satellite is contained in the navigation message and can be computed by the receiver. It changes relatively slowly during the 20 ms or less required by the receiver to make four time-of-arrival measurements, so we shall assume that for a particular set of four time-of-arrival measurements, Atsvj are constants.
To assist in making clear the distinctions among the various times which occur in the
R following discussion, we shall use a subscript for times measured on the receiver clock, and SV for times measured on the transmitter clock. Both of these may differ from GPS
system times for events, due to noise, drifts, offsets, and relativistic effects.

18

tp on SV clock

GPS system time

tQ = chosen GPS system time
Figure 1. Illustration of SV and receiver clock offsets.

4.2 Phase Time

The "phase time" r is transmitted from the satellite at SV clock time tsv and propa-

gates to the receiver where it arrives at receiver clock time Ir. In a particular measurement,

r is some number. If this number is an integral multiple of 1 ms so that it lies on a code

boundary, then it may be particularly convenient for the receiver to implement the mea-

We surement of the arrival time.

shall assume in the following that r lies on a code

boundary. Thus for the channel corresponding to the jth satelhte (with j = 1, 2, 3, 4), the

sequence of events may be described as follows.

(a) A common system time tc is chosen for the measurements, which are to
be time tagged at the transmitter at this system time tc- This system time tc
is some integral multiple of 1 ms and thus there is a code boundary in the phase
time which is numerically equal, t = tc-

(b) While the phase code is numerically equal to tc and is transmitted when
the SV clock time tsv3 = t = tc-, this does not usually occur at the GPS system time t = tc- Because of clock and relativistic corrections, transmission of the true GPS system time tc occurs ^tsvj later.
(c) The wave train propagates to the receiver where the code boundary r =
tc is detected and the receiver clock time tjij of this event is recorded.

This sequence of events is illustrated in Figure 1.

The arrival of a point in the code corresponding to tru(> GPS system tinu^ ic ^vill
occur an amount A^;?j later. Because of the relative motion of the riHcMVtn and tiansniitlcr

19

through the ECI, Atjij will in general be different from Atsvj due to the first-order Doppler effect. The receiver needs to estimate this time difference in order to find the arrival time, on the local clock, of the point in the code corresponding to the true time tag tc-
The small interval Atjij can be very accurately estimated by using the Doppler formula, eq (D.6) in Appendix D, which relates the coordinate time interval dtn between signal arrivals at the receiver, and the corresponding time interval dtsv at the transmitter

/ N-v\

/ N-V\

dtRil

= dtsv 1

,

{D.6)

V where and v are the transmitter and receiver velocities through the ECI, respectively,
at the instants of transmission and reception. This gives the relation between the time interval separating transmission of two nearly equal phase times, and the time interval separating reception of the two phase times. The time intervals dtsvj are less than 1 ms. The first-order Doppler corrections are typically of order 10~^. Therefore we need good estimates of the velocities of receiver and transmitter, so that the Doppler corrections are known within 10~^ in order to calculate dtn to better than 1 ns.
Because the Doppler corrections in eq (D.6) are small, we can use eq (D.6) for small time intervals without making appreciable error. Thus, to sufficient accuracy for GPS, eq
(D.6) can be rewritten as

Atn,=Atsvj^r^-^^.

(16)

N The receiver needs estimates of the receiver velocity v and the unit vector pointing
from transmitter to receiver in the ECI. These can be obtained from estimates of position and velocity which have been built up over time during the computation of successive
navigation solutions. The transmitter velocity Vj can be computed from the prescriptions
given in Table 20-IV in the ICD-GPS-200, by taking derivatives of the satellite position with respect to coordinate time. Since all these vectors must be represented in the chosen
ECI, some rotations may be necessary. It is useful to keep in mind, however, that dot products may be computed in any convenient ECI frame.
Thus, in brief, if the phase time corresponding to true GPS system time tc is later than tsvj — tc by Atsvj, then the time of reception of the phase time corresponding to
tc at the receiver is later by Atjij given by eq (16) above. These time relationships are illustrated in Figure 1.

4.3 Accounting for Receiver Motion In the following discussion we will assume that the code boundary measurements have
been made and corrected as described above. After making the necessary corrections,
denote the values of the arrival times of the code corresponding to the true GPS time
20

tc on the receiver clock by tcj = tRj + ^tjij. These will still not be GPS system times
because in general there will be a bias on the local clock. Denote the bias by be, and
denote the GPS system times of arrival by tj; then the relation

tj=tcj+bc

(17)

R holds. Note that there is no subscript on the time on the left side of eq (17) because

this symbol now refers to the GPS system time of an event.

Further, because of the receiver's motion through the ECI, the receiver will not be at

a fixed position during the measurements, which occur from 67 to 86 ms after tc if the

receiver is near the earth. This means it is necessary to have available to the receiver an

estimate of receiver velocity in the chosen ECI frame, so that the receiver positions at the

various reception times can be calculated. Let the receiver position and velocity in the ECI

vc at the system time tc be re and

respectively. We assume vc is known to sufficient

accuracy, but re and be are to be determined by the measurement. Thus the positions of

the receiver at the measurement times are given by

rRj = TRitj) = rc + ^citcj +bc- tc)

(18)

The sateUite positions at the instant of transmission of the signals at the common GPS
system time tc wiU be denoted by r^. These positions are obtained from the navigation
messages and are specified in the WGS-84 reference frame. This means that all of these position vectors are specified in the same reference frame. Therefore it is natural when setting up the propagation delay equations to choose an ECI reference frame which coin-
cides with the ECEF frame at the instant tc- Thus no rotations are necessary in order to express these SV positions in a common inertial frame. If we solve for re, this will be in the ECIF (and the ECEF frame) at tc, so no rotation will be required for that vector. However, if we want trj in the ECEF frame, rotations will be necessary.
The propagation delay equations are

i^Rj -^jf = c^(^j -tcf

or

(re + ^citcj +bc- tc) - Fjf - c^itcj +bc- tc)^ . j = 1.2.3.4.

(19)

These four equations can be solved iteratively for the unknown bias be and the tlu-ee
position components re of the receiver in the ECEF frame at the chosen system time tcOne iterative method of solution is described in Appendix C. We briefly summai-ize that
process here.
We write the position and bias in terms of estimates, plus small corrections:

re = r^c + ^r- ^c = i^c + ^^

(20)

21

Substituting these expressions into eq (19) and expanding, keeping only linear terms in the
small quantities Ar and At, leads to a linear system of four equations in four unknowns:

R^;^ -Ar-c^i(T7^^)-v..Ry))A.4 c^(lj^O'-K)

(21)

where

T^/^=tcj+b)^^-tc,

(22)

and

Rf =r^)+vcTf)-r,.

(23)

This system of equations is solved by iteration in the manner described in Appendix C.

— 4.4 Time Tagging at the Transmitter General Prescription:
We now describe the specific procedure in parallel with the steps described above for
time tagging at the receiver.

(1) Measure the times tjij on the local clock at the reception of a given code
epoch r, corresponding to a single GPS phase time r in each code train, from

each of four satellites.
It may be convenient for a receiver to measure the time on its local

clock for the arrival of a particular point of the received code train,
such as one of the code boundaries. A code boundary corresponds to a

particular millisecond of the week. This millisecond is our chosen time

tag. However, each SV clock, and hence the clock's transmitted code,

is generally offset from from GPS time, and the offset will be different

We for each clock. The differences will be less than 1 ms.

will need to

account for this in step 2.

(2) Apply the prescription for corrections as given in the "User Algorithm for SV
Clock Correction" of the ICD-200, section 20.3.3.3.3.1, to obtain Atsvj For a given satellite, Atsv is the difference between the instant the
code for phase time r was transmitted and the value of GPS time r.
(See eqs (14-15).) Atsv includes both the relativistic correction Atr and the broadcast second-order polynomial correction for the SV clock
offset from GPS time.

(3) Compute each satellite's position in the ECEF frame at the time of transmission of GPS time t.
This is a straightforward application of the equations of Table 20-
GPS IV, for determining the z .x, t/, coordinates of the satellite at the

time tc-

22

(4) Choose an ECIF frame for computation of the one-way hght times.
In this case the natural choice is to freeze the ECEF frame at the GPS transmission time tc(5) If the ECI is defined by freezing the ECEF frame at GPS transmission time
tc, the coordinates of each sateUite will not need to be rotated. If we estimate receiver position at tc we do not need to rotate the receiver position coordinates
either. We need to adjust the local clock measurements, however, for the time
interval Atsv and the motion of the receiver over the time Atsv- This is because
the code boundary numerically corresponding to tc is not necessarily transmitted
at GPS time tcWe may use a deterministic estimate of the change in receiver po-
sition over the reception interval. Since this model is usually a velocity, we will refer to it that way here for simplicity. This velocity estimate remains fixed during the solution of user position and time. It is neces-
sary to obtain the velocity estimate in the coordinates of the chosen ECI.
Then we may estimate receiver coordinates in the ECEF frame which coincides at tc with the ECI. We then use the velocity to estimate receiver positions at the other reception times. We assume the user clock offset from GPS time is constant over the interval of reception.
First apply the corrections represented in eq (16) to compute Atfij. Then compute the times of arrival on the local clock tcj, of the signal
transmitted at GPS time r. Iftjij is the local clock time for the received
epoch from satellite j, then

tcj = tRj + AtRj

(24)

is the local clock time of reception of the signal transmitted at GPS time
t from satellite j. The correction AtRj is calculated from eq (16) using estimates of vrj, Tj, and Vj, and Vj in the ECI which have been built up from a filtering algorithm or from prior solutions.
The time relationships are illustrated in Figure 1 . For a given satellite vehicle SV, Atsv is the diff'erence between the instant the code for
GPS time r was transmitted and the value of GPS time r. Atsv includes
both the relativistic correction Atr and the broadcast second-order poly-
nomial correction for the SV clock offset from GPS time. Atjij is related
to Atsv by first-order Doppler corrections given in eq (16).

(6) Solve the equations for the geometric path delays to find the receiver's position

and time offset at tc-

Again we linearize the propagation delay equations and solve Ihcm

We iteratively.

may iterate to find the position at the first reception

23

.

time, using the velocity estimate to extend to the other reception times. This iteration is similar to step (6) in the case of time tagging at the
receiver.

(7) If we have solved for the user position at tc there is no final rotation. In
this case the ECIF and ECEF frame coincide at tc- If we want the ECEF frame
position at a different reception time, a rotation will be required.

It would be possible to update the velocity estimate using positions obtained from

two transmission time tags. This leads into other concerns and options associated with

using GPS: filtering estimates over time. Techniques for estimating position, velocity, and

We acceleration could be coupled with strategies for filtering these estimates over time.

We will not discuss these options here.

will mention, however, that if there are infrequent

measurements of GPS signals from individual satellites, the local clock may be used to

flywheel between them and find a solution. This may require careful filtering algorithms

for estimation of the local clock frequency. In particular, relativistic effects on the local

clock frequency due to velocity or gravitational potential may have to be considered.

4.5 EXAMPLES

Example 4.51 Time Tagging at Transmitter
We give here an example of the analysis of data when the measurements are time
tagged at the transmitter. This means that the receiver measures the arrival times of at
least four signals which were transmitted simultaneously, at GPS system time tc, from four satellites which are labelled with index j = 1, 2, 3, 4. This method of time tagging has
the great advantage that the satellite positions at the instant of transmission, as specified
in the navigation message, are all given in the same coordinate system. Thus no rotations
are necessary to express these positions in a common inertial frame. For this example we
choose a code boundary exactly corresponding to
tc = 240.000 000 000 s
The following steps are in the order of the seven steps listed in Sections 2 and 3 above.
(1) The times tfij on the local clock of reception of the code epoch r = 240
sec exactly are recorded. These are given in the first column of Table 7. (2) The satellite clock corrections Atsv, including the relativity corrections,
are obtained from the navigation messages; these are given in the second column
of Table 7. In this simulation these offsets are randomly generated within a 1 ms
range using a uniform distribution.
(3) Compute each satellite's position in the ECEF frame at the time of transmission of GPS system time r.

24

This is a straightforward appHcation of the equations of Table 20-IV in the ICD-200, for determining the x, y, z coordinates of the satelhtes at
the GPS time tc- The results are given in Table 7 for the four satelUtes.

(4) Choose an ECI frame for computation of the one-way Hght times.
In this case the natural choice is to freeze the ECEF frame at the GPS transmission time tc-

(5) Compute the received time offsets At/^j and the times tcj of arrival on
the local clock, of the signals transmitted at GPS system time tc = t. We must make the first-order Doppler corrections, then compute

the time intervals ^tnj. This requires estimating the Doppler correc-

tion. However since the Doppler corrections are very small, in a first

approximation we might neglect the Doppler correction factor entirely
and assume Ai/?j ^ Atsvj- We might then use the resultant positions
to estimate Nj to within a few percent. If the estimate of position in the
chosen ECIF is within 890 km of the true position and the receiver is
near the earth, then the error made in estimating the direction cosines of the line-of-sight unit vector Nj will then be less than 10""^. In a

second approximation, estimates of the Doppler correction factor which
use velocities and a value of Nj accurate within a few percent should be

adequate.

Estimating receiver velocity. In the present example we consider
a receiver which has an estimated ground velocity of Vg = 300 m/s in

a direction 30° north of east, parallel to earth's surface. Suppose the

receiver at the instant tc — 240 s exactly is at the position

i? COS 35°

m 5 224 663.389

X,iyG584 =

—

i? sin 35°

m 3 658 348.690

(25)

which has been used in our earlier examples. The net velocity will be

a linear combination of the ground velocity and velocity of the point

XivG584 due to earth rotation. From geometric considerations, the net

receiver velocity at time tc in this example will be

-Vy sin 30° sin 35°

vc

t;^ cos 30° + nei?cos35°

Vg sin 30° cos 35°

(26)

If the position and ground velocity were accurately known, the net re-

ceiver velocity in the ECIF would then be

-86.036 465 453"

Vc = 640.796 091481 m/s.

(2-

122.872 806 643

25

Therefore in the present example, we take the velocity estimate in the

ECIF to be

"-86"

vc = 641 m/s.

123 J

We Estimating receiver position.

illustrate the calculation of the

first-order Doppler correction with satellite 1 in Table 7. The unit vector
N is computed from the satellite and receiver positions, as in eq (D.3):

(28)

N= ' \rR{tR) -ftI^t)!

{D.3)

where rR{tji) and ftI^t) are the ECIF position vectors of receiver and
transmitter at reception and transmission times tR and tx- At the receive time, with bias of correction 423 ns, the receiver position vector is
actually

rRi = rc+vc(240.069 124 140 3 - 0.000 000 423

5 224 656.660

50.117

m

3 658 358.300

- tc)

(29)

and we therefore take as an estimate of this position the vector

6 000 000

ri?i

23 000 m,

3 100 000

(30)

and the receiver velocity is given in eq (28), which was discussed above.
The bias correction was not significant in arriving at this rough estimate
of receiver position.
Estimating satellite velocity. The satellite position at transmit time tc is given in the first row of Table 7. The satellite velocity is also needed. It may be obtained by differentiating the satellite positions, which are derivable from the prescription of Table 20-IV in the ICD-GPS-200, with respect to time, and then applying a rotation into the ECIF. For the present example the satellite velocity in the ECIF is

Vi =

-304.267 760 3 853.051 168 -257.220 637

m/s

(31)

26

The position estimates result in the following calculated values of the components of Ni:

-0.437 588

Ni

-0.114 125

-0.891 904

while the "true" value of Ni obtained from the simulation is

(32)

-0.477 846 083

Ni

-0.115 974 357

-0.870 754 311

(33)

For the present case, Atsv = -0.000 056 554 783 4 s. Then eq (16)
gives the estimate

Atm ~ -0.000 056 554 784 s

This estimate agrees extremely well with the "true" value of Atjn obtained from the simulation. For details of this comparison, see Appendix
E.
Thus the Doppler correction factor given in eq (16) can be computed by the receiver and added to tjij to predict, with excellent accuracy, the
arrival time tjij + dtjij of the code which corresponds to transmission
of tc from the satellite. Thus only times of arrival of code boundaries would actually be measured by the receiver.
Three additional calculations using eq (16) are given in the third column of Table 7 for the satellites used in the present example. Here
the time tc equals 240.000 000 000 s exactly. The first column of Table 7
= gives local clock arrival times of the phase time r tc- The second column lists the values of the SV clock and relativistic corrections derived
from the navigation messages. The third column lists the corrections AtRj obtained from eq (16). Note that the largest difference between
Atsvj and Atnj is little more than 0.5 ns. The maximum possible such
difference is about 4 ns and could occur for a clock in a satellite in low earth orbit. The resulting values of local clock arrival times tcj are listed in the last column of Table 7 and are repeated in the last column
of Table 8. There may still be a local clock bias in these data. In this
example, the receiver clock has a bias of —423 ns.
27

Table 7. Computation of first-order Doppler corrections.
= (units s)

ID #

tRj

AtsVj

^tRj eq (16)

tcj = tRj + AtRj

1

240.069 180 695 048 -0.000 056 554 796 2 -0.000 056 554 784 240.069 124 140 263

2

240.078 615 674 081

0.000 493 792 700 6

0.000 493 792 194 240.079 109 466 275

3

240.075 360 518 523

0.000 259 725 929 9

0.000 259 725 401 240.075 620 243 923

4

240.072 939 857 576

0.000 305 714 376 7

0.000 305 713 988 240.073 245 571 564

Table 8. Satellite positions r^ and mecisiured arrival times tcj-

ID #
I
2 3 4

Satellite Xj (in)
15 126 951.488 20 996 573.361 23 391 384.147 10 656 983.742

Satellite IJj (m)
2 403 354.115
-15 349 939.783
10 956 612.394
-11 270 652.822

Satellite Zj (m)
21 702 797.709
-5 178 329.880 —4 333 140.777
21 703 578.299

Local clock time tcj (s)
240.069 124 140 3 240.079 109 466 3 240.075 620 243 9 240.073 245 571 6

Table 8 gives the data from the all four measurements: the derived
positions of the satellites at GPS system transmission times and the

corresponding arrival times tcj on the receiver clock.

(6) The iterative solution, which has been carried out as described in Ap-

pendix C is given at successive stages in Table 9. At the initial stage, estimates

ECEF for the

receiver position are inserted. For the sake of illustration, we as-

sume the initial estimate of the receiver position places it at earth's center. The value of the local clock bias is initially estimated to be 300 ns. The iteration pro-
cess converges rapidly, to within a millimeter or so of the correct ECEF receiver

position.

(7) No rotations are necessary in the final stage since the position has been
determined in the ECEF frame.

The receiver positions at the times tcj of reception might be needed for other purposes. These could be obtained using a best estimate of the velocity of the receiver in the ECIF
(the ECEF frame at tc) using eq (18). The components of re are in this reference frame. A second series of measurements at a slightly later time t^ will give a new value of tq in
the new reference frame. These two values could be differenced or filtered to provide an update to the ground velocity. Such an updated value could be transformed into an ECIF
28

Table 9. Iteration solution for re and receiver clock bias be-

# Trial
1
2 3 4 5

Receiver Xq (m)

Receiver yc (m)

Receiver Zc (m)

Local clock bicus be (s)

0.0
5 051 460.568 5 226 991.162 5 224 663.799 5 224 663.390 5 224 663.390

0.0
-9 835.265
132.182
-0.023 -0.000 -0.000

0.0 3 545 131.610
3 659 870.280 3 658 349.573 3 658 348.690 3 658 348.690

0.000 000 300 -0.002 453 425 2
0.000 032 544 3 -0.000 000 417 2 -0.000 000 423 -0.000 000 423

with the help of a rotation such as in eq (28), to be used in further position measurements.
Example 4.52 Time Tagging at Transmitter Neglecting Doppler Corrections The first-order Doppler corrections to Atsv are very small, as seen in Example 4.51,
above. Here we repeat Example 4.51 with one change: the Doppler corrections are completely neglected. Since including these first-order Doppler corrections usually amount to
changing the arrival times on the local clock by less than 1 ns, only a few centimeters of
error would be expected to occur when such corrections are neglected. Table 10 gives the satellite positions Vj and the measured arrival times tcj when this approximation is made. The table is identical to Table 8 except in the fourth column.
Solutions of Eqs. (21) by iteration are given in Table 11. Thus neglecting the first-
order Doppler correction causes a position error of about 75 cm in this case.

Table 10. Satellite positions Vj and arrival times tcj neglecting Doppler correction.

ID #
1
2 3 4

satellite Xj (m)
15 126 951.488 20 996 573.361 23 391 384.147 10 656 983.742

satellite yj (m)
2 403 354.115
-15 349 939.783
10 956 612.394
— 11 270 652.822

satellite Zj (m)
21 702 797.709
-5 178 329.880 -4 333 140.777
21 703 578.299

~ tcj tjij -\- /\tsvj (s)
240.069 124 140 3 240.079 109 466 8 240.075 620 244 5 240.073 245 572

5. TIME TRANSFER RECEIVER
In this case we may use each pseudorange measurement from each satellite sepai^ately
to estimate our clock offset from GPS time. We either time tag measurements at a trans-
mission time and measure time of reception, or time tag measurements at a reception time
and determine time of transmission from the code lock of the receiver. We then use the
29

Table 11. Solution for re and receiver clock bias be neglecting Doppler correction.

# Trial
1
3 4 5 6

Receiver Xc (m) Receiver yc (m) Receiver Zq (m) Local clock bias be (s)

0.0
5 051 461.207 5 226 991.866 5 224 664.501 5 224 664.092 5 224 664.092

0.0
-9 835.359
132.086 -0.072
-0.096 -0.096

0.0
3 545 131.862 3 659 870.566 3 658 349.243 3 658 348.976 3 658 348.976

0.000 000 300 -0.002 453 427 6
0.000 032 542 4 -0.000 000 419 1 -0.000 000 424 9 -0.000 000 424.9

prescriptions from the ICD-200, Table 20-IV, to obtain the satellite position in the ECEF frame. We must estimate the geometric path delay t^ in addition to the ionospheric
and tropospheric corrections. Since we know the receiver position we may compute t£) as
follows.
If ECI position vectors are referenced to the time of signal transmission, then

\v{tn)-Ii{tT)\ ^Jr{tT)-mT)\ [r{tT)-n{tT)]-v

._

,

c

c

c^

where

R is the SV position vector,

r is the receiver position vector,
tT and tn are GPS time at respectively transmit and receive times, and
V is the receiver velocity vector in the ECIF.

The last term in eq (34) is sometimes called the Sagnac correction. In the case of an

earth-fixed user,

= V i7e X r.

(35)

If ECIF position vectors are referenced to the time of signal reception (instead of
satellite transmission time), then

to =tR — tT = \v{tR)-K{tT)\ «^^ \r{tR)-K{tn)\

c

c

V - [r(t;,) R(^h)] •
5
c^

.... (3b)

where:

V is the SV velocity vector in the ECI.

A single- channel receiver with antenna placed at a known location can provide a means for accurate transfer of GPS time to a local clock. Let the antenna be at the position
^ECEF whose earth-fixed coordinates are accurately known and do not change with time.

30

The appropriate expression for the total propagation delay for a signal from one satellite is eq (34). In eq (34) the positions are ECI positions but are referred to the time of transmission t^. The first term on the right side of eq (34) is proportional to the distance between receiver and transmitter at the transmit time tx', if we call this distance L, then

L = \r{tT)-Ti{tT)\.

(37)

Although the position vectors r{t'r) and R(tT) are ECI position vectors, the distance L is invariant under spatial rotation and can therefore be computed in any convenient
coordinate system. The ECEF position of the satellite at transmit time is directly derivable
from the prescription given in the ICD-200, Table 20-IV, and therefore no rotations are necessary to compute L:

L = IrECEFitr) —'R'ECEF{tT)\-

(38)

The second term in eq (34), the Sagnac correction term, can be further reduced for this case since it also is a scalar with respect to rotations. The ECI velocity of the earth-fixed receiver is due to earth rotation and is

= V 17e X r{tT).

(39)

The Sagnac correction term, which we denote by Atsagnaci is therefore

[r{tT)-R{tT)]-^e^r{tT) At Sagnac
C2

tie r{tT) X [r{tT)-K{tT)]

C2
^e r{tT)x [-R(iT)]

C2
He • K{tT) XX Tr{tT)

(40)

c2

^e • ^ECEPitT) X ^ECEFitr)

c2
l^el {XEC E F{iT)y EC EFitT) — yECEF{tT)x ECEFitr))
C2

where standard vector operations have been used to simplify the scalai* product. In the last step, the only component of the cross product which contributes is parallel to the earth's angular velocity vector; this product is invariant under rotations about this axis,
so the computation can be performed in the ECEF frame at the instant of transmission.
Thus for a timing receiver no rotations need be applied to compute the local time.

31

The final result in eq (40) above can be put in a form which is easy to interpret and can be found in many places in the literature:

"^^S AtSagnac =

,

(41)

where Ae — | [R(iT) x ^{tT)]z i^ the equatorial projection of the area swept out by a
vector from the center of the earth to the path of the timing pulse as it propagates from
transmitter at R(iT) to receiver at r{tT)- Clearly this area can be calculated in either
ECIF or ECEF coordinates.

Example 5.1 Earth-fixed Single- Channel Timing Receiver
We give four examples here based on the data from four different satellites given in Table 1. For all these examples the ECEF position of the receiver is that specified in eq
(7).
In Table 12 the first column labels different satellites corresponding to the data entries in Table 1. The second column gives the transmission epochs. The third column gives the transmission delay contributions X/c, the fourth column gives the Sagnac corrections, and the fifth column sums the entries in the previous three columns to give the time at the receiver. The agreement of the entries in the last column refiects consistency of the procedure across the satellite constellation.

Table 12. Examples of time transfer to a single-channel earth-fixed receiver.

# SV

Transmission epoch ij- (s)

L/C (s)

^^Sagnac (s)

tjij (s)

1

37 239.924 422 365 6

0.075 577 714 8 -0.000 000 080 5 37 240.000 000 000

2

37 239.920 713 391 8

0.079 286 677 5 -0.000 000 085 5 37 240.000 000 000

3

37 239.925 307 870

0.074 692 197 4 -0.000 000 067 4 37 240.000 000 000

4

37 239.929 346 353 9

0.070 653 609 2

0.000 000 036 9 37 240.000 000 000

6. REFERENCES
[1] "Navstar GPS Space Segment/Navigation User Interfaces," lCD-GPS-200, Revision C (10 October 1993), Code Went. No. 29562, ARINC Research Corporation, 4055 Hancock Street, San Diego, CA 92110.

[2] Ashby, N.; J.J. Spilker, Jr., "Introduction to Relativistic Effects in the Global Position-

ing System," Ch. 18 in Global Positioning System: Theory and Applications, Vol. /,

eds.,

B.

W.

Parkinson

and

J.

J.

Spilker,

Jr.,

Progress

in

Astronautics

and

Aeronautics ,

32

Vol. 163, AIAA, Inc., 370 L'Enfant Promenade, SW, Washington, DC 20024-2518
(1996).

Ashby, N. "Relativistic Effects in the Global Positioning System," Chinese J. Systems Eng. Elect. 6, 199-237 (1995).

Gibson, L. Ralph, "A Derivation of Relativistic Effects in Satellite Tracking," Naval

TR Surface Weapons Center

83-55, (April 1983).

Deines, S., "Missing Relativity Effects in GPS," paper presented at ION-GPS-90, Colorado Springs, CO, September 1990. See also "Missing Relativity Terms in the Current
GPS Algorithm," S. Deines, Contract F04701-87-C-0027, (June 1990), and "Missing
Relativity Terms in the GPS," Navigation 39, 111-131 (1992).
FUegel, H.; R. DiEsposti (1997), "GPS and Relativity: an Engineering Overview,"
Aerospace Report ATR-97(3389-l).

"Department of Defense World Geodetic System 1984, Its Definition and Relation-

DMA ships with Local Geodetic Systems,"

Technical Report 8350.2, Stock Number

DMATR83502WGS84, U.S. Geological Survey, Denver, CO, 80225.

Malys, S.; J. Slater, "Maintenance and Enhancement of the World Geodetic System 1984" (1994), Proc. ION-GPS-94, Salt Lake City, Sept. 20-23, 17-24 (1994).

"1992 lERS Annual Report" (July 1993), Central Bureau of lERS, Observatoire de
Paris, 61 avenue de I'Observatoire, F-75014 Paris, France.

10] "lERS Standards (1989)," D.D. McCarthy, ed., lERS Technical Note No. 3, Observa-
toire de Paris, 61 avenue de I'Observatoire, F-75014 Paris, France.

11] Moritz, H.; I.I. Meuller (1987), Earth Rotation, Theory and Observation, Ungar, New
York.
12] Hofmann-Wellenhof, B.; H. Lichtenegger, J. Collins (1993), GPS Theory and Practice, Second Edition, Springer- Verlag, New York.

13] Goldstein, H., Classical Mechanics, 2nd Ed. (1980), Addison-Wesley Pubhshing.
Reading, MA, Sect. 4-9, p. 174.
14] "Relativistic Effects in the Global Positioning System," D. Eardley et al., JASON,
MITRE Corporation, 1820 DoUey Madison Blvd, McLean, VA 22102, May 1985.

15] "Accuracy of Time Transfer in SateUite Systems," Committee on Accuracy of Time Transfer in Satelhte Systems, Air Force Studies Board, National Academy Press, Wash-
ington, D.C. (1986).

16] See eq (2) of "Lectures on Frequency Stabihty and Clocks," by R.F.C. Vessot. in

33

Experimental Gravitation, ed. B. Bertotti, Academic Press, Inc., 155 (1974). [17] Boucher, C; Z. Altamimi, M. Feissel, P. Sillard, "Results and Analysis of the ITRF94,"
lERS Technical Note 20, Central Bureau of lERS (March 1996), Observatoire de Paris, 61, avenue de I'Observatoire, F-75014 PARIS, Prance.
34

Appendix A.

The Eccentricity Correction

The epoch of the transmission event from each satelhte must be corrected to account

for variations in the proper time of each satellite clock due to orbital eccentricity. This
aspect of system design of the GPS was fixed very early when satellites did not have much
computing power. The correction compensates the atomic clock in the satellite for periodic gravitational frequency shifts and second-order Doppler shifts due to the satellites' periodic change of radius and velocity in its eccentric orbit. The correction converts the proper time
— kept on the atomic clock, to coordinate time that is, to GPS time. To put it another

way, a satellite in an eccentric orbit samples different values of earth's gravitational field as the satellite moves toward and away from the earth. Consequently it suffers a periodic
gravitational frequency shift. At the same time, conservation of angular momentum causes
the satellite to speed up or slow down, giving rise to a periodic second-order Doppler shift.
Thus SV clocks must be corrected for these effects so that the corrected clocks effectively

run at a constant rate, the rate of coordinate time. It is important to realize that this correction is the same for all receivers, no matter
where they are or what their state of motion is. The correction relates solely to the
transmitter. In particular it is correct for users of GPS anywhere in space. It is incorrect

to claim, as is done in Ref. [6], that the correction applies only to users near earth's
surface.
Let us denote the time that the user estimates from the information encoded in the signal as tsv a^ in the ICD-GPS-200. This must be corrected according to the section
describing the user algorithm to obtain GPS system time t (Sect. 20.3.3.3.3.1 in the ICD-
GPS-200C). The correction has several different, but equivalent forms. Using the notation of ICD-GPS-200C, Sect. 20.3.3.3.3.1, the relativistic correction to the transmission epoch

tsv is

Atr

=

-

^VGMA

2

^

es. mE^k

c^

= Fe^sinEk

2R V

(^-1)

where

F = -4.442 807 633 X 10"^° s//l^,

(^-2)

GM and where

is the product of the Newtonian gravitational constant and the mass of

A the earth,

is the semi-major axis of the satellite in meters, e is the eccentricity, Ek

R is the satellite's eccentric anomaly,

is the instantaneous position vector of the SV at

V transmission time, and is the instantaneous velocity vector of the SV. We need GPS

35

time t before we can obtain Ek as prescribed in Table 20- IV, so the equations are coupled.
It suffices to use tsv to determine E^, then obtain GPS time t. Also, since the scalar
product in the last form of the correction is invariant under spatial rotations, it can be evaluated in any inertial frame. This correction Atr is subtractedhom the transmitted time
tsv to obtain GPS system time t. The time t may then be applied to the user algorithm
for ephemeris (Table 20-IV in the ICD-GPS-200C), to obtain the x, y, z coordinates of each
satellite at transmission time.
Derivations of the above results can be found in references [2], [3], and [4].
36

Appendix B.

Linearization Algorithm for
Time Tagging at Receiver

We now explain one possible method of linearization of the one-way light time equa-
tions which are valid in an inertial coordinate system.
Consider four GPS satellites which at GPS times ti, ^2, ^3, ^4 send out signals from the ECIF locations ri, r2, r3, r^. We assume that these data have been expressed in the
chosen ECI reference frame. These four signals are received simultaneously at time t by a receiver at position r. The problem is to determine t and r at the receiver. The velocity of the receiver does not have to be considered. The receiver's position r at the time of
reception t will be determined by the solution of four equations which express the condition
that the speed of propagation is c. The four equations are

{r-rjf-c\t-t-,f=0; i = 1,2,3,4.

(B.l)

Such equations are clumsy to solve in general. However, if the position r and time t are known approximately, the equations can be reduced to a system of four linear equations which can be solved by standard matrix inversion techniques.
Suppose that for user position and time we write

r = r(*)+Ar; = + t t^'^ At,

{B.2)

where the quantities r^*^, t^'^^ represent guesses or trial values for the position and time, and Ar, At are assumed to be small quantities which are to be determined by solving eqs (B.l). The index i in eq (B.2) labels the zth trial value in a repeated or iterative process
of solving Eqs. (B.l).
We substitute eqs (B.2) into eqs (B.l), expand and neglect squares of the small
quantities Ar and At. Then the equations take the following form

= - (r(*) - r,) • Ar - ^{t^'^ - tj)At

- ((P{t^'^ - tjf

- (r^'' r,)^) ,

{B.3)

— for j 1, 2, 3, 4. In eqs (B.3), the quantity c(t^*^ ~^j) is the ith estimate of the pseudorange
from the receiver to the jfth satellite. Eqs (B.3) are a system of linear inhomogeneous equations which can be solved by matrix inversion. Numerical examples of the solution
are discussed in the main text. The matrix of coefficients of the unknowns Ar and At
will usually be nonsingular, unless the configuration of satellites is so unfavorable that eqs
(B.l) do not have a solution. The solutions obtained will be approximate, but can be used to obtain new trial values; the process of iteration can then be repeated.

37

Thus, suppose at stage i in the iteration, eqs (B.3) have been solved for corrections
Ar and At; then better approximations can be found by setting

^ = + r(*+i)

j.{r)
_^ ^j.;

i(^+i)

t(^)

At.

(BA)

The iteration converges rapidly, particularly if the initial trial guesses were close. The solutions will be found at some stage when the inhomogeneous terms on the right sides of eqs (B.3) become negligible. Then the solutions satisfy

(rW_r^.)2_c2(i(0_t^,)2^0; i = 1,2,3,4,

(R5)

which are equivalent to eqs (B.l).
In summary the computation scheme proceeds as follows: Step 1. Choose approximate values for user position and time. Step 2. Solve eqs (B.3) to obtain corrections Ar and At for user position
and time. Step 3. Test whether the corrections are sufficiently small that computa-
tion can be halted; if not, combine the corrections obtained in Step 2 with the
approximate values of position to obtain a new approximate position. Then go
back to Step 2.

38

—

Appendix C.

Linearization Algorithm for
Time Tagging at Transmitters

When the measurements are time tagged at the transmitters, then the measurements
at the receiver will occur at unequal times. This entails a different method for linearization

of the propagation delay equations than that discussed in Appendix B.
Suppose that we decide to time tag the signals from four GPS satellite transmitters at the common GPS system time tc- Let the arrival times of these signals on the receiver
clock be tjij, and let the positions of the sateUites at time tc from which these signals
were transmitted be r^, with j = 1, 2, 3, 4. These positions are derived from the navigation messages; they will be expressed in the ECEF frame at the instant tc- Therefore it is highly advantageous to work with an ECI reference frame which coincides with the ECEF

frame as though it were frozen at tc- Then no rotations will be necessary to bring the
— satellite positions into a common ECI frame they are already expressed in the ECIF
defined at tc- If we solve the propagation delay equations for the receiver position at tc
we have the receiver position in this same ECIF. Thus time tagging at the transmitters

.has the great advantage that no rotations are needed.

In general the signal propagation times from satellites to receiver will not be equal
they can vary roughly from 67 to 86 ms for a receiver on earth's surface. This means that

it will be necessary to have an estimate of the receiver velocity at the time tc in order that

the different receiver positions at the different signal arrival times can be accounted for.

We We may also need the receiver's acceleration.

assume these are available, to sufficient

accuracy, from previous position solutions or from a filtering algorithm.
We will illustrate this here by assuming the receiver position is expressible as a Taylor

series

r/?(t)-rc + vc(i-ic) + -,

(C.l)

where re and vc are the receiver position and velocity in the chosen ECIF, at time tc- An
additional acceleration term in this expansion might be needed for a receiver in a highly
dynamic trajectory. We assume here that vc is known but r^ is unknown.
The measured arrival times tjij of the signals on the receiver clock will in general not
be true GPS system times because of a possible bias be on the receiver clock. Thus the GPS system times of arrival will be

tj = tRj + be.

(C.2)

and as part of the data analysis, the bias will be found. The receiver positions at times tj

will therefore be

mitj) = TRj =rc + ycit'Rj +bc- tc)

{C.2)

39

In this Appendix we do not discuss code boundary measurements and corrections. This

We important topic is discussed in detail in Sect. 3.

simply assume that the appropriate

places in the code have been identified by the receiver.

The four propagation delay equations are then

(VRj -Tjf ^ C^itRj +bc -tcf

(C.4)

or

(re + y^citRj +bc- tc) - Yjf = c^{tRj + bc- tcf .

(C.5)

The unknown position vc and bias be may be found from these four equations cis follows. As in Appendix B, we assume we can write the unknown quantities in terms of good
estimates, plus small corrections. Thus we write

re

c + Ar, be = b^fi + Ai

(C.6)

where the quantities Yq\ b^ represent trial estimates for the position vector and clock bias, and Ar, At are small quantities. The index (i) in eq (C6) labels the i*'* trial value
in a repeated or iterative process of solving eqs (C.5).
We substitute {C.Q) into (C.5), expand and neglect squares of small quantities Ar
and At. Then the equations reduce to

R{^) Ar-c' ^(.) _ ^c R;

At = 1

R (')
T.

-

M)

(C.7)

where

Tf^ = tRj + 6^' - tc ,

(C.8)

and
Rf = r^;) + ycT^;^

(C.9)

This system of equations is of the same form as that discussed in Appendix B and can be

solved by iteration in a similar manner as that discussed at the end of Appendix B. The

main difference from eq (B.3) is that the velocity estimate vc appears in eq (C.7) in three

places:

both on the left and the right

sides in

the definition of R^(i) ,

and in the coefficient

of At on the left side.

Thus the iteration scheme proceeds as follows.

Step 1. Choose tc; estimate vq in the ECIF. This will be a combination of
velocity relative to the ground, and velocity of the ground relative to the chosen ECIF, due to earth rotation.

40

Step 2. Choose approximate values for the user position re and clock bias he. Compute ¥if and T^/K
Step 3. Solve eqs (C7) to obtain corrections Ar and At for user position
and clock bias. Step 4. Test whether the corrections are sufficiently small that computa-
tion can be halted; if not combine the corrections obtained in Step 3 with the
approximate values of position and bias, to obtain a new approximate position and bias.
Step 5. Recompute R^*^ and T-*^ and return to Step 3.
41

Appendix D.

Equivalence of Integrated Doppler
and One-Way Light Time

In this Appendix we show the complete equivalence of the one-way light time mea-
surement to integrated Doppler measurements. Neglecting gravitational path bending and gravitational time delay is well justified;
electromagnetic signals travel with speed c in straight lines in the ECI reference frame.
N The direction of propagation may then be described by a constant unit vector along the
path from transmitter to receiver. Let the coordinate time of transmission of a signal be denoted by tT and the coordi-
nate time of reception be denoted by tji. The positions of the transmission and reception events in ECI coordinates are denoted by rT{tT) and rR{tR). Also, the velocities of the
V transmitter and receiver at the transmission and reception epochs are denoted by and v,
respectively; this notation is introduced in eqs (1-2) where the Doppler effect was stated. In this Appendix we derive the Doppler formula including gravitational frequency shifts.
The one-way propagation of a signal from transmitter to receiver gives

K tR = + tT + -IrnitR) - ftI^t)! ,

{D.l)

K where represents a possible systematic error or bias in the coordinate time of the receiver
clock. The term \rji{tfi) — rT{tT)\/c is the propagation delay in the equation, which must
be accounted for in any use of GPS signals for navigation or time transfer. We show here

that if a receiver properly accounts for this term, then it has automatically accounted for

the relativistic Doppler effect, integrated over time.

Let there be a second signal sent from transmitter to receiver after one period of

oscillation of the oscillator in the transmitter. Let the times of transmission and reception

of this second signal be tT + dtx and tn -\- dtn. The frequencies will be very high, so the
K small differentials may be calculated by taking differentials of eq (D.l). The term in

drops out, and

— dtR = dtT+ --. 7-—

-—- -{drR-drT)

,

c\rR{tR) -rritTn

{D.2)

where dvR and drx are the changes in positions of the receiver and transmitter during the
N time intervals dtR and dix- But the unit vector is exactly

\rR{tR) -rT{tT)\
and during the small time intervals dtR, dtx we have
dvR = \dtR, dvT = VdtT .
42

{DA)

Therefore

dtR = dtr + -N • {\dtR - VdtT) •

•

(D.5)

Rearranging gives

— —j /
,
dtR ( 1

N-v\
j

=

,
dtT

/
1
(

N-V\
.

{D.6)

Next we relate the coordinate time intervals dtR, dtT to the proper frequencies of

received signal and transmitted signal. This part of the calculation is only carried to terms

of order c~^ . Let $(r) be the Newtonian gravitational potential at position r in ECI

coordinates. Let

A$ = $-$o,

(^-7)

where the constant $o is the effective gravitational potential on the rotating geoid, in-
cluding the centripetal contribution arising from earth's spin. Then the proper time dr

that elapses on an atomic clock as it moves through the position increment dr during the

coordinate time increment dt is given approximately according to general relativity by [2,3]

^ c^dT^ =[l + t^\c^dt^- {dvf = h +

i^^^ -

) c^dt^ .

{D.8)

This may be solved for the coordinate time increment dt in terms of the proper time:

— dr

dr

dt =

,
v/(l

+

2A$/c2-(drM)2/c2)

~

l+

^-liJrLr-MJliTlTT

_ .

.

{D.9)

where the square root has been expanded using the first term in a binomial expansion,
y/l + x/c^ ~ 1 + x/{2c^). This may be applied to the coordinate time intervals dtR and
dtx between receptions and transmissions of the two successive signals. Accurate to terms
of order 1/c^,

dtR^

rx7-T

T-

(D.IO)

dtT =

dTT T^n-A

^ "•"

c2

T-
2 c2

We substitute these results into eq (D.6) and obtain

(^11)

^) dra (1 -

. A^(rfl) _ IyI

^ "•"

C2

2 C2

^) drr (l -

_ A^(rT)

1

I

1 V^

^ '^

C^

2 C2

{D.12)

The last step is to relate the proper time intervals to the proper frequencies. Using
F and / for the proper transmitted and received frequencies, respectively, we find

1
.
dTT

f= 1
dTR

{D.13)

43

.

And substitution into eq (D.12) gives the final relation between proper frequencies:

1_N^ / (1 + A$(ri?)/c2 - v2/2c2) - F (1 + A$(rT)/c2 - Y'^ I2c^)

1_N:X

•

c

(^-14)

Upon neglecting the gravitational frequency shifts, this form agrees with eq (1) in the main
text to the required order of approximation. Several observations about the result in eq (D.14) are appropriate. First, all quanti-
ties, except the proper frequencies, are measured in ECI coordinates. Second, there are no
terms of the form v • V/c^. Third, to obtain the ratio f I F or F// correct to order 1/c^, we
N N must take care to expand the aberration factors (1 — • \/c^)~^ and (1 — • V/c^)~^ to
the required order. To illustrate this last point, let us calculate the ratio F//. Expanding
eq (D.14) to second order gives

P/,

^(rH)-^(rT)
,
,

v2

V2

N- v-V N-vN-V // N-v^x^2

+

^

+'

(^-15)

^

^

Conservation of coordinate frequency. Consider the special case when both the transmitter and the receiver are at rest in the ECI reference frame. Then eq (D.6) gives
dtji = dtT
This means that the coordinate time between reception of successive wavefronts is the same as the coordinate time between the emission of successive wavefronts. The coordinate
frequency is the reciprocal of this time, so this result states that for observers at rest in the
ECI reference frame, coordinate frequency is conserved. This is a well-known fact that can
be proved by separate arguments whenever the metric tensor coefficients do not depend on the time. This is the case here.
We may readily reverse the argument given above and pass from eq (D.14) back to
eq (D.l) by integration rather than differentiation. Equation (D.6) is an intermediate step in such an argument, which is obtained by combining eqs (D.IO, D.ll, and D.12) so as
to eliminate the proper time intervals drji and drT- Then eq (D.2) follows from (D.3-5). Integration then yields eq (D.l) with an unknown constant of integration.
Thus the change in the one-way light time, as used in GPS for pseudorange, is equal to
the integrated Doppler frequency shift over the interval of the change using the frequency
relation (D.13).

44

.

.

Appendix E.

Numerical Check on the First-Order Doppler Correction, Equation (16)

In this Appendix we use data provided in Example 5.1 to compare the "true" value

of Atjij from the simulation with values obtained using eq (16). For the case of satellite 1

in Example 5.1,

Atsv = -0.000 056 554 796 2 s.

Then eq (16) gives

Atm = -0.000 056 554 784 s

This estimate using first-order Doppler on the coordinate time may be compared with
the simulation which gives for the two arrival times of the position in the code train corresponding to tsv and tc the following time difference:

Atm =240.069 124 140 263 - 240.069 180 695 048 s = - 0.000 056 554 785 s
— The extremely good agreement between these numbers one calculated using the — Doppler formula and one using propagation delay equations suggests the following strat-
egy for measuring the coordinate time of arrival of the code corresponding to the chosen
measurement time tc- If we choose a time tc as some multiple of 1 ms so that it lies on a code boundary, the SV clock correction Atsvj can be computed from the navigation message. The receiver can then measure the arrival time tjij of the code boundary
tsvJ which is equal numerically to tc T (arrival of phase time = tc)- Then the Doppler
correction factor given in eq (16) can be computed by the receiver and added to tRj to
+ predict, with excellent accuracy, the arrival time tjij dtjij of the code which corresponds
to transmission of tc from the satellite. Thus only times of arrival of code boundaries would actually be measured by the receiver.
Three additional calculations based on eq (16) are given in Table E.l, for the four satellites used in Example 5.1. The simulated intervals Atsvj are given in the first column
of Table E.l. The second column lists the values of At^j obtained using eq (16). For comparison, the last column gives the "true" values of Atjij obtained from the simulation.
In all cases the first-order Doppler correction, eq (16), agrees extremely well with the
"true" values. The corrections have been incorporated into the data listed in Table 7.

45

Table E.l. Computation of first-order Doppler corrections.
= (units s)

ID #
1
2 3 4

AtsVj
-0.000 056 554 796 2 0.000 493 792 700 6 0.000 259 725 929 9 0.000 305 714 376 7

^tRj eq (16)
-0.000 056 554 784 0.000 493 792 194 0.000 259 725 401 0.000 305 713 988

AtRj
-0.000 056 554 785 0.000 493 792 194 0.000 259 725 400 0.000 305 713 988

46 U.S. GOVERNMENT PRINTING OFFICE I 999-0-773-022/40 1 3 I

NIST Technical Publications
Periodical
Journal of Research of the National Institute of Standards and Technology— Reports NIST research and development in those disciplines of the physical and engineering sciences in which the Institute is active. These include physics, chemistry, engineering, mathematics, and computer sciences. Papers cover a broad range of subjects, with major emphasis on measurement methodology and the basic technology underlying standardization. Also included from time to time are survey articles on topics
closely related to the Institute's technical and scientific programs. Issued six times a year.
Nonperiodicals
Monographs— Major contributions to the technical literature on various subjects related to the Institute's
scientific and technical activities.
Handbooks— Recommended codes of engineering and industrial practice (including safety codes)
developed in cooperation with interested industries, professional organizations, and regulatory bodies.
Special Publications— Include proceedings of conferences sponsored by NIST, NIST annual reports, and
other special publications appropriate to this grouping such as wall charts, pocket cards, and
bibliographies.
Applied Mathematics Series— Mathematical tables, manuals, and studies of special interest to physicists, engineers, chemists, biologists, mathematicians, computer programmers, and others engaged in scientific and technical work.
National Standard Reference Data Series— Provides quantitative data on the physical and chemical properties of materials, compiled from the world's literature and critically evaluated. Developed under a worldwide program coordinated by NIST under the authority of the National Standard Data Act (Public
Law 90-396). NOTE: The Journal of Physical and Chemical Reference Data (JPCRD) is published bimonthly for NIST by the American Chemical Society (ACS) and the American Institute of Physics (AIP). Subscriptions, reprints, and supplements are available from ACS, 1155 Sixteenth St., NW, Washington,
DC 20056.
Building Science Series— Disseminates technical information developed at the Institute on building materials, components, systems, and whole structures. The series presents research results, test methods, and performance criteria related to the structural and environmental functions and the durability and safety characteristics of building elements and systems.
Technical Notes— Studies or reports which are complete in themselves but restrictive in their treatment of a subject. Analogous to monographs but not so comprehensive in scope or definitive in treatment of the subject area. Often serve as a vehicle for final reports of work performed at NIST under the sponsorship of other government agencies.
Voluntary Product Standards— Developed under procedures published by the Department of Commerce in Part 10, Title 15, of the Code of Federal Regulations. The standards establish nationally recognized requirements for products, and provide all concerned interests with a basis for common understanding of the characteristics of the products. NIST administers this program in support of the efforts of private-
sector standardizing organizations.
Consumer Information Series— Practical information, based on NIST research and experience, covering
areas of interest to the consumer. Easily understandable language and illustrations provide useful background knowledge for shopping in today's technological marketplace.
Order the above NIST publications from: Superintendent of Documents, Government Printing Office,
Washington, DC 20402.
Order the following NIST publications—FIPS and NISTIRs—from the National Technical Information Service, Springfield, VA 22161.
Federal Information Processing Standards Publications (FIPS PUB)— Publications in this series
collectively constitute the Federal Information Processing Standards Register. The Register serves as the official source of information in the Federal Government regarding standards issued by NIST pursuant to
the Federal Property and Administrative Services Act of 1949 as amended, Public Law 89-306 (79 Stat. 1127), and as implemented by Executive Order 11717 (38 FR 12315, dated May 11, 1973) and Part 6 of Title 15 CFR (Code of Federal Regulations). NIST Interagency Reports (NISTIR)— A special series of interim or final reports on work performed by NIST for outside sponsors (both government and non-government). In general, initial distribution is
handled by the sponsor; public distribution is by the National Technical Information Service, Springfield,
VA 22161, in paper copy or microfiche form.

U.S. Department of Commerce
National Institute of Standards and Technology 325 Broadway Boulder, Colorado 80303-3328
Official Business Penalty for Private Use, $300