The Science of Where – Code Based Positioning with GPS and Galileo¶

Linear model GNSS SPS

$~\color{blue}{Input:}$

  • $\color{blue}{\text{Pseudoranges}} ~\color{blue}{\rho_j}$

Navigation message:

  • $\color{blue}{\text{Satellite position}}$ when transmitting signal:$ ~\color{blue}{s_j =(x_{sat}^j, y_{sat}^j, z_{sat}^j)}$
  • Satellite clock $\color{blue}{\text{offset}}:\color{blue}{dt^j}$ (j= 1,2,3...n) (n>=4)

$~\color{red}{Unknowns:}$

  • $\color{red}{\text{Receiver position}}$:$~\color{red}{(r =x_{rec}, y_{rec}, z_{rec})}$
  • $\color{red}{\text{Receiver clock offset}}$:$~\color{red}{\delta_{rec}}$

Pseudo Code of Standard Positioning Service (SPS) solution¶

Input: Raw pseudoranges from Rinex, satellite ephemeris information from broadscast ephemeris.

Initialise user clock bias = 0, user position at earth center ECEF = 0,0,0, alternativly use previous user position or 'Bancrofts method'

Initialise estimated satellite clock bias and pseudoranges for each satellite:

  1. Calculate satellite clock bias at receiver time tag
  2. Calculate satellite signal emission time using measured pseudorange, receiver time tag, satellite clock bias and speed of light.
  3. Estimate ionospherec delay using Ionosphere-free Combination for dual frequency receivers.
  4. Calculate satellite position at satellite signal emission time in two steps:
    1. Perifocal coordinate system
    2. Earth-fixed coordinate system (ECEF)
  5. Transform satellite coordinates from the system tied to the earth at signal emission time to the system tied to the earth at signal reception time. Correct the ECEF coordinates of satellite position for earth rotation during signal propagation time.

Do until change in user user position and receiver clock < threshold

For each satellite:

  1. Update the modeled psudorange with previous estimate of user position
  2. Calculate the prefit residual which is the difference between measured and modeled pseudorange, update troposheric slant correction from new receiver position.
  3. Form the linearized GPS measurement matrix G linearising around 'a priori' receiver position $(X_0, Y_0, Z_0)$
  4. Use the Least square method or Kalman filter method to solve y = Gx for correction in receiver position and receiver clock bias.
  5. Update user position.

End For

End Do

Output: User position in ECEF and user clock bias.

Sources of GNSS pseudorange measurement errors¶

GNSS pseudorange measurement errors

$PR_{rec}^{sat} = \rho + c (\delta t_{rec} - \delta t^{sat}) + \sum \delta_k + \epsilon$

$\sum \delta_k = Iono_{rec}^{sat} + Trop_{rec}^{sat} + TGD^{sat} $

Ionospheric Delay¶

The propagation speed of GNSS electromagnetic signals in the ionosphere depends on its electron density and frequncy of signal. Codes and phases propagate at different speeds in the ionosphere, resulting in a delay in the code and an advance in the carrier measurements. Phase measurements are advanced on crossing the ionosphere, that is a negative delay, while the code measurements undergo a positive delay. The difference between the measured range (with frequency f signal) and the euclidean distance between the satellite and receiver is given by:

Phase ionospheric refraction $\Delta_{phase, f}^{iono} =-\frac{\text{40.3}}{{f}^2}\int_{}^{} N_e dl$

Code ionospheric refraction $\Delta_{group, f}^{iono} =+\frac{\text{40.3}}{{f}^2}\int_{}^{} N_e dl$

Where $N_e$ is electron density.

Ionosphere-Free Combination for Dual-Frequency Receivers¶

The ionospheric effects on code $R$ and carrier phase $\Phi$ measurements depend (99.9%) on the inverse of squared signal frequency f. Therefore, the dual-frequency receivers can eliminate these effect through a linear combination of code or carrier measurements:

Ionosphere-Free Combination:

Phase $ \Phi_{_{\mbox{iono-free}}} = \frac{f_1^2\;\Phi_{_{L_1}}-f_2^2\;\Phi_{_{L_2}}}{f_1^2-f_2^2}$

Pseudorange $ R_{_{\mbox{iono-free}}}=\frac{f_1^2\;R_{_{P_1}}-f_2^2\;R_{_{P_2}}}{f_1^2-f_2^2}$

Timing Group Delay or Total Group Delay (TGD) of the associated codes cancel in this combination and are not needed (since the satellite clocks are referred to the ionospheric-free combination of such codes).

In [74]:
# GPS frequencies
F1 = 1575.42e6
F2 = 1227.60e6
F5 = 1176.60e6
# Galileo frequncies
E1 = 1575.420
E6 = 1278.750
E5 = 1191.795
E5a = 1176.450
E5b = 1207.140
C1C = 6
C2L = 9
C5Q =10
GPSTime = "2022 06 15 13 59 50.0000000"
SAT = "G01"
P1 = 22009709.740
P2 = 22009701.580
P5 = 22009698.680
ionosfreeP1P2=(F1**2*P1-F2**2*P2)/(F1**2-F2**2)
print('Ionosfree P1 P2: '+ str(ionosfreeP1P2))

Ionosfree P1 P2: 22009722.35313868
In [75]:
ionosfreeP1P5=(F1**2*P1-F5**2*P5)/(F1**2-F5**2)
print('Ionosfree P1 P5: '+ str(ionosfreeP1P5))

Ionosfree P1 P5: 22009723.690324184
In [76]:
TEC_P1P2 =1/40.3/(1/F1**2-1/F2**2)*(P1-P2)
print('TEC P1 P2: '+str(TEC_P1P2))

TEC P1 P2: -7.768028923398451e+17
In [77]:
P1_TECcoor =  P1-(40.3/F1**2)* TEC_P1P2
print('P1 TEC coorection: '+str(P1_TECcoor))

P1 TEC coorection: 22009722.353138685
In [78]:
TEC_P1P5 =1/40.3/(1/F1**2-1/F5**2)*(P1-P5)
print('TEC P1 P5: '+str(TEC_P1P5))

TEC P1 P5: -8.591558727486689e+17
In [79]:
P2_TECcoor =  P2-(40.3/F2**2)* TEC_P1P2
print('P2_TEC coorection: '+str(P2_TECcoor))

P2_TEC coorection: 22009722.353138685
In [80]:
P5_TECcoor =  P5-(40.3/F5**2)* TEC_P1P5
print('P5 TEC coorection: '+str(P5_TECcoor))

P5 TEC coorection: 22009723.690324184

Elliptic satellite orbit and perifocal coordinate system¶

Elliptic satellite orbit

In [81]:
# External dependencies and constants
import math
import datetime
C = 2.99792458E+08 #Speed of light in vacuum 
WE = 7.2921151467E-5  #Earth’s angular velocity rad/sec
A_WGS84 = 6378137        #WGS84 ellipsoid semi major axis
B_WGS84 = 6356752.31425  #WGS84 ellipsoid semi minor axis

Orbital elements¶

Orbital elements

In [82]:
class Eph():
    """ class to define GPS/GAL/QZS/CMP ephemeris """                                       #Rinex NAV RECORD    Cell
    sat = 0     # Satellite system (G), sat number (PRN)                                    #SV / EPOCH / SV CLK [0]
    iode = 0    # Issue of data, ephemeris                                                  #BROADCAST ORBIT - 1 [10]
    iodc = 0.0  # IODC issue of data, clock                                                 #BROADCAST ORBIT - 6 [33]
    af0 = 0.0   # SV clock correction (seconds)                                             #SV / EPOCH / SV CLK [7]
    af1 = 0.0   # SV clock clock drift correction coefficient (sec/sec)                     #SV / EPOCH / SV CLK [8]
    af2 = 0.0   # SV clock clock drift rate (sec/sec2)                                      #SV / EPOCH / SV CLK [9]
    toc = 0     # Epoch: Toc - Time of clock (GPS) year, month, day, hour, minute, second   #SV / EPOCH / SV CLK [1-6]
    toe = 0     # Toe Time of Ephemeris (sec of GPS week)                                   #BROADCAST ORBIT - 3 [18]
    tot = 0     # Transmission time of message                                              #BROADCAST ORBIT - 7 [34]
    week = 0    # GPS Week # (to go with TOE) continuous number, not mod(1024)!             #BROADCAST ORBIT - 5 [28]
    crs = 0.0   # radius sine harmonic correction (meters)                                  #BROADCAST ORBIT - 1 [11]
    crc = 0.0   # radius cosine harmonic correction (meters)                                #BROADCAST ORBIT - 4 [23]
    cus = 0.0   # argument of latitude correction sine harmonic correction (radians)        #BROADCAST ORBIT - 2 [16]
    cuc = 0.0   # argument of latitude cosine harmonic correction (radians)                 #BROADCAST ORBIT - 2 [14]
    cis = 0.0   # inclination sine harmonic correction (radians)                            #BROADCAST ORBIT - 3 [21]
    cic = 0.0   # inclination cosine harmonic correction  (radians)                         #BROADCAST ORBIT - 2 [19]
    e = 0.0     # eccentrisity                                                              #BROADCAST ORBIT - 3 [15]
    i0 = 0.0    # inclination angle at reference time (radians)                             #BROADCAST ORBIT - 4 [22]
    sqrtA = 0.0 # sqrt of semi-major axis (sqrt(m))                                         #BROADCAST ORBIT - 2 [17]
    deln = 0.0  # Delta n, mean motion difference from computed value                       #BROADCAST ORBIT - 1 [12]
    M0 = 0.0    # M0 mean anomaly at reference time(radians)                                #BROADCAST ORBIT - 1 [13]
    OMG0 = 0.0  # OMEGA0, longitude of ascending node (radians)                             #BROADCAST ORBIT - 3 [20]
    OMGd = 0.0  # OMEGA DOT, rate of right ascension (radians/sec)                          #BROADCAST ORBIT - 3 [25]
    omg = 0.0   # Argument of perigee (radians)                                             #BROADCAST ORBIT - 4 [24]
    idot = 0.0  # Rate of inclination angle T (radians/sec)                                 #BROADCAST ORBIT - 5 [26]
    tgd = 0.0   # Group delay differential (only to be used with single freq. recvr.) (sec) #BROADCAST ORBIT - 6 [32]
    sva = 0     # SV accuracy (meters)                                                      #BROADCAST ORBIT - 6 [30]
    health = 0  # SV health, URA Index                                                      #BROADCAST ORBIT - 6 [31]
    fit = 0     # Fit Interval in hours see section 6.11.(BNK)                              #BROADCAST ORBIT - 7 
    toes = 0

    def __init__(self, sat=0):
        self.sat = sat

