Calculating Position from Raw GPS Data

We are all familiar with GPS (Global Positioning System) and its myriad applications. From getting directions using Google maps to hailing a ride using a ride sharing app, countless individuals and businesses rely on accurate position estimation using GPS. Position estimation using GPS is now so accurate that GPS is being used for measuring plate tectonics and continental drift. Indeed, GPS is so much a part of our lives that few of us stop and wonder how it actually works. The purpose of this post is to show you how.  More specifically, we’ll first learn how position is calculated from range measurements to GPS satellites and then consider a concrete example where we’ll process raw data collected by a commercial GPS receiver to obtain user position estimates.

GPS is owned by the US government and operated by the US air force. As such, access to GPS data can be degraded or denied to other nations due to geopolitical concerns. For example, according to this wikipedia article, access to GPS data was denied to the Indian military during the Kargil conflict. Other countries therefore decided to implement their own GPS like systems. Examples are GLONASS by Russia, Galileo by the European Union and BeiDou by China. This may seem like a waste of resources, but many of these systems operate using the same standards as GPS and are interoperable. For example, latest GPS receivers can receive signals from both GPS and GLONASS satellites, improving positioning accuracy. The wikipedia article referenced above also provides an interesting history of the origins of the GPS system dating back to Sputnik, the first man made satellite launched in space.

The post is divided into 3 parts – the first part describes the various coordinate systems used to express position. An understanding of these coordinate systems is essential for understanding how GPS works and is useful in its own right. For example, the concept of height above the earth surface is something that people intuitively understand, but turns out to be tricky idea to precisely define. The second part describes the basic principle of triangulation used by GPS to calculate position and develops the mathematics for calculating the user position from satellite measurements. A casual user, interested in a general knowledge of how GPS operates but not in the details can read part 1 and the first few sections of part 2 before the math starts.  Finally, the appendix provides the Matlab code for implementing various ideas discussed in this post.

Most of the information in this post derives from the book “Global Positioning Systems: Signals, Measurements and Performance” 1, an excellent book about GPS systems. The matlab code presented in this post was written as solutions to the various end of chapter problems in the book. I also make use of material in the GPS Interface Specification document2. My contribution is to present the keys ideas related to GPS in a (hopefully) clear and concise manner and show the Matlab code that implements the positioning algorithms on real world data collected using commercial GPS receivers that can be purchased from the internet. After reading this post, a casual reader should have a better understanding of GPS, so the next time they use Google maps on their phone, they are aware of the amazing technical infrastructure that makes positioning and navigation applications possible. A more technical reader should be able to replicate my experiments by collecting their own data and running the Matlab code.

Part 1

The fundamental task of a GPS system is to calculate position. To do so, we must first define what position means. Most people think about latitude/longitude and height when they think about position. While latitude/longitude are useful to represent a position on the earth surface, they are not suitable for mathematical calculations as they don’t provide a cartesian coordinate system. The physical distance represented by a unit difference between two longitudes is not constant, and depends on the position. For example, the distance between two longitudes one degree apart is greatest at the equator and approaches 0 at the poles.

To be useful for mathematical calculations, a coordinate system where a unit difference between coordinates represents a constant physical distance is required. A family of such coordinate systems can be created by defining a set of perpendicular axes intersecting at an origin that is rigidly attached to the earth. Such coordinate systems are called “ECEF” (Earth Centered, Earth Fixed). These coordinate systems work well to express the position of a user on earth as they rotate with the earth and the position of a stationary user on the earth surface is constant.

The most commonly used ECEF coordinate system is called the WGS 1984 system developed by the department of defense. WGS 84 is an ECEF frame, defined as follows:

  • Origin at the center of mass of the earth
  • z axis passing through the CTP (Conventional Terrestrial Pole). CTP is the average of the earth’s pole’s position between the years 1900 and 1905. An average needs to be used as the position of the pole is not fixed and wanders around in a circle of radius ~15m.
  • x axis passing through the intersection of the CTP’s equatorial plane and a reference meridian, often called the Mean Greenwich Meridian

ECEF frames are convenient to represent positions with respect to the earth, but they are not inertial as they are rigidly attached to a spinning earth. To formulate the problem of satellite motion around the earth in accordance with Newton’s laws of motion, we need an inertial coordinate system in which to express acceleration, position and velocity vectors.

An inertial frame can be defined in a manner similar to the ECEF frame except that the x-axis points to the vernal equinox (the direction of intersection of the equatorial plane of the earth with the plane of the earth’s orbit around the sun). Strictly speaking, this frame isn’t inertial either since the earth is moving around the sun. However, it can be considered inertial over short periods (just as the ECEF frame can be considered inertial for analyzing the movement of a body on earth over a short period).

In the rest of this post, we’ll be using the WGS 84 ECEF frame to express the user and satellite position. This is possible because the GPS interface specification document (referred to as GPS-IS in the rest of this post) provides a step by step procedure to calculate the satellite position in the ECEF frame at a given time instant that accounts for the non-inertial nature of the ECEF frame. This is fortunate because ECEF frames are ideal to express the position of a user on the surface of the earth as they rotate with the earth and thus the position of a stationary user is constant, as one would expect. This would not be the case if an inertial reference frame is used to express position as such a frame would be fixed in space, and thus would appear rotating with respect to the earth.

