fcT'L INST. OF STAND & TEf"

f.:]f.

IIIDS blfiT17

NIST

PUBLICATIONS

NisrUnited 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 implement receiver designs consistently with the requirements of relativity. We 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 synchronization). 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. For satellites in GPS 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.

(3) Gravitational frequency shift. A 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

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

straight line in an ECEF 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 measurement of the difference between the time associated with the phase of a GPS signal

and the time on the local clock at the arrival of that phase. We 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 earthcentered 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^ riMciviMsystem 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

Euclidean straight lines. The time for a signal to travel from transmitter at R 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
in Appendix D with accuracy sufficient for GPS. This relationship is [2,3,4]

/ (1 + A^{rR)/c^ - v2/2c2) F{1 + A^{rT)/c^ - y^/2c^)

. Nv ~

1 NV ' {^-^^J

c 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-^)'

where

(1)

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

position solution) and integrated Doppler frequency (to obtain changes in position). A

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 WGS84 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

^ECI =

xeci

Veci

zeci

(3)
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

XeCEF = R {t)^ECI —

cos0(t) sin0(t)

-sin0(t) cos0(i)

1

XeCI

Veci

ZECI

(4)

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

larger than 0.05 mm 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

can determine the offset of the locally generated PRN 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 eccentricity, and E 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

satellite positions at the instant of transmission in ECEF 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 corrections may be small. If the local clock is sufficiently close to GPS 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 geometric 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

process. We 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 system and the correct GPS system time, the receiver's ECI position co- ordinates are rotated into the ECEF reference frame at the instant cor- responding to the measured reception time. If the ECI frame chosen

in step (4) is sufficiently close to the ECEF 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 approximately 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 Multi- 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

ECEF latitude and longitude 35° N, 0° E, and distance R = Q 378 136.3 m 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 comparison with the navigation solution, will be

WGS84

i? cos 35°

E sin 35°

5 224 663.389 m

3 658 348.690 m

(7)

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.

(2) We apply the prescription according to Sect. 20.3.3.3.3.1 from the ICDGPS-200 describing the user algorithm. We then obtain the system time t^, the GPS 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# Transmission epoch ti {s) X^ (m) y^ (m) Zj (m)

1 37 239.924 422 365 6 13 005 878.255 18 996 947.213 13 246 718.721

2 37 239.920 713 391 8 20 451 225.952 16 359 086.310 -4 436 309.875

3 37 239.925 307 870 20 983 704.633 15 906 974.416 3 486 595.546

4 37 239.929 346 353 9 13 798 849.321 -8 706 113.822 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

of each satellite are given in the ECEF 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

tc = 37 239.000 000 000 s. We 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

{R^{t3-tc)y

cosCte{t3 -tc) -sinpeitd -tc)

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

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

matrix is different for each transmission epoch. We use a subscript C to indicate

that the inertial system has been arbitrarily chosen. The results from transforming 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

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

13

(8)
Table 2. Transmitted data transformed to the chosen ECI system

SV # Transmission Epoch t^ (s) Xa (m) ya (m) ZCi (ni)

1 37 239.924 422 365 6 13 004 597.642 18 997 823.895 13 246 718.721

2 37 239.920 713 391 8 20 450 127.566 16 360 459.358 -4 436 309.875

3 37 239.925 307 870 20 982 631.270 15 908 390.245 3 486 595.546

4 37 239.929 346 353 9 13 799 439.294 -8 705 178.668 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

(r(^) _ r^.) . Ar - 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.

The calculation converges to within 1 m 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

R^tR-tc) =

COSfteitR -tc) - sin ft,,{tR - tc]

sniVt^{tR-tc)

COS^e{tR-tc)

1

(11)

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

WGSBA

5 224 663.388 m

• 0.000 m

3 658 348.689 m

(12)

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

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

Trial # User's X position User's y position User's z position User's clock time

(start) m m m 37 239.999 422 365 6 s

1 5 057 363.392 m 2 355.126 m 3 541 092.792 m 37 239.997 532 305 8 s

2 5 226 931.552 m 354.224 m 3 659 938.391 m 37 240.000 033 455 9 s

3 5 224 663.780 m 380.983 m 3 658 348.973 ui 37 240.000 000 006 s

4 5 224 663.374 m 380.988 m 3 658 348.689 ni 37 240.000 000 000 s