Example data¶

Example measurements TOC 2022-6-15 13:59:50 and 2022-6-15 14:00:30

In [83]:
XYZapprox = [1962040.2281, 844038.2429, 5989768.7110]          # APPROX POSITION XYZ from Rinex

# prn, time  og measurement,    C1C,          C2W/C2L
satobs1 = [['G01', 2022, 6, 15, 13, 59, 50, 22009709.740, 22009701.580],
         ['G08', 2022, 6, 15, 13, 59, 50, 21982005.060, 21981997.620],
         ['G10', 2022, 6, 15, 13, 59, 50, 21597421.060, 21597412.480],
         ['G14', 2022, 6, 15, 13, 59, 50, 22853592.900, 22853583.060],
         ['G21', 2022, 6, 15, 13, 59, 50, 21228529.140, 21228516.720],
         ['G22', 2022, 6, 15, 13, 59, 50, 24111518.280, 24111510.480],
         ['G24', 2022, 6, 15, 13, 59, 50, 23071483.120, 23071476.140],
         ['G27', 2022, 6, 15, 13, 59, 50, 24314409.040, 24314403.480]]

satobs = [['G01',2022,6,15,14,0,30,21985760.860,21985752.700],
         ['G08',2022,6,15,14,0,30,22000879.460,22000872.020],
         ['G10',2022,6,15,14,0,30,21611138.380,21611129.860],
         ['G14',2022,6,15,14,0,30,22855354.620,22855344.820],
         ['G21',2022,6,15,14,0,30,21222574.420,21222561.980],
         ['G22',2022,6,15,14,0,30,24085793.080,24085785.080],
         ['G24',2022,6,15,14,0,30,23063882.340,23063875.380],
         ['G27',2022,6,15,14,0,30,24344764.260,24344756.740]]
        
ephs1 = [['G01', 2022, 6, 15, 14, 0, 0, .340794678777E-03, -.784439180279E-11, .000000000000E+00,
         .520000000000E+02, -.842500000000E+02, .416731644260E-08, .264235247830E-01,
         -.429712235928E-05, .120090576820E-01, .735744833946E-06, .515366528702E+04,
         .309600000000E+06, -.279396772385E-07, .224388442451E+01, .875443220139E-07,
         .987856820887E+00, .374968750000E+03, .921269443462E+00, -.838463496799E-08,
         -.380730144653E-09, .000000000000E+00, .221400000000E+04, .000000000000E+00,
         .000000000000E+00, .000000000000E+00, .512227416039E-08, .520000000000E+02,
         .309335000000E+06, .000000000000E+00, .000000000000E+00, .000000000000E+00],
        ['G08', 2022, 6, 15, 14, 0, 0, -.722706317902E-04, -.159161572810E-11, .000000000000E+00,
         .820000000000E+02, .140875000000E+03, .459733435457E-08, .219640037703E+01,
         .740587711334E-05, .744616868906E-02, .344030559063E-05, .515378232765E+04,
         .309600000000E+06, -.184401869774E-06, .116075602777E+01, .195577740669E-06,
         .963246002285E+00, .315625000000E+03, .173596929569E+00, -.824212903204E-08,
         .228938107620E-09, .000000000000E+00, .221400000000E+04, .000000000000E+00,
         .000000000000E+00, .000000000000E+00, .512227416039E-08, .820000000000E+02,
         .309335000000E+06, .000000000000E+00, .000000000000E+00, .000000000000E+00],
        ['G10', 2022, 6, 15, 14, 0, 0, -.455882865936E-03, -.147792889038E-10, .000000000000E+00,
         .580000000000E+02, -.136031250000E+03, .365872382910E-08, -.198612952341E+01,
         -.702962279320E-05, .765492708888E-02, .103227794170E-04, .515368768883E+04,
         .309600000000E+06, -.279396772385E-07, -.300850082132E+01, -.137835741043E-06,
         .974342578239E+00, .189125000000E+03, -.249896124194E+01, -.737780731529E-08,
         .143934566881E-09, .000000000000E+00, .221400000000E+04, .000000000000E+00,
         .000000000000E+00, .000000000000E+00, .232830643654E-08, .580000000000E+02,
         .309336000000E+06, .000000000000E+00, .000000000000E+00, .000000000000E+00],
        ['G14', 2022, 6, 15, 14, 0, 0, -.111524946988E-03, -.682121026330E-12, .000000000000E+00,
         .150000000000E+02, .985625000000E+02, .447661504041E-08, -.130199759818E+01,
         .498257577419E-05, .199141632766E-02, .110100954771E-04, .515366280174E+04,
         .309600000000E+06, .707805156708E-07, .173366849128E+00, -.931322574615E-08,
         .952816756894E+00, .163218750000E+03, -.313134934801E+01, -.778175271266E-08,
         -.286440502825E-09, .000000000000E+00, .221400000000E+04, .000000000000E+00,
         .000000000000E+00, .000000000000E+00, -.791624188423E-08, .783000000000E+03,
         .309337000000E+06, .000000000000E+00, .000000000000E+00, .000000000000E+00],
        ['G21', 2022, 6, 15, 14, 0, 0, .161311123520E-03, .000000000000E+00, .000000000000E+00,
         .480000000000E+02, -.690000000000E+02, .473019703169E-08, .235557719582E+01,
         -.333972275257E-05, .242386998143E-01, -.119768083096E-05, .515374702072E+04,
         .309600000000E+06, .441446900368E-06, .214969297667E+01, .117346644402E-06,
         .960199874762E+00, .397781250000E+03, -.937445853228E+00, -.824212903204E-08,
         -.492163357722E-09, .000000000000E+00, .221400000000E+04, .000000000000E+00,
         .000000000000E+00, .000000000000E+00, -.102445483208E-07, .480000000000E+02,
         .309338000000E+06, .000000000000E+00, .000000000000E+00, .000000000000E+00],
        ['G22', 2022, 6, 15, 14, 0, 0, .276219565421E-03, .636646291241E-11, .000000000000E+00,
         .790000000000E+02, .637500000000E+01, .447268630534E-08, .209109460950E+01,
         .281259417534E-06, .134111947846E-01, .758469104767E-05, .515369627380E+04,
         .309600000000E+06, .134110450745E-06, -.187625357347E+01, -.329688191414E-06,
         .961125351849E+00, .230781250000E+03, -.184820560104E+01, -.769603485646E-08,
         .706100840506E-09, .000000000000E+00, .221400000000E+04, .000000000000E+00,
         .000000000000E+00, .000000000000E+00, -.838190317154E-08, .790000000000E+02,
         .309485000000E+06, .000000000000E+00, .000000000000E+00, .000000000000E+00],
        ['G24', 2022, 6, 15, 14, 0, 0, .220152083784E-03, -.841282599140E-11, .000000000000E+00,
         .840000000000E+02, -.106875000000E+02, .558523264736E-08, .346157288052E+00,
         -.622123479843E-06, .125489170896E-01, .565126538277E-05, .515367090797E+04,
         .309600000000E+06, -.782310962677E-07, -.991197407488E+00, .219792127609E-06,
         .933901572841E+00, .259656250000E+03, .817978952381E+00, -.853106963901E-08,
         -.304655547269E-09, .000000000000E+00, .221400000000E+04, .000000000000E+00,
         .000000000000E+00, .000000000000E+00, .279396772385E-08, .840000000000E+02,
         .309339000000E+06, .000000000000E+00, .000000000000E+00, .000000000000E+00],
        ['G27', 2022, 6, 15, 14, 0, 0, .215651467443E-03, .329691829393E-11, .000000000000E+00,
         .700000000000E+02, .148187500000E+03, .440161191623E-08, .225401738364E+01,
         .754743814468E-05, .105409963289E-01, .368244946003E-05, .515365896797E+04,
         .309600000000E+06, -.240281224251E-06, .118118244570E+01, -.150874257088E-06,
         .971982136512E+00, .313812500000E+03, .688193783108E+00, -.803997775449E-08,
         .273225666660E-09, .000000000000E+00, .221400000000E+04, .000000000000E+00,
         .100000000000E+01, .000000000000E+00, .186264514923E-08, .700000000000E+02,
         .309340000000E+06, .000000000000E+00, .000000000000E+00, .000000000000E+00]] 