Let’s now turn to the issue of defining height that we alluded to earlier. The key question that needs to be answered when defining height is “with respect to what”? Consider for example the height of a point in a desert. If we define height to be the distance of the point from the “ground”, then such a height would change constantly as the level of the ground itself changes due to natural phenomenon (for example, the wind depositing more sand).  Thus, to define height, we must first define “ground” (note that strictly speaking, this issue doesn’t arise in a ECEF reference frame, as all distances are measured with respect to the center of mass of the earth. However generally people are not interested in distance from the center of the earth, but from the local earth surface). To define “ground” in a consistent manner, we need a model for the surface of the earth. Two models are generally used:

Model 1: Reference Ellipsoid

The surface of the earth is modeled as an oblate ellipse called the “reference ellipsoid”, centered at the earth’s center with the axis of revolution coincident with the z axis of an ECEF frame. The lengths of the semi-major and semi-minor axes (denoted as a and b) are a = 6378137m and b = 6357002m. As expected, the value generally used for earth’s radius (with the earth modeled as a sphere) is 6371Km which lies between the lengths of the semi-major and minor-axis. The reference ellipsoid is merely an abstraction for the shape of the earth. It doesn’t have any physical significance. An actual point on the earth surface will usually lie above or below the reference ellipsoid.

Model 2: Geoid

Another surface commonly used to measure height that does have physical significance is called the geoid. This surface is defined as the locus of all points with the same value of the gravitational potential. If the earth was a sphere with a uniform composition, this surface would be a regular surface that can be parameterized as a mathematical function. However since the earth is neither spherical nor uniform, the geoid is usually specified as a series of heights above the reference ellipsoid. Height relative to the geoid is called the orthometric height, or height above the mean sea level (MSL). This makes sense as if the oceans covered the surface of the earth, the shape of the oceans would be a close approximation to the geoid. Strictly speaking, the shape of the geoid itself is not fixed as the surface of the earth is constantly changing due to human activities and natural phenomenon. However given the bulk of the earth, these only have a negligible impact on the earth’s gravitational field and thus on the shape of the geoid.

Note that the height of the geoid is specified with respect to the reference ellipsoid and the height of a given point can be specified with respect to either. The relationship between the ellipsoid, geoid and the actual surface of the earth is shown below.

Latitude and Longitude

Armed with an ellipsoidal representation of the shape of the earth and with an ECEF coordinate system, we can define the position of a point P in ellipsoidal coordinates (commonly known as latitude, longitude and height) as follows:

  • Geodetic Latitude (\phi): angle measured in the meridian plane through the point P between the equatorial plane of the ellipsoid and the line perpendicular to the surface of the ellipsoid at P (positive north from the equator, negative south)
  • Geodetic Longitude (\lambda): angle measured in the equitorial plane between the reference meridian and the meridian plane through P (positive east from the zero meridian)
  • Geodetic height (h): measured along the normal to the ellipsoid through P.

Note that the latitude is the angle between the equatorial plane and the line perpendicular to the surface of the ellipse at point P, not the line joining point P and the earth center (origin of the ECEF frame). Angle between the equatorial plane and the line joining point P and the earth center is called the geocentric (as opposed to geodetic) latitude. If the earth were a perfect sphere, the normal to a point would pass through the earth center and the geocentric and geodetic latitude would coincide.

Conversion between Geodetic (Ellipsoidal) and Cartesion Coordinates

Conversion from ellipsoidal to cartesian coordinates is straightforward and can be implemented in one step.

\begin{bmatrix}x \\ y \\ z \end{bmatrix} = \begin{bmatrix}(N+h)cos(\phi)cos(\lambda) \\ (N+h)cos(\phi)sin(\lambda) \\ (N(1-e^2)+h)sin(\phi) \end{bmatrix}

Going from cartesian to ellipsoidal is trickier and involves an iterative procedure that quickly converges. See Appendix 4.A of 1  for details. The corresponding Matlab code is listed in section 1.e of the appendix. The more common use case is to convert ECEF coordinates to ellipsoidal as the input and output of GPS processing algorithms are usually in ECEF coordinates.

Local Coordinate Systems

So far, we have mostly talked about “global” coordinate systems, centered at the earth center. In some applications, it is more convenient to use local coordinate systems, centered on the user’s position. These coordinate systems are called East-North-Up (ENU) systems. Given the user’s position in elliptical coordinates (\phi: latitude, \lambda: longitude), user coordinates can be easily converted from a ECEF frame to a local ENU frame using the following matrix multiplication.

R_{ENU} = \begin{bmatrix}-sin(\lambda) & cos(\lambda) & 0 \\ -sin(\phi)cos(\lambda) & -sin(\phi)sin(\lambda) & cos(\phi) \\ -cos(\phi)cos(\lambda) & cos(\phi)sin(\lambda) & sin(\phi)\end{bmatrix}R_{ECEF}

