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

LAB 4

# From GPS ephemerides to ECEF satellite positions

GPS satellite ephemerides are broadcast to users as part of the GPS signal in the navigation message. They are calculated by predicting orbits forward on the basis of a curve-fit to 4 to 6 hours of actual orbit data. GPS broadcast ephemerides are accurate to ~10 m. They include the following parameters:

• toe: epoch of reference of the ephemerides in seconds of the current GPS week. It is the center of the interval over which the ephemeris is valid.
• m0, Dn, e, sqrt(a), W0, i0, w0: orbital elements (= Keplerian parameters) at toe.
• Wdot, idot: rate of right ascension and inclination angle.
• Cwc Cws Crc Crs Cic Cis: corrections coefficients to the argument of perigee, orbit radius and inclination for the perturbations caused by variations in the Earth’s gravity field, solar radiation pressure, and attraction from Sun and Moon.

Assignment:

1.      Write a program to convert the ephemerides given in a RINEX navigation file into an ECEF coordinate system. Input file = epgga2.010:

·        Write a function to read the broadcast ephemeris file. This function should return a matrix with 17 rows (the parameters needed later on in the calculations) and as many columns as there are ephemerides. You will need the 17 following parameters:

 svprn satellite PRN number Mo mean anomaly roota sqrt(semi-major axis) deltan variation of mean angular velocity ecc eccentricity omega0 argument of perigee cuc, cus, crc, crs, cic, cis correction coefficients i0 inclination idot rate of inclination Omega0 right ascension Omegadot rate of right ascension toe time of ephemeris

The call to that function could look like:

eph = read_rinexn('epgga2.010');

·        Write a function to convert the ephemerides matrix (output of previous function) for a given satellite at a given time (i.e. 3 input arguments). This function should return the corresponding X,Y,Z coordinates in ECEF frame, for instance:

[X,Y,Z] = get_satpos(t,sv,eph);

A possible algorithm is given below, together with the formulas that give the values for W, w, i, E, r, and v needed in the equations above:

1. Extract needed parameters from the GPS broadcast ephemerides file. Note that the RINEX navigation files provide information every 2 hours only. Therefore, you first need to find the correct time of ephemeris (toe) in the navigation file given the time for which you need to extrapolate the satellite position (tdata).

2. Compute basic parameters at request time tdata:

Time elapsed since toe:

Mean anomaly at t:

Iterative solution for E:              (use iterative solution)

True anomaly v:

3. Correct for orbital perturbations:

Argument of perigee:

Radial distance:

Inclination:

4. Compute the right ascension, accounting for Earth's rotation (we = mean angular velocity of the Earth) and variations of W (= Wdot):

5. Convert satellite position from orbital frame to ECEF frame:

Write the satellite position in the orbital frame:

Build the rotation matrix to ECEF frame:

Apply rotation:

6. Check your program. Use values from the sp3 file igs11484.sp3 as “ground truth”:

t=346500 and sv = 31 Ž X=11660.379642, Y=11313.211213, Z=-21326.822815

2.      Compute ECEF coordinates for satellite 31 every 15 minutes and compare with the coordinates given in the corresponding sp3 file (precise IGS orbits).

·        Write a function that extract the XYZ position of a satellite ‘sv’ from an sp3 file and returns that position:

[Xs,Ys,Zs,Ts] = read_sp3(sp3file,sv);

·        Compute the 3-D residual position = vector difference between sp3 and broadcast XYZ position. Plot the norm of this vector as a function of time for satellite 31.

3.      Comment on the differences between the 2 orbits.

Short description of the RINEX navigation file format. The format of fully described in http://www.ngs.noaa.gov/CORS/Rinex2.html.

+----------------------------------------------------------------------------+

|                                  TABLE A4                                  |

|           GPS NAVIGATION MESSAGE FILE - DATA RECORD DESCRIPTION            |

+--------------------+------------------------------------------+------------+

|    OBS. RECORD     | DESCRIPTION                              |   FORMAT   |

+--------------------+------------------------------------------+------------+

|PRN / EPOCH / SV CLK| - Satellite PRN number                   |     I2,    |

|                    | - Epoch: Toc - Time of Clock             |            |

|                    |          year (2 digits, padded with 0   |            |

|                    |                if necessary)             |  1X,I2.2,  |

|                    |          month                           |   1X,I2,   |

|                    |          day                             |   1X,I2,   |

|                    |          hour                            |   1X,I2,   |

|                    |          minute                          |   1X,I2,   |

|                    |          second                          |    F5.1,   |

|                    | - SV clock bias       (seconds)          |  3D19.12   |

|                    | - SV clock drift      (sec/sec)          |            |

|                    | - SV clock drift rate (sec/sec2)         |            |

+--------------------+------------------------------------------+------------+

| BROADCAST ORBIT - 1| - IODE Issue of Data, Ephemeris          | 3X,4D19.12 |

|                    | - Crs                 (meters)           |            |

|                    | - Delta n             (radians/sec)      |            |

|                    | - M0                  (radians)          |            |

+--------------------+------------------------------------------+------------+

| BROADCAST ORBIT - 2| - Cuc                 (radians)          | 3X,4D19.12 |

|                    | - e Eccentricity                         |            |

|                    | - Cus                 (radians)          |            |

|                    | - sqrt(A)             (sqrt(m))          |            |

+--------------------+------------------------------------------+------------+

| BROADCAST ORBIT - 3| - Toe Time of Ephemeris                  | 3X,4D19.12 |

|                    |                       (sec of GPS week)  |            |

|                    | - Cic                 (radians)          |            |

|                    | - OMEGA               (radians)          |            |

|                    | - CIS                 (radians)          |            |

+--------------------+------------------------------------------+------------+

| BROADCAST ORBIT - 4| - i0                  (radians)          | 3X,4D19.12 |

|                    | - Crc                 (meters)           |            |

|                    | - omega               (radians)          |            |

|                    | - OMEGA DOT           (radians/sec)      |            |

+--------------------+------------------------------------------+------------+

| BROADCAST ORBIT - 5| - IDOT                (radians/sec)      | 3X,4D19.12 |

|                    | - Codes on L2 channel                    |            |

|                    | - GPS Week # (to go with TOE)            |            |

|                    |   Continuous number, not mod(1024)!      |            |

|                    | - L2 P data flag                         |            |

+--------------------+------------------------------------------+------------+

| BROADCAST ORBIT - 6| - SV accuracy         (meters)           | 3X,4D19.12 |

|                    | - SV health        (bits 17-22 w 3 sf 1) |            |

|                    | - TGD                 (seconds)          |            |

|                    | - IODC Issue of Data, Clock              |            |

+--------------------+------------------------------------------+------------+

| BROADCAST ORBIT - 7| - Transmission time of message        *) | 3X,4D19.12 |

|                    |         (sec of GPS week, derived e.g.   |            |

|                    |    from Z-count in Hand Over Word (HOW)  |            |

|                    | - Fit interval        (hours)            |            |

|                    |         (see ICD-GPS-200, 20.3.4.4)      |            |

|                    |   Zero if not known                      |            |

|                    | - spare                                  |            |

|                    | - spare                                  |            |

+--------------------+------------------------------------------+------------+