ephs =[['G01', 2022,6,15,16,0,0,.340738333762E-03,-.784439180279E-11,.000000000000E+00,
        .540000000000E+02,-.773437500000E+02,.429303596504E-08,.107670373452E+01,
        -.395812094212E-05,.120082932990E-01,.612810254097E-06,.515366287804E+04,
        .316800000000E+06,-.912696123123E-07,.224382395479E+01,.147148966789E-06,
        .987854126192E+00,.385250000000E+03,.921173225877E+00,-.836713423901E-08,
        -.345014371233E-09,.000000000000E+00,.221400000000E+04,.000000000000E+00,
        .000000000000E+00,.000000000000E+00,.512227416039E-08,.540000000000E+02,
        .309628000000E+06,.000000000000E+00,.000000000000E+00,.000000000000E+00],
        ['G08',2022,6,15,16,0,0,-.722822733223E-04,-.159161572810E-11,.000000000000E+00,
        .830000000000E+02,.150718750000E+03,.456054710795E-08,-.303652790620E+01,
        .783987343311E-05,.744804507121E-02,.411830842495E-05,.515377832794E+04,
        .316800000000E+06,-.152736902237E-06,.116069700634E+01,-.447034835815E-07,
        .963248680888E+00,.303375000000E+03,.173456537710E+00,-.807605068564E-08,
        .300012496725E-09,.000000000000E+00,.221400000000E+04,.000000000000E+00,
        .000000000000E+00,.000000000000E+00,.512227416039E-08,.830000000000E+02,
        .309618000000E+06,.000000000000E+00,.000000000000E+00,.000000000000E+00],
        ['G10',2022,6,15,16,0,0,-.455989036709E-03,-.147792889038E-10,.000000000000E+00,
        .840000000000E+02,-.154406250000E+03,.376515683390E-08,-.935626010933E+00,
        -.817887485027E-05,.765509542543E-02,.104885548353E-04,.515369258881E+04,
        .316800000000E+06,.203028321266E-06,-.300855339128E+01,.912696123123E-07,
        .974342123272E+00,.186187500000E+03,-.249930062284E+01,-.761103131572E-08,
        -.496449250533E-10,.000000000000E+00,.221400000000E+04,.000000000000E+00,
        .100000000000E+01,.000000000000E+00,.232830643654E-08,.840000000000E+02,
        .309619000000E+06,.000000000000E+00,.000000000000E+00,.000000000000E+00],
        ['G14',2022,6,15,16,0,0,-.111529603601E-03,-.682121026330E-12,.000000000000E+00,
        .160000000000E+02,.826562500000E+02,.437018203562E-08,-.251014454845E+00,
        .437907874584E-05,.199296663050E-02,.110045075417E-04,.515366699409E+04,
        .316800000000E+06,.391155481339E-07,.173310839846E+00,.372529029846E-07,
        .952814682476E+00,.163093750000E+03,-.313214701434E+01,-.779961059937E-08,
        -.283226083217E-09,.000000000000E+00,.221400000000E+04,.000000000000E+00,
        .000000000000E+00,.000000000000E+00,-.791624188423E-08,.784000000000E+03,
        .309620000000E+06,.000000000000E+00,.000000000000E+00,.000000000000E+00],
        ['G21',2022,6,15,16,0,0,.161311123520E-03,.000000000000E+00,.000000000000E+00,
        .520000000000E+02,-.630625000000E+02,.467519474063E-08,-.287745316399E+01,
        -.290945172310E-05,.242399923736E-01,-.149756669998E-05,.515374453163E+04,
        .316800000000E+06,.204890966415E-07,.214963212220E+01,.389292836189E-06,
        .960198189480E+00,.405531250000E+03,-.937461246052E+00,-.809105131048E-08,
        -.133219834855E-09,.000000000000E+00,.221400000000E+04,.000000000000E+00,
        .000000000000E+00,.000000000000E+00,-.102445483208E-07,.520000000000E+02,
        .309621000000E+06,.000000000000E+00,.000000000000E+00,.000000000000E+00],
        ['G22',2022,6,15,16,0,0,.276264734566E-03,.636646291241E-11,.000000000000E+00,
        .820000000000E+02,.171875000000E+01,.443554190098E-08,.314129354829E+01,
        .184401869774E-06,.134117096895E-01,.768713653088E-05,.515369504356E+04,
        .316800000000E+06,.145286321640E-06,-.187630757124E+01,-.875443220139E-07,
        .961130650538E+00,.228937500000E+03,-.184823787448E+01,-.753888545341E-08,
        .637526555540E-09,.000000000000E+00,.221400000000E+04,.000000000000E+00,
        .000000000000E+00,.000000000000E+00,-.838190317154E-08,.820000000000E+02,
        .309622000000E+06,.000000000000E+00,.000000000000E+00,.000000000000E+00],
        ['G24',2022,6,15,16,0,0,.220091082156E-03,-.841282599140E-11,.000000000000E+00,
        .850000000000E+02,-.113750000000E+02,.552022993973E-08,.139631594392E+01,
        -.678002834320E-06,.125492629595E-01,.570528209209E-05,.515367191124E+04,
        .316800000000E+06,-.189989805222E-06,-.991258418486E+00,.838190317154E-07,
        .933898960069E+00,.256125000000E+03,.818009805323E+00,-.849606818106E-08,
        -.373229832235E-09,.000000000000E+00,.221400000000E+04,.000000000000E+00,
        .000000000000E+00,.000000000000E+00,.279396772385E-08,.850000000000E+02,
        .309623000000E+06,.000000000000E+00,.000000000000E+00,.000000000000E+00],
        ['G27',2022,6,15,16,0,0,.215674750507E-03,.329691829393E-11,.000000000000E+00,
        .710000000000E+02,.162375000000E+03,.431732269096E-08,-.297898170629E+01,
        .817701220512E-05,.105392733822E-01,.352598726749E-05,.515366197205E+04,
        .316800000000E+06,-.726431608200E-07,.118112470578E+01,-.257045030594E-06,
        .971984339666E+00,.317718750000E+03,.688196155961E+00,-.801533387083E-08,
        .314655963827E-09,.000000000000E+00,.221400000000E+04,.000000000000E+00,
        .100000000000E+01,.000000000000E+00,.186264514923E-08,.710000000000E+02,
        .309624000000E+06,.000000000000E+00,.000000000000E+00,.000000000000E+00]]

Step 1. Calculate satellite clock bias¶

  1. Calculate GPST second of week from receiver time tag GPST = "2022 06 15 13 59 50.0000000"
  2. Select the last transmitted navigation message block before instant GPST of week = 309590 (wsec)
In [84]:
# Example calculating sec of GPS week from Time of clock (GPS) year, month, day, hour, minute, second 
yy = 2022
M = 6
dd = 15
hh = 13
mm = 59
ss = 50
obstime = datetime.datetime(yy, M, dd, hh, mm, ss)
gps_day_of_week = obstime.weekday() + 1 #GPS week starts midnight between saturday and sunday
wsec = gps_day_of_week *86400 + hh*3600 + mm*60 + ss
print(wsec)
# Epoch: Toc - Time of clock (GPS) year, month, day, hour, minute, second   #SV / EPOCH / SV CLK
yy = 2022
M = 6
dd = 15
hh = 14
mm = 0
ss = 0
obstime = datetime.datetime(yy, M, dd, hh, mm, ss)
gps_day_of_week = obstime.weekday() + 1 #GPS week starts midnight between saturday and sunday
wsec = gps_day_of_week *86400 + hh*3600 + mm*60 + ss
print(wsec) # Toe Time of Ephemeris (sec of GPS week)   

309590
309600
In [85]:
def gps_sec_of_week(toc):
    yy = toc[0]
    M = toc[1]
    dd = toc[2]
    hh = toc[3]
    mm = toc[4]
    ss = toc[5]
    obstime = datetime.datetime(yy, M, dd, hh, mm, ss)
    gps_day_of_week = obstime.weekday() + 1 #GPS week starts midnight between saturday and sunday
    gps_sec = gps_day_of_week *86400 + hh*3600 + mm*60 + ss
    return gps_sec

Clock modelling¶

The clock offsets are due to clock synchronisation errors referring to the GNSS time scale. The modelling of such offsets, as well as its effects on the navigation solution, are described as follows:

Clock modelling

Receiver clock offset¶

$\delta t_{rcv}$ is estimated together with the receiver coordinates, hence no modelling is needed in this case.

Satellite clock offset¶

This can be split into two terms:

$\delta t_{sat} = \delta t_{clockdrift}+\Delta_{rel}$

The first term $\delta t_{sat}$ can be calculated from values broadcast in the navigation messages (for SPP) or from the precise products available from IGS centres or other providers (for PPP). The broadcast navigation message provides the clock information as the coefficients of a polynomial, in a given reference epoch (t0), to compute the satellite clock offset as

${\color{black}{\delta t_{clockdrift} = a_0 + a_1(t − t_0) + a_2(t − t_0)^2} }$

The relativistic effect is a function of Eccentic anomali:

$\Delta_{rel} = eph.e \cdot\sqrt{eph.A}\cdot sinE_k$

In [86]:
def calculate_clock_bias(eph, t):
    F = -4.442807633E-10              # -2*sqrt(MU)/c**2
    MU = 3.986005E+14                  # value of Earth's universal gravitational parameters
    A = eph.sqrtA**2                  # semi-major axis
    cmm = math.sqrt(MU/A**3)          # computed mean motion
    tk = t - eph.toe                  # time from ephemeris reference epoch
    if (tk > 302400):                 # account for beginning or end of week crossover
        tk = tk-604800
    if (tk < -302400):
        tk = tk+604800
    # apply mean motion correction
    n = cmm + eph.deln
    # Mean anomaly
    Mk = eph.M0 + n*tk
    # Eccentric anaomaly
    Ek = Mk                             # initial value    
    for x in range(3):                  # iterate 2 times
        Ek=Ek+(Mk-Ek + eph.e*math.sin(Ek))/(1-eph.e*math.cos(Ek))
    dt_rel = F*eph.e*eph.sqrtA*math.sin(Ek)
    dt_sat = eph.af0 + eph.af1*(t-eph.toc) + eph.af2*(t-eph.toc)**2 + dt_rel
    print(eph.sat," SV clock drift bias(m): ", dt_sat *C, " SV relativity clock bias(m): ", dt_rel * C)
    return [dt_sat, dt_rel]
In [87]:
def make_eph_list():
    eph_list = []
    for eph in ephs:  
        sv = Eph(eph[0])
        sv.iode = eph[10]
        sv.iodc =eph[33]
        sv.af0 = eph[7]
        sv.af1 = eph[8]
        sv.af2 = eph[9]
        sv.toc = gps_sec_of_week(eph[1:7])
        sv.toe = eph[18]
        sv.tot = eph[34]
        sv.week = eph[28]
        sv.crs = eph[11]
        sv.crc = eph[23]
        sv.cus = eph[16]
        sv.cuc = eph[14]
        sv.cis = eph[21]
        sv.cic = eph[19]
        sv.e = eph[15]
        sv.i0 = eph[22]
        sv.sqrtA = eph[17]
        sv.deln = eph[12]
        sv.M0 = eph[13]
        sv.OMG0 = eph[20]
        sv.OMGd = eph[25]
        sv.omg = eph[24]
        sv.idot = eph[26]
        sv.tdg = eph[32]
        sv.sva = eph[30]
        sv.health = eph[31]
        eph_list.append(sv)
    return eph_list
In [88]:
sv_ephs = make_eph_list()
satobs_eph = [] # list of observations:
                # time tag of observation                 rec 0
                # pseudorange measurement C1C,  C2W/C2L   record 1 and 2
                # ephemeris of sv linked to observation   record 3
                # computed clock bias of sv               record 4
                # computed emmision time                  record 5
                # computed Xsat at emmision time          record 6
                # computed Ysat at emmision time          record 7
                # computed Zsat at emmision time          record 8
                # computed Xsat at reception time         record 9
                # computed Ysat at reception time         record 10
                # computed Zsat at reception time         record 11
                # computed troposphere correction         record 12
                # computed sv relativistic clock bias     record 13
                