We’ll see an application of ECEF to ENU conversion when we compute the azimuth and elevation of one of the GPS satellites later.

Part 2: Using GPS to Calculate User Position

Determining the position of a user using GPS is essentially a triangulation problem. If the distance of the user from three or more satellites is known, then the 3D position can be calculated using triangulation. The basic idea in 2 dimensions is shown below.

The problem of locating the user can thus can be divided into two subproblems:

  1. Finding the distance of the user from each satellite
  2. Determining the position of this satellite in the user’s coordinate system

Let’s look at each of these problems step by step.

Step 1: Determining the Position of a Satellite

The ideal satellite orbit is an elliptical orbit and specified completely by 5 Keplerian parameters five of which determine the size and shape of the ellipse and the orientation of the orbital plane relative to the fixed stars (i.e., an inertial reference frame). The sixth parameter specifies the position of the satellite at a particular time instant of epoch. Given the six elements, the satellite position and velocity can be computed at any other epoch.

The orbit of a satellite is not an exact ellipse however as the earth is not uniform in composition and the movement of a satellite is perturbed by the gravitational forces of the sun and the moon. The resulting perturbations are small, but must be accounted for to obtain accurate position. GPS accounts for these perturbations by transmitting an expanded set of 16 orbital parameters that can be used to accurately compute the position of a satellite at a given time instant. We won’t go into the details of how these corrections are applied. The step-by-step procedure to compute the satellite position (including applying the orbital corrections) is described in table 20-IV in the GPS IS. The Matlab code that implements this procedure is shown in section 1.b of the appendix.

Rotating the Satellite Reference Frame

The user position is being calculated at time t. The GPS signal left the satellite at time t-\tau and arrives at the user \tau seconds later. While calculating the user’s position, we’ll calculate the satellite positions at t-\tau, the time of signal transmission. In time interval \tau however, the earth has rotated by \omega_E\tau where \omega_E is the rotation rate of the earth. We must rotate the satellite position vectors by the same amount, so that they are expressed in the user’s reference frame, i.e., ECEF frame at time t. This is done by multiplying the satellite position vector x^k(t-\tau) by the following rotation matrix:

x^k(t) = \begin{bmatrix} cos(\omega_E\tau) & sin(\omega_E\tau) & 0 \\ -sin(\omega_E\tau) & cos(\omega_E\tau) & 0 \\ 0 & 0 & 1 \end{bmatrix}x^k(t-\tau)

Note that this is not equivalent to calculating the position of the satellites at time t!

Step 2: Computing the Distance of the User from the Satellite

The satellite signal received by a GPS receiver bears the time stamp at which the signal was sent from the satellite. By taking the difference between its own time and the timestamp of the GPS signal and multiplying by the speed of light, the receiver calculates a rough measure of the distance between the receiver and the satellite. This measure is called the pseudorange. If the clock on the satellite and the receiver were perfectly in sync and the GPS signal traveled in a straight line at the speed of light, then this measure would be the true distance between the satellite and the receiver. This however is not the case and the difference between the satellite and user clocks as well as the delays caused by the atmosphere must be modeled to calculate an accurate position of the user.

Let’s consider first the effect of the offset between the user and satellite clocks. Let’s denote the transmission time from the satellite to the user by \tau and a common time reference (referred to as GPS Time (GPST)) by t. Denoting the receiver clock bias at time t by \delta{t}_u(t) and satellite clock bias at time t-\tau by \delta{t}_s(t-\tau), the pseudorange \rho(t) measured by the GPS receiver is given as:

(1)   \begin{equation*} \rho(t) = c[t+\delta{t}_u(t)-(t-\tau+\delta{t}_s(t))] + \epsilon_\rho(t) = c\tau + c[\delta{t}_u(t)-\delta{t}_s(t)] + \epsilon_\rho(t) \end{equation*}

If there wasn’t a layer of atmosphere between the satellite and the user, the satellite signal would travel along a straight line with the speed of light. Thus, the distance to the satellite at time t denoted by r(t, t-\tau) would be equal to c\tau.

In reality (thankfully), there is thick layer of atmosphere in the path of the satellite signal which slows it down and bends its path. The change in the path length due to bending of the signal ray is generally not significant. However the change in the speed of propagation can be significant and can result in positioning errors of several meters or more. This change in speed has two components – I_\rho arising from the propagation of the signal through the ionosphere and T_\rho arising from the propagation of the signal through the troposphere (see appendix for more information about these delays and how they can be estimated). Thus,

c\tau = r(t, t-\tau) + I_\rho + I_\rho

Combining the above two equations and dropping the reference to measurement epoch t,

(2)   \begin{equation*} \rho(t)  = r(t, t-\tau) + I_\rho + I_\rho + c[\delta{t}_u(t)-\delta{t}_s(t)] + \epsilon_\rho(t) \end{equation*}

Thus, the true range (actual distance of the user from a satellite) can be calculated from the raw pseudo ranges output by the GPS device by first subtracting the user’s estimate of the Tropo and Iono delays and then adjusting for the satellite clock offset. The Iono and the Tropo delay generally lead to a position error of ~25m and ~2m respectively. More information about how to estimate these delays are provided in the appendix. In the rest of our analysis in this blog, we’ll ignore these delays.

