University of Alaska Fairbanks
Department of Geosciences

GEOS F493 / F693 - Geodetic Methods

Lab 3: Pseudorange Position Estimation

"I like my crust deformed."
UNAVCO bumper sticker

Note that you don't have to work on the GIseis network today (but you can, if you want to)!

Introduction

We will start making our way to calculating a position from pseudorange data as recorded by the receiver. We'll start out looking at rinex and orbit files and then you'll get to write up the inversion we talked through in lecture today. This lab is to some degree inspired by a lab Eric Calais posted on his website when he was at Purdue.

Rinex Files

First look at this ...

		     2.11           OBSERVATION DATA    G (GPS)             RINEX VERSION / TYPE
		teqc  2015Jun23     UNAVCO Archive Ops  20150720 00:13:21UTCPGM / RUN BY / DATE
		Solaris x86 5.10|AMD64|cc SC5.8 -xarch=amd64|=+|=+          COMMENT
		BIT 2 OF LLI FLAGS DATA COLLECTED UNDER A/S CONDITION       COMMENT
		ABBZ                                                        MARKER NAME
		                                                            MARKER NUMBER
		Kyle, Philip        NMT                                     OBSERVER / AGENCY
		5115K74978          TRIMBLE NETR9       4.85                REC # / TYPE / VERS
		12561303            TRM41249.00     NONE                    ANT # / TYPE
		 -1353856.8945   314830.6876 -6205742.1059                  APPROX POSITION XYZ
		        0.0000        0.0000        0.0000                  ANTENNA: DELTA H/E/N
		     1     1                                                WAVELENGTH FACT L1/2
		     7    L1    L2    C1    P2    P1    S1    S2            # / TYPES OF OBSERV
		    15.0000                                                 INTERVAL
		    17                                                      LEAP SECONDS
		input file: abbz201507190000a.tgd                           COMMENT
		RINEX file created by UNAVCO GPS Archive.                   COMMENT
		For more information contact archive@unavco.org             COMMENT
		Monument ID: 14840                                          COMMENT
		UNAVCO 4-char name:   ABBZ                                  COMMENT
		4-char name from Log or data file: ABBZ                     COMMENT
		Monument location: -77.4569 166.9089 1734.2656              COMMENT
		Visit ID: 112562                                            COMMENT
		End of DB comments                                          COMMENT
		 SNR is mapped to RINEX snr flag value [0-9]                COMMENT
		  L1 & L2: min(max(int(snr_dBHz/6), 0), 9)                  COMMENT
		  2015     7    19     0     0    0.0000000     GPS         TIME OF FIRST OBS
		                                                            END OF HEADER
		 15  7 19  0  0  0.0000000  0 10G27G22G14G28G24G04G11G15G18G19
		 125212588.537 6  97568267.12843  23827164.883    23827169.719
		        40.800          22.600
		 109422791.431 8  85264519.77646  20822467.477    20822467.918
		        51.200          40.200
		 127128337.81515                  24191714.477
		        35.700
		 115450099.390 7  89961166.65045  21969433.359    21969435.793
		        47.600          30.500
		 120380962.929 7  93803351.19444  22907727.289    22907732.887
		        43.200          28.900
		 122440368.761 8  95408123.50345  23299637.109    23299640.129
		        48.000          30.200
		 126295265.169 7  98411944.66444  24033202.547    24033204.746
		        43.500          27.000
		 118875481.704 7  92630241.51045  22621249.406    22621251.781
		        47.000          32.400
		 109981293.845 8  85699702.35446  20928744.125    20928745.914
		        51.000          36.800
		 117698935.043 8  91713448.20746  22397361.672    22397362.699
		        49.400          36.200
		 15  7 19  0  0 15.0000000  0 10G27G22G14G28G24G04G11G15G18G19

This is a rinex file (receiver independent exchange format). There's a lot going on here, but basically this separates into a header section until END OF HEADER and then we the observation for 1 epoch 15 7 19 0 0 0.0000000, which is 2015-July-19 at midnight. The line 10G27G22G14G28G24G04G11G15G18G19 gives you the satellites that were in view (G refers to GPS), so here we have PRN 27, 22, 14... The first number (10, here) gives the number of satellites in view at that epoch.

The next 20 lines are the observations for each satellite given in the order as defined in the 10G27... string:

	 125212588.537 6  97568267.12843  23827164.883    23827169.719
	        40.800          22.600

... is one observation for PRN10. These 2 lines are composed of the individual observables. The order and character of the observables is defined in the header line # / TYPES OF OBSERV. If you go and look you'll find that this line starts with a 7, meaning we have 7 observables, which is followed by 7 strings. Here L1, L2 for the phase observations on L1, L2; C1 for the pseudorange inferred from C/A on L1; and P1,P2 which are pseudoranges inferred from the P-code on L1, L2; and lastly S1,S2, which are indicate the signal strengths on L1, L2.

In this lab we will use the pseudorange on L1, so for each satellite you will always be looking at the 3rd column in the first line of its 2 line observation!

The rinex file you'll be using is from station Abbot Peak on Ross Island, Antarctica: abbz2000.15o

Orbit Files

In class I talked a bit about satellite orbit files. The IGS provides precise, post-processed orbits about 2 weeks after the fact.

		#cP2015  7 19  0  0  0.00000000      96 ORBIT IGb08 HLM  IGS
		## 1854      0.00000000   900.00000000 57222 0.0000000000000
		+   31   G01G02G03G04G05G06G07G09G10G11G12G13G14G15G16G17G18
		+        G19G20G21G22G23G24G25G26G27G28G29G30G31G32  0  0  0
		+          0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
		+          0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
		+          0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
		++         2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
		++         2  2  2  2  2  2  2  2  2  2  2  2  2  2  0  0  0
		++         0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
		++         0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
		++         0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
		%c G  cc GPS ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc
		%c cc cc ccc ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc
		%f  1.2500000  1.025000000  0.00000000000  0.000000000000000
		%f  0.0000000  0.000000000  0.00000000000  0.000000000000000
		%i    0    0    0    0      0      0      0      0         0
		%i    0    0    0    0      0      0      0      0         0
		/* FINAL ORBIT COMBINATION FROM WEIGHTED AVERAGE OF:        
		/* cod emr esa gfz grg jpl mit ngs sio                      
		/* REFERENCED TO IGS TIME (IGST) AND TO WEIGHTED MEAN POLE: 
		/* PCV:IGS08_1854 OL/AL:FES2004  NONE     Y  ORB:CMB CLK:CMB
		*  2015  7 19  0  0  0.00000000
		PG01  18829.174425  13507.523320 -13215.249692     -2.670727  8  8  9 109       
		PG02   3708.172047 -15380.790615  21695.248921    579.065558  9  3  8  88       
		PG03  13652.481176  20791.848200   9285.060119      8.499972  8  8  6  80       
		PG04  10011.699639  16084.389751 -19050.244266     -2.469151  7  9  6  87       
		PG05   -939.523866 -24616.626155   9624.904794   -216.443730  7  4  8  92       
		PG06  17807.401446  -8812.348945  17624.772591     19.301639  8  6  7  97       
		PG07  26053.972176   6163.904952     -3.860242    475.220308  8  7 10 105       
		PG09  17947.907290    578.189652  19572.892136     -7.756942  7  6  8  81       
		PG10  12231.245947 -10800.089940  20489.425535   -183.089361  6  7  7  70       
		PG11  14697.108223  12954.318873 -18522.373638   -592.243492  9 10 12 116       
		PG12 -11059.392698 -23259.388628   6832.355549    311.126639 10  8  9 108       
		PG13   9094.018662 -20832.138975 -13973.496637   -133.344970  9  8  7 100       
		PG14 -21395.252372  13872.727013  -6882.375231     47.983976 10  6  9  95       
		PG15  -3750.820337 -17956.885654 -19189.766833   -261.224294  9  9  7 131       
		PG16   -617.107855  24330.865403  10088.693768   -107.080805  8  8  8 121       
		PG17  20671.227040 -15283.721033  -5935.131157   -180.682624  7  9  9 115       
		PG18 -17051.363484  -4161.527113 -19503.448336    418.339886  8  9  7 107       
		PG19    -75.524186  18164.222089 -19415.548695   -516.232497  6  7  9  96       
		PG20 -13621.803708 -22088.770198  -5575.510729    339.555062 10 10  7 115       
		PG21 -26365.185049  -3888.470674  -1564.859572   -464.856887 12  8  8 130       
		PG22 -12933.794843  10205.759827 -20569.279215    372.435926  8  5  7 117       
		PG23  10235.841275  11819.396635  21604.947116   -121.830762  8  6  6  99       
		PG24 -16548.513468 -14236.584786 -15267.729561     -1.557331  9  9  9 110       
		PG25 -16057.516800 -13045.207716  16646.609876    -14.652893  8 10  7  78       
		PG26  -7070.428060  18357.038062  17834.006702    -10.811951  9  6  8  90       
		PG27  -9678.067077  22140.496825 -10913.108494      5.315444  6  6  6 129       
		PG28  13506.960104  -4894.537400 -21719.611795    454.532750  9 10  8 120       
		PG29 -16439.453511  -3356.044382  20594.646988    623.146711  9  9 10 101       
		PG30  24254.540458   -996.882647 -10878.584342    -11.110350 10  9  9 142       
		PG31 -16263.550371   9327.107306  19072.314740    309.359520  9  7  6 126       
		PG32   4490.858537  26368.590945  -2451.536040    -96.311874  9  6  9 101       
		*  2015  7 19  0 15  0.00000000

Again, this file has a lot of header information. It gives satellite positions in 15 minute intervals. Each line starting with a PG** defines the X,Y,Z position of a satellite rotated into ECEF (here ITRF2008b) in kilometers. In the fifth column, we have the satellite clock bias in microseconds and then the X,Y,Z standard deviations in mm/sec, and a clock std-dev in psec/sec. The full file for that day is igs18540.sp3

Your Task!

I want you to take the C/A pseudoranges and satellite positions for the epoch given above and compute the receiver position in ECEF reference frame. Once you've done that you will use the script from the last lab and convert this position into the WGS84 reference frame

Preparation (this requires some manual labor):

Assemble an XYZ apriori position, use the one given in the rinex header, set the apriori receiver clock bias to 0. You should have a column vector with 4 rows.
Assemble a vector of the n pseudorange values (C1) for the n satellites visible at epoch 2015 7 19 0 0 0.00000000 (first column with big numbers that start with 2).
Assemble a n x 4 matrix (n rows, 4 columns) of satellite positions for epoch 2015 7 19 0 0 0.00000000, which is the one listed above. These positions will be X,Y,Z in km and the clock bias in microseconds, i.e. straight from the sp3 file above (or really, from the table). Make sure the satellite rows are ordered as in the rinex file: G27G22G14G28G24G04G11G15G18G19.

Now, write a function, solve_pseudo_r, that takes the 3 vectors you assembled above as input and returns the updated/adjusted position (as a column vector with values (Xa, Ya, Za, Ta), which are the adjusted position in m, and clock bias in nanoseconds). The signature of the function should be:

Matlab:

			function updated_postion_vector = solve_pseudo_r(apriori_pos, sat_pos, pseudoranges)
			   ...

Python:

			def solve_pseudo_r(apriori_pos, sat_pos, pseudoranges):
			   ...
			   return updated_postion_vector;

The algorithm is as follows:

convert the satellite positions to meters
convert the satellite clock bias to seconds
correct pseudoranges for satellite clock bias by adding it to the pseudoranges (you probably want to convert the clock bias into a distance first by multiplying by the speed of light,c , in m/s!)
compute the modeled observables rho_o^k = sqrt((X^k-Xo)^2 + (Y^k-Yo)^2 + (Z^k-Zo)^2), i.e. the geometric distance between satellite k and a priori receiver position.
compute the observation vector deltaRHO, which is the difference of pseudoranges and modeled observations above.
compute partial derivatives of Taylor Series (dx = -(X-Xo)/rho_o; ..., dt=c*1e-9). The multiplication of the speed of light will give us the receiver bias in nanoseconds and circumvent numerical instabilities
assemble the n x 4 design matrix A = [dx dy dz dt] (note that dx,... are column vectors!)
solve the inverse problem using the pseudoinverse: Matlab: X = pinv(A) * deltaRHO (pinv computes generalized inverse via SVD); Python: X = numpy.linalg.pinv(A)*deltaRHO
Compute the new position by correcting Xo, Yo, Zo, To for the results of the inversion

Now write complete the initial script you started to call this function and iterate 4-5 times over the results to find a stable position solution. Convert the final position to geodetic coordinates using last week's results.

Solution

Your code should approximate these ECEF and WGS84 coordinates for a GPS station on Ross Island, Antarctica (5 iterations):

	ECEF		WGS84
X (m)	-1353875.822	lat (deg)	-77.453051
Y (m)	314824.972	lon (deg)	166.909306
Z (m)	-6205811.52	height (m)	1805.735
T (nsec)	-795.52

Real Test

Set your apriori position to [0 0 0 0] and your iterations to, maybe, 10. Collect the results from each run and plot them in 4 subplots showing the convergence of each of the 4 parameters we are solving for. Astonishing, ey?

Deliverables: (submit via canvas!)

the function solve_pseudo_r
script that calls solve_pseudo_r, contains the vectors and matrices that you assembled and iterates over the position solution using solve_pseudo_r
plot of the apriori position [0 0 0 0] convergence

rg <at> nmt <dot> edu | Last modified: September 18 2019 23:12.

University of Alaska Fairbanks Department of Geosciences GEOS F493 / F693 - Geodetic Methods