for obs in satobs:
    rcv = gps_sec_of_week(obs[1:7]) # Receiver reception time of signal
    sv = obs[0]
    for eph in sv_ephs: 
        if sv == eph.sat:
            satobs_eph.append([rcv, obs[7], obs[8], eph ,0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
i = 0
for mea in satobs_eph:
    dt_sat = calculate_clock_bias(mea[3], mea[0]) #
    satobs_eph[i][4] = dt_sat[0] # dt sat including realtivistic bias
    satobs_eph[i][13] = dt_sat[1]  # dt sat realtivistic bias
    i = i + 1

G01  SV clock drift bias(m):  102167.38652414233  SV relativity clock bias(m):  -0.2577029075724685
G08  SV clock drift bias(m):  -21670.371872938864  SV relativity clock bias(m):  -4.1126808044774945
G10  SV clock drift bias(m):  -136665.50278015406  SV relativity clock bias(m):  4.803098164396179
G14  SV clock drift bias(m):  -33432.9500639449  SV relativity clock bias(m):  1.3177120849272912
G21  SV clock drift bias(m):  48348.336859136216  SV relativity clock bias(m):  -11.521363666201287
G22  SV clock drift bias(m):  82800.48503129893  SV relativity clock bias(m):  -7.914015016737841
G24  SV clock drift bias(m):  65996.73687962558  SV relativity clock bias(m):  -2.993093583476232
G27  SV clock drift bias(m):  64645.02356264339  SV relativity clock bias(m):  -5.553255202087061

Step 2. Compute emission time of satellite signal¶

We need to know the emission time of the satellite signal to calculate the position of the satellite in the reference frame. As the measurements are linked to the signal reception time, which is given by the receiver time tags (i.e. in the receiver clock), an algorithm is needed to obtain the signal emission time in the GNSS (GPS, Glonass, Galileo, etc.) time scale. The emission time can be directly obtained from the reception time by taking into account that the pseudorange R is a direct measurement of the time difference between both epochs, each one being measured in the corresponding clock.

Emission time computation from receiver time-tag and code pseudorange:

$T[emission]$ in this system time scale can be computed from the receiver measurement time tags $t_{rcv}$ receiver as

${\color{black}{T_{[emission]} = t_{sat[emission]} − \delta t_{sat} = t_{rcv[reception]} − R/C −\delta t_{sat}}}$

where $R$ is the pseudorange and $\delta t_{sat}$ satellite clock bias

In [89]:
def emission_time(r, t_rec, dt_sat):
    tems = t_rec-(r/C+dt_sat)
    return tems
In [90]:
i = 0
for mea in satobs_eph:
    p1 = mea[1] 
    p2 = mea[2]
    ionosfreeP1P2=(F1**2*p1-F2**2*p2)/(F1**2-F2**2) # Using inosphere free pseudrorange
    reception_time_rec = mea[0]
    dt_sv = mea[4]
    satobs_eph[i][5] = emission_time(ionosfreeP1P2, reception_time_rec, dt_sv)
    print('SV ', mea[3].sat)
    print('Ionos-free pseudorange ', ionosfreeP1P2)
    print('Computed emission time of SV signal: ',satobs_eph[i][5])
    i = i + 1

SV  G01
Ionos-free pseudorange  21985773.473138683
Computed emission time of SV signal:  309629.92632255994
SV  G08
Ionos-free pseudorange  22000890.960214686
Computed emission time of SV signal:  309629.9266852117
SV  G10
Ionos-free pseudorange  21611151.549600683
Computed emission time of SV signal:  309629.92836882494
SV  G14
Ionos-free pseudorange  22855369.76813225
Computed emission time of SV signal:  309629.92387421295
SV  G21
Ionos-free pseudorange  21222593.648853585
Computed emission time of SV signal:  309629.9290477748
SV  G22
Ionos-free pseudorange  24085805.445822243
Computed emission time of SV signal:  309629.9193822083
SV  G24
Ionos-free pseudorange  23063893.09826535
Computed emission time of SV signal:  309629.92284699227
SV  G27
Ionos-free pseudorange  24344775.883872915
Computed emission time of SV signal:  309629.918578936

Step 3. Solve satellite coordinates ECEF at emission time of pseudorange¶

Orbital calculation of satellite coordinates in Conventional Terrestrial System - CTS, Earth Centered Earth Fixed

Compute the time $t_k$ (time of data) from the ephemerides reference epoch $toe$ ($t$ and $toe$ are expressed in seconds in the GPS week):

$toe$ - Time of Ephemeris (sec of GPS week)
$t$ - emission time

$t_k = t − toe$

If $t_k$ > 302 400 s, subtract 604 800 s from tk. If $tk$ < −302 400 s, add 604 800 s.( week rollover)

Compute the mean anomaly for tk:

$M_k = eph.M0 + eph.deln\cdot t_k$

Solve (iteratively) the Kepler equation for the eccentric anomaly Ek:

$M_k = E_k − e\cdot sinE_k$

$E_k = M_k$ initial value

Iterate:

$ E_k=E_k+ \frac{M_k-E_k-eph.e\cdot sinE_k}{1-eph.e\cdot cosE_k}$

Compute the true anomaly $v_k$ (use atan2):

$v_k = atan2(\frac{(1-eph.e^2)\cdot sinE_k}{1-eph.e \cdot cosE_k}, \frac{cosE_k-eph.e}{1-eph.e\cdot cosE_k})$

Compute the argument of latitude $u_k$ from the argument of perigee $\omega$, true anomaly $v_k$ and corrections $cuc$ and $cus$:

$u_k = \omega + v_k + eph.cuc \cdot cos2(\omega + v_k) + eph.cus \cdot sin2(\omega + v_k)$

Compute the radial distance $r_k$, considering corrections $crc$ and $crs$

$r_k = eph.sqrtA^2(1 − eph.e\cdot cosE_k) + eph.crc\cdot cos2(\omega + v_k) + eph.crs\cdot sin2(\omega + v_k)$

Compute the inclination $i_k$ of the orbital plane from the inclination $io$ at reference time $toe$, and corrections $cic$ and $cis$

$\begin{align*} i_k = eph.io+ eph.idot t_k + eph.cic\cdot cos2(\omega + v_k) + eph.cis\cdot sin2(\omega + v_k) \end{align*}$

Compute the longitude of the ascending node $\lambda_k$ (with respect to Greenwich). This calculation uses the right ascension at the beginning of the current week ($\Omega_0$), the correction from the apparent sidereal time variation in Greenwich between the beginning of the week and reference time $t_k = t − toe$, and the change in longitude of the ascending node from the reference time $toe$:

$\lambda_k = \Omega_0 + t_k(\dot\Omega − ωE) − ωE\cdot toe$

Compute the satellite coordinates in the orbital plane:

Perifocal coordinate system

$ x_k = r_k\cdot cos(u_k)$

$ y_k = r_k\cdot sin(u_k)$

Compute the satellite coordinates in the CTS frame, applying three rotations (around $u_k, i_k$ and $\lambda_k$):

$ x_{CTS} = x_k\cdot cos(\omega) - y_k\cdot cos(i_k)\cdot sin(\omega)$

$ y_{CTS} = x_k\cdot sin(\omega) + yk\cdot cos(i_k)\cdot cos(\omega)$

$ z_{CTS} = y_k\cdot sin(i_k)$

In [91]:
def orbit_calculation(eph, t):
    F = -4.442807633e-10              # -2*sqrt(MU)/c**2
    MU = 3.986005e14                  # value of Earth's universal gravitational parameters
    A = eph.sqrtA**2                  # semi-major axis
    cmm = math.sqrt(MU/A**3)          # computed mean motion
    tk = t - eph.toe                  # time from ephemeris reference epoch
    if (tk > 302400):                 # account for beginning or end of week crossover
        tk = tk-604800
    if (tk < -302400):
        tk = tk+604800
    # apply mean motion correction
    n = cmm + eph.deln
    # Mean anomaly
    Mk = eph.M0 + n*tk
    # Eccentric anomaly
    Ek = Mk                             # initial value    
    for x in range(3):                  # iterate 2 times
        Ek=Ek+(Mk-Ek+eph.e*math.sin(Ek))/(1-eph.e*math.cos(Ek))
    # True anomaly
    vk = math.atan2((math.sqrt(1-eph.e**2))*math.sin(Ek)/(1-eph.e*math.cos(Ek)), (math.cos(Ek)-eph.e)/(1-eph.e*math.cos(Ek)))
    
    # vk = 2*math.atan((math.sqrt((1+eph.e)/(1-eph.e))*math.tan(Ek/2)))

    # argument of latitude
    phi = vk + eph.omg
    
    # compute harmonic corrections
    du = eph.cus*math.sin(2*phi) + eph.cuc*math.cos(2*phi); #argument of lattitude correction
    dr = eph.crs*math.sin(2*phi) + eph.crc*math.cos(2*phi); #radius correction
    di = eph.cis*math.sin(2*phi) + eph.cic*math.cos(2*phi); #inclination correction
    
    # corrected argument of lattitude
    uk = phi + du
    
    # corrected radius
    rk =  A*(1-eph.e*math.cos(Ek)) + dr
    
    # corrected inclination angle at reference time
    ik = eph.i0 + eph.idot*tk + di; 
    
    # corrected longitude of ascending node
    omega = eph.OMG0 + (eph.OMGd - WE)*tk - WE*eph.toe 
    
    xk = rk*math.cos(uk) #x posistion in orbital plane
    yk = rk*math.sin(uk) #y posistion in orbital plane
    
    x = xk*math.cos(omega) - yk*math.cos(ik)*math.sin(omega); #Earth-fixed coordinates ECEF
    y = xk*math.sin(omega) + yk*math.cos(ik)*math.cos(omega);
    z = yk*math.sin(ik);
    print('SV',eph.sat)

    return[x,y,z]
In [92]:
i = 0
print('Satellite coordinates in ECEF at signal emission time')
for mea in satobs_eph:
    svXYZ = orbit_calculation(mea[3],mea[5])  #input eph and sv emission time, ouput sv coordinates at emmision time in CTS
    satobs_eph[i][6] = svXYZ[0]
    satobs_eph[i][7] = svXYZ[1]
    satobs_eph[i][8] = svXYZ[2]
    print('Emission time (s): ',mea[5])
    print('X ', satobs_eph[i][6],'Y ', satobs_eph[i][7],'Z ', satobs_eph[i][8])
    i = i + 1
    

Satellite coordinates in ECEF at signal emission time
SV G01
Emission time (s):  309629.92632255994
X  13031293.310108224 Y  -14140924.611738503 Z  17855617.049962882
SV G08
Emission time (s):  309629.9266852117
X  21981431.907177202 Y  1766152.6526481416 Z  15015223.581840554
SV G10
Emission time (s):  309629.92836882494
X  1242509.956379449 Y  15655321.7354242 Z  21522353.842483167
SV G14
Emission time (s):  309629.92387421295
X  758181.5523897447 Y  -16481033.125235895 Z  20796258.11599192
SV G21
Emission time (s):  309629.9290477748
X  15365982.640044274 Y  -3228753.5219289847 Z  21995975.30020505
SV G22
Emission time (s):  309629.9193822083
X  17509064.222957592 Y  19347881.04037465 Z  5853841.371709985
SV G24
Emission time (s):  309629.92284699227
X  -14337055.003536966 Y  10177482.993712874 Z  19488564.12592963
SV G27
Emission time (s):  309629.918578936
X  23309804.967000924 Y  12499352.057187416 Z  3929408.615245241

Step 4. Transform satellite coordinates to signal reception time¶

Account for earth rotation $\alpha$ during signal propagation time from emission to reception $ \delta t$. Earth fixed reference system (Conventional Terrestrial System - CTS) at reception time is used to build the navigation equations.

Looking down to the North Pole from space: Sagnac effect

$ \alpha = \omega E\cdot \delta t_{prop} = ~ 0.07 sec $

Methodology:

  • Solve the position of satellite at signal emission in previous step $(x_E,y_E,x_E)$
  • Rotate into ECEF reference frame at time og reception $(x_R,y_R,x_R)$

Neclecting atmospheric delay, the signal propagation time is given by:

Signal propagation time $\delta t_{prop} = \frac{|P_{sat} - P_{rec}|}{C} = \frac{p0_{sat -rec}}{C}$

$P_{sat}$ = satellite ECEF position vector

$P_{rec}$ = receiver ECEF position Vector

approximate value of $ p0_{sat-rec}$ can be used in this equation

$ [Xsat_R,Ysat_R,Zsat_R]_{CTS-reception} = R3(\alpha)\cdot [Xsat_E,Ysat_E,Zsat_E]_{CTS-emission}$

$\begin{bmatrix}x_R \\y_R \\z_R \end{bmatrix}=\begin{bmatrix}cos\alpha & \sin\alpha & 0 \\-sin\alpha & \cos\alpha & 0\\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix}x_E \\y_E\\z_E \end{bmatrix}$

The geometric range between satellite coordinates at emission time and the “approximate position of the receiver” at reception time (both coordinates given in the same reference system CTS system) is computed by:

$ p0_{sat-rec} = \sqrt{(x_{sat}-x_{0rec})^2+(y_{sat}-y_{0rec})^2+(z_{sat}-z_{0rec})^2}$

In [93]:
def earth_rotation_coor(Xsat, Ysat, Zsat, transmission_time):
    rot=WE*transmission_time 
    rsin = math.sin(rot)
    rcos = math.cos(rot)
    #Rotating the Satellite Reference Frame
    Xsat_reception = rcos*Xsat + rsin*Ysat 
    Ysat_reception = -rsin*Xsat + rcos*Ysat
    Zsat_reception = Zsat
    return [Xsat_reception, Ysat_reception, Zsat_reception]
    
def geometric_range(X0rec, Y0rec, Z0rec, Xsat, Ysat, Zsat):
    p0sat_rec =  math.sqrt((Xsat-X0rec)**2+(Ysat-Y0rec)**2+(Zsat-Z0rec)**2)
    return p0sat_rec

def transmission_time(p0sat_rec):
    transmission_time = p0sat_rec/C # negative time but positiv definition in rotation
    return transmission_time
    
def reception_time_satcoords(X0rec, Y0rec, Z0rec, Xsat, Ysat, Zsat):
    p0 = geometric_range(X0rec, Y0rec, Z0rec, Xsat, Ysat, Zsat)
    dt = transmission_time(p0)
    svCTS_reception = earth_rotation_coor(Xsat, Ysat, Zsat, dt)
    return svCTS_reception
In [94]:
i = 0
for mea in satobs_eph:
        print('SV',mea[3].sat, 'ECEF-coordinates of sattelite at signal reception time')
        svCTS_reception = reception_time_satcoords(XYZapprox[0], XYZapprox[1], XYZapprox[2], mea[6], mea[7], mea[8])
        for x in range(1):  # check difference p0sat_rec
            p0 = geometric_range(XYZapprox[0], XYZapprox[1], XYZapprox[2], svCTS_reception[0], svCTS_reception[1] ,svCTS_reception[2])
            dt = transmission_time(p0)
            print('Propagation time of signal:', dt)
            svCTS_reception = earth_rotation_coor(mea[6], mea[7], mea[8], dt)
            print('X ', svCTS_reception[0], 'Y ', svCTS_reception[1],'Z ', svCTS_reception[2]) # Reception time sat 
            satobs_eph[i][9] = svCTS_reception[0]
            satobs_eph[i][10] = svCTS_reception[1]
            satobs_eph[i][11] = svCTS_reception[2]
        i = i + 1 
for mea in satobs_eph:
    print('SV ',mea[3].sat, ': Impact of earth rotaion on geometric range:')
    print(math.sqrt((mea[6]-mea[9])**2 + (mea[7]-mea[10])**2 +(mea[8]-mea[11])**2)) # earth rotation impact (distance)
        

SV G01 ECEF-coordinates of sattelite at signal reception time
Propagation time of signal: 0.0736773734213029
X  13031217.335838398 Y  -14140994.623967856 Z  17855617.049962882
SV G08 ECEF-coordinates of sattelite at signal reception time
Propagation time of signal: 0.07331472738948755
X  21981441.349058382 Y  1766035.135616039 Z  15015223.581840554
SV G10 ECEF-coordinates of sattelite at signal reception time
Propagation time of signal: 0.07163111158659244
X  1242591.7307322798 Y  15655315.245055374 Z  21522353.842483167
SV G14 ECEF-coordinates of sattelite at signal reception time
Propagation time of signal: 0.07612571895154371
X  758090.0632776495 Y  -16481037.3337805 Z  20796258.11599192
SV G21 ECEF-coordinates of sattelite at signal reception time
Propagation time of signal: 0.07095216096756103
X  15365965.93454789 Y  -3228833.024147361 Z  21995975.30020505
SV G22 ECEF-coordinates of sattelite at signal reception time
Propagation time of signal: 0.08061771194983812
X  17509177.963747263 Y  19347778.10886743 Z  5853841.371709985
SV G24 ECEF-coordinates of sattelite at signal reception time
Propagation time of signal: 0.07715293631283941
X  -14336997.743966822 Y  10177563.654983908 Z  19488564.12592963
SV G27 ECEF-coordinates of sattelite at signal reception time
Propagation time of signal: 0.08142096809997065
X  23309879.179127373 Y  12499213.659411553 Z  3929408.615245241
SV  G01 : Impact of earth rotaion on geometric range:
103.31409359130132
SV  G08 : Impact of earth rotaion on geometric range:
117.89572492010201
SV  G10 : Impact of earth rotaion on geometric range:
82.03151631174106
SV  G14 : Impact of earth rotaion on geometric range:
91.585858513515
SV  G21 : Impact of earth rotaion on geometric range:
81.23839200892466
SV  G22 : Impact of earth rotaion on geometric range:
153.4009857006593
SV  G24 : Impact of earth rotaion on geometric range:
98.91864848388583
SV  G27 : Impact of earth rotaion on geometric range:
157.03943477895933

Step 5. Calculate elevation and azimut from receiver to satellite¶

From ECEF(x,y,s) to Local Geodetic(e,n,u) coordinates:

  1. Compute vector from receiver to satellite $V_{rec-sat}$ at reception time in ECEF (start with approx values)
  2. Compute geodetic coordinates of receiver ($ \lambda, \phi, h$) (start with approx values)
  3. Transform $V_{rec-sat}$ from ECEF to Local Geodetic System ($e,n,u$)

$\begin{bmatrix}\Delta e\\\Delta n\\\Delta u\end{bmatrix}=\bf{R_1 }\begin{bmatrix}\pi/2-\phi \end{bmatrix}\bf{R_3}\begin{bmatrix}\pi/2+\lambda \end{bmatrix}\begin{bmatrix}\Delta x\\\Delta y\\\Delta z\end{bmatrix} =\begin{bmatrix}-sin\lambda & cos\lambda & 0 \\-cos\lambda sin\phi & -sin\lambda sin\phi & cos\phi\\ cos\lambda cos\phi & sin\lambda cos\phi & sin\phi \end{bmatrix}\begin{bmatrix}\Delta x \\\Delta y\\\Delta z \end{bmatrix}$

  1. Calculate asimut $Az$ and elevation $E$ of satellite in Local Geodetic System.

$~~Az = atan2(\frac{\Delta e}{\Delta n})$ $~~E = atan(\frac{\Delta u}{\sqrt{\Delta e^2 + \Delta n^2 }})$

  1. Use $E_{sat}^j$ to calculate tropospehric slant correction for each satellite $j$.

Local geodetic system

In [95]:
def dXYZ_receiver_sat(xrec, yrec, zrec, xsat, ysat, zsat): #ECEF rec-sat vector 
    dxyz =[xsat-xrec, ysat-yrec, zsat-zrec]
    return dxyz
In [96]:
def rad_deg(radians):
    degrees = 180 * radians/math.pi
    return degrees
In [97]:
def deg_rad(degrees):
    radians = math.pi * degrees/180
    return radians
In [98]:
 def cartesian_geographic(da, db, dx, dy, dz):
    # da major axis
    # db minor axis
    # d_la latitude of point          *Changing arg
    # d_lg longitude of point         *Changing arg
    # dh ellipsoidal height           *Changing arg
    # dx x-axis
    # dy y-axis
    # dz z-axis

    de1 = (da ** 2 - db ** 2) / (da ** 2)
    d_e2 = (da ** 2 - db ** 2) / (db ** 2)
    dp = math.sqrt(dx ** 2 + dy ** 2)
    d_t = (dz * da) / (dp * db)
    d_teta = math.atan(d_t)
    d_c2 = math.cos(d_teta)
    d_s2 = math.sin(d_teta)
    d_t1 = (dz + (d_e2 * db * d_s2 ** 3)) / (dp - de1 * da * d_c2 ** 3)
    d_t2 = dy / dx
    d_la = math.atan(d_t1)
    d_lg = math.atan(d_t2)
    d_c = math.cos(d_la)
    d_s = math.sin(d_la)
    d_n = da ** 2 / (math.sqrt((da ** 2 * d_c ** 2) + (db ** 2 * d_s ** 2)))
    dh = (dp / math.cos(d_la)) - d_n
    return d_lg, d_la, dh
In [99]:
# test cartesian to geodetic

point = cartesian_geographic(A_WGS84, B_WGS84, 2660520.2488, 589219.3944, 5747982.8939)
print('Longitude: ',rad_deg(point[0])) # longitude
print('Latitude: ',rad_deg(point[1])) # latitude
print('Ellipsoidal height: ',point[2]) # ellipsoidal height

Longitude:  12.48760652265354
Latitude:  64.78409006861396
Ellipsoidal height:  538.7430106811225
In [100]:
def dxyz_denu(long, lat, vector): # Transform ECEF rec-sat vector to local geodetic system enu
    dx = vector[0]
    dy = vector[1]
    dz = vector[2]
    sinl = math.sin(long)
    sinb = math.sin(lat)
    cosl = math.cos(long)
    cosb = math.cos(lat)
    de = -sinl*dx + cosl*dy
    dn = -cosl*sinb*dx - sinl*sinb*dy + cosb*dz
    du = cosl*cosb*dx + sinl*cosb*dy + sinb*dz
    return de,dn,du
In [101]:
def azimut(de,dn):
    azi = math.atan2(de, dn)
    return azi
In [102]:
def elevation(de,dn,du):
    elev = math.atan2(du, math.sqrt(de**2 + dn**2))
    return elev
In [103]:
i = 0
xrec = XYZapprox[0]
yrec = XYZapprox[1]
zrec = XYZapprox[2]
point = cartesian_geographic(A_WGS84, B_WGS84, xrec, yrec, zrec)
print('Position of receiver: ')
print('Longitude of revceiver: ',rad_deg(point[0])) # longitude
print('Latitude of receiver: ',rad_deg(point[1])) # latitude
print('Ellipsoidal height of receiver: ',point[2]) # ellipsoidal height
for mea in satobs_eph:
    print('SV',mea[3].sat)
    dxyz = dXYZ_receiver_sat(xrec, yrec, zrec, mea[9], mea[10], mea[11])
    print ('Geometric range SV-receiver: ', math.sqrt(dxyz[0]**2+dxyz[1]**2+dxyz[2]**2))
    dneu = dxyz_denu(point[0], point[1], dxyz)
    print('Azimut SV:',rad_deg(azimut(dneu[0], dneu[1])))
    print('Elevation angle SV:' ,rad_deg(elevation(dneu[0], dneu[1], dneu[2])))

Position of receiver: 
Longitude of revceiver:  23.276599787253613
Latitude of receiver:  70.49576914993929
Ellipsoidal height of receiver:  38.70248999912292
SV G01
Geometric range SV-receiver:  22087920.87696028
Azimut SV: -90.13007647727463
Elevation angle SV: 34.79023770350955
SV G08
Geometric range SV-receiver:  21979202.331695005
Azimut SV: -154.27954135049987
Elevation angle SV: 42.218620172251455
SV G10
Geometric range SV-receiver:  21474467.011819255
Azimut SV: 88.79645013548918
Elevation angle SV: 49.687197212120154
SV G14
Geometric range SV-receiver:  22821916.40150329
Azimut SV: -51.14121857614808
Elevation angle SV: 29.683929402086285
SV G21
Geometric range SV-receiver:  21270922.736877814
Azimut SV: -117.70007173418381
Elevation angle SV: 61.320476707944806
SV G22
Geometric range SV-receiver:  24168582.02377926
Azimut SV: 151.98471741159088
Elevation angle SV: 17.038037029769498
SV G24
Geometric range SV-receiver:  23129868.419146217
Azimut SV: 44.755381600680195
Elevation angle SV: 22.782574605077336
SV G27
Geometric range SV-receiver:  24409392.159429844
Azimut SV: 174.4847310279382
Elevation angle SV: 14.576198068402974

Satellite positions in space¶

In [104]:
# Import libraries

from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt

%matplotlib notebook
 
# Define Data
x = []
y = []
z = []
recx = []
recy = []
recz = []
sat = []
for mea in satobs_eph:
    x.append(mea[9])
    y.append(mea[10])
    z.append(mea[11])
    sat.append(mea[3].sat)
recx.append(XYZapprox[0])
recy.append(XYZapprox[1])
recz.append(XYZapprox[2])

 
# Create Figure

fig = plt.figure(figsize = (10, 7))
ax = plt.axes(projection ="3d")

 
# Create Plot and change color

ax.scatter3D(x, y, z, color='red')
ax.scatter3D(recx, recy, recz, color='green')
for mea in satobs_eph: 
    ax.text(mea[9], mea[10], mea[11], mea[3].sat)
ax.text(recx[0], recy[0], recz[0], 'REC')

for mea in satobs_eph: 
    x, y, z = [recx[0], mea[9]], [recy[0], mea[10]], [recz[0], mea[11]]
    ax.plot(x, y, z, color='blue', linewidth=0.2)

u = np.linspace(0, 2 * np.pi, 25)
v = np.linspace(0, np.pi, 25)
x = 6371000 * np.outer(np.cos(u), np.sin(v))
y = 6371000 * np.outer(np.sin(u), np.sin(v))
z = 6371000 * np.outer(np.ones(np.size(u)), np.cos(v))

ax.plot_surface(x, y, z, color='grey' )

ax.view_init(10, -75)
   
# Show plot

plt.show()

Step 6. Calculate Troposheric slant correction¶

The tropospheric delay does not depend on frequency and affects both the code and carrier phases in the same way. It can be modeled (about 90%) as:

  • dry delay corresponds to the vertical delay of the dry atmosphere (basically oxygen and nitrogen in hydrostatical equilibrium) It can be modeled as an ideal gas.
  • wet delay corresponds to the vertical delay of the wet component (water vapor) - difficult to model.

Troposheric delay

A simple model is model presented here is from [Collins, 1999] and it is the model adopted by the SBAS (WAAS, EGNOS...)

$ Trop_{sat-rec} = (d_{dry} + d_{wet})\cdot M_{E_{sat}} $

$ d_{dry} = 2.3exp(-0.116 \cdot 10^{-3}\cdot H) m $

$ d_{wet} = 0.1m$

$ M_{E_{sat}} = \frac{\text{0.01}}{\sqrt{0.002001 + sin^2(elev)}} $

The obliquity factor $M(E)$ as function of elevation angle

In [105]:
def tropospheric_delay(h_sealevel, sat_elev):
    d_dry = 2.3* math.exp(-0.116E-03 * h_sealevel)
    d_wet = 0.1
    m_elev = 1.001/math.sqrt(0.002001 + math.sin(sat_elev)**2)
    troposd_sat_rec = (d_dry + d_wet) * m_elev
    return troposd_sat_rec
In [106]:
print('Tropsheric delay at 42m and elev 14.865201084274346 :' ,tropospheric_delay(42, deg_rad(14.865201084274346)))

Tropsheric delay at 42m and elev 14.865201084274346 : 9.18225409265146
In [107]:
xrec = XYZapprox[0]
yrec = XYZapprox[1]
zrec = XYZapprox[2]
point = cartesian_geographic(A_WGS84, B_WGS84, xrec, yrec, zrec)
print('Receiver position:')
print('Longitude: ',rad_deg(point[0])) # longitude
print('Latitude: ',rad_deg(point[1])) # latitude
print('Height: ', point[2]) # ellipsoidal height
print("Tropospheric delay:")
i = 0
for mea in satobs_eph:
    dxyz = dXYZ_receiver_sat(xrec, yrec, zrec, mea[9], mea[10], mea[11])
    dneu = dxyz_denu(point[0], point[1], dxyz)
    tropo_d = tropospheric_delay(point[2], elevation(dneu[0], dneu[1], dneu[2]))
    satobs_eph[i][12] = tropo_d
    print('SV',mea[3].sat, satobs_eph[i][12])
    i = i + 1 

Receiver position:
Longitude:  23.276599787253613
Latitude:  70.49576914993929
Height:  38.70248999912292
Tropospheric delay:
SV G01 4.179599477413308
SV G08 3.551996515901868
SV G10 3.1316813259152916
SV G14 4.810814002878507
SV G21 2.7230540901157143
SV G22 8.070429923387122
SV G24 6.136540694026958
SV G27 9.35827713799435

Step 7. Calculate prefit residual¶

Prefit residual is the difference between measured and modeled pseudorange

$PR_{rec}^{sat} = \rho + c (\delta t_{rec} - \delta t^{sat}) + \sum \delta_k + \epsilon$

Prefit residual:

$ Prefit_{rec}^{sat}= PR_{rec}^{sat} - C_{modeled}$

Modeled range includes estimated "true range" $\rho0_{rec}^{sat}$:

$ C_{modeled} = \rho0_{rec}^{sat} - C(\delta t_{sat} + \delta_{rel}) + Trop_{rec}^{sat} + Iono_{rec}^{sat} + TGD$

In [108]:
xrec = XYZapprox[0]
yrec = XYZapprox[1]
zrec = XYZapprox[2]
for mea in satobs_eph:
    print('SV: ', mea[3].sat)
    p0_rs = geometric_range(xrec, yrec, zrec, mea[9], mea[10], mea[11])
    print('P0 geometric range: ',p0_rs )
    cdt_sat = C*mea[4]
    print('Sv clock bias(s): ',mea[4])
    print('Sv clock bias(m):', cdt_sat)
    tropo = mea[12] 
    print('Tropospehric delay(m): ', tropo)
    c_modeled = p0_rs-cdt_sat+tropo
    print("C modeled :" ,c_modeled)
    p1 = mea[1] 
    p2 = mea[2]
    ionosfreeP1P2=(F1**2*p1-F2**2*p2)/(F1**2-F2**2) # Using inosphere free pseudrorange 
    prefit  = ionosfreeP1P2 - c_modeled
    prefitL1  = p1 - c_modeled
    print('Prefit residual ionos: ',prefit) 
    print('Prefit L1: ',prefitL1) 
    

SV:  G01
P0 geometric range:  22087920.87696028
Sv clock bias(s):  0.0003407937184468541
Sv clock bias(m): 102167.38652414233
Tropospehric delay(m):  4.179599477413308
C modeled : 21985757.67003562
Prefit residual ionos:  15.803103063255548
Prefit L1:  3.1899643801152706
SV:  G08
P0 geometric range:  21979202.331695005
Sv clock bias(s):  -7.228457986404335e-05
Sv clock bias(m): -21670.371872938864
Tropospehric delay(m):  3.551996515901868
C modeled : 22000876.25556446
Prefit residual ionos:  14.704650226980448
Prefit L1:  3.2044355422258377
SV:  G10
P0 geometric range:  21474467.011819255
Sv clock bias(s):  -0.0004558670477966262
Sv clock bias(m): -136665.50278015406
Tropospehric delay(m):  3.1316813259152916
C modeled : 21611135.646280736
Prefit residual ionos:  15.90331994742155
Prefit L1:  2.733719263225794
SV:  G14
P0 geometric range:  22821916.40150329
Sv clock bias(s):  -0.00011152031737884782
Sv clock bias(m): -33432.9500639449
Tropospehric delay(m):  4.810814002878507
C modeled : 22855354.16238124
Prefit residual ionos:  15.605751011520624
Prefit L1:  0.45761876180768013
SV:  G21
P0 geometric range:  21270922.736877814
Sv clock bias(s):  0.00016127269238753236
Sv clock bias(m): 48348.336859136216
Tropospehric delay(m):  2.7230540901157143
C modeled : 21222577.123072766
Prefit residual ionos:  16.52578081935644
Prefit L1:  -2.703072763979435
SV:  G22
P0 geometric range:  24168582.02377926
Sv clock bias(s):  0.0002761926887143336
Sv clock bias(m): 82800.48503129893
Tropospehric delay(m):  8.070429923387122
C modeled : 24085789.609177887
Prefit residual ionos:  15.836644355207682
Prefit L1:  3.470822110772133
SV:  G24
P0 geometric range:  23129868.419146217
Sv clock bias(s):  0.00022014141823282818
Sv clock bias(m): 65996.73687962558
Tropospehric delay(m):  6.136540694026958
C modeled : 23063877.818807285
Prefit residual ionos:  15.279458064585924
Prefit L1:  4.521192714571953
SV:  G27
P0 geometric range:  24409392.159429844
Sv clock bias(s):  0.00021563258793736364
Sv clock bias(m): 64645.02356264339
Tropospehric delay(m):  9.35827713799435
C modeled : 24344756.49414434
Prefit residual ionos:  19.3897285759449
Prefit L1:  7.7658556625247

Step 8. Least square solution of SPS position equation¶

Solving SPS position equation and receiver position

  • Pseudorange $PR_{rec}^{sat}$
    • raw measurement of receiver
    • includes user clock bias, satellite clock bias, ionospeherec delay, tropospheric delay, multipath and other signal errors
  • Modeled range $C_{modeled}$
    • constructed from asumed location of receiver, navigation message and satellite position
    • includes estimates of satellite clock bias, troposheric and ionospheric signal delays
  • Difference of Pseudorange an Modeled range (theorange) is a linear function of unknowns
  • Solve system of linear unknowns for position increment and clock bias

$\color{blue} {color~blue = measured}$

$\color{green} {color~green = computed}$

$\color{red} {color~red = estimated}$

$\color{blue}{PR_{rec}^{sat}} = \color{red}{\rho_{rec}^{sat}} + \color{green}c \cdot (\color{red}{\delta t_{rec}} - \color{green}{\delta t^{sat}}) + \color{green}{\sum \delta_k} + \epsilon$

Form the ionos-free combination of L1 and L2 or L5 pseudoranges:

$\color{blue}{PR_{_{\mbox{iono-free}}}^{sat-rec}}=\color{blue}{\frac{f_1^2\;R_{_{P_1}}-f_2^2\;R_{_{P_2}}}{f_1^2-f_2^2}}$

Prefit residual:

$ ResPrefit_{rec}^{sat}= PR_{rec}^{sat} - C_{modeled}$

Modeled range includes estimated "true range" $\rho0_{rec}^{sat}$:

$\color{green} { C_{modeled} = \rho0_{rec}^{sat} - c(\delta t_{sat} + \delta_{rel}) + Trop_{rec}^{sat} + Iono_{rec}^{sat} + TGD}$

The following equation can be written for geometric range, neglecting the multipath and receiver noise terms:

$ PR^j- D^j\simeq\sqrt{(x^j-x)^2+(y^j-y)^2+(z^j-z)^2} +c(\delta t)~~$

$ D^j = - c(\delta t_{sat} + \delta_{rel}) + Trop_{rec}^{j} + Iono_{rec}^{j} + TGD$

$ j=1,2,...,n~~~~ (n \geq 4)~~satellite $

Linearising $\color{red}{\rho_{rec}^{sat}}$ around "an apriori" receiver position $(x_0, y_0, z_0)$

$ \color{blue}{PR_{rec}^{j}}=\color{red}{\rho_0^j+\frac{x_0-x^j}{\rho_0^j} dx +\frac{y_0-y^j}{\rho_0^j} dy+\frac{z_0-z^j}{\rho_0^j}dz} +c (\delta t_{rec} - \delta t^{sat}) + \sum \delta_k \qquad with \qquad dx=x-x_0; dy=y-y_0;dz=z-z_0 \qquad$

We can rearrange equation by separating the known and unknown terms of each side:

$\color{blue}{PR_{rec}^{j}}+\color{green}{{\rho_0^j + c\cdot\delta t^{sat}}-\sum \delta_k} =\color{green}{\frac{x_0-x^j}{\rho_0^j}} dx +\color{green}{\frac{y_0-y^j}{\rho_0^j}} dy+\color{green}{\frac{z_0-z^j}{\rho_0^j}}dz +\color{green}c\cdot\delta t_{rec})\qquad with \qquad dx=x-x_0; dy=y-y_0;dz=z-z_0 \qquad$