Estimating the Satellite Clock Bias

The satellite clock bias is important to estimate as it can result in a position error of thousands of meters. The procedure to calculate the satellite clock bias is described on page 96 of the GPS Interface Specification document. The procedure consists of evaluating a polynomial whose coefficients are provided in the GPS ephemeris message and adding a relativistic term. The polynomial provides most of the correction, with the relativistic effect contributing about 1-10 m depending on the position of the satellite. The code for calculating the satellite clock bias is given in section 1.d of the appendix.

User Clock Bias

The user clock bias is an unknown quantity, just like the user position. The user clock bias will be estimated along with the user position.

Step 3: User Position and Clock Bias Estimation

We are now ready to look at the algorithm used to compute the user position. A quick point about the notation before we look at the math. Generally, vector quantities are expressed using boldface (\boldsymbol{x}) and scalar quantities without emphasis. However in the math shown below, (mostly due to laziness), I’m skipping making this distinction. I think it is clear from the context which quantities are vectors and which ones are not.

Once the satellite clock bias has been accounted for and all available corrections have been applied, the corrected pseudorange measurement for satellite k can be written as

(3)   \begin{equation*} \rho^k = r^k + c\delta{t}_u + \epsilon^k \end{equation*}

Here r^k is the true distance between the user and the satellite and \epsilon^k denotes the combined effect of the unmodeled errors. Note that all elements in this equation are distances. We have converted the receiver clock bias \delta{t}_u into a distance by multiplying by the speed of light c. Also note that all of these calculations are carried out at time t and all position vectors are expressed in the ECEF frame at time t. We drop explicit reference to time for convenience.

Let the position of the user in the ECEF frame be denoted by x, and the position of satellite k by x^k. Then,

r^k = \lVert x^k - x\rVert

and

\rho_k = \lVert x^k - x\rVert + b + \epsilon^k

We use b = c\delta{t}_u to represent the user clock bias in units of distance

We wish to determine x and b that minimize the difference between the measured and the estimated pseudorange for each satellite. In other words, our task is to determine x and b such that

(4)   \begin{equation*}\delta{\rho}^k = \rho^k -  (\lVert x^k - x\rVert + b) \end{equation*}

is minimized for all satellites k.

This problem is solved using an iterative procedure by starting with an estimate of x and b and finding the corrections that minimize the equation above. Let our initial estimates be denoted by x_0 and b_0. We wish to find \delta{x} and \delta{b} such that the true user position x = x_0 + \delta{x} and the true clock bias b = b_0 + \delta{b} minimize

(5)   \begin{equation*} \begin{split} \delta{\rho}^k &= \rho^k -  \rho_0^k  \\ &= \lVert x^k - x\rVert + b - \lVert x^k - x_0\rVert - b_0 + \epsilon^k\\ &= \lVert x^k - x_0 - \delta{x}\rVert + b - \lVert x^k - x_0\rVert - b_0 + \epsilon^k \\ &\approx \lVert x^k - x_0\rVert -\frac{x^k -x_0}{\lVert x^k - x_0\rVert}\cdot\delta{x}+ b - \lVert x^k - x_0\rVert - b_0 + \epsilon^k \\ &= -\frac{x^k -x_0}{\lVert x^k - x_0\rVert}\cdot\delta{x} + \delta{b} + \epsilon^k \\ &= -\hat{x}_{uk}\cdot\delta{x} + \delta{b} + \epsilon^k \\ &= \begin{bmatrix}-\hat{x}_{uk} &  1\end{bmatrix}\begin{bmatrix}\delta{x} \\ \delta{b}\end{bmatrix} + \epsilon^k \end{split} \end{equation*}

Here \hat{x}_{uk} is the unit vector from the user to satellite k and we used the Taylor series expansion of the vector norm (see appendix) to express the delta pseudorange as a matrix multiplication between known and unknown quantities. Concatenating the linear equation developed above for K satellites,

\boldsymbol{\delta{\rho}} = \begin{bmatrix}\delta{\rho}^1 \\ \delta{\rho}^2 \\ \vdots \\ \delta{\rho}^K \end{bmatrix}  = \begin{bmatrix} -\hat{x}_{u1} &1 \\ -\hat{x}_{u2} &1 \\ \vdots \\ -\hat{x}_{uK} &1\end{bmatrix} \begin{bmatrix}\boldsymbol{\delta{x}} \\ \delta{b} \end{bmatrix} + \boldsymbol{\epsilon}

Setting

 \boldsymbol{G}  = \begin{bmatrix} -\hat{x}_{u1} &1 \\ -\hat{x}_{u2} &1 \\ \vdots \\ -\hat{x}_{uK} &1\end{bmatrix}

The expression above can be written more compactly as

\boldsymbol{\delta{\rho}} = \boldsymbol{G}\begin{bmatrix}\boldsymbol{\delta{x}} \\ \delta{b} \end{bmatrix} + \boldsymbol{\epsilon}

