EAS-591T – Space Geodetic Measurements of Active Crustal Motions



From GPS data to a position: processing pseudorange data


The pseudorange measurements jRi(t) can be modeled as:


t = time of epoch

jRi = pseudorange measurement

jri = satellite-receiver geometric distance

c = speed of light

jd = satellite clock bias

di = receiver clock bias

DI = ionospheric propagation error

DT = tropospheric propagation error

MP = multipath

e = receiver noise

(ranges in meters, time in seconds)


Neglecting the propagation, multipath, and receiver errors, eq.(1) becomes:


The geometric distance between satellite j and receiver i is given by:


with [jX, jY, jZ] = satellite position, [Xi, Yi, Zi] = receiver position in an ECEF coordinate system.


Our mission, if we accept it, is to solve for [Xi, Yi, Zi, di], assuming that we know [jX, jY, jZ, jd]. A major problem here is that the unknowns [Xi, Yi, Zi] are not linearly related to the observables


Assuming that we now the approximate coordinates of the receiver [Xo, Yo, Zo], one can write that the actual coordinates equal the approximate coordinates plus a slight adjustment:


DXi, DYi, DZi are our new unknowns. We can now write:


Since we know the approximate point [Xo, Yo, Zo], we can now expand f(Xo+DXi, Yo+DYi , Zo+DZi) using a Taylor’s series with respect to that point:


We intentionally truncate the Taylor’s expansion after the linear terms. Recall from eq.(3) that:



The partial derivatives in eq.(6) are therefore given by:



We can now substitute eq.(8) into eq.(6):



We now have an equation that is linear with respect to the unknowns DXi, DYi, DZi.


Now let us go back to our pseudorange measurements jRi(t) and rewrite eq.(2):


We can rearrange eq.(10) by separating the known and unknown terms of each side (recall that the satellite clock correction  jd(t) is provided in the navigation message):



We can simplify the notation by assigning:



Let us assume that we have 4 satellites visible simultaneously. We use eq.(11) and write it for the 4 satellites::



Tired of carrying along all these terms, subscripts, and superscripts? Me too. Let us introduce:


L = vector of n observations. Must have at least 4 elements (i.e. 4 satellites), but in reality will have from 4 to 12 elements depending on the satellite constellation geometry.

X = vector of u unknowns. Four elements in our case.

A = matrix of linear functions of the unknowns (= design matrix), n rows by u columns.


Now we can write our problem (eq.13) in a matrix-vector form:


In general, n > u, leading to an overdetermined system. Because actual data contain observational errors and noise, this system is in non-consistent. In order to make it consistent, one must account for a noise vector r. Eq.(15) becomes:


The “noise vector” r represents residuals, i.e. observations (L) minus model (AX). The least squares solution to eq.(16) is:


P is the weight matrix, defined by:


so2 = a priori variance

SL = covariance matrix of the observations.


The law of covariance propagation gives the covariance matrix of the unknowns SX:


In the case of pseudoranges, the observations are independant and have equal variance so2. Therefore SL is the diagonal matrix:


Assuming that the weight matrix is I, eq.(17) can be simplified to:


Now that DXi, DYi, DZi are found, the antenna coordinates [Xi, Yi, Zi] are obtained using eq.(4).


The associated covariance matrix of the unknowns Sx is:



We can transform Sx from an ECEF frame to a local topocentric frame using the law of variance propagation (disregarding the time-correlated components of Sx:



where R is the rotation matrix (cf. lab 1):



with j = geodetic latitude of the site, l = geodetic longitude of the site.


The DOP factors (Dilution Of Precision) are given by:





Write a MATLAB program to compute the position and clock bias of a GPS receiver and the GDOP using:

Compare solutions using C1, P1, and P2


The a priori position of the receiver in ECEF frame (in meters) is:





Possible program structure:

  1. Define constants (c) and a priori GPS receiver position and clock bias;
  2. Read satellite positions and clock biases, convert to meters and seconds;
  3. Read pseudorange data;
  4. Correct pseudorange for satellite clock bias. Trick: satellite clock biases must be added to the measured pseudoranges;
  5. Compute modeled observables jri (eq.(7));
  6. Compute observation vector L (eq.(12)). Trick: discard satellite 1 because of its large clock bias (it must be a flag);
  7. Compute partial derivatives of Taylor’s series (eq.(8));
  8. Form design matrix (eq.(14)). Trick: multiply the c by 1e-9 in the design matrix in order to avoid numerical instabilities in the inversion. The receiver clock bias will be output in nanoseconds.
  9. Invert the design matrix (using inv) and find the vector of unknowns (eq.(21)), or solve the least squares problem directly using pinv or lscov);
  10. Compute adjusted parameters (eq.(4));
  11. Compute covariance in ECEF frame (eq.(22));
  12. Compute site ellipsiodal coordinates. Trick: use the xyz2wgs.m routine that you wrote for lab 1 or get it from the class web site;
  13. Form the ECEF to topocentric rotation matrix (eq.(24));
  14. Compute covariance in topocentric frame (eq.(23));
  15. Compute DOPs (eq.(25)).
  16. Go to bed.



I find:

















DT (nsec)