For all satellites in view calculate prefit residual:

$\color{green}{Prefit_{rec}^{j}}= \color{green}{PR_{rec}^{j} - C_{modeled}^j} = \color{green}{PR_{rec}^{j}-{\rho_0^j + c\cdot dt^{sat}}-\sum \delta_k} $

$\color{green}{Prefit_{rec}^{j}} = \color{green}{\frac{x_0-x^j}{\rho_0^j}} dx -\color{green}{\frac{y_0-y^j}{\rho_0^j}} dy+ \color{green}{\frac{z_0-z^j}{\rho_0^j}}dz + \color{green}c\cdot dt_{rec}$

We can simplify notation:

$ a_x^j = \frac{x_0-x^j}{\rho_0^j}~~~$ $ a_y^j = \frac{y_0-y^j}{\rho_0^j}~~~$ $ a_z^j = \frac{z_0-z^j}{\rho_0^j}$

$Y = Prefit_{rec}^{j}= {PR_{rec}^{sat} - C_{modeled}}$

Let us assume that we have n satellites visible simultaneously. We can form the design matrix of unit vectors and the matrix of unknowns:

$\vec{Gx} = \begin{bmatrix}a_x^1 & a_y^1 & a_z^1 & c\\a_x^2 & a_y^2 & a_z^2 & c\\a_x^3 & a_y^3 & a_z^3 & c\\a_x^4 & a_y^4 & a_z^4 & c\\ \vdots&\vdots&\vdots&\vdots\\a_x^n & a_y^n & a_z^n & c\end{bmatrix}\begin{bmatrix}\Delta x \\\Delta y\\\Delta z\\\Delta t \end{bmatrix}$