For a non-degenerate (\hat{x}_{uk} are non-coplanar) configuration of K=4 satellites, the equation above can be solved directly. In general, if the sky is not obstructed, many more satellites are visible, and a exact solution is not available. We can use linear algebra techniques to find a least squares solution:

\begin{bmatrix}\boldsymbol{\delta{\hat{x}}} \\ \delta{\hat{b}} \end{bmatrix}  = (\boldsymbol{G^T}\boldsymbol{G})^{-1}\boldsymbol{G}^T\boldsymbol{\delta{\rho}}

The hat above the estimated quantities indicates that they are a least squares solution, not an exact solution.

The complete procedure to calculate the user position and clock bias is as follows.

Note about the Matlab Code

It is worth noting that while most of the equations in the Matlab code shown in the Appendix involve evaluation of expressions where the known quantities are on the right hand side and unknown quantities on the left hand side (and thus can be evaluated in one step), there are a few that requires invocation of a solver. One example is the calculation of the eccentric anomaly (E) from the mean anomaly (M) required during the calculation of the satellite position from ephemeris parameters. The relation between the two is expressed as:

M = E-eSin(E)

where e is the eccentricity of the earth’s orbit. Typically, M is calculated first, and calculating E from M can’t be done in closed form and requires the use of a solver. Part of the code is shown below.

Experimental Setup

In this section, I’ll describe the hardware setup I used to collect raw GPS data which was then processed using the algorithms described above to calculate user position and clock bias. An interested user should be able to easily replicate my setup using cheap and commercially available GPS receivers and open source software.

The figure above shows the main steps involved in my setup. For collecting raw GPS data, special GPS units that output “timing” information consisting of raw pseudoranges and satellite ephemeris information must be used. Regular GPS units calculate the user position internally and don’t output the raw data we need to calculate the user position using the algorithm described above. The NEO-M8T and 6T chips from u-blox fit our needs. Complete hardware assemblies with the GPS unit, antennae and serial output port can be purchased from Amazon for ~40 dollars.

