GPS Timing-Code-Based Positioning in Mathematica

Thomas H. Meyer

Department of Natural Resources and the Environment
University of Connecticut

Background and Big Picture

Name: U.S. NAVSTAR Global Positioning System (GPS)

Distances from known points to a point of interest: multilateration

Determine ranges to the satellites

Determine SVs’ positions

ranges + satellite positions => receiver position

More specifically

Space-based timing and ranging:

constellation of 24-32 satellites in medium Earth orbits

broadcasting radio transmissions L-Band

timing codes and user information (the navigation message)

timing codes used to determine ranges to the satellites

Navigation message used to determine SVs’ positions

multilateration requires a nonlinear least squares computation

The GPS Constellation

These are pictures of a Block II a GPS satellite that is on display at the San Diego Air and Space Museum in Balboa Park.

TabView[{,,,}]

The satellites (abbreviated SV for “space vehicle”) are arranged in orbital planes nominally inclined 55° to the plane of the equator at a nominal orbital radius of around 26560 km (or altitude of around 20200 km). The orbits have very little eccentricity, they are nearly circular.

By design, there are six orbital planes and four SV’s per plane = 24 SVs. The SVs are not necessarily spaced evenly around the orbits and, in fact, there are 29 SVs available to users at the moment (4 Oct 2011). This SV configuration is called the “bird cage.”

({{, }})

Radio transmissions in the L-Band cannot penetrate through the Earth, so ground-based receivers can receive transmissions only from SVs above the local horizon. Therefore, receivers cannot usually acquire more than a dozen SVs or so, and often are limited to much fewer due to interference from buildings, trees, canyon walls, and other obstructions.

({{, }})

The Positioning Geometry Problem

Trilateration: given three points in space -- called control points -- and distances from those points to a point-of-interest (POI), determine the coordinates of the POI.

Any Cartesian coordinate system will do. For GPS, the coordinate system is called a reference frame and its name is the World Geodetic System 1984 or WGS 84. This reference system has been realized more than once. The realizations are named according to the GPS week the realization was put into use. The current realization is G1150, and a new one is due out next year.

The coordinate-system origin is at Earth’s center-of-mass. Z-axis is parallel with a mean rotational Pole and the X-Y axes form the equatorial plane. The X-Z plane is the plane of the Prime Meridian. Y-axis forms a right-handed system.

The (errorless) observation equation is Pythagoras’s formula in three dimensions. The distance from SV i to receiver j is

(Subscript[ρ, j]^i)=(Sqrt[((ΔX)^2)+((ΔY)^2)+((ΔZ)^2)])

where ΔX is the difference of the SV’s X coordinate (known) and the receiver’s X coordinate (unknown, sought), etc.

Satellite Positions

Although the SVs are in nearly circular, medium Earth orbits, they are far too close to treat Earth as a point mass. Its varying gravitational field, and other environmental factors such as the solar wind, prevent the orbits from being strictly Keplerian. Therefore, GPS SVs broadcast a navigation message that includes the six Keplerian elements plus time-dependent element correction terms, among other things. (See the Interface Control Document for full details).

With the broadcast ephemerides, a receiver can compute a SV’s position simply by knowing the transmission time.

Ranging With Timing Codes

GPS SVs broadcast a civilian timing code called the “coarse/acquisition code” (C/A code) on the L1 frequency. Receivers use these codes to roughly synchronize their internal clocks with GPS time and to measure ranges to the SVs.

The codes are comprised of “chips.” Each (entire) C/A code has 1023 chips, and each (entire) code is transmitted in 1μs, so each chip is 1×(10^(-3))/1023=9.77517×(10^(-7)) sec multiplied by the speed of light gives 293.052 m per chip. Timing codes are pseudo-random noise codes (PRN). The least expensive receivers measure ranges to around 300 m, better receivers measure ranges to around 30 m.

The codes allow the receiver to determine when the transmission left the SV. The receiver has its own internal clock. It reads its own clock at the moment of observation and subtracts the time of transmission yielding the time it took for the transmission to travel from the SV to the receiver. Geometric separation between the receiver and the SV (called the range ρ to the SV) equals time-of-flight of the transmission multiplied by the speed of light.

ρ=Δt c

The receiver’s clock is not synchronized with GPS time. This prevents ranges from being measured without bias. This unknown time bias (Subscript[ϵ, Δt]) introduces a bias into the observed range equal to c (Subscript[ϵ, Δt]), where c is the speed of light. The bias can be large: 60 km and more is possible. Because the range is bias, a range determined using timing codes is called a pseudorange and denoted r.

r=(Δt+(Subscript[ϵ, Δt])) c

The Positioning Problem

The observation equation is

(Subscript[r, j]^i)=(Sqrt[((((Subscript[X, j])-(X^i)))^2)+((((Subscript[Y, j])-(Y^i)))^2)+((((Subscript[Z, j])-(Z^i)))^2)])+c (Subscript[Δt, r])

There are four unknowns -- (Subscript[X, j]), (Subscript[Y, j]), (Subscript[Z, j]),(Subscript[Δt, r]) -- so four or more SVs must be acquired. The unknowns are found using non-linear least squares. The iteration typically converges in three to six steps even if the initial guess for the receiver’s position is, say, at (0,0,0). Data can be collected over time. SV positions and ranges are computed at each observation epoch resulting in a single least-squares position estimate. Expected accuracy is better than 3 m in the absence of environmental errors.

Solve[{ρ1==(Sqrt[(((x1-x))^2)+(((y1-y))^2)+(((z1-z))^2)]),ρ2==(Sqrt[(((x2-x))^2)+(((y2-y))^2)+(((z2-z))^2)]),ρ3==(Sqrt[(((x3-x))^2)+(((y3-y))^2)+(((z3-z))^2)])},{x,y,z}]//Simplify

D[(Sqrt[(((x1-x))^2)+(((y1-y))^2)+(((z1-z))^2)]),x]

Input Data

Receiver observations are downloaded off a receiver in manufacturer-specific, proprietary format. These data are converted to

RINEX: Receiver-INdependent EXchange format

Mathematica does not import RINEX files (yet), so I wrote my own Import routines. RINEX files are pretty complicated in structure and full of highly technical, detailed content. My code currently reads the data into RAM but I want to have it populate a relational database instead. Unfortunately, the database interface in Mathematica does not function properly with Access databases and I have not yet installed an enterprise-level alternative, so this is work-in-progress.

End of year problem at NGS

Sometimes it just doesn’t work. That SV is actually on the other side of the Earth

Thomas H. Meyer
Department of Natural Resources and the Environment
University of Connecticut
thomas.meyer@uconn.edu