5 5 224 663.374 m 380.988 m 3 658 348.689 m 37 240.000 000 000 (1 s

15
Table 4. Iterative solution ignoring earth rotation.

1(staxt)

1 5 057 359.485 2 012.579 3 541 090.386 m

2 5 226 927.086 00.194 3 659 935.704 m

3 5 224 659.325 27.107 3 658 346.291 m

4 5 224 658.919 27.112 3 658 346.008 m

5 5 224 658.919 27.112 3 658 346.008 m

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

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 propagation 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

in Table 1. We 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 efficient 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)

to obtain inertial coordinates of the transmission events. We will first make a reasoned

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

estimate for Ir. We 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)

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.

(start)

1 5 057 363.550 1 986.551 3 541 092.792

2 5 226 931.564 -25.990 3 659 938.391

3 5 224 663.794 -0.018 3 658 348.973

4 5 224 663.388 0.000 3 658 348.689

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

(start) 3 313 469.280 -2 090 568.570 5 032 997.819 37 240.004 346 353 9

1 5 190 759.859 400.888 3 634 586.629 37 239.999 499 918 6

2 5 224 663.389 -0.892 3 658 412.593 37 240.000 000 000

3 5 224 663.389 -0.000 3 658 348.689 37 240.000 000 000

4 5 224 663.389 0.000 3 658 348.689 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 approach 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

following discussion, we shall use a subscript R 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 propagates 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 measurement of the arrival time. We 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 tcThe 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)

where V 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)

The receiver needs estimates of the receiver velocity v and the unit vector N 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)

holds. Note that there is no subscript R 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

at the system time tc be re and vc 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 coincides 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^iT(7^^)-v..Ry))A.4 c^(lj^O'-K) (21)

where

and

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

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

for each clock. The differences will be less than 1 ms. We 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-

IV, for determining the .x, t/, z coordinates of the satellite at the GPS

time tc22
(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 position 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 re- ceiver 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 polynomial 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 tcAgain we linearize the propagation delay equations and solve Ihcm

iteratively. We 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

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

will not discuss these options here. We 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 correction. 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

X,

i? COS 35° 5 224 663.389 m

iyG584 =

i? sin 35° —

3 658 348.690 m

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

vc

-Vy sin 30° sin 35°

t;^ cos 30° + nei?cos35°

Vg sin 30° cos 35°

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.

122.872 806 643

(25)

(26)

(2-

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

ECIF to be

"-86"

vc = 641 m/s.

123 J

Estimating receiver position. We 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):

N =

\rR{tR) -ftI^t)!

'

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 re- ceive time, with bias of correction 423 ns, the receiver position vector is actually

(28)

{D.3)

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

5 224 656.660

50.117

3 658 358.300

m

(29)

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

ri?i

6 000 000

23 000

3 100 000

m, (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:

Ni

-0.437 588

-0.114 125

-0.891 904

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

(32)

Ni

-0.477 846 083

-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 tcjI 15 126 951.488

2 20 996 573.361

3 23 391 384.147

4 10 656 983.742

ID # Satellite Xj (in) Satellite IJj (m) Satellite Zj (m) Local clock time tcj (s)

2 403 354.115 21 702 797.709 240.069 124 140 3 -15 349 939.783 -5 178 329.880 240.079 109 466 3

10 956 612.394 —4 333 140.777 240.075 620 243 9

-11 270 652.822 21 703 578.299 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 Appendix C is given at successive stages in Table 9. At the initial stage, estimates

for the ECEF 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 process 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 beTrial # Receiver Xq (m) Receiver yc (m) Receiver Zc (m) Local clock bicus be (s)

0.0 0.0 0.0 0.000 000 300

1 5 051 460.568 -9 835.265 3 545 131.610 -0.002 453 425 2

2 5 226 991.162 132.182 3 659 870.280 0.000 032 544 3

3 5 224 663.799 -0.023 3 658 349.573 -0.000 000 417 2

4 5 224 663.390 -0.000 3 658 348.690 -0.000 000 423

5 5 224 663.390 -0.000 3 658 348.690 -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 # satellite Xj (m) satellite yj (m) satellite Zj (m) tcj ~ tjij -\- /\tsvj (s)

1 15 126 951.488 2 403 354.115 21 702 797.709 240.069 124 140 3 2 20 996 573.361 -15 349 939.783 -5 178 329.880 240.079 109 466 8

3 23 391 384.147 10 956 612.394 -4 333 140.777 240.075 620 244 5

4 10 656 983.742 — 11 270 652.822 21 703 578.299 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 transmission 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 # Receiver Xc (m) Receiver yc (m) Receiver Zq (m) Local clock bias be (s)

0.0 0.0 0.0 0.000 000 300

1 5 051 461.207 -9 835.359 3 545 131.862 -0.002 453 427 6

3 5 226 991.866 132.086 3 659 870.566 0.000 032 542 4

4 5 224 664.501 -0.072 3 658 349.243 -0.000 000 419 1

5 5 224 664.092 -0.096 3 658 348.976 -0.000 000 424 9

6 5 224 664.092 -0.096 3 658 348.976 -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

\v{tR)-K{tT)\ ^^ \r{tR)-K{tn)\ [r(t;,) - R(^h)] • V .... to =tR — tT = « 5 (3b)

c c c^

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) AtSagnac

tie r{tT) X

C2

[r{tT)-K{tT)]

^e r{tT)x

C2

[-R(iT)]

He

C2

• K{tT) X T{tT)

(40)

X r{tT)

c2

^e • ^ECEPitT) X ^ECEFitr)