In order to receive and save the raw GPS signal, I use the STRSVR utility in RTKLib (RTKLIB is an open source program package for standard and precise positioning that supports all common Global Navigation Satellite Systems (GNSS) – GPS, Glonass, Galileo, Baidu etc. See: http://www.rtklib.com). STRSVR takes the custom u-blox formatted output from the u-blox receiver and converts it into RTCM standard format. RTCM (Radio Technical Commission for Maritime Services, see: http://www.rtcm.org/differential-global-navigation-satellite–dgnss–standards.html) is a standard for GNSS data and content designed to support a variety of GNSS applications in air/land/sea navigation, surveying, radio navigation/location etc.

RTCM standard provides the definition of various messages that provide specific GPS information. The information we are interested in are the raw pseudoranges and satellite ephemeris. This information is contained in messages 1002 and 1019.

Part of GPS 1019 message content

We must request the u-blox receiver to send 1002 and 1019 message information as part of the data it sends to STRSVR. This is done by using the following commands in the Cmd window:

The first command sets the update rate to 1Hz and the other two enable RAW messages. The details on how this happens are a bit murky. I tried going through parts of the UBX documentation to understand this better but soon got lost. Suffices to say that the above commands work and cause the receiver to send the data we need.

I configure the STRSVR utility to receive data from the serial port at 9600 Baud and save the data to a file in RTCM 3 format.

To collect the data, I went to the roof of my apartment building and placed the GPS receiver at a location where it could get an unobstructed view of the sky. I used the u-blox configuration software (https://www.u-blox.com/en/product/u-center-windows) to verify that a sufficient number of satellites were visible and a good position fix could be obtained. I then used STRSVR to collect about an hour of raw GPS data and saved the data to a file.

Processing Raw GPS Data

STRSVR saves the raw GPS data into binary RTCM3 format. We must decode the RTCM3 data to create Matlab data structures out of it. I considered writing my own RTCM decoder, but then found an excellent Matlab library called goGPS (http://www.gogps-project.org/about/) that provides many useful routines to read and process GPS data. I used the load_stream function that reads an RTCM formatted file and extracts the RTCM messages (line 95 in my version of goGPS). I then saved the extracted data to a .mat file to serve as input to my the position calculation algorithm that implements the equations above. The Matlab code is shown in the section 1.a of the appendix. I’m also attaching the rtcm_data since many people have asked for it. I changed the file extension from .mat to .txt due to wordpress security restrictions. Rename the file back to .mat when you download it.

Analysis of Results

In this section, I’ll analyze the positions and clock bias computed by my algorithms. We will look at the variation of the north/east/up components of the position with time, variation of the clock bias with time and introduce the concept of Dilution of Precision (DOP), a common metric used to measure the quality of GPS position estimates. Since the GPS receiver was stationary during data collection, the variation of the calculated position reflects the true performance of the algorithm used to calculate the position.

East North Up
Std-Dev 14.00 39.88 47.35

The plot of E/N/U components of the user position in a user centered ENU frame along with the position scatter plot are shown above. As expected, the variation in position is around 30m in the East and North direction and slightly higher along the Up direction. We’ll explain the reason for this during the discussion about DOP (Dilution of Precision) below.

Clock Bias Drift

The variation of the receiver clock bias against time is shown below. Note that the clock bias b used in the algorithm above is in the unit of distance. In the plot below, the bias has been converted into unit of time by dividing by the speed of light.

We can see that the clock bias is not a constant but drifts linearly with time. The amount of drift is 4.27 e-7sec/sec.

Computing the Satellite Azimuth/Elevation

As an exercise, let’s compute the azimuth and elevation of all the satellites visible at a given time instant. This will neatly tie in many of the ideas discussed above and also provide a good segue to discussing DOP.

Azimuth and elevation angles are calculated from the perspective of the user. Thus, they are expressed in an ENU frame centered about the user position.

If the position of the satellite in the user centered ENU frame is (x_s, y_s, z_s), then the azimuth (az) and elevation (el) is given by:

tan(az) = \frac{x_E}{x_N}

sin(el) = \frac{x_U}{\sqrt{x_E^2 + x_N^2 + x_U^2}}

The procedure to extract azimuth and elevation from satellite and user positions (calculated in the ECEF frame) is as follows:

  1. Calculate the position vector from the user to the satellite (in ECEF frame)
  2. Calculate the user position in Ellipsoidal coordinates (lat/lng)
  3. Rotate the position vector to the ENU frame centered about the user position
  4. Calculate azimuth/elevation using equations above)

The Matlab code is shown below:

The computed azimuth and elevations at a given epoch are shown in the table and the corresponding 3D positions in the ENU frame are shown in the plot below

Azimuth Elevation
69.68 62.89
-20.46 81.63
96.38 16.91
34.71 5.21
-145.06 12.59
-165.22 7.44
-110.29 39.03
-49.30 43.11

 

As can be seen the azimuth angles are both positive and negative, while the elevations are only positive. This makes sense as the user can’t see the satellites below the horizon.

You may have seen satellite track charts that many GPS processing software output as a visualization aid. Those charts are created by calculating the satellite positions using the procedure described here at multiple epochs.

Dilution of Precision (DOP)

DOP helps us answer the question, “how good are my position estimates”? There are two components of the errors in our position estimate. The first one is the measurement noise, which is obvious. The noisier our measurements (pseudorange, satellite position etc), higher are the expected position errors. However, there is another, less obvious component of the positioning errors. This component is the user-satellite geometry. As explained in chapter 6 of 1, it the easiest to see this in a 2D example.

A user measures his distance from a pair of signal sources S1 and S2 at known locations. If the range measurements were perfect, the user could determine his position exactly as lying on the intersection of two circles centered on S1 and S2 and radii equal to the measured ranges. The measurements however are imperfect and have a maximum uncertainty of \pm\epsilon. The figure below shows how the user-source geometry affects the amount of uncertainty in the user position.

The uncertainty in the position and clock bias is measured by the covariance of the estimated position and clock bias error. Representing the true position and clock bias by x and b, the position and clock bias estimation errors are given as:

\delta{x} = \hat{x}-x

\delta{b} = \hat{b}-b

Where \hat{x} and \hat{b} are our estimates of the user position and clock bias (calculated using the algorithm above). As shown in chapter 6 of 1, the covariance matrix of position and clock bias is given as:

Cov\begin{bmatrix}\delta{x} \\ \delta{b} \end{bmatrix} = \sigma^2(G^{T}G)^{-1}

Here \sigma is the “user range error” which captures the uncertainty in the pseudoranges and satellite positions. The G matrix is the matrix consisting of the unit vectors from the user to the satellite that we saw earlier. This matrix depends entirely on the user-satellite geometry. Thus, the covariance can be neatly factored into a product of measurement uncertainty and a function of the user-satellite geometry matrix. The elements of the G matrix are expressed in the ECEF frame. It is more convenient to rotate it into the user’s local ENU frame. Let R_L be the required rotation matrix and \tilde{G} be the resulting geometry matrix.

\tilde{G} = [G(1:N, 1:3)R_L^{T}, ones(N,1)]

Here N is the number of satellites. The rotation is applied to the position components of the geometry matrix. it can be easily shown that

Cov\begin{bmatrix}\delta{x_{ENU}} \\ \delta{b}\end{bmatrix} = \sigma^2(\tilde{G^{T}}\tilde{G})^{-1}

Setting

H =  (\tilde{G^{T}}\tilde{G})^{-1},

The East, North and Up components of DOP are defined as

\sigma_E^2 = H_{11}, \sigma_N^2 = H_{22}, \sigma_U^2 = H_{33}

We can further define a horizontal DOP (HDOP) by combining the E and N term (since the East and North directions define the horizontal plane at the user’s location) and vertical DOP (VDOP) consisting of the U term.

HDOP = \sqrt{\sigma_E^2 + \sigma_N^2}

VDOP = \sqrt{\sigma_U^2}

DOP provides a simple characterization of the user-satellite geometry. The more favourable the geometry, the lower the DOP. More favourable means the satellites are spread apart in azimuth and elevation. The lower the DOP and user range error (\sigma), the better the quality of the position estimate. To understand the relationship between DOP and satellite geometry better, lets consider two simple examples. The first example looks at the variation between DOP and the angle between the satellite vectors in a three satellite configuration shown in the figure below. The first satellite is located along the y axis and the other two are located symmetrically about the y axis, at an angle \theta from the x-axis. The G matrix and the plot of DOP against angle theta is shown below.

The minimum DOP occurs at \theta = -30^o. This makes intuitive sense as at this angle, the satellites are maximally separated from each other. Note that the DOP does not vary symmetrically with positive and negative \theta. This is because for a given \theta, the separation between the satellite vectors is higher for a negative value than for a positive value.

Next let’s look at a more realistic example of four satellites located in 3D. First satellite is located at the zenith and the others three are located 120 degrees apart in azimuth and at an elevation of \theta (same for the three satellites). We’ll look at the variation of VDOP with the elevation angle.

 

As expected, the VDOP decreases with increasing elevation angle, as the satellites get farther apart from each other. It continues to decrease as the elevation goes past the horizon (elevation = 90 degrees). Since a user located on the earth surface can’t observe satellites below the horizon (in fact, signal from satellites < 10 degrees above the horizon are too noisy and generally not used), the minimum value of VDOP is not achievable. This example demonstrates why VDOP is generally higher than HDOP.

Let’s now look at the variation of HDOP and VDOP with time calculated by our position algorithm on real data.

No surprises here. The HDOP/VDOP values are generally < 2.5, which is considered adequate and as expected, VDOP is higher than HDOP.

This concludes the post! I hope that next time you use Google maps to get directions, you’ll think about the incredible scientists, engineers and policy makers who have given us this incredible system that makes possible so many applications that we now take for granted. According to 1, the GPS constellation cost about 30 billion dollars to put in place and costs the US government about 1B annually to maintain. The valuation of Uber alone, which wouldn’t exist without GPS is upwards of 70B dollars. When you include all of the other amazing applications made possible by GPS, the public investment in GPS appears to be one of the best public investments ever made! 🙂

Appendix

1.a Code for Calculating User Position and Clock Bias

1.b Code for Calculating Satellite Position

1.c Code for Computing the Least Squares Solution for User Position and Clock Bias

1.d Code for Calculating Satellite Clock Bias

1.e Code for converting ECEF (WGS84) to Ellipsoidal Coordinates

1.f format_ephemeris3.m

2.a Taylor Series Expansion of Vector Modulus

In the derivation of the position estimation formula, we used the taylor series expansion of the modulus function. The proof is as follows:

\lVert x \rVert = (x\cdot{x})^{1/2}

First two terms of the taylor series expansion of \lVert x \rVert about x are:

\lVert x + \delta{x}\rVert = \lVert x \rVert + \frac{\partial{\lVert x \rVert }}{\partial{x}}\cdot{\delta{x}}

 \frac{\partial{\lVert x \rVert }}{\partial{x}} = \frac{1}{2}(x\cdot{x})^{-\frac{1}{2}}2x

 \frac{\partial{\lVert x \rVert }}{\partial{x}} = \frac{x}{\lVert x \rVert }

2.b Ionospheric and Tropospheric Delay

The passage of the satellite signal through the earth’s atmosphere changes the speed of the signal and bends the signal path. These effects result in the distance to the satellite calculated by multiplying the time delay between signal transmission and reception with the speed of light to be different from the true distance. Estimating the effect of the atmosphere is important to calculate accurate pseudoranges and thus accurate user position. In the treatment above, we neglected the atmospheric effects. In this section, we’ll provide some information about atmospheric effects and how they are typically modeled. The interested reader is referred to chapter 5 of 1 for details.

For the purpose of analyzing its interaction with GPS signals, the earth’s atmosphere can be split into two parts – troposphere and the ionosphere. The ionosphere extends from a height of about 50km to about 1000km above the earth is a region of ionized gases. The ionization is caused by the sun’s radiation and the state of the ionosphere is determined primarily by the intensity of solar activity. Electron density generally builds up during the day as the sun rises, peaking at around 2 PM local time and then starts declining. There can be considerable variability from day to day depending upon the solar activity. The speed of propagation of the signal in the ionosphere depends upon the number of free electrons in the path of the signal, defined as the total electron count (TEC). The number of electrons in a tube of 1m^2 cross section extending from the receiver to the satellite is given as:

TEC = \int_S^R n_e(l)dl

where n_e(l) is the electron density along the signal path and the integral is along the signal path from the satellite to the receiver. The path length through the ionospheric is shortest in the zenith direction (when the satellite is directly overhead) and longest when the satellite is close to the horizon.

In general, ionosphere delays are hard to model as solar flares and magnetic storms can large and rapidly varying fluctuations in electron densities. Two methods are primarily used to estimate ionospheric delays. The first method uses dual frequency GPS measurements and the fact that the ionospheric delay is inversely proportional to the square of the frequency. This enables us to set up a system of equations to estimate the ionosphere free pseudoranges from the noisy pseudoranges.

Receivers limited to single frequency measurements can use an empirical model called the Klobuchar model parameterized by the user’s latitude, longitude, satellite azimuth/elevation, local time and a set of parameters broadcast by the satellites. A step by step implementation of this model is detailed in the GPS Interface Specification and can be easily implemented in Matlab or C.

The troposphere extends to about 16 Km above the equator and contains roughly three quarters of the gaseous mass of the atmosphere. The troposphere is primarily composed of the dry gasses – nitrogen and oxygen, and water vapour. The primary effect of the troposphere on GPS signal propagation is to delay it slightly, causing the apparent range to the satellite to appear longer by about 2.5m, depending on the satellite elevation angle. Unlike the ionosphere, the troposphere is non-dispersive for GPS frequencies, and thus the tropospheric effect can’t be estimated by dual frequency measurements. The amount of delay caused by the troposphere depends on the pressure, temperature and humidity. Various models exist to calculate the delay from meteorological measurements. Refer to chapter 5 of 1 for details.

2.c Plotting code

Some people have asked for the plotting code – here it is.

 

1.
Misra P, Enge P. Global Positioning System. Vol 0. 0th ed. Ganga-Jamuna Press; 2011.
2.
None S. Global Positioning System Interface Specification. Global Position Systems Directorate; 2013:213.

27 Comments

      • Hi sir,
        can you suggest the python code?
        It would be helpful if you provide the code.
        Also please suggest Matlab code to plot receiver position of the satellite
        Thank you

        • I’m unable to support requests such as these. I’m not actively working on GPS anymore, and can’t provide support for custom code requests. The Matlab code provided with this post is complete and people have been able to use the Matlab code provided to recreate the results shown. Please refer to Matlab documentation and stackoverflow for custom code requests.

  1. Hi ankur,

    you did a really good work here!
    Thanks for sharing this with us.

    Using your code I was able to recreate your results you have shown here.
    So for every other user I can say that this code works as it should be.

    Currently the function “format_ephemeris3(eph_)” in line 89 of the main file is missing (Appendix 1.a) but I contacted ankur and he said that he will upload this soon. But it is also really easy to write this on your own since it just formats the input data into a struct. Using https://github.com/goGPS-Project/goGPS_MATLAB/blob/master/goGPS/input_output/RINEX/RINEX_get_nav.m I was able to write it myself.

  2. Hi Ankur,
    Thanks for sharing, it’s really useful for me to learn the GPS data processing.
    I have a question, how do you convert the RTCM file data into .mat file?
    Thank you.

  3. hi, Ankur

    Thanks for the nice blog. It is right the knowledge I am seeking for to start my own project. I would like to try to understand this blog by reproducing your work. However it will take some time before my GPS receiver kit is ready. Is it possible to also provide a sample rtcm_data.mat and the missing format_ephemeris3.m, so that I can start play with it before my kit arrive to sample my own data.

    A executable example is always better to help understand the codes.

    • The main block loop is over 1002 messages, while the inner loop is the optimization loop for position and clock bias. In other words, the position and clock bias is calculated for every 1002 message.

  4. Hi. What is the best way to plot the satellite postions like you have shown in the one figure. I can’t seem to figure out the best way to do a sky plot like that using your sample data

  5. Hi Ankur,

    Thanks for sharing. The Position Scatter Plot is not available in the make_position_scatter_plot.m script. Could you please share that part with us?

  6. Hello ankur6ue,

    You did a very good job. But I have one question, why do you use the RTCM format for calculating a gnss position?
    Is it possible to do the same thing with any conversion (u-blox raw data to RTCM) ?

  7. Hello Ankur,

    Are you aware of any other GPS devices that support RTCM 3.0 output of the required messages? I have a u-blox EVK-M8T and it does not appear to support that output. Also, NEO-M8T appears to be unavailable at present.

    For my application I cannot use the LEA-6T/M8T because I need an external antenna.

    Thanks,
    Frank Crow

  8. Really nice coding, and excellent writing! The demo data file works well for me.

    Has anyone looked at porting your code to C? I don’t have Matlab Coder (expensive!) in my home matlab license, so i may just do the conversion myself.

    Thx for posting! Don

    • I don’t think anyone has! I’d think porting the matlab code to Python would make the most sense. The code is relatively simply, so shouldn’t be too hard to do a port. If you do it, please share a link! 🙂

  9. For all who still finds this incredibly good readable article via internet search like me, pay attention that the implementation presented here is not enough to use in production. Because it is lacking the logic to process some essential pieces of data broadcast by satellites e.g. group delay parameter available in navigation frames (only ephemeris and clock bias of satellite are used in this article/code but it is not the full set of User Algorithm defined by GPS/etc ICD). So as I also ran the code from here, it results in wrong coordinates for some epochs (wrong continent) and obviously less accurate coordinates for the rest epochs. Also e.g. calculation of residual pseudoranges, typically being part of GNNS calculation, is not presented here.

    Nevertheless this article and its runnable code is super valuable for those who enters this topic. The reduced complexity of code here helped me better understand the main steps of positioning engines running in GNSS chips – thank you ankur6ue!

    • You are absolutely correct, the code is most definitely not production-ready. It is only meant to convey the basics of how GPS works

Leave a Reply

Your email address will not be published.


*