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

LAB 7

From GPS data to a position: processing pseudorange data

The pseudorange measurements ^jR_i(t) can be modeled as:

(1)

t = time of epoch

^jR_i = pseudorange measurement

^jr_i = satellite-receiver geometric distance

c = speed of light

^jd = satellite clock bias

d_i = 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:

(2)

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

(3)

with [^jX, ^jY, ^jZ] = satellite position, [X_i, Y_i, Z_i] = receiver position in an ECEF coordinate system.

Our mission, if we accept it, is to solve for [X_i, Y_i, Z_i, d_i], assuming that we know [^jX, ^jY, ^jZ, ^jd]. A major problem here is that the unknowns [X_i, Y_i, Z_i] are not linearly related to the observables…

Assuming that we now the approximate coordinates of the receiver [X_o, Y_o, Z_o], one can write that the actual coordinates equal the approximate coordinates plus a slight adjustment:

(4)

DX_i, DY_i, DZ_i are our new unknowns. We can now write:

(5)

Since we know the approximate point [X_o, Y_o, Z_o], we can now expand f(X_o+DX_i, Y_o+DY_i , Z_o+DZ_i) using a Taylor’s series with respect to that point:

(6)

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

(7)

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

(8)

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

(9)

We now have an equation that is linear with respect to the unknowns DX_i, DY_i, DZ_i.

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

(10)

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

(11)

We can simplify the notation by assigning:

(12)

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

(13)

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

(14)

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:

(15)

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:

(16)

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

(17)

P is the weight matrix, defined by:

(18)

s_o² = a priori variance

S_L = covariance matrix of the observations.

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

(19)

In the case of pseudoranges, the observations are independant and have equal variance s_o². Therefore S_L is the diagonal matrix:

(20)

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

(21)

Now that DX_i, DY_i, DZ_i are found, the antenna coordinates [X_i, Y_i, Z_i] are obtained using eq.(4).

The associated covariance matrix of the unknowns S_x is:

(22)

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

(23)

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

(24)

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

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

(25)

Assignment:

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

• Pseudorange data at epoch 00:15:00.0 from rinex observation file sjdv0100.02o
• Satellite position and clock bias from orbit file igs1484.sp3 (satellite positions in kilometers in ECEF frame, clock biases in microseconds). Note that the satellite clock bias should be added to the pseudorange.

Compare solutions using C1, P1, and P2

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

Xo=4433470.0

Yo=362670.0

Zo=4556210.0

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 ^jr_i (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:

	C1	P1	P2
DX	-37.448	-36.926	-36.136
DY	52.132	52.521	53.487
DZ	-60.883	-60.628	-59.357
DT (nsec)	187.16	186.30	172.53
Xa	4433432.552	4433433.074	4433433.864
Ya	362722.131	362722.521	362723.487
Za	4556149.117	4556149.372	4556150.643
GDOP	5.2	5.2	5.2