c2

l^el {XECEF{iT)yECEFitT) — 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:

AtSagnac = "^^S , (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 Positioning 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

Surface Weapons Center TR 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 Relationships with Local Geodetic Systems," DMA 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, Observatoire 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, Washington, 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 ICDGPS-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 ^VGMA . ^

Atr = - 2 ^ esmEk

c^ = Fe^sinEk (^-1)

2R V

where

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

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

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

is the satellite's eccentric anomaly, R is the instantaneous position vector of the SV at

transmission time, and V 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 equations 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 computation 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 may also need the receiver's acceleration. We 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

important topic is discussed in detail in Sect. 3. We 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 T.

(') - RM) (C.7)

where

and

Tf^ = tRj + 6^' - tc ,

Rf = r^;) + ycT^;^

(C.8)

(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

(i) places: both on the left and the right sides in the definition of R^ , 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 computation 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 measurement 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.

The direction of propagation may then be described by a constant unit vector N along the

path from transmitter to receiver.

Let the coordinate time of transmission of a signal be denoted by tT and the coordinate 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

transmitter and receiver at the transmission and reception epochs are denoted by V 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

tR = K + tT + -IrnitR) - ftI^t)! , {D.l)

where K 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

small differentials may be calculated by taking differentials of eq (D.l). The term in K

drops out, and

dtR = dtT+ --.—7-— -—- -{drR-drT) , {D.2) c\rR{tR) -rritTn

where dvR and drx are the changes in positions of the receiver and transmitter during the

time intervals dtR and dix- But the unit vector N is exactly

\rR{tR) -rT{tT)\

and during the small time intervals dtR, dtx we have

dvR = \dtR, dvT = VdtT . {DA)

42
Therefore

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

• (D.5)

Rearranging gives

, / N-v\ , / N-V\

dtR ( 1 —j = dtT ( 1 —j . {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, including 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 = , ~ — rr-JiTTT {D.9)

v/(l + 2A$/c2-(drM)2/c2) l + ^-liJLMll

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 T- (^11)

^ "•" c2 2 c2 We substitute these results into eq (D.6) and obtain

dra (1 - ^) drr (l - ^)

. A^(rfl) _ IyI 1 I A^(rT) _ 1 V^

^ "•" C2 2 C2 ^ '^ 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

{D.13)

1 dTT

. f =

43

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

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

1_N:X - 1_N^ • (^-14)

c 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 quantities, 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

must take care to expand the aberration factors (1 — N • \/c^)~^ and (1 — N • 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

/ x2 (^-15) N- v-V N-vN-V /N-v^^ ^

^ + ^ +

'

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 (arrival of phase time T = 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 # AtsVj ^tRj eq (16) AtRj

1 -0.000 056 554 796 2 -0.000 056 554 784 -0.000 056 554 785

2 0.000 493 792 700 6 0.000 493 792 194 0.000 493 792 194

3 0.000 259 725 929 9 0.000 259 725 401 0.000 259 725 400

4 0.000 305 714 376 7 0.000 305 713 988 0.000 305 713 988

46

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


NISTTechnical 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 bi- monthly 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 privatesector 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