Matrix of observations (Prefit residuals):

$\vec{F} = \begin{bmatrix}Prefit_{rec}^1\\Prefit_{rec}^2\\Prefit_{rec}^3\\Prefit_{rec}^4\\\vdots\\Prefit_{rec}^n \end{bmatrix}$

Estimating the unknowns by LS method:

$ F~= Gx$

$\hat x~= (G^TG)^{-1} \cdot G^T \cdot F$

Iterativ least square solution:

Update $(x_{rec},~y_{rec},~z_{rec})$ and use it to calculate geometric range $\rho0_{rec}^{sat}$ and fill in $(x0_{rec},~y0_{rec},~z0_{rec})$ in the next step of calculation of $C_{modeled}$.

$\begin{bmatrix}x_{rec}\\y_{rec}\\z_{rec}\end{bmatrix} = \begin{bmatrix}x0_{rec}\\y0_{rec}\\z0_{rec}\end{bmatrix} + \begin{bmatrix}dx\\dy\\dz\end{bmatrix}$

In [109]:
xrec = XYZapprox[0]
yrec = XYZapprox[1]
zrec = XYZapprox[2]
xrec = 0
yrec = 0
zrec = 0
cdt_rec = 0
residuals_prev = 0
i = 1
for x in range(11):
    y = []
    G = []
    for mea in satobs_eph:
        p0_rs = geometric_range(xrec, yrec, zrec, mea[9], mea[10], mea[11])
        cdt_sat = C*mea[4]
        tropo = mea[12] 
        c_modeled = p0_rs-cdt_sat+tropo # + cdt_rec
        p1 = mea[1]
        p2 = mea[2]
        ionosfreeP1P2=(F1**2*p1-F2**2*p2)/(F1**2-F2**2)
        prefit  = ionosfreeP1P2 - c_modeled 
        y.append([prefit])
        x1 = xrec-mea[9]/p0_rs
        y1 = yrec-mea[10]/p0_rs
        z1 = zrec-mea[11]/p0_rs
        G.append([x1,y1,z1,C])
    ymat = np.array(y)
    Gmat = np.array(G)
    res = np.linalg.inv(np.dot(Gmat.transpose(),Gmat))
    res = np.dot(res,Gmat.transpose()) 
    res = np.dot(res,ymat)
    xrec = xrec + res[0,0]
    yrec = yrec + res[1,0]
    zrec = zrec + res[2,0]
    cdt_rec = res[3,0] # + cdt_rec 
    residuals = math.sqrt(res[0,0]**2+res[1,0]**2+res[2,0]**2)
    print('iter:', i, math.sqrt(res[0,0]**2+res[1,0]**2+res[2,0]**2) )
    print('residual change: ', residuals-residuals_prev)
    residuals_prev = residuals
    i = i+1
print('Residuals:')
print('dx:', XYZapprox[0]-xrec)
print('dy:', XYZapprox[1]-yrec)
print('dz:', XYZapprox[2]-zrec)
print('Receiver clock bias: ', cdt_rec)
print('X_ECEF:', xrec)
print('Y_ECEF:', yrec)
print('Z_ECEF:', zrec)
point = cartesian_geographic(A_WGS84, B_WGS84, xrec, yrec, zrec)
print('Latitude:',rad_deg(point[1]))
print('Longitude:',rad_deg(point[0]))
print('Ellipsoidal height:',point[2])

iter: 1 7527002.431962806
residual change:  7527002.431962806
iter: 2 940255.84193763
residual change:  -6586746.590025175
iter: 3 220827.24033477184
residual change:  -719428.6016028582
iter: 4 6696.084249741599
residual change:  -214131.15608503023
iter: 5 2677.3390988890087
residual change:  -4018.74515085259
iter: 6 833.2741860400477
residual change:  -1844.064912848961
iter: 7 6.985661685345771
residual change:  -826.2885243547018
iter: 8 0.7385928389945531
residual change:  -6.247068846351217
iter: 9 0.06951584022928747
residual change:  -0.6690769987652656
iter: 10 0.01568315619266222
residual change:  -0.05383268403662525
iter: 11 0.0019362998563241288
residual change:  -0.013746856336338091
Residuals:
dx: 0.30062739574350417
dy: -0.005998189211823046
dz: -2.135871988721192
Receiver clock bias:  -4.097940710301451e-05
X_ECEF: 1962039.9274726042
Y_ECEF: 844038.2488981892
Z_ECEF: 5989770.846871989
Latitude: 70.49577785482565
Longitude: 23.276603121883067
Ellipsoidal height: 40.624387479387224